Estudo computacional da dependência dos parâmetros termodinâmicos em sistemas de bósons

Wenderson Rodrigues Fialho da Silva - Viçosa, julho de 2022.

    Os bósons são partículas cujo spin pode assumir somente valores inteiros, os quais não obedecem ao princípio de exclusão de Pauli, portanto, pode-se encontrar mais de um bóson no mesmo estado quântico. O potencial químico dos bósons apresentam valores negativos ou nulo, sendo esse ultimo relacionado ao condensado de Bose-Einstein. Esse trabalho teve como objetivo estudar por meio de gráficos as dependência dos parâmetros termodinâmicos de bósons, como a fugacidade $z$, energia interna U e capacidade térmica $C_{v}$, ambos em função de $T/T_{c}$, onde $T$ é a temperatura e $T_{c}$ a temperatura crítica do sistema. 

    Desenvolveu-se um programa de computador em linguagem Python que gera os dados referentes as curvas $z$, $U$ e $C_{v}$, ambas em função de $T/T_{c}$. As equações (1, 2 e 3) foram utilizadas para obtenção dessas curvas. A demostração da equação (3) é apresentada no apêndice A.

$$ \frac{T}{T_{c}} = \left ( \frac{\zeta(\frac{3}{2})}{Li_{3/2}(z)} \right)\qquad(1)$$

$$U = \frac{3}{2}NK_{B}T\frac{Li_{5/2}(z)}{Li_{3/2}(z)}\qquad(2)$$

$$ C_{v} = \begin{cases}\frac{15}{4}\frac{\zeta(\frac{5}{2})}{\zeta(\frac{3}{2})}NK_{B}\left ( \frac{T}{T_{c}} \right)^{3/2} se\qquad T < T_{c}\qquad\qquad\qquad\qquad(3) \\\frac{3}{2}NK_{B}\left ( \frac{5}{2}\frac{Li_{5/2}(z)}{Li_{3/2}(z)} - \frac{3}{2}\frac{Li_{3/2}(z)}{Li_{1/2}(z)} \right)\qquad se\qquad T > T_{c}\end{cases}$$

    A equação (1) e (2) são apresentadas no capítulo 30 na referencia [1]. Com auxílio do programa Origin pro 8.5, fez-se os gráficos com os dados gerados. Desenvolveu-se uma rotina que realiza-se os cálculos das equações das equações (1), (2) e (3) e salva-se tais dados em arquivos $.txt$, os quais foram utilizados para fazer os gráficos. O código de programação em linguagem Python é apresentado no apêndice B.

    O resultado com as curvas obtidas é apresentado nos gráficos das Figuras (1) a seguir:

Figura 1. (A) fugacidade, (B) energia interna e em (C) capacidade térmica de bósons em função da temperatura sobre a temperatura crítica.

    A energia interna está normalizada pelo número de partículas e pela energia térmica. A capacidade térmica pelo número de partículas e pela constante de Boltzmann. Para obtenção das curvas para o potencial químico de bósons e férmions, utilizou-se as equações para $v_{B}$ e  $v_{F}$ abaixo, obtidas da referência [2]. O código de programação em Python utilizado é apresentado no apêndice C.

    No gráfico (A) da Figura 1, como descrito na referencia [1] e em sala, a fugacidade para $T \leq T_{c}$ é igual a um. Para $T > T_{c}$ ela cai continuamente. Em (B), para a energia interna, a linha pontilhada representa a previsão dada pelo teorema da equipartição de energia. Para altas temperaturas, ou seja, $T >> T_{c}$ as curvas tendem a se encontrar, onde o caráter quântico do sistema não tem grandes contribuições, e a abordagem clássica já descreve o problema. Em (C), abaixo da temperatura crítica $T_{c}$, o comportamento quântico relacionado a capacidade térmica tem que ser levando em consideração para a descrição do sistema. Já para $T >> T_{c}$, o sistema se comporta como previsto pelo teorema da equipartição de energia.

    Para obtenção das curvas dos potenciais químicos para bósons e para férmions foi utilizado as equações da referência [2].

    Para bósons, foi utilizado a equação (4) abaixo, sendo $v_{B}=\frac{u_{B}}{T_{0}}$.

$$ v_{B} = \begin{cases}0, \qquad T\leq 1 \\\sum_{k=1}^{4} a_{k}(t-1)^k, \qquad 1 \leq T \leq 3.68 \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(4) \\T\ln(\zeta(3/2)+T^{-3/2}) + T\ln(1+\zeta(3/2)(2T)^{-3/2}), \qquad T \ge 3.68\end{cases}$$

    Para férmions, foi utilizado a equação (5) abaixo, sendo $v_{F}=\frac{u_{F}}{T_{F}}$

$$ v_{F} = \begin{cases}1+\sum_{k=1}^{4} a_{k}T^k, \qquad 0 \leq T \leq 1.36 \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(5) \\T\ln \left ( \frac{2}{3\Gamma(3/2)} \right ) -T^{-3/2}) + T\ln\left ( 1-\frac{2}{3\Gamma(3/2)}(2T)^{-3/2} \right ), \qquad T \ge 1.36\end{cases}$$

    Os coeficientes dos 4 primeiros temos da soma  do cálculo do potencial químico são apresentados na Tabela 1 abaixo.

Tabela 1. Valores dos coeficiente utilizado para o calculo dos potenciais quimicos para bósons e fémions [2].

    Fazendo $\mu_{B}=v_{B}T_{0}$, pode-se as curvas do potencial químico para bósons $\mu_{B}$ e para $\mu_{F}=v_{F}T_{0}$, obteve-se as curvas para férmios $\mu_{F}$. Os gráficos obtidos são apresentados na Figura (2) abaixo.

