Fator de compressibilidade de um fluido bidimensional

Compressibility factor of a two-dimensional fluid

Resumos

As simulações computacionais têm assumido um papel cada vez mais importante no cenário de ciência e tecnologia e também têm se tornado uma grande ferramenta no ensino de física. De forma a contribuir com o ensino de técnicas de simulação aplicadas à matéria condensada, este trabalho visa demonstrar o desenvolvimento teórico do cálculo do fator de compressibilidade e seu respectivo resultado numérico através de simulações de dinâmica molecular clássica para um sistema formado por discos rígidos.

Palavras-chave:
Simulação Computacional; Dinâmica Molecular; Discos Rígidos; Teoria cinética; Teorema do Virial

Computer simulations play an important role in science and technology and have become a great teaching tool. In order to contribute to the field of computer simulation teaching applied to condensed matter, this work aims to demonstrate the theoretical development of compressibility factor calculation and its respective numerical result by classical molecular dynamics simulations of a system formed by hard disks.

Keywords
Computer Simulation; Molecular Dynamics; Hard-disks; Kinetic theory; Virial theorem

1. Introdução

As simulações computacionais têm assumido um papel cada vez mais importante no cenário internacional de ciência e tecnologia. O impacto desse tipo de estudo é mais diretamente percebido quando é oferecida uma alternativa viável aos experimentos de laboratório que, por vezes, requerem um dispendioso orçamento ou até mesmo quando as simulações possibilitam inferir certas propriedades de um sistema em condições extremas, muitas vezes inacessíveis experimentalmente.

Do ponto de vista de um estudante de ciências físicas, a literatura disponível para a aprendizagem de modelos e técnicas computacionais é vasta. Particularmente, a modelagem computacional aplicada à física da matéria condensada é bem documentada em livros [¹, ², ³] e artigos [⁴, ⁵]. Também são encontradas iniciativas lúdicas para o aprendizado de simulações físicas de propósito variado [⁶, ⁷]. Por outro lado, ainda que técnicas computacionais, como a dinâmica molecular, tenham sido desenvolvidas a partir das décadas de 1950 e 1960 e alguns trabalhos visando o seu aprendizado tenham sido publicados desde então, observamos, a partir de trabalhos de conclusão de curso em nível de graduação e projetos de iniciação científica, que os aspectos técnicos sobre linguagens de programação e modelos computacionais nem sempre são os únicos desafios a serem superados pelos estudantes. Em geral, a habilidade em lidar com modelos físicos discutidos desde o ensino médio, manipulando-os de forma conveniente com o objetivo de extrair equações adequadas para serem utilizadas em simulações computacionais é algo não-trivial e que merece atenção nesse nível de ensino.

Nesse contexto, este trabalho aborda um assunto tradicionalmente presente nas ementas de disciplinas de física e química tanto no ensino médio quanto no ensino superior em ciências exatas: o gás ideal e o gás real [⁸, ⁹, ¹⁰]. O desvio da idealidade em experimentos com gases reais é tradicionalmente estudado pelo comportamento do fator de compressibilidade Z≡(PV)/(nRT) em função da pressão P. Como a equação de estado do gás ideal é dada por PV = nRT, o sistema idealizado exibe o valor unitário para Z no eixo das ordenadas para qualquer valor de pressão P no eixo das abscissas. Em gases reais, o fator de compressibilidade pode assumir valores menores ou maiores que a unidade e assume Z = 1 apenas na situação em que se comporta como um gás ideal em uma determinada condição. Isto ocorre pois, diferentemente do gás ideal, o gás real é submetido às forças de interação atrativas e/ou repulsivas entre as partículas influindo, portanto, na magnitude da pressão do sistema em relação àquela experimentada por um gás ideal nas mesmas condições, já que interações atrativas/repulsivas tendem a aproximar/afastar as partículas uma das outras, eventualmente modificando a taxa de colisões das moléculas com as paredes da caixa, consecutivamente alterando a pressão, diminuindo-a ou aumentando-a.

O artigo visa demonstrar as principais equações que viabilizam um algoritmo que permite calcular o fator de compressibilidade de um sistema formado por discos rígidos. Essa grandeza é relevante quando o aspecto do tamanho relativo das moléculas torna-se importante e tem influência na equação de estado do sistema, afastando-se da conhecida equação de estado do gás ideal. Simulações de dinâmica molecular clássica são utilizadas com essa finalidade e dessa forma, o leitor é apresentado às técnicas teóricas e numéricas da área. O modelo computacional utiliza a teoria da colisão bidimensional, discutido com certo rigor matemático nos primeiros anos de um curso superior em ciências exatas e, portanto, apenas a compreensão das leis de Newton e um pouco de termodinâmica são suficientes para o acompanhamento das ideias, sem a necessidade de um profundo conhecimento em mecânica estatística, necessário para uma completa e rigorosa descrição teórica do sistema em questão.

2. Modelo e Metodologia

A implementação computacional de sistemas formados por discos rígidos foi fruto do trabalho pioneiro de Alder e Wainwright [¹¹, ¹²] ao longo das décadas de 1950 e 1960, quando a técnica conhecida por Dinâmica Molecular apresentava os primeiros resultados. Embora o principal código computacional tenha sido reescrito em FORTRAN 90 especialmente para este trabalho, é possível encontrar uma versão em linguagem Java no excelente livro introdutório às técnicas de física computacional escrito por Harvey Gould, Jan Tobochnik e Wolfgang Christian [³].

O modelo físico-químico faz alusão a um fluido confinado em uma caixa, modelado por discos rígidos idênticos de massa m e diâmetro σ colidindo elasticamente entre si e com uma caixa quadrada no papel de recipiente. Por julgarmos mais intuitivo para o estudante conceber o sistema confinado através da abordagem explícita de colisões das partículas com a caixa, suprimimos o emprego de condições periódicas de contorno (do inglês, Periodic Boundary Conditions (PBC)) que usualmente são utilizadas nesse tipo de simulação visando minimizar efeitos de superfície indesejados no estudo do interior de materiais. Portanto, a parede da caixa é concebida como um conjunto de quatro anteparos ideais de comprimento L, situação na qual os partículas, ao se chocarem com ela, são rebatidas com o mesmo ângulo de incidência à colisão, respeitando às leis de conservação de momento linear e energia cinética. A dinâmica é descrita em relação às posições ${\mathbf{\text{r}}}_i(t)=x_i(t)\widehat{\mathbf{\text{i}}}+y_i(t)\widehat{\mathbf{\text{j}}}$ e às velocidades ${\mathbf{\text{v}}}_i(t)=v_{i_x}(t)\widehat{\mathbf{\text{i}}}+v_{i_y}\widehat{\mathbf{\text{j}}}(t)$ de cada disco rígido i ao longo do tempo t.

As simulações preveem dois tipos de colisões: partícula-partícula e partícula-parede. A Hamiltoniana do sistema de N discos rígidos é dada por:

(1) $$H=\underset{\text{termo cinético}}{\underbrace{\sum _{i=1}^N\frac{p_i^2}{2m}}}+\underset{\text{partícula-partícula}}{\underbrace{\sum _{i=1}^{N-1}\sum _{j=i+1}^{N}u(r_{ij})}}+\underset{\text{partícula-parede}}{\underbrace{\sum _{i=1}^N\sum _{k=1}^{4}w(r_{ik})}},$$

em que

$$u(r_{ij})=\{\begin{array}{cc}\mathrm{\infty },\hfill & r_{ij}\le \sigma \hfill \\ 0,\hfill & r_{ij}>\sigma \hfill \end{array}$$

e p_i representa a magnitude do momento linear da partícula i. A distância entre os centros de massa das partículas i e j é dada por r_ij = |r_i−r_j|. O terceiro termo da equação (1) refere-se às interações entre um disco e os quatro anteparos:

$$w(r_{ik})=\{\begin{array}{cc}\mathrm{\infty },\hfill & r_{ik}\le \frac{\sigma }{2}\hfill \\ 0,\hfill & r_{ik}>\frac{\sigma }{2}\hfill \end{array}$$

sendo r_ik = |r_i−r_k|, onde r_i é a posição do centro de massa da partícula i e r_k representa a posição do k-ésimo anteparo da parede da caixa de simulação.

2.1. Dinâmica do sistema

Dada a Hamiltoniana escrita na forma da equação (1), – que exprime apenas interações instantâneas entre partícula-partícula e partícula-parede – a dinâmica a ser simulada dependerá apenas dos eventos de colisão entre estes entes. Diferentemente de simulações em que o tempo é discretizado em intervalos Δt regulares e suficientemente pequenos, os sistemas formados por discos rígidos permitem que o lapso temporal da simulação ocorra por meio de intervalos irregulares de tempo, uma vez que a ausência de forças externas entre colisões propicia um movimento retilíneo e uniforme (MRU) das partículas entre sucessivos choques.

Com o objetivo de deixarmos o texto o mais auto-contido quanto possível, de forma a deixá-lo conveniente ao uso pedagógico e instrucional, as próximas seções visam explicitar a dedução das principais equações que possibilitam a escrita de um código computacional descrito pelo fluxograma da Figura 1.

Figura 1
Fluxograma da simulação de discos rígidos.

2.1.1. Determinação do instante de colisão entre duas partículas

Duas partículas i e j irão colidir de acordo com a geometria ilustrada na Figura 2(b). O instante de colisão t_ij entre elas é determinado de tal modo que seus vetores-posição r_i(t) e r_j(t) satisfaçam a condição

(2) $$\left|{\mathbf{\text{r}}}_i(t_{ij})-{\mathbf{\text{r}}}_j(t_{ij})\right|=\sigma .$$

Uma vez que a força de interação entre as partículas atua apenas no instante de contato entre suas superfícies, nenhuma outra força externa é aplicada no decorrer das trajetórias entre duas colisões sucessivas, de forma que o movimento de todos os discos é um MRU entre colisões. Em particular, as equações horárias dos discos que irão colidir entre si até o instante do choque são dadas por:

$${\mathbf{\text{r}}}_i(t_{ij})={\mathbf{\text{r}}}_i(t_o)+{\mathbf{\text{v}}}_i(t_o)t_{ij}$$ (3) $${\mathbf{\text{r}}}_j(t_{ij})={\mathbf{\text{r}}}_j(t_o)+{\mathbf{\text{v}}}_j(t_o)t_{ij},$$

em que v_i(t_o) e v_j(t_o) são os vetores-velocidade das partículas i e j em um instante inicial t_o, respectivamente.

Substituindo as equações (3) na equação (2), obtêm-se

$$\sigma =\left|{\mathbf{\text{r}}}_i(t_{ij})-{\mathbf{\text{r}}}_j(t_{ij})\right|$$ $$=\left|\left({\mathbf{\text{r}}}_i(t_o)-{\mathbf{\text{r}}}_j(t_o)\right)+\left({\mathbf{\text{v}}}_i(t_o)-{\mathbf{\text{v}}}_j(t_o)\right)t_{ij}\right|$$ (4) $$=\left|\mathrm{\Delta }{\mathbf{\text{r}}}_o+\mathrm{\Delta }{\mathbf{\text{v}}}_o{t}_{ij}\right|,$$

em que Δr_o≡r_i(t_o)−r_j(t_o) e Δv_o≡v_i(t_o)−v_j(t_o).

Desenvolvendo a equação (4), chega-se a uma equação de segundo grau em t_ij

(5) $${|\mathrm{\Delta }{\mathbf{\text{v}}}_o|}^2{t}_{ij}^2+\left(2\mathrm{\Delta }{\mathbf{\text{r}}}_o\cdot \mathrm{\Delta }{\mathbf{\text{v}}}_o\right)t_{ij}+{|\mathrm{\Delta }{\mathbf{\text{r}}}_o|}^2-{\sigma }^2=0,$$

cuja solução

$$t_{ij}=-{\displaystyle \frac{\mathrm{\Delta }{\mathbf{\text{r}}}_o\cdot \mathrm{\Delta }{\mathbf{\text{v}}}_o}{{|\mathrm{\Delta }{\mathbf{\text{v}}}_o|}^2}}$$ (6) $$\mathrm{}\pm {\displaystyle \frac{\sqrt{{(\mathrm{\Delta }{\mathbf{\text{r}}}_o\cdot \mathrm{\Delta }{\mathbf{\text{v}}}_o)}^2-{|\mathrm{\Delta }{\mathbf{\text{v}}}_o|}^2\left({|\mathrm{\Delta }{\mathbf{\text{r}}}_o|}^2-{\sigma }^2\right)}}{{|\mathrm{\Delta }{\mathbf{\text{v}}}_o|}^2}}$$

admite, em princípio, até dois valores reais distintos. No contexto do problema, apenas o menor valor positivo possui significado físico, que é o tempo necessário para se observar a condição da equação (2) pela primeira vez.

Do ponto de vista da simulação, é necessário calcular e armazenar os tempos de colisões entre todos os pares de partículas que formam o sistema para que seja identificado o par que colidirá primeiro de acordo com

(7) $$t_p=\text{min}\{t_{12},t_{13},\mathrm{\dots },t_{ij},\mathrm{\dots },t_{N-1N}\},$$

em que t_ij é o instante de tempo em que ocorrerá a colisão entre a i e j-ésima partículas. Certamente, este é o trecho do programa mais crítico e que demanda mais tempo de processamento.

2.1.2. Determinação do instante de colisão entre uma partícula e a caixa

Em relação à origem do sistema de coordenadas mostrado na Figura 2(c), os quatro anteparos que compõem a caixa quadrada são representados pelas retas y = 0 e y = L (anteparos horizontais) e x = 0 e x = L (anteparos verticais). Considerando a posição ${\mathbf{\text{r}}}_i(t_o)=x_i^o\widehat{\mathbf{\text{i}}}+y_i^o\widehat{\mathbf{\text{j}}}$ e velocidade ${\mathbf{\text{v}}}_i(t_o)=v_{i_x}^o\widehat{\mathbf{\text{i}}}+v_{i_y}^o\widehat{\mathbf{\text{j}}}$ do i-ésimo disco em um instante t_o, pode-se determinar o instante de tempo da colisão desse disco com cada um dos quatro anteparos que compõem a caixa:

(8) $$\{\begin{array}{cc}{t}_{w_1}^i=\frac{1}{v_{i_y}^o}\left(\frac{\sigma }{2}-y_i^o\right),\hfill & \text{lado}y=0.\hfill \\ t_{w_2}^i=\frac{1}{v_{i_y}^o}\left(L-y_i^o-\frac{\sigma }{2}\right),\hfill & \text{lado}y=L.\hfill \\ t_{w_3}^i=\frac{1}{v_{i_x}^o}\left(\frac{\sigma }{2}-x_i^o\right),\hfill & \text{lado}x=0.\hfill \\ t_{w_4}^i=\frac{1}{v_{i_x}^o}\left(L-x_i^o-\frac{\sigma }{2}\right),\hfill & \text{lado}x=L.\hfill \end{array}$$

Dessa forma, a determinação do disco que colidirá primeiro com a caixa é feita a partir do menor valor calculado pelas equações (8) para os N discos:

(9) $$t_w=\text{min}\{t_{w_1}^i,t_{w_2}^i,t_{w_3}^i,t_{w_4}^i,\mathrm{\dots },t_{w_1}^N,t_{w_2}^N,t_{w_3}^N,t_{w_4}^N\}.$$

Por fim, a identificação do próximo evento de colisão é determinado pela comparação entre t_p (equação (7): colisão com outra partícula) e t_w (equação (9): colisão com a parede). O menor valor entre os instantes t_p e t_w, obtidos entre todos os pares de partículas e partícula-parede, fornece o lapso temporal da simulação:

(10) $$\mathrm{\Delta }t=\text{min}\{t_w,t_p\}.$$

2.1.3. Cálculo das velocidades após uma colisão

Um disco pode colidir com a caixa ou com outro disco. Após a identificação do tipo de alvo e do instante de tempo do próximo evento de colisão, desloca-se em MRU todos os discos durante um intervalo de tempo Δt. Em seguida, o vetor-velocidade da(s) partícula(s) envolvida(s) de acordo com o tipo de choque é atualizado.

2.1.4. Colisão partícula-partícula

A partir do teorema impulso-variação do momento linear, a grandeza impulso J pode ser escrita como:

(11) $$\mathbf{\text{J}}=\mathrm{\Delta }\mathbf{\text{p}},$$

em que p = m v é o momento linear de uma partícula de massa m e velocidade v.

Denota-se por ${\mathbf{\text{p}}}_i^o=m_i{\mathbf{\text{v}}}_i^o$ o momento inicial do disco i e por ${\mathbf{\text{p}}}_i^f=m_i{\mathbf{\text{v}}}_i^f$ seu momento posterior à uma colisão com outro disco. Deste modo, os vetores-velocidade pós-colisão entre as partículas i e j são escritos como

(12) $${\mathbf{\text{v}}}_i^f={\mathbf{\text{v}}}_i^o+{\displaystyle \frac{1}{m_i}}{\mathbf{\text{J}}}_{ij},$$ $${\mathbf{\text{v}}}_j^f={\mathbf{\text{v}}}_j^o+{\displaystyle \frac{1}{m_j}}{\mathbf{\text{J}}}_{ji},$$

em que J_ij é o impulso aplicado ao disco i devido ao choque com o disco j e consequentemente, pela terceira lei de Newton, J_ji = −J_ij. Contudo, para que as equações (12) sejam úteis do ponto de vista algorítmico na atualização das velocidades das partículas após um choque, deve-se encontrar J_ij em termos das variáveis conhecidas antes da colisão entre os discos, isto é, das posições ${\mathbf{\text{r}}}_i^o$ e ${\mathbf{\text{r}}}_j^o$, das velocidades ${\mathbf{\text{v}}}_i^o$ e ${\mathbf{\text{v}}}_j^o$ e das massas m_i, m_j.

Sem perda de generalidade, o impulso J_ij pode ser decomposto ao longo do eixo “normal” $\widehat{\mathbf{\text{n}}}$ cuja direção é dada pela reta que interliga os centros de massa dos discos e ao longo do eixo “tangencial” $\widehat{\mathbf{\text{t}}}$, perpendicular ao anterior e que é tangente às circunferências dos dois discos no instante de colisão:

(13) $${\mathbf{\text{J}}}_{ij}=({\mathbf{\text{J}}}_{ij}\cdot \widehat{\mathbf{\text{n}}})\widehat{\mathbf{\text{n}}}+({\mathbf{\text{J}}}_{ij}\cdot \widehat{\mathbf{\text{t}}})\widehat{\mathbf{\text{t}}}.$$

Todavia, um modelo simplificado da dinâmica das colisões considera que as superfícies não apresentam rugosidade. Isto significa que, após uma colisão, os discos não adquirem rotação em torno de seu próprio eixo, por não sofrerem um torque promovido por uma eventual força de atrito tangencial ao disco; situação inexistente para discos perfeitamente lisos. Consequentemente, o termo $({\mathbf{\text{J}}}_{ij}\cdot \widehat{\mathbf{\text{t}}})\widehat{\mathbf{\text{t}}}$ da equação (13) é nulo:

(14) $${\mathbf{\text{J}}}_{ij}\cdot \widehat{\mathbf{\text{t}}}={\mathbf{\text{J}}}_{ji}\cdot \widehat{\mathbf{\text{t}}}=0\mathit{}(\text{sem atrito/rugosidade}).$$

Logo, o vetor J_ij deve estar posicionado integralmente ao longo do eixo normal à colisão, conforme ilustrado na Figura 2(b). Além disto, como não há forças externas atuando em um subsistema isolado formado pelos discos i e j, é garantida a conservação do momento linear deste par:

Figura 2
Caixa de simulação e momento do impacto de uma partícula com outra partícula e com a parede. (a) Representação de um instante da simulação para N partículas. As setas representam os respectivos vetores-velocidade. (b) Colisão partícula-partícula. (c) Colisão partícula-parede.

(15) $${\mathbf{\text{p}}}_i^o+{\mathbf{\text{p}}}_j^o={\mathbf{\text{p}}}_i^f+{\mathbf{\text{p}}}_j^f,$$

cujo caráter vetorial pode ser rearranjado em termos escalares em função das componentes normal e tangencial à colisão:

(16) $$\{\begin{array}{cc}{p}_{i_n}^o+p_{j_n}^o=p_{i_n}^f+p_{j_n}^f\hfill & \\ p_{i_t}^o+p_{j_t}^o=p_{i_t}^f+p_{j_t}^f.\hfill \end{array}$$

Como apenas colisões elásticas são consideradas nas simulações, a energia cinética dessse subsistema também é conservada:

(17) $${\displaystyle \frac{p_{i_n}^o{}^2+p_{i_t}^o{}^2}{2m_i}}+{\displaystyle \frac{p_{j_n}^o{}^2+p_{j_t}^o{}^2}{2m_j}}={\displaystyle \frac{p_{i_n}^f{}^2+p_{i_t}^f{}^2}{2m_i}}+{\displaystyle \frac{p_{j_n}^f{}^2+p_{j_t}^f{}^2}{2m_j}}.$$

Neste procedimento, surgem quatro incógnitas $(p_{i_n}^f,p_{i_t}^f,p_{j_n}^f{p}_{j_t}^f)$ que podem ser resolvidas com as equações (14), (16) e (17). Particularmente, a partir da equação (14) e do teorema impulso-variação do momento linear (equação (11)) é deduzido que os momentos lineares na direção tangencial à colisão não se alteram. Com isto, as incógnitas podem ser determinadas através do sistema:

(18) $$\{\begin{array}{c}{p}_{i_n}^o+p_{j_n}^o=p_{i_n}^f+p_{j_n}^f,\hfill \\ \frac{p_{i_n}^o{}^2}{2m_i}+\frac{p_{j_n}^o{}^2}{2m_j}=\frac{p_{i_n}^f{}^2}{2m_i}+\frac{p_{j_n}^f{}^2}{2m_j},\hfill \\ p_{i_t}^o=p_{i_t}^f,\hfill \\ p_{j_t}^o=p_{j_t}^f,\hfill \end{array}$$

que após alguma álgebra, resolve-se para $p_{i_n}^f$ e $p_{j_n}^f$:

(19) $$p_{i_n}^f={\displaystyle \frac{2m_i}{m_i+m_j}}{p}_{j_n}^o+{\displaystyle \frac{m_i-m_j}{m_i+m_j}}{p}_{i_n}^o,$$ $$p_{j_n}^f={\displaystyle \frac{2m_j}{m_i+m_j}}{p}_{i_n}^o+{\displaystyle \frac{m_j-m_i}{m_i+m_j}}{p}_{j_n}^o.$$

A partir deste resultado e verificando que o vetor-impulso possui apenas componente normal à colisão, a equação (13) pode ser reescrita levando em conta apenas os respectivos momentos nesta direção. Uma vez que ${\mathbf{\text{J}}}_{ij}\cdot \widehat{\mathbf{\text{n}}}\equiv J_{i{j}_n}=p_{i_n}^f-p_{i_n}^o$, pode-se resolver para o disco i:

(20) $$J_{i{j}_n}={\displaystyle \frac{2m_i{m}_j}{m_i+m_j}}\mathit{}(v_{j_n}^o-v_{i_n}^o),$$ $${\mathbf{\text{J}}}_{ij}\cdot \widehat{\mathbf{\text{n}}}=-{\displaystyle \frac{2m_i{m}_j}{m_i+m_j}}\mathit{}({\mathbf{\text{v}}}_i^o-{\mathbf{\text{v}}}_j^o)\cdot \widehat{\mathbf{\text{n}}},$$

sendo o vetor unitário $\widehat{\mathbf{\text{n}}}$ definido por

(21) $$\widehat{\mathbf{\text{n}}}=\frac{{\mathbf{\text{r}}}_i(t_{ij})-{\mathbf{\text{r}}}_j(t_{ij})}{|{\mathbf{\text{r}}}_i(t_{ij})-{\mathbf{\text{r}}}_j(t_{ij})|}\equiv \frac{\mathrm{\Delta }{\mathbf{\text{r}}}_f}{\sigma },$$

que pode ser obtido pelas equações (3) do MRU no instante de colisão calculado pela equação (6). Neste trabalho, as massas dos discos são consideradas idênticas, logo m_i = m_j = m. Assim, o vetor-impulso escrito em sua forma final é dado por

(22) $${\mathbf{\text{J}}}_{ij}=-m\left(\mathrm{\Delta }{\mathbf{\text{V}}}_o\cdot \mathrm{\Delta }{\mathbf{\text{r}}}_f\right)\frac{\mathrm{\Delta }{\mathbf{\text{r}}}_f}{{\sigma }^2}.$$

Apesar da equação (22) ter sido obtida a partir dos eixos normal e tangencial à colisão, ela também é aplicável ao sistema canônico quando os vetores-posição e velocidade são escritos em termos das coordenadas retangulares x e y, como é o caso das simulações realizadas. Em vista disso, os vetores-velocidade dos discos i e j após colidirem um com o outro podem ser atualizados empregando-se as equações (12) e (22).

2.1.5. Colisão partícula-parede

A caixa pode ser considerada um alvo de massa infinitamente maior que a de um disco, de forma que o momento linear de uma partícula i logo após a colisão com a parede j pode ser obtido analisando o limite m_j→∞ na equação (19). Pode-se concluir, portanto, que $p_{i_n}^f\approx -p_{i_n}^o$, isto é, a colisão de uma partícula com uma parede plana resulta simplesmente na modificação do sinal do momento da componente normal à ela e consequentemente modificará apenas a componente da velocidade perpendicular à parede da caixa onde ocorre o choque.

3. Simulações

O esquema geral da simulação está ilustrado no fluxograma da Figura 1. Os parâmetros previamente definidos para a simulação são: o número de discos rígidos (N) de mesma massa, o diâmetro (σ), o tamanho do lado da caixa de simulação (L) e a quantidade máxima de colisões permitidas.

3.1. Posições iniciais

O conjunto das posições iniciais das partículas $\{{\mathbf{\text{r}}}_i^o\}$ poderia ser obtido de forma aleatória no interior da caixa de simulação tomando-se o cuidado para que os discos não se sobreponham uns nos outros, sendo uma abordagem satisfatória para o estudo de gases mas inadequada para o estudo de cristais, por exemplo. Por outro lado, essa abordagem introduziria uma dificuldade adicional no estudo de sistemas suficientemente densos, onde o espaço entre discos é escasso e um arranjo possível dificilmente seria obtido de maneira aleatória. Para contornar essa dificuldade, a melhor prática em simulações de dinâmica molecular é iniciar o sistema em uma configuração tipo rede cristalina perfeita e aguardar um certo número de passos até que o sistema atinja o equilíbrio.

3.2. Velocidades iniciais

O conjunto de velocidades iniciais dos discos $\{{\mathbf{\text{v}}}_i^o\}$ também poderia ser escolhido de forma aleatória semelhantemente às posições iniciais. Após algum tempo de evolução do sistema, a distribuição de velocidades atinge o equilíbrio e permite que medidas de grandezas físicas representativas deste estado possam ser calculadas. Entretanto, um método mais eficiente para ajustar as velocidades iniciais é empregar o algoritmo de Box-Muller [¹⁹] que amostra as velocidades iniciais de forma a satisfazer a distribuição de Maxwell-Boltzmann em uma dada temperatura T, ou equivalentente, em uma dada velocidade mais provável v_P da distribuição de velocidades, já no ponto de partida da simulação [¹⁴].

4. Propriedades de Interesse

Pretende-se avaliar o comportamento do sistema sobretudo em relação à fração de ocupação ϕ bem como ante a abordagem de colisões explícitas das partículas com as paredes da caixa que introduz uma situação diferente em relação às simulações tradicionais que empregam PBC.

A grandeza mecânica diretamente relacionada às colisões com a caixa é a pressão. Outra grandeza associada à pressão mas que também está ligada à termodinâmica é o fator de compressibilidade. Ambas propriedades são discutidas nas próximas seções.

4.1. Fração de ocupação

Dado o conjunto de parâmetros iniciais da simulação, a densidade

(23) $$\rho =\frac{N}{A}$$

é calculada pela razão entre o número de partículas e a área A = L ² da caixa. Porém, no caso de discos com tamanho não-desprezível, a grandeza denominada “fração de ocupação” é mais relevante para caracterizar a compacidade do sistema:

(24) $$\varphi =\frac{A_{\sigma }}{A}⟹\varphi =\frac{N\pi }{4}{\left(\frac{\sigma }{L}\right)}^2=\frac{\pi \rho {\sigma }^2}{4},$$

e mensura a área A_σ ocupada por todos os discos em relação à área da caixa. Para um sistema de esferas rígidas, por exemplo – análogo tridimensional dos discos rígidos – o valor ϕ_3D≈0,50 define uma transição de fase de primeira ordem líquido-sólido [¹⁵, ¹⁶] enquanto que para o caso bidimensional, o valor encontrado na literatura é ϕ_2D≈0,70 [¹⁷, ¹⁸].

Como apenas a fase fluida é de interesse neste trabalho, o intervalo escolhido para a fração de ocupação é 0 < ϕ < 0,6.

4.2. Pressão

A “pressão” ($\tilde{P}$), interpretada aqui como força por unidade de comprimento e com origem puramente mecânica, pode ser escrita para o sistema em questão como sendo

(25) $$\tilde{P}=\frac{|{\mathbf{\text{F}}}_{\text{total}}^{\text{ext}}|}{4L},$$

em que 4L é o perímetro da caixa e $|{\mathbf{\text{F}}}_{\text{total}}^{\text{ext}}|$ é o módulo da força externa total que o sistema de partículas está submetido devido à caixa que, por sua vez, é constantemente colidida pelos discos. Devido ao confinamento, o sistema de partículas possui velocidade do centro de massa nula e, portanto, a soma vetorial das forças externas oriundas dos lados da caixa será zero, uma vez que a própria caixa não se move. Por isso, consideramos uma parede de comprimento igual ao perímetro da caixa e somamos o módulo da força que cada um dos quatro lados exerce nos discos:

(26) $$|{\mathbf{\text{F}}}_{\text{total}}^{\text{ext}}|=|{\mathbf{\text{F}}}_{w_1}|+|{\mathbf{\text{F}}}_{w_2}|+|{\mathbf{\text{F}}}_{w_3}|+|{\mathbf{\text{F}}}_{w_4}|.$$

Evocando-se a segunda lei de Newton, na sua forma da variação do momento linear em um intervalo de tempo, e levando-se em conta que as colisões dos discos com a caixa são instantâneas, a pressão média neste intervalo advém da força média associada às colisões ocorridas com a caixa em um intervalo de tempo τ:

(27) $$\tilde{P}=\frac{1}{4L\tau }\sum _{i=1}^N|\mathrm{\Delta }{\mathbf{\text{p}}}_i|,$$

em que as variações dos momentos Δp_i são aquelas relacionadas aos choques dos discos com a caixa no intervalo de tempo τ.

4.3. Fator de compressibilidade

O teorema do Virial [²⁰] é bastante útil na análise de sistemas de partículas, pois permite estabelecer uma relação direta entre a energia cinética e a energia potencial do sistema. Em particular, simulações de Dinâmica Molecular e de Monte Carlo utilizam o resultado deste teorema para derivar um modo de cálculo da equação de estado a partir do estudo da grandeza G, definida como

(28) $$G\equiv \sum _{i=1}^N{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{p}}}_i,$$

onde r_i e p_i representam o vetor-posição e o vetor-momento linear da partícula i e o somatório corre para todas as N partículas do sistema.

A variação temporal instantânea da equação (28) é dada por

(29) $$\frac{dG}{dt}=\sum _{i=1}^N\frac{d{\mathbf{\text{r}}}_i}{dt}\cdot {\mathbf{\text{p}}}_i+\sum _{i=1}^N{\mathbf{\text{r}}}_i\cdot \frac{d{\mathbf{\text{p}}}_i}{dt},$$

que pode ser analisada por meio de sua média temporal em um intervalo de tempo τ:

(30) $$\frac{1}{\tau }{\int }_0^{\tau }\frac{dG}{dt}\mathit{d}t=\frac{1}{\tau }\left(G(\tau )-G(0)\right).$$

Em um sistema confinado em uma caixa, como é o caso do modelo proposto, as coordenadas x e y das partículas estão restritas ao intervalo [0,L]. No equilíbrio, as velocidades são limitadas pela distribuição de Maxwell-Boltzmann e, consequentemente, a grandeza G tem uma cota superior. Portanto, para intervalos τ suficientemente longos, a média temporal de $\frac{dG}{dt}$, dada pela equação (30), converge a zero. A partir deste resultado e em conjunto com as relações dinâmicas convencionais

(31) $${\mathbf{\text{v}}}_i={\displaystyle \frac{d{\mathbf{\text{r}}}_i}{dt}},$$ $${\mathbf{\text{p}}}_i=m_i{\mathbf{\text{v}}}_i,$$ $${\mathbf{\text{F}}}_i={\displaystyle \frac{d{\mathbf{\text{p}}}_i}{dt}},$$

faz-se a média temporal da equação (29):

(32) $${⟨{\displaystyle \frac{dG}{dt}}⟩}_{\tau }={⟨{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{d{\mathbf{\text{r}}}_i}{dt}}\cdot {\mathbf{\text{p}}}_i+{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\displaystyle \frac{d{\mathbf{\text{p}}}_i}{dt}}⟩}_{\tau }$$ $$0={⟨{\displaystyle \sum _{i=1}^N}{\mathbf{\text{v}}}_i\cdot m_i{\mathbf{\text{v}}}_i+{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i⟩}_{\tau }$$ $$0=2{<K>}_{\tau }+{⟨{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i⟩}_{\tau },$$

sendo K a energia cinética do sistema de partículas, F_i a força resultante aplicada à i-ésima partícula e o símbolo “ < _{> τ}” denotando a média temporal em um intervalo de tempo τ→∞.

Apesar do sistema de N discos rígidos trocar constantemente momento linear nas colisões, os choques elásticos fazem com que a energia cinética total não varie, uma vez que esse tipo de energia não é armezenada em nenhum tipo de energia potencial, logo _⟨K⟩τ = K. Além disto, o sistema apresenta 2N graus de liberdade e, de acordo com o teorema de equipartição de energia, a energia cinética do sistema é associada à temperatura T na forma:

(33) $$K=2N\left(\frac{1}{2}{k}_{B}T\right)⟹2K=2N{k}_{B}T,$$

cujo resultado pode ser substituído nas equações (32).

O termo relativo à força sob a partícula i pode ser decomposto em mais dois termos: um deles refere-se às forças externas que são aplicadas pelos anteparos da caixa às partículas, e o outro às forças internas aplicadas pelas demais partículas sobre a partícula i levando as equações (32) a

(34) $$0=2N{k}_{B}T+{⟨{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i^{\text{ext}}⟩}_{\tau }+{⟨{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i^{\text{int}}⟩}_{\tau },$$

em que o segundo termo dessa equação pode ser investigado no instante de colisão dos discos com a caixa, estando relacionado à pressão $\tilde{P}$:

(35) $${\mathbf{\text{F}}}_i^{\text{ext}}=\{\begin{array}{cc}\tilde{P}L\widehat{\mathbf{\text{j}}},\hfill & \text{lado}y=0.\hfill \\ -\tilde{P}L\widehat{\mathbf{\text{j}}},\hfill & \text{lado}y=L.\hfill \\ \tilde{P}L\widehat{\mathbf{\text{i}}},\hfill & \text{lado}x=0.\hfill \\ -\tilde{P}L\widehat{\mathbf{\text{i}}},\hfill & \text{lado}x=L.\hfill \end{array}$$

No instante de impacto com a caixa, os discos podem assumir os seguintes vetores-posição:

(36) $${\mathbf{\text{r}}}_i=\{\begin{array}{cc}{x}_i\widehat{\mathbf{\text{i}}}+\left(\frac{\sigma }{2}\right)\widehat{\mathbf{\text{j}}},\hfill & \text{lado}y=0.\hfill \\ x_i\widehat{\mathbf{\text{i}}}+\left(L-\frac{\sigma }{2}\right)\widehat{\mathbf{\text{j}}},\hfill & \text{lado}y=L.\hfill \\ \left(\frac{\sigma }{2}\right)\widehat{\mathbf{\text{i}}}+y_i\widehat{\mathbf{\text{j}}},\hfill & \text{lado}x=0.\hfill \\ \left(L-\frac{\sigma }{2}\right)\widehat{\mathbf{\text{i}}}+y_i\widehat{\mathbf{\text{j}}},\hfill & \text{lado}x=L.\hfill \end{array}$$

Com isso, o produto escalar no segundo termo da equação (34) pode ser realizado a partir dos vetores da equação (35) e equação (36):

(37) $${⟨{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i^{\text{ext}}⟩}_{\tau }=2\times \left[\tilde{P}L{\displaystyle \frac{\sigma }{2}}-\tilde{P}L\left(L-{\displaystyle \frac{\sigma }{2}}\right)\right]$$ $$=-2\tilde{P}{L}^2\stackrel{\equiv \eta }{\overbrace{\left(1-{\displaystyle \frac{\sigma }{L}}\right)}}$$ $$=-2\tilde{P}A\eta ,$$

sendo A = L ² a área da caixa e η um fator de correção devido ao tamanho relativo do diâmetro dos discos em relação ao lado L.

Associando as equações (34) e (37), pode-se exprimir a equação de estado desse sistema bidimensional de partículas interagentes na forma

(38) $$\tilde{P}=\frac{1}{\eta A}\left(N{k}_{B}T+\frac{1}{2}{⟨\sum _{i=1}^N{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i^{\text{int}}⟩}_{\tau }\right).$$

Quando σ≪L, o fator de correção η tende à unidade e aproxima-se da situação em que partículas infinitamente diminutas estão no interior de uma caixa macroscópica. Além disto, como já era esperado, se as partículas não interagem, isto é, se não há forças internas no sistema, a equação de estado de um gás ideal é recuperada, ou seja, $\tilde{P}A=N{k}_{B}T$, que é a equação análoga a de Clapeyron na forma bidimensional. Portanto, o surgimento de interações entre as partículas favorece o desvio da equação de estado do gás ideal e seu respectivo impacto pode ser analisado pelo comportamento do fator de compressibilidade:

(39) $$Z\equiv \frac{\tilde{P}A}{N{k}_{B}T}=B_1+B_2\varphi +B_3{\varphi }^2+B_4{\varphi }^3+\mathrm{\cdots }$$

O fator Z é estudado em função da fração de ocupação ϕ (equação (24)) em uma caixa de área constante A e em uma temperatura T previamente fixada conforme discutido na seção (^3.2). O conjunto {B_i} é conhecido como conjunto dos coeficientes do Virial, sendo que B₁ = 1. Se B_{i > 1} = 0, a equação (39) se torna a equação de estado do gás ideal com Z = 1. Particularmente para discos rígidos, alguns elementos do conjunto {B_i} são determinados teoricamente a partir de cálculos sofisticados utilizando teoria de grafos e integração Monte Carlo [²¹].

Como a Hamiltoniana do sistema de discos rígidos admite apenas forças repulsivas de ação instantânea entre duas partículas, é conveniente utilizar a segunda lei de Newton na forma de momento linear e produzir a média do termo que representa as interações entre partículas na equação (38), a partir da variação do momento linear na colisão entre uma partícula i e outra j:

(40) $${⟨\sum _{i=1}^N{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i^{\text{int}}⟩}_{\tau }=\frac{1}{\tau }\sum _{t=0}^{\tau }\sum _{i=1}^N{\mathbf{\text{r}}}_i(t)\cdot \mathrm{\Delta }{\mathbf{\text{p}}}_i(t).$$

Assumindo que há apenas um evento de colisão a cada instante t e admitindo-se a conservação do momento linear, conclui-se que Δp_i(t) = −Δp_j(t). Decorre disto, que a equação (40) pode ser reescrita como

(41) $${⟨{\displaystyle \sum _{i=1}^N}{\mathbf{\text{r}}}_i\cdot {\mathbf{\text{F}}}_i^{\text{int}}⟩}_{\tau }={\displaystyle \frac{1}{\tau }}{\displaystyle \sum _{t=0}^{\tau }}\left({\mathbf{\text{r}}}_i(t)-{\mathbf{\text{r}}}_j(t)\right)\cdot \mathrm{\Delta }{\mathbf{\text{p}}}_i(t).$$

Sendo assim, o fator de compressibilidade para um sistema de discos rígidos confinados em uma caixa é dado por

(42) $$Z=\frac{1}{\eta }\left[1+\frac{1}{2K}\frac{1}{\tau }\sum _{t=0}^{\tau }{\mathbf{\text{r}}}_{ij}\cdot \mathrm{\Delta }{\mathbf{\text{p}}}_i\right],$$

em que $\eta =(1-\frac{\sigma }{L})$ é o fator de correção devido ao tamanho relativo dos discos em relação à caixa, K é a energia cinética do sistema, r_ij = r_i−r_j é a diferença das posições dos centros de massas das partículas i e j no instante da colisão entre ambas e Δp_i é a respectiva variação do momento linear da partícula i nessa colisão.

4.4. Incertezas das medidas

A estratégia de cálculo da incerteza do valor médio da pressão (equação (27)) e do fator de compressibilidade (equação (42)) consiste em se analisar a estatística de blocos de configurações. Nesta abordagem, um bloco é formado por um conjunto de configurações {r_i,v_i} geradas a partir de n_c colisões sucessivas de onde se extrai a média da grandeza de interesse denominando-a “medida” da grandeza x. Ao longo de uma simulação de n_max colisões, há n_b blocos tal que n_max = n_b×n_c. Ao término da simulação, o erro padrão das n_b medidas da grandeza x é avaliado por

(43) $$\mathrm{\Delta }x=\sqrt{\frac{{<x^2>}_{n_b}-{<x>}_{n_b}^2}{n_b}}$$

em que “ < >_nb” representa a média nos n_b blocos. Portanto, o valor médio da grandeza em questão junto à sua incerteza é expresso na forma

(44) $$x={<x>}_{n_b}\pm \mathrm{\Delta }x.$$

5. Resultados

Em todas as simulações a aresta da caixa foi fixada em L = 100 em unidades arbitrárias [u. a.]. As posições iniciais do centro de massa de cada disco foram iniciadas em uma rede quadrada com determinado número de partículas em cada simulação: N = k ², em que k ∈ ℕ com 2≤k≤80. As velocidades iniciais foram amostradas pelo algoritmo de Box-Muller com v_p = 10 [u. a.] e que é equivalente a todas as simulações estarem virtualmente na mesma temperatura. Em geral, um grande número de colisões na etapa de equilibração e de tomada das medidas é necessário, tipicamente da ordem de 10⁴ a 10⁶ colisões, dependendo do tamanho do sistema.

A seguir, os resultados para a pressão e fator de compressibilidade são apresentados e discutidos.

5.1. Pressão

Com o objetivo de testar o código computacional e estabelecer uma referência de comparação para simulações com σ≠0, avaliou-se o comportamento do sistema quando os diâmetros σ são nulos. O resultado é exibido na Figura 3 em que se observa um excelente ajuste linear entre a pressão $\tilde{P}$ e o número de partículas N. Ainda que longe do limite termodinâmico (no sentido de N→∞) e na ausência de PBC, as simulações indicam o comportamento esperado de um gás ideal, uma vez que $\tilde{P}A=N{k}_{B}T$ implica que $\tilde{P}\propto N$ com A e T fixados.

Figura 3
Pressão em função do número de discos, mantidas a área (L ²) da caixa e temperatura constantes. O ajuste linear é realizado para simulações que desprezam colisões entre discos (“gás ideal”) e um polinômio do terceiro grau ajusta a curva para simulações de discos com diâmetros de tamanho não-desprezível (σ/L = 0,01). As barras de erros são menores que o tamanho dos símbolos utilizados e por isso não são perceptíveis na escala do gráfico.

Todavia, quando a dimensão σ dos discos é uma fração da aresta da caixa, o comportamento da pressão deixa de ser uma função linear do número de partículas, mesmo quando a razão σ/L é relativamente pequena. Como ainda se observa na Figura 3, a característica não-linear é mais pronunciada em sistemas suficientemente densos, relacionados a uma fração de ocupação ϕ significativa. Essa fração desempenha um papel importante no resultado pois, quando σ/L é mantida fixa em 0, 01, por exemplo, ϕ cresce tão somente em função do número N de discos. Nesta condição, quando são utilizados 900 discos a fração ϕ é de aproximadamente 7 %, enquanto que para N = 6400, a fração de ocupação é da ordem de 50 %. Portanto, a fração de ocupação causa o efeito de área “excluída” promovendo uma diminuição virtual da área total da caixa propiciando a elevação da pressão decorrente de uma maior taxa de colisões dos discos com as paredes em relação a um gás ideal nas mesmas condições de densidade ρ e temperatura T.

5.2. Fator de compressibilidade

A partir do perfil da pressão do sistema em função do número de discos, verificou-se que o efeito da área “excluída” decorrente das interações repulsivas entre as partículas contribui ao desvio do comportamento esperado da equação de estado do gás ideal. Outra forma de investigar a desvio da idealidade é por meio do fator de compressibilidade, cujo resultado é mostrado na Figura 4. Nesse resultado, nota-se um excelente acordo entre os dados obtidos por nossas simulações e os valores teóricos derivados por E. J. Janse van Rensburg [²¹] que consideram os oito primeiros coeficientes do Virial da equação de estado. Como esperado, o fator de compressibilidade Z converge à unidade para sistemas suficientemente rarefeitos, como é o caso do gás ideal, mas à medida que o sistema de discos rígidos se torna mais denso no interior da caixa, o efeito relacionado ao diâmetro se torna cada vez mais importante, como pode ser verificado pelo crescimento do valor de Z enquanto ϕ aumenta.

Figura 4
Fator de compressibilidade $Z=\frac{\tilde{P}A}{N{k}_{B}T}$ em função da fração de ocupação ϕ. A curva teórica considera os primeiros oito coeficientes do Virial descritos na literatura [²¹]. As setas indicam desvios sistemáticos da curva, sugerindo que efeitos relacionados ao tamanho do sistema estejam presentes nas simulações. As barras de erros são menores que o tamanho dos símbolos utilizados e por isso não são perceptíveis na escala do gráfico.

Contudo, também é observado na Figura 4um pequeno desacordo entre os valores de Z obtidos por simulação e pela curva teórica. Particularmente, quando a razão σ/L cresce de 0, 01 para 0, 05, há uma tendência do fator de compressibilidade ser maior, cujo efeito é mais perceptível para os valores de Z indicados pelas setas em ϕ aproximadamente igual a 0, 28, 0, 39 e 0, 50. Há dois aspectos que estão relacionados a esse pequeno desacordo entre as curvas. Primeiro, é preciso salientar que o modelo computacional mais adequado para emular um sistema microscópico “infinito” – porém implicitamente confinado macroscopicamente – utiliza PBC e considera o fator η = 1 na equação (37) uma vez que, neste limite, os diâmetros dos discos são infinitamente menores que as dimensões da caixa. No entanto, o sistema simulado é explicitamente confinado e o fator de correção η = 1−σ/L abarca o denominador da equação (42) fazendo, portanto, que razões σ/L maiores contribuam ao incremento de Z. O segundo aspecto está relacionado a presença de um número relativamente pequeno de discos nas simulações com o pretenso objetivo em representar o limite termodinâmico, situação em que o número de partículas é da ordem do número de Avogadro sendo, portanto, intratável do ponto de vista computacional. A propósito, os valores de Z destacados na Figura 4 para ϕ≈0,50 foram obtidos com N = 256 e σ/L = 0,05 e também com N = 6400 discos e σ/L = 0,01 já que, da própria definição da equação (24), dois valores distintos de σ/L estão associados a outros dois valores de N para uma mesma fração ϕ.

As três regiões destacadas por setas na Figura 4 são analisadas em maiores detalhes na Figura 5. Nesta figura, nota-se uma forte correlação entre Z e o número de partículas N nas três frações de ocupação mostradas. Quanto maior é a quantidade de discos nas simulações mais sugestiva é a convergência a um determinado valor de Z – por extrapolação do ajuste linear dos dados quando $1/\sqrt{N}\to 0$. A convergência é provável não apenas devido ao aumento de N mas também pela simultânea diminuição da razão σ/L em uma dada fração de ocupação ϕ, indicando que tal procedimento pode ser usado como uma boa estratégia para estimar o fator de compressibilidade no limite termodinâmico N→∞.

Figura 5
Fator de compreessibilidade Z em função do número de discos para três valores da fração de ocupação ϕ. A quantidade de partículas utilizada em cada simulação está ao lado de cada ponto. A extrapolação do ajuste linear estima Z no limite termodinâmico em que $1/\sqrt{N}\to 0$ quando N→∞.

Por fim, esse procedimento foi aplicado a todas as frações de ocupação simuladas e o resultado é mostrado na Tabela 1 em que, dentro das incertezas do ajuste, os valores de Z obtidos por extrapolação são indistintos dos valores teóricos preditos para um fluido de discos rígidos no limite termodinâmico.

Tabela 1
Valores de Z no limite N→∞.

ϕ	Z N→∞	$Z_{\text{teórico}}^{\mathrm{a}}$	desvio (%)b
0,0079	1,0160 ± 0,0001	1,0159	0,01
0,0314	1,0660 ± 0,0003	1,0661	0,01
0,0707	1,1583 ± 0,0008	1,1587	0,03
0,1257	1,3146 ± 0,0031	1,3107	0,30
0,1963	1,5561 ± 0,0008	1,5558	0,02
0,2827	1,9632 ± 0,0036	1,9624	0,04
0,3848	2,6866 ± 0,0039	2,6807	0,22
0,5027	4,1441 ± 0,0164	4,0674	1,89

^aExpansão de Z até o oitavo coeficiente do virial, cujos valores B₁ = 1, B₂ = 2, B₃ = 3,12801775, B₄ = 4,25785446, B₅ = 5,336897, B₆ = 6,3626, B₇ = 7,351 e B₈ = 8,338 na equação (39) foram obtidos por E. J. Janse van Rensburg [²¹] no limite termodinâmico.^bDesvio percentual em relação ao valor teórico, calculado por |ZN→∞-Zteórico |Zteórico ×100%.

6. Considerações finais

Os fundamentos, principais equações e procedimentos que viabilizam a simulação computacional de um sistema formado por discos rígidos foram discutidos. Acreditamos que o texto pode servir como guia para uma discussão a posteriori de um algoritmo que pode ser trabalhado em várias disciplinas de graduação em ciências físicas instigando os estudantes na escrita do próprio código computacional na linguagem de programação preferida.

No regime simulado, o procedimento de cálculo do fator de compressibilidade para o sistema de discos rígidos – através da extrapolação do ajuste linear dos valores obtidos pelas simulações – mostrou-se em excelente acordo com os valores teóricos preditos por técnicas independentes no limite termodinâmico. Porém, devido aos discos não sofrerem qualquer sobreposição com os demais, há um limite máximo possível de ocupação da caixa de simulação e transições de fase em frações ϕ superiores as simuladas neste trabalho podem ser observadas. Diante disso, novas simulações e discussões podem ser aprofundadas nesse sentido, em que limitações do tratamento explícito das colisões dos discos com a caixa podem dificultar a descrição apropriada do bulk de um líquido ou sólido tendo maior ou menor impacto em certas propriedades do sistema.

Não obstante, resultados interessantes sobre a equação de estado para sistemas bidimensionais formados por discos rígidos comparadas à famosa equação de Clapeyron dos gases ideais puderam ser obtidos a partir da abordagem computacional apresentada, tornando-se uma excelente oportunidade para professores e estudantes discutí-los à luz de bases físicas elementares.

Agradecimentos

Os autores agradecem ao Grupo de Estudos em Modelagem Computacional e Aplicações (GEMCA) e a equipe do Núcleo de Tecnologia da Informação (NTI) do Instituto Federal de Educação, Ciência e Tecnologia do Sul de Minas Gerais – Campus Inconfidentes, pelo suporte computacional às simulações realizadas neste trabalho.

Referênces

• ^[1]
  M.P. Allen e D.J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1996).
• ^[2]
  M.E.J. Newman e G.T. Barkema, Monte Carlo Methods in Statistical Physics (Clarendon Press, Oxford, 1999).
• ^[3]
  H. Gold, J. Tobochnik, W. Christian, An Introduction to Computer Simulation Methods: Applications to Physical Systems (Addison-Wesley, Boston, 2006), 3^a ed.
• ^[4]
  M.A. dos Reis e S.A. Vitiello, Rev. Bras. Ens. Fis. 28 , 45 (2006).
• ^[5]
  L. Madeira e S.A. Vitiello, Rev. Bras. Ens. Fis. 34 , 4303 (2012).
• ^[6]
  OPS (The Open Source Physics Project), disponível em https://www.compadre.org/osp/, acessado em 18/06/2021.
  » https://www.compadre.org/osp/
• ^[7]
  PhET (Physics Education Technology Project), disponível em https://phet.colorado.edu/, acessado em 18/06/2021.
  » https://phet.colorado.edu/
• ^[8]
  H.M. Nussenzveig, Curso de Física Básica 1 (Blucher, São Paulo, 2014).
• ^[9]
  H.M. Nussenzveig, Curso de Física Básica 2 (Blucher, São Paulo, 2014).
• ^[10]
  G. Castellan, Fundamentos de Físico-Química (LTC, Rio de Janeiro, 2014).
• ^[11]
  B.J. Alder e T.E. Wainwright, J. Chem. Phys. 27 , 1208 (1957).
• ^[12]
  B.J. Alder, W.G. Hoover e D.A Young, J. Chem. Phys. 49 , 3688 (1968).
• ^[13]
  A. Stukowski, Modeling Simul Mater. Sci. Eng. 18 , 015012 (2010).
• ^[14]
  M.A. dos Reis e G.H. Batista, A Física na Escola 18 , 41 (2020).
• ^[15]
  M.D. Rintoul e S. Torquato, Physical Review E 58 , 532 (1998).
• ^[16]
  M.N. Bannerman, L. Lue e L.V. Woodcock, The Journal of Chemical Physics 132 , 084507 (2010).
• ^[17]
  M.A.A. Silva, A. Caliri e B.J. Mokross, Physical Review Letters 58 , 2312 (1987).
• ^[18]
  A. Huerta, D. Henderson e A. Trokhymchuk, Physical Review E 74 , 061106 (2006).
• ^[19]
  G.E.P. Box e M.E. Muller, Ann. math. Stat. 29 , 610 (1958).
• ^[20]
  J.L. Safko, H. Goldstein e C.P. Poole, Classical Mechanics (Pearson, Edinburgh, 2014), 3ª ed.
• ^[21]
  E.J. Janse van Rensburg, J. Phys. A: Math. Gen. 26 , 4805 (1993).

Datas de Publicação

• Publicação nesta coleção
  22 Nov 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  18 Jun 2021
• Revisado
  03 Set 2021
• Aceito
  21 Out 2021