Figura 2. Potenciais químicos para bósons (A) e férmions (B) em função da temperatura.

    Para $T \leq T_{0}$, o gráfico (A) para bósons da Figura 2 nos mostra que o potencial químico é zero. Para valores de $T > T_{0}$, o potencial químico assume valores negativos. Já para os férmions, o potencial químico apresentado no gráfico (B) assume valores positivos para $T < T_{F}$, e negativos para $T > T_{F}$.

    Afim de poder-se comprar os valores de $U$ para $T < T_{c}$ com a curva obtidas pela equação (2), sendo ela uma aproximação da equação (6), fez-se o gráfico da Figura (3) a seguir, onde, na equação (6), $n_{k} = \frac{1}{e^{\frac{T_{c}}{K_{B}T}(E_{k}-\mu)}-1}$, $E_{k} = \frac{\hbar^{2}k^{2}}{2m}$ e $\mu = 0$ e fazendo as constantes igual a um. O código de programação em Python utilizado é apresentado no apêndice D.

$$ U = \sum_{k}n_{k}E_{k}\qquad(6) $$

Figura 3. Curvas para Energia interna $U$ obtidas da equação (1), a qual representa a aproximação do somatório da equação (6) para a integral, bem como a curva da para o somatório da equação (6).

    Nota-se que a uma semelhança entre as curvas obtidas para energia interna $U$ para $T < T_{c}$, entretanto, quando $T$ se aproxima de $T_{c}$, observa-se que o valor previsto pelo somatório começa a adquirir valores menores em relação a integral. Para valores de $T < T_{c}$, os bósons tendem a se acumularem no nível de energia mais baixo, formando o condensado de Bose-Einstein. A aproximação feita do somatório para integral não é boa quando energia térmica $(K_{B}T)$ é da ordem da energia característica do sistema, e caráter quântico associado deve ser levado em consideração. Quando $(K_{B}T)$ é  muito maior que a energia característica do sistema quântico, podemos fazer a aproximação e tratar o sistema como um contínuo de energia, de modo que o uso da integral não trará grandes desvios. 

Apêndice A - Dedução da equação (3) utilizada para obtenção da curva da capacidade térmica $C_{v}$

Para $T < T_{c}$:

$$U=\frac{3}{2}NK_{B}T_{c}\frac{\zeta(\frac{5}{2})}{\zeta(\frac{3}{2})} \left ( \frac{T}{T_{c}} \right )^\frac{5}{2}, \qquad(7) $$

logo:

$$C_{v}=\left ( \frac{\partial U}{\partial T} \right )_{v}=\frac{15}{4}NK_{B}\frac{\zeta(\frac{5}{2})}{\zeta(\frac{3}{2})}T^\frac{3}{2}\qquad(8)$$

Para $T > T_{c}$, usando a equação 30.51 e $ \lambda _{T}^3=KT^\frac{3}{2} $ podemos escrever $U$ como:

$$U=\frac{3}{2}K_{B} \left ( 2S+1 \right )KVT^\frac{5}{2}Li_{5/2}(z),\qquad(9) $$

logo:

$$C_{v}=\left ( \frac{\partial U}{\partial T} \right )_{v}=\frac{15}{4}K_{B} \left ( 2S+1 \right ) KVT^\frac{5}{2}Li_{5/2}(z) + \frac{3}{2}K_{B} \left ( 2S+1 \right ) KVT^\frac{5}{2}\frac{dLi_{5/2}(z)}{dz}.\qquad(10)$$

Derivando a equação C.32 do apendice C da referencia [1], temos que $ \frac{dLi_{n}}{dz}=\frac{L_{i_{n-1}}}{z}\frac{dz}{dT}$. Portanto,

$$C_{v}=\frac{15}{4}K_{B}\left (2S+1 \right )KVT^\frac{3}{2}Li_{5/2}(z) + \frac{3}{2}K_{B}\left (2S+1 \right )KVT^\frac{5}{2}Li_{3/2}(z)\frac{1}{z}\frac{dz}{dT}.\qquad(11)$$

Como: $Li_{5/2}(z) \propto T^{-3/2}$,

$$\frac{d}{dT}(Li_{3/2}(z)) = Li_{1/2}(z)\frac{dz}{dT} = -\frac{3}{2}\beta T^{-5/2}.\qquad(12)$$

Portanto:

$$Li_{1/2}(z)\frac{dz}{dT} = -\frac{3}{2T}Li_{3/2}(z) \longrightarrow \frac{dz}{dT} = -\frac{3z}{2T}\frac{Li_{3/2}(z)}{Li_{1/2}(z)}.\qquad(13)$$

Substituindo em $C_{v}$:

$$C_{v} = \left(\dfrac{15}{4}K_{B}(2S+1)KVT^{3/2}Li_{5/2}(z) + \frac{3}{2}K_{B}(2S+1)KVT^{3/2}\left(\dfrac{-3}{2T} \right)\right)\frac{Li^{2}_{5/2}(z)}{Li_{1/2}(z)}$$

$$C_{v} = \frac{3}{2}K_{B}(2S+1)KVT^{3/2}\Bigg(Li_{5/2}(z)\frac{5}{2} - \frac{3}{2}\frac{Li^{2}_{5/2}(z)}{Li_{1/2}(z)}\Bigg).\qquad(14)$$

E como:

$$Li_{3/2}(z) = \frac{N}{V(2S+1)KT^{3/2}} = \frac{N\lambda^{3}_{T}}{V(2S+1)} \longrightarrow \frac{N}{Li_{3/2}(z)} = \frac{3}{2}K_{B}(2S+1)KT^{3/2}V,\qquad(15) $$

finalmente:

$$C_{v} = N\Bigg(\frac{5}{2}\frac{Li_{5/2}(z)}{Li_{3/2}(z)} - \frac{3}{2}\frac{Li_{3/2}(z)}{Li_{1/2}(z)}\Bigg).\qquad(16)$$

Apêndice B - Código para obtenção das curvas de $z$, $U$ e $C_{v}$

\begin{lstlisting}[language=Python]

# -*- coding: utf-8 -*-

"""

Created on Sun Jul 17 11:24:21 2022

@author: Wenderson R. F. da Silva

"""

from mpmath import*

import numpy  as np

pontos = 1000

n1 = 5/2

n2 = 3/2

n3 = 1/2

z = 0.0000001

delta = 0.001

Q = T = U = Cv = 0

deltaT = 0.002612

with open('C_termica.txt','w') as C_termica:

    with open('E_interna.txt','w') as E_interna:  

        with open('fugacidade.txt','w') as fugacidade:

            for j in range(0, pontos):

                Li1 = Li2 = Li3 = 0

                for i in range(1, 10000):

                    Li1 = Li1 + z**i/i**n1

                    Li2 = Li2 + z**i/i**n2

                    Li3 = Li3 + z**i/i**n3

                if T < 1: # T<Tc                               

                    U = 0.77*T**(5/2)                   

                    cv = 1.92*T**(3/2)

                    fugacidade.write('%.4f %7.5f\n' % (T, 1))

                    E_interna.write('%.4f %7.5f\n' % (T, U))

                    C_termica.write('%.4f %7.5f\n' % (T, cv))

                else:

                    Q = (2.61/Li2)**(2/3)   # Q = T/Tc para zeta(3/2)= 2.612

                    U = (3/2)*Q*(Li1/Li2)                       

                    cv = 3.75*(Li1/Li2) - (2.25*(Li2/Li3))    

                    fugacidade.write('%.4f %7.5f\n' % (Q, z))

                    E_interna.write('%.4f %7.5f\n' % (Q, U))

                    C_termica.write('%.4f %7.5f\n' % (Q, cv))

                z = z + delta

                T = T + deltaT

        fugacidade.close()

    E_interna.close()

C_termica.close()

\end{lstlisting}

Apêndice C - Código para obtenção das curvas de $\mu_{bóson}$ e $\mu_{férmion}$

\begin{lstlisting}[language=Python]

import numpy  as np

from numpy import log as ln

from numpy import pi

pontos = 1001

T = 0.0001

deltaT = 0.004

with open('vb.txt','w') as vb: 

    for j in range(1, pontos):

        

        if T <= 1:

            vb.write('%.4f %7.5f\n' % (T, 0))            

        elif T > 1 and T <= 3.68:

            b = - 0.016*(T-1) - 1.064*(T-1)**2 + 0.299*(T-1)**3 - 0.043*(T-1)**4          

            vb.write('%.4f %7.5f\n' % (T, b))   

        else:         

            b = T*ln(2.612*T**(-3/2)) + T*ln(1 - 2.612*(2*T)**(-3/2))         

            vb.write('%.4f %7.5f\n' % (T, b))          

        T = T + deltaT

vb.close()

T=0

G = 0.88622692

with open('vf.txt','w') as vf: 

    for j in range(1, pontos):     

        if T <= 1.36:         

            f = 1 + 0.016*T - 0.957*T**2 - 0.293*T**3 + 0.209*T**4            

            vf.write('%.4f %7.5f\n' % (T, f))   

        else:           

            f = T*ln(2/(3*G)*T**(-3/2)) + T*ln(1 + 2/(3*G)*(2*T)**(-3/2))           

            vf.write('%.4f %7.5f\n' % (T, f))          

        T = T + deltaT

vf.close()

\end{lstlisting}

Apêndice D - Código para obtenção das curvas de $U$ pela equação (6)

\begin{lstlisting}[language=Python]

import numpy  as np

from numpy import log as ln

from numpy import exp

pontos = 1000

T = 0.0001

deltaT = 0.001

with open('Ub.txt','w') as Ub:

    for i in range (1, pontos):

        U = 0

        for k in range(1, 100):    

            U += (k**2) / (exp((k**2)/T) - 1)            # u = 0           

        Ub.write('%.4f %7.5f\n' % (T, U))             

        T = T + deltaT        

Ub.close()

\end{lstlisting}

Referências

[1] S. J. Blundell, K. M. Blundell. Concepts in Thermal Physics, 2ª ed. Oxford University Press, 2010.

\newline

[2] A. G. Sotnikov. Chemical potentials and thermodynamic characteristics of ideal Bose- and Fermi-gases in the region of quantum degeneracy, Low Temperature Physics 43, 144, 2017.

Construção de curvas isoentálpicas e curva de inversão para o gás de van der Waals em uma expansão de Joule-Kelvin

                                                                                      Wenderson Rodrigues F. da Silva - Viçosa, 07 de abril de 2022. 

Objetivo

        Esse trabalho tem como objetivo construir curvas isoentálpicas para o gás de van der Waals em uma expansão de Joule-Kelvin, bem como a curva de inversão curva de inversão. Esse trabalho foi desenvolvido para disciplina de termodinâmica estatística I da UFV.

Introdução

        Em processos onde se deseja liquefazer gases a expansão de joule-kelvin é um mecanismo de grande aplicabilidade. Uma vez que o gás, ao passar por uma constrição, pode resfriar, com o mecanismo adequado, pode-se projetar dispositivos eficientes destinados a essa função. Para tanto, faz-se necessário entender acerca dos conhecimentos teóricos envolvidos no processo de resfriamento e que pode levar a liquefação do gás, e o coeficiente de joule-kelvin (equação 1 abaixo), que descreve como o gás muda de temperatura quando reduzimos sua pressão, a entalpia constante, é o ferramental matemático que nos auxilia nesse entendimento.

$$\mu_{JK} = \left ( \frac{\partial P}{\partial T} \right)_{H} = \frac{1}{C_{p}}\left [ T\left ( \frac{\partial V}{\partial T} \right)_{P} -V\right ].  \qquad(1)$$

        Como se trata de um processo que ocorre a entalpia constante, num gráfico $T^{*} x P^{*}$, uma vez que o coeficiente de joule-kelvin nos dá a variação da temperatura com a pressão para entalpia constante, quando esse é zero (um ponto crítico), tal coeficiente dará informações de interesse tanto teórico como prático, pois, para diferentes valores de entalpia, esses pontos estarão contidos numa curva denominada curva de inversão, por meio da qual pode-se conhecer a região de resfriamento do gás. Aqui, utilizaremos o modelo do gás de Van der Waals, com objetivo de construir curvas isoentálpicas para a expansão de Joule-Kelvin.

Metodologia

        Desenvolveu-se um programa de computador em linguagem Python que gera os dados reverentes as curvas isoentálpicas, bem como a curva de inversão. As curvas isoentálpicas foram obtidas usando as equações paramétricas (2) e (3), onde $\textit{T*}$ e $\textit{P*}$ são a temperatura reduzida e a pressão reduzida, respectivamente, e $\textit{V*}$ é o parâmetro. A demostração da equação (3) é apresentada no apêndice A. A equação (2) foi resolvida no exemplo 26.3.

$$ p^{*} = \frac{8T^{*}}{3V^{*} - 1} - \frac{3}{{(V^{*}})^{2}}, \qquad(2)$$

$$ T^{*} = \frac{3V^{*} - 1}{4(5V^{*} - 1)}\left ( \frac{H}{P_{c}V_{c}} + \frac{6}{V^{*}} \right ).\qquad(3)$$

        Para diferentes valores de entalpia, criou-se 10 curvas, com $\frac{H}{p_{c}V_{c}}$ variando de 4 a 40, uma vez que, por se tratar de um mesmo gás em análise, $p_{c}$ e $V_{c}$ são também constantes. Um loop com o comando $\textit{while}$ foi usado para isso. A curva de inversão foi obtida analiticamente e, com a equação resolvida, implementou-se uma rotina no código para gerar os pontos correspondente. Com base nas equações (2) de estado para o gás de van der Waals, fazendo $\left ( \frac{\partial P}{\partial T} \right)_{H} = 0$ e $h = \frac{H}{P_{C} V_{C}}$, temos que:

$$\mu_{JK} = \left ( \frac{\partial P}{\partial T} \right)_{H} = 0 \qquad \rightarrow \qquad \left ( \frac{\partial T}{\partial V} \right)_{P} = \frac{T^{*}}{V^{*}},  \qquad(4) $$

da equação (2), isolando $T^{*}$ e derivando:

$$\left ( \frac{\partial T}{\partial V} \right)_{P} = \frac{1}{8} \left(3P^{*} +  \frac{6}{(V^{*})^{3}} - \frac{9}{(V^{*})^{2}} \right) = \frac{T^{*}}{V^{*}},\qquad(5)$$

substituindo na equação (5) a equação (2), obtemos:

$$T^{*} = \frac{3(V^{*}-1)^{2}}{4(V^{*})^{2}},\qquad(6)$$

daí, isolando V:

$$V^{*} = \frac{1}{3-2\sqrt{\frac{T^{*}}{3}}}.\qquad(7)$$

        Por fim, substituindo a equação(7) na equação(2), obtemos:

$$P^{*} = 9\left ({3-2\sqrt{\frac{T^{*}}{3}}}\right)\left ({2\sqrt{\frac{T^{*}}{3}}-1}\right),\qquad(8)$$

que representa a equação para a curva de inversão num grafico $\textit{T*}$ x $\textit{P*}$. O valores obtidos para cada uma das curvas estudadas e para a curva de inversão foram escritos em arquivos $\textit{.txt}$ separados, os quais foram utilizados para construção dos gráficos com auxílio do software $\textit{Origin}$.

Resultados e Discussões

        Com base nos dados gerados fez-se o grafico da figura (1) e (2) abaixo, que são referentes as curvas isoentálpicas e a curva de inversão do gás de van der Waals em uma expansão de Joule-Kelvin. 

Figura 1. Curvas isoentálpicas juntamente com a curva de inversão (curva pontilhada verde) para o gás de Van der Waals. Os pontos em verde sobre a intersecção das curva isoentálpicas e a curva de inversão correspondem as pontos de inversão.

Figura 2. Curva de inversão para dois valores de $h$. Em (A), $h = 4$  com $T \simeq 1,2T_{c}$ e, em (B), $h = 40$, com $T \simeq 6T_{c}$. Plotando as curvas isoentálpicas separadamente, a visualização do ponto de inversão (ponto verde) fica mais evidente.

    Se a pressão do gás de entrada for maior que a pressão no ponto de inversão, ao expandir, o gás irá aquecer. Por outro lado, se a pressão for menor do que a pressão no ponto de inversão, ao expandir, o gás se resfriará. Em outras palavras, o resfriamento ocorre à esquerda da curva de inversão quando $\left ( \frac{\partial T}{\partial P} \right)_{H} > 0$ e o aquecimento à direita, onde $\left ( \frac{\partial T}{\partial P} \right)_{H} < 0$ . Portanto, com o conhecimento acerca do fenômeno e tendo em vista que poderá ocorrer aquecimento do gás, se objetivo é que ele resfrie, é necessário que o gás a pressão mais alta que entra na máquina tenha uma pressão menor do que pressão do ponto de inversão. 

Apêndice A - Dedução da equação (3)

Partida da equação para energia interna para o gás de van der Waals:

$$ U = F + TS = \frac{3NK_{B}T}{2} - \frac{N^{2}a}{V}$$

Como $ N = \frac{8P_{c}V_{c}}{3T_{c}}$ , $b = \frac{T_{c}}{8P_{c}}$

e $a = \frac{27T_{c}^{2}}{64P_{c}}$, e temos também que: $T^{*} = \frac{T}{T_{c}}$ , $P^{*} = \frac{P}{P_{c}}$ e $V^{*} = \frac{V}{V_{c}}$.

Logo:

$$ U = P_{c}V_{c} \left (4T^{*} - \frac{3}{V^{*}} \right).\qquad(A1)$$

Como $H = U + PV$, podemos escrever $H$ como:

$$ H = P_{c}V_{c} \left (4T^{*} - \frac{3}{V^{*}} \right) + P_{c}P^{*}V_{c}V^{*},\qquad(A2)$$

o que resulta em:

$$ H = P_{c}V_{c} \left (\left (4T^{*} - \frac{3}{V*} \right) + P^{*}V^{*}\right).\qquad(A3)$$

Logo, da equação (2), podemos expressar $H = H(T^{*}, V^{*})$:

$$ \frac{H}{P_{c}V_{c}} =\left (\left (4T^{*} - \frac{3}{V^{*}} \right) + \left (\frac{8T^{*}}{(3V^{*}-1)} - \frac{3}{V*^{2}} \right)V^{*}\right),\qquad(A4)$$

assim, obtemos:

$$ \frac{H}{P_{c}V_{c}} = 4T^{*}-\frac{3}{V^{*}} + \frac{8T^{*}V^{*}}{\left (3V^{*}-1\right)} -\frac{3}{V^{*}} = 4T^{*} + \frac{8T^{*}V^{*}}{\left (3V^{*}-1\right)} -\frac{6}{V^{*}}.\qquad(A5) $$

Simplificando a expressão, obtemos:

$$ \frac{H}{P_{c}V_{c}} = \frac{12T^{*}V^{*} + 8T^{*}V^{*} - 4T^{*}}{\left (3V^{*}-1\right)} - \frac{6}{V^{*}} \qquad \rightarrow \qquad \frac{H}{P_{c}V_{c}} = \frac{4T^{*}(5V^{*} - 1)}{\left (3V^{*}-1\right)} - \frac{6}{V^{*}}.\qquad(A6)$$

Isolando $T^{*}$, obtemos a equação (3)

$$ T^{*} = \frac{3V^{*} - 1}{4(5V^{*} - 1)}\left ( \frac{H}{P_{c}V_{c}} + \frac{6}{V^{*}} \right ).\qquad(3)$$

(obs: A dedução da origem de U para o gás de VDW está relacionada com a dedução do exercício 26.5, o qual não foi exigida demostração por envolver a função de partição que não estudamos nesse curso. Por esse motivo, iniciei a dedução partindo do resultado de U. As deduções foram baseadas em [1]).

Apêndice B - Código-fonte

#import matplotlib.pyplot as plt

from math import sqrt

h = 4.0

#h = float(input('H / (p_c V_c) = '))

V = 0.4 #float(input('Valor inicial de V* = '))

delta = 0.01 #float(input('Tamanho dos subintervalos = ')

n = 1000

while h <= 40:

    # Gás de Van der Waals

    with open(f'T2_VDW_H={h}.txt', 'w') as file:

        for i in range(0, n):

            T = (3 * V - 1)*(h + 6 / V) / (4 * (5 * V - 1))

            p = (8 * T) / (3 * V - 1) - 3 / (V ** 2)

            V = V + delta

            file.write('%.4f %7.4f\n' % (p, T))

    file.close()

    h = h + 4

# Escreve a curva de inversão

T = 0.4 # valor escolhido para melhor ajustar as isoentalpicas

with open(f'T2_inv.txt', 'w') as file:

    for i in range(0, n):

        p = 9 * (3 - 2 * sqrt(T / 3)) * (2 * sqrt(T / 3) - 1)

        T = T + delta

        file.write('%.4f %7.4f\n' % (p, T))

file.close()

#plt.plot(p, T)

#plt.show()

Referências

[1] JOHNSTON, D. C.. Thermodynamic Properties of the van der Waals Fluid. Department of Physics and Astronomy, Iowa State University, Ames, Iowa, USA. (2014). Disponível em: \https://arxiv.org/pdf/1402.1205.pdf.

Construção de isotermas dos gases de van der Waals e de Dieterici

(Construction of the van der Waals and Dieterici gas isotherms)

Wenderson Rodrigues F. da Silva - Viçosa, 30 de setembro de 2021.

Objetivos

        Esse trabalho tem como objetivo construir os gráficos $PxV$ das isotermas e da função livre de Gibbs $G$ em função da pressão $P$ para os gases de van der Waals e de Dieterici, bem como identificar neles a temperatura crítica e a linha de coexistência de fases.

Introdução

        Quando se deseja estudar as propriedades termodinâmicas de um gás a escolha de um modelo matemático adequando é determinante para obter boa precisão na descrição do sistema. No estudo de gases muito rarefeitos (baixa pressão) e a temperaturas elevadas, um modelo simples é o de um gás ideal, mostrado na equação (1):

$$PV_m = RT,\qquad(1)$$

onde $P$ é a pressão, $V_m$ é o volume molar, $T$ a temperatura e $R$ a constante real dos gases. Tal modelo não leva em consideração o volume ocupado pela molécula, bem como a interação entre elas.

        Um modelo mais completo que consegue descrever com mais exatidão as propriedades termodinâmicas dos gases é apresentado na equação (2) abaixo, denominado gás de van der Waals, em que é levado em consideração as interações intermoleculares (parâmetro $a$), assim como o volume ocupado pela molécula (parâmetro $b$, que leva em conta o tamanho da molécula do gás, limitando sua liberdade em mover-se em todo volume do recipiente que o contem).

$$\left(P + \dfrac{a}{V_m^2}\right) (V_m - b) = RT.\qquad(2)$$

        Uma reformulação da equação acima, leva ao modelo de e Dieterici, proposta por Conrad Dieterici em 1899.

$$P (V_m - b) = RT e^{-a / RTV_m}.\qquad(3)$$

        Aqui será construir gráficos da função livre de Gibbs $G$ em função da pressão $P$ para o gás de van der Waals e para o gás de Dieterici, identificando a temperatura crítica (temperatura acima da qual não pode coexistir, em equilíbrio, as fases liquidas e vapor para uma substância) no gráfico e a linha de coexistência de fases.

Metodologia

        Implementou-se um programa em Python 3.7 para obter as isotermas e a função livre de Gibbs para o gás de van der Waals e para o gás de Dieterici. Para isso, utilizou-se das equações (2) e (3) com as coordenadas reduzidas $T^* = \dfrac{T}{T_c}$, $P^* = \dfrac{P}{P_c}$ e $V^* = \dfrac{V}{V_c}$, onde $T_c$, $P_c$ e $V_c$ são a temperatura, pressão e volume reduzidos, respectivamente. Então, para cada valor de $T$ escolhido para análise, fez-se também o cálculo da integral $\Delta G = \int Vdp$ via integração numérica utilizando o método dos trapézios.

        Desse modo, fez-se com que os valores obtidos de $p$, $V$ e $\Delta G$ fossem escritos em três colunas em arquivos separados, um para o gás de van der Waals e outro para o gás de Dieterici.

Resultados e Discussões

        Com os valores obtidos de $p$, $V$ e $\Delta G$, um para o gás de van der Waals e outro para o gás de Dieterici, e com o auxílio do software Origin, plotou-se os gráficos das Figuras (1), (2) e (3) abaixo.

        Na Figura (1), pode-se observar que, para certas regiões do diagrama $\Delta G \times p$, a energia livre de Gibbs assume mais de um valor. Nesses casos, os estados estáveis são aqueles para os quais a energia livre de Gibbs é mínima. 

Figura 1 - Isotermas para os gases de (A) van der Waals e (C) Dieterici. Em (B) e (D), vemos os respectivos comportamentos de $\Delta G$ em função de $p$.

        Na Figura (2), os dois pontos destacados na isoterma no diagrama $p \times V$, ligados pela reta horizontal (construção de Maxwell), possuem um mesmo valor de $\Delta G$, indicado pelo ponto B, e representam duas fases em equilíbrio. Ao repetirmos este procedimento para os diferentes valores de $T*$, podemos determinar a região de coexistência das fases.

Figura 2 - Isoterma de van der Waals para a temperatura reduzida $T^* = 0,9$. A reta tracejada, que liga os pontos em equilíbrio, corresponde a $p = 0,64$.

        Os valores de $p$ utilizados para traçar as retas horizontais nos gráficos da Figura (3) abaixo foram obtidos dos gráficos de $\Delta G$, na Figura (2) (B e D), levando o cursor até o ponto de equilíbrio de fases. Se traçarmos uma curva ligando os pontos destacados no gráfico, teremos uma região delimitada de coexistência das fases líquida e gasosa.

Figura 3 - Isoterma de van der Waals para a temperatura reduzida $T^* = 0,9$. A reta tracejada, que liga os pontos em equilíbrio, corresponde a $p = 0,64$.

Código em Python 

from math import exp

T = float(input('T* = '))

V = 0.51                                 #float(input('Valor inicial de V* = '))

delta = 0.0005                        #float(input('Tamanho dos subintervalos = '))

n = 10000

# Gás de Van der Waals

V1 = V

area = 0

with open(f'T3_VDW_T={T}.txt', 'w') as file:

    for i in range(0, n):

        p1 = (8 * T) / (3 * V1 - 1) - 3 / (V1 ** 2)

        V2 = V1 + delta

        p2 = (8 * T) / (3 * V2 - 1) - 3 / (V2 ** 2)

        h = p2 - p1

        area += (V1 + V2) * h / 2

        file.write('%.4f %7.4f %7.4f %7.4f \n' % (V1, p1, p1, area))

        V1 = V2

file.close()

# Gás de Dieterici

V1 = V

area = 0

with open(f'T3_Dieterici_T={T}.txt', 'w') as file:

    for i in range(0, n):

        p1 = (T * exp(2 * (1 - 1 / (T * V1)))) / (2 * V1 - 1)

        V2 = V1 + delta

        p2 = (T * exp(2 * (1 - 1 / (T * V2)))) / (2 * V2 - 1)

        h = p2 - p1

        area += (V1 + V2) * h / 2

        file.write('%.4f %7.4f %7.4f %7.4f \n' % (V1, p1, p1, area))

        V1 = V2

file.close()

Modelagem matemática para predição de casos da Covid-19 no Brasil com e sem a adesão da população à Quarentena.

Mathematical modeling for the prediction of Covid-19 cases in Brazil with and without population adhesion to Quarantine.

Wenderson Rodrigues F. da Silva - Viçosa, 13 de abril de 2020.

OBJETIVO         

         Estudar um modelo matemático de uso epidemiológico para propagação de doenças infecciosas denominado SIR, afim de reforçar a importância da quarentena, bem como prever o número de casos de Covid-19 no Brasil para os próximos 10 dias, por meio do ajuste de uma curva aos dados de pessoas infectadas no país e pelo modelo SIR.

INTRODUÇÃO           

      Modelos matemáticos são primordiais para se estudar a evolução de sistemas. As ciências da natureza, como a  Física, Química, Meteorologia e Biologia, assim como todos os ramos de estudos que se deseje modelar, ou seja, desenvolver um conjunto de regras que ordene a forma como a dinâmica de um sistema evolui no tempo e/ou no espaço, se fundamenta no ferramental matemático para tal. E a teoria que serve de base para se investigar comportamentos de casos onde se observa a mudança de algo no tempo é a descrito pelas Equações Diferenciais, sejam elas ordinárias (EDOs) quando o sistema evoluí apenas temporalmente, ou as Equação Diferencias Parciais (EDPs) quando se observa evolução no tempo e no espaço.

      A epidemiologia é a ciência que lida com a propagação de doenças infecciosas e faz uso das Equações Diferencias para conhecer e prever o estágio de evolução da disseminação de um dado agente infeccioso, afim de possibilitar o controle da doença bem como nortear a respeito das melhores tomadas de decisões para conter seu avanço. Aqui abordaremos um modelo para propagação de doenças infecciosas chamado SIR que faz uso de um sistema de EDOs para modelagem dos casos, para prever a evolução de casos confirmados da Covid-19 no Brasil, bem como prever qual será esse número dentro de 10 dias.

      Um modelo epidêmico simples e bem conhecido no meio dos epidemiologistas é o chamado Modelo SIR (susceptible, infectious, recovered), o qual classifica os integrantes de uma população como indivíduo susceptível (S) ao contágio com patógeno, indivíduo infectado (I) pelo patógeno e aqueles recuperados (R) com êxito.

Figura 1- Fluxo epidêmico adotado no modelo SIR para propagação de uma doença infecciosa. Um indivíduo S interage com um indivíduo I e entra no grupo dos infectados. Por fim, o indivíduo do grupo I passa para o grupo R.

Tal modelo é descrito pelo sistema de Equações Diferencias Ordinárias abaixo.²

$$ \frac{dS}{dt}=-\frac{{\beta}SI}{N} \qquad (1)$$

$$ \frac{dI}{dt}=\frac{{\beta}SI}{N}-{\gamma}I \qquad (2)$$

$$ \frac{dR}{dt}={\gamma}I \qquad (3)$$

Nas expressões (1), (2) e (3), $dS/dt$, $dI/dt$ e $dR/dt$ indicam como (com que velocidade) o número de pessoas susceptíveis, infectadas e recuperadas estão mudando no tempo, respectivamente. O sinal de negativo (-) na equação (1) indica que essa quantidade de indivíduos susceptíveis diminuí ao longo do tempo. Portanto, o primeiro termo da equação (2) representa o crescimento do número de infectados. O segundo termo dessa mesma equação está relacionado com o número de indivíduos recuperados, por isso tem um peso negativo sobre os indivíduos infectados. Já a equação (3) descreve o aumento de indivíduos recuperados.

No caso mais simples do modelo SIR (ignorando a demografia da população, nascimentos, mortes e migrações e assumimos uma mistura homogênea), temos apenas as transições S → I e I → R. Tal modelo assume uma taxa de remoção do patógeno ou recuperação do indivíduo infectado $\gamma$, que depende, dentre outras coisas, de como o sistema de saúde está comportando a demanda de pacientes. Para o modelo SIR, usa-se $\frac{1}{{\gamma}}$, que determina o período médio de infecção (ou tempo de incubação do vírus, que no caso do Sars-Cov-2 é algo entorno de 12 dias)³, que é o período médio durante o qual um indivíduo infectado pode transmitir o patógeno. Tal valor é estimado a partir de dados epidemiológicos.

Outro fator que o modelo leva em consideração é a taxa de contato efetiva entre indivíduos infectados (I) com susceptíveis (S), denominada $\beta$. Tal fator está correlacionado com o número de indivíduos (I) e (S) interagindo uns com os outros em um dado intervalor de tempo (por dia, por exemplo). $N$ é o número de indivíduos de uma dada população (N = S + I + R). É fácil perceber que, se essa taxa aumenta, o número de infectados também aumenta. Medidas como a QUARENTENA são extremamente importantes para contenção de epidemias e estão relacionadas com o fator $\beta$ do modelo SIR. Verifica-se na simulação abaixo que, modificando o valor desse fator, muda-se toda a dinâmica da epidemia e, consequentemente, os cenários nos quais teríamos que lidar com a doença.

RESULTADOS E DISCUSSÕES

Com base nos dados fornecidos pelo ministério da saúde¹, até o dia 12/04/2020, o número total de casos da doença no país era de 20.727 infectados, e a evolução desses casos se deu como mostrado no gráfico da Figura (2) abaixo.

Figura 2 - Evolução temporal do número de casos da COVID-19 no Brasil. Os pontos representam os dados do Ministério da Saúde. A curva em azul é um ajuste com uma função SLogistica (sigmoide) para predição de casos futuros.
O ajuste gerou a seguinte curva regida pela Equação (4) abaixo, cujo coeficiente de correlação r² =  0,998.
$$Nº{\quad}de{\quad}infectados=\frac{37583,25}{1+ e^{-0,17(Dia-45,74)}}\qquad(4)$$

Com base na equação (4), pode-se fazer previsões a respeito do número de casos da doença. A tabela (1) abaixo mostra a previsão para os próximos 10 dias com o modelo e com os dados reais a partir da data de hoje, 13/04/2020.¹

Tabela 1 - Dados do Ministério da Saúde sobre os infectados e previsão teórica para o modelo matemático (Equação (4)) e o modelo SIR, em função dos dias desde do primeiro contágio. Os valores em vermelho são os divulgados pelo governo após essa publicação.

################################

Atualização 12/05/2020.

Estudo publicado pela Universidade de São Paulo (USP) mostra que o numero de casos pode ser até 15 vezes maior do que aquele divulgado pelo Ministério da Saúde.

Fonte: https://oglobo.globo.com/sociedade/coronavirus/coronavirus-numero-de-infectados-quinze-vezes-maior-aponta-estudo-24369657. 

Nesse site da USP, pode-se acompanhar o número de casos nacionais e estaduais: https://ciis.fmrp.usp.br/covid19-subnotificacao/

################################

     O modelo SIR representado pelas EDOs (equações (1), (2) e (3)) não pode ser resolvido explicitamente, ou seja, não podemos obter uma expressão analítica exata para a dinâmica de S e I no tempo, portanto, o modelo deve ser resolvido numericamente. Preparou-se uma simulação computacional utilizando a linguagem de programação Python, cujo código será disponibilizado no final desse artigo.

       Fixando os valores de $\frac{1}{{\gamma}}=$ 12 dias e N = 210 milhões de habitantes, verifica-se como a taxa de contato efetiva entre indivíduos infectados (I) com susceptíveis (S) $\beta$ influenciará no número de casos de infectados e no tempo no qual esses irão ocorrer. No gráfico da Figura (2) abaixo pode-se verificar que para $\beta=$1 (que ocorreria caso o país não estivesse aderido a quarentena) teríamos um pico da doença em menos que 30 dias, atingindo mais que metade da população, o que seria catastrófico para qualquer sistema de saúde mundial e acarretaria em milhões de mortes.

Figura 3 - Cenário do contágio da Covid-19 no Brasil sem adesão à quarentena.

Já o gráfico da Figura (4) trata-se do caso em que estamos vivenciando atualmente, no estado de quarentena. O valor de $\beta=$0,290 foi ajustado aos dados do Ministério da Saúde do Brasil e estão de acordo com o número de infectados, que no dia 12/04/2020, 48 dias após o início da pandemia no Brasil, era de 22.169 infectados. Os valores para os próximos 10 dias estão apresentados na tabela (1) acima. O aumento de casos de pessoas infectadas é distribuído num intervalo de tempo maior, quanto menor for o valor de $\beta$, o que torna possível uma organização para enfrentar a pandemia. Observa-se também que, se a quarentena for mantida, o pico dos números de casos será por volta de 100 dias do início da transmissão comunitária (25/02/2020) e ocorrerá no início do mês de junho totalizando um número de casos de superior a 74 milhões de pessoas. Tal modelo, comparado com os dados do Ministério da saúde, parece estar superestimando esse número, entretanto, há cientistas que confirmam essa previsão⁴, alegando que o número de casos confirmados no Brasil é muito maior que o divulgado pelo ministério, principalmente pela falta de testes massivos na população não estarem ocorrendo.

Figura 4 – Cenário atual do contágio da Covid-19 no Brasil com adesão à quarentena.

     Existem outros modelos matemático mais abrangente, que leva em consideração mais fatores a partir dos quais pode-se obter mais precisão nas estimativas. O modelo SEIR é um deles, considerando o número de indivíduos expostos a doença denominado. Tal modelo pode ser simulado online na plataforma da USP "Epcalc", pelo link: https://bit.ly/3aaPMF1. Modelo SIR também pode ser simulado online, nesse link: https://bit.ly/3cwgutB. Todavia, não se consegue salvar os dados da simulação em ambos os sites citados acima. 

           CONCLUSÃO 

Com base nos dados e nos modelos (regressão e SIR) apresentados, pode-se concluir que a matemática é uma importante aliada para descrição de dinâmicas de sistemas como aqueles encontrados na propagação de doenças na sociedade. Com base em modelos matemáticos os agentes públicos podem nortear suas decisões a respeito de como conter o avanço da doença. No Brasil, tais previsões reforçam a necessidade das políticas de quarentena, como ressaltado pelo Ministro da Saúde no programa Fantástico da Rede Globo, exibido no dia 12/04/2020. Com base no ajuste dos dados do Ministério pode-se prever o número de casos no país nos próximos dias. Com o modelo SIR pode-se verificar a importância da quarentena, sem a qual qualquer país do mundo exposto ao vírus colapsaria seu sistema de saúde em poucos dias.

REFERÊNCIAS

1.      Painel Coronavírus - Ministério da Saúde. Acessado: 13/04/2020. Disponível em: https://covid.saude.gov.br/.

2.      M. J. Keeling and P. Rohani, Modeling Infectious Diseases in Humans and Animals, Princeton (2007).

3.      Saude de A a Z. Novo Coronavírus - Ministério da Saúde. Acessado: 13/04/2020. Disponível em: https://www.saude.gov.br/o-ministro/746-saude-de-a-a-z/46490-novo-coronavirus-o-que-e-causas-sintomas-tratamento-e-prevencao-3.

4.      Covid-19 Brasil - estimativa de população infectada.  Acessado: 15/04/2020. Disponível em: https://ciis.fmrp.usp.br/covid19/estimativa-de-populacao-infectada/.

CÓDIGO PYTHON DA SIMULAÇÃO DO MODELO SIR

"""

Created on Fri Apr 10 20:27:17 2020

@author: Wenderson Rodrigues F. da Silva

"""

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

import csv

import pandas as pd

# População total, N.

N = 210000000

time = 70

# Número inicial de indivíduos infectados e recuperados, I0 e R0.

I0, R0 = 1, 0

# Todos os outros, S0, são suscetíveis à infecção inicialmente.

S0 = N - I0 - R0

# Taxa de contato, beta e taxa de recuperação média, gama, (em 1 / dias).

beta, gamma = 0.290, 1./12

# intervalo de pontos no tempo (em dias).

t = np.linspace(0, time, time)

# As equações diferenciais do modelo SIR.

def deriv(y, t, N, beta, gamma):

    S, I, R = y

    dSdt = -beta * S * I / N

    dIdt = beta * S * I / N - gamma * I

    dRdt = gamma * I

    return dSdt, dIdt, dRdt

# Condições iniciais.

y0 = S0, I0, R0

# Integração das equações do modelo SIR ao longo da série de tempo, t.

ret = odeint(deriv, y0, t, args=(N, beta, gamma))

S, I, R = ret.T

###################################################################

t = t.astype(int)

S = S.astype(int)

I = I.astype(int)

R = R.astype(int)

#criando o dataframe

df = pd.DataFrame({'tempo (dias)': t, 'Susceptveis': S, 'Infectados': I, 'Recuperados': R})

# salvando a tabela

df.to_csv('test.csv', sep=',')

###################################################################

# Plotagem os dados em três curvas separadas para S (t), I (t) e R (t).

fig = plt.figure(facecolor='w')

ax = fig.add_subplot(111,axisbelow=True)

ax.plot(t, S/N, 'b', alpha=0.6, lw=2, label='Susceptveis')

ax.plot(t, I/N, 'r', alpha=0.6, lw=2, label='Infectados')

ax.plot(t, R/N, 'g', alpha=0.6, lw=2, label='Recuperado')

ax.set_xlabel('tempo (dias)')

ax.set_ylabel('Fração da população em cada classe')

ax.set_ylim(0,0.000009)

ax.yaxis.set_tick_params(length=0)

ax.xaxis.set_tick_params(length=0)

ax.grid(b=True, which='major', c='w', lw=2, ls='-')

legend = ax.legend()

legend.get_frame().set_alpha(0.5)

for spine in ('top', 'right', 'bottom', 'left'):

    ax.spines[spine].set_visible(False)

plt.savefig('teste.png',dpi=300, format='png')

plt.show()