MBA: Gás simples

Anterior: Por que usar e o que são modelos baseados em indivíduos | Índice: Ecologia | Próximo: MBA: Caminhante aleatório

Colisão entre partículas

Considerando um modelo simples de gás ideal, não há força atuando sob as partículas, então a interação que ocorre entre as partículas se dá apenas por meio de colisões. Assim é necessário calcular a variação na velocidade de cada partícula após a colisão. Começando em uma dimensão, precisamos garantir a conservação do momento: $$m_1{v}_{1^{\prime }}+m_2{v}_{2^{\prime }}=m_1{v}_{1^{″}}+m_2{v}_{2^{″}}$$E da energia cinética: $$\frac{1}{2}{m}_1{v}_1{\text{'}}^2+\frac{1}{2}{m}_2{v}_2{\text{'}}^2=\frac{1}{2}{m}_1{v}_{1^{\prime }}{\text{'}}^2+\frac{1}{2}{m}_2{v}_{2^{\prime }}{\text{'}}^2$$Colocando o referencial de forma que ${v_2}^{\prime }=0$, então as velocidades no novo referencial podem ser escritas como $u_i=v_i-{v_2}^{\prime }$, de forma que da conservação do momento ficamos com: $$u_{1^{″}}=u_{1^{\prime }}-\frac{m_2}{m_1}{u}_{2^{″}},(1)$$ Elevando ao quadrado: $$u_{1^{\prime }}{\text{'}}^2=u_1{\text{'}}^2-{\left(\frac{m_2}{m_1}\right)}^2{u}_{2^{\prime }}{\text{'}}^2-2\frac{m_2}{m_1}{u}_{1^{\prime }}{u}_{2^{″}}$$ Substituindo na equação de conservação de energia, podemos encontrar ${u_2}^{″}$: $$\begin{array}{c}{m}_1{u}_1{\text{'}}^2=m_1{u}_{1^{\prime }}{\text{'}}^2+m_2{u}_{2^{\prime }}{\text{'}}^2\\ m_1{u}_1{\text{'}}^2=m_1[u_{1^{\prime }}{}^2+{\left(\frac{m_2}{m_1}\right)}^2{u}_{2^{\prime }}{\text{'}}^2-2\frac{m_2}{m_1}{u}_{1^{\prime }}{u}_{2^{″}}]+m_2{u}_{2^{\prime }}{\text{'}}^2\\ m_1{u}_1{\text{'}}^2=m_1{u}_{1^{\prime }}{}^2+\frac{m_2^2}{m_1}{u}_{2^{\prime }}{\text{'}}^2-2m_2{u}_{1^{\prime }}{u}_{2^{″}}+m_2{u}_{2^{\prime }}{\text{'}}^2\\ 0=(\frac{m_2^2}{m_1}+m_2)u_{2^{\prime }}{\text{'}}^2+(-2m_2{u}_{1^{\prime }})u_{2^{″}}\\ 0=\left(\frac{m_2+m_1}{m_1}\right)u_{2^{\prime }}{\text{'}}^2+(-2u_{1^{\prime }})u_{2^{″}}\end{array}$$ Calculando as raízes do segundo grau, temos: $$u_{2^{″}}=\frac{2u_{1^{\prime }}\pm \sqrt{4u_1{\text{'}}^2}}{2\left(\frac{m_2+m_1}{m_1}\right)}=\frac{m_1}{m_2+m_1}(u_{1^{\prime }}\pm u_{1^{\prime }})$$E como uma solução é $u_{2^{\prime }}=0$, mas queremos a situação em que $u_{2^{\prime }}\ne 0$, logo: $$u_{2^{″}}=\frac{2m_1}{m_2+m_1}{u}_{1^{\prime }}$$Substituindo em (1), temos: $$\begin{array}{c}{u}_{1^{″}}=u_{1^{\prime }}-\frac{m_2}{m_1}\left(\frac{2m_1}{m_2+m_1}{u}_{1^{\prime }}\right)\\ u_{1^{″}}=(1-\frac{2m_2}{m_2+m_1})u_1\\ u_{1^{″}}=\left(\frac{m_2+m_1-2m_2}{m_2+m_1}\right)u_1\end{array}$$ Ou seja: $$u_{1^{″}}=\frac{m_1-m_2}{m_2+m_1}{u}_{1^{\prime }}$$Retornando ao referencial original, sendo $v_i=u_i+{v_2}^{\prime }$, temos então para ${v_1}^{\prime }$: $$\begin{array}{c}{u}_{1^{″}}=\frac{m_1-m_2}{m_2+m_1}{u}_{1^{\prime }}\\ v_{1^{″}}-v_{2^{\prime }}=\frac{m_1-m_2}{m_2+m_1}(v_{1^{\prime }}-v_{2^{\prime }})\\ v_{1^{″}}=\frac{m_1-m_2}{m_2+m_1}{v}_{1^{\prime }}+(1-\frac{m_1-m_2}{m_2+m_1})v_{2^{\prime }}\\ v_{1^{″}}=\frac{m_1-m_2}{m_2+m_1}{v}_{1^{\prime }}+\frac{2m_2}{m_2+m_1}{v}_{2^{\prime }}\end{array}$$E de maneira análoga para ${v_2}^{\prime }$: $$\begin{array}{c}{u}_{2^{″}}=\frac{2m_1}{m_2+m_1}{u}_{1^{\prime }}\\ v_{2^{″}}-v_{2^{\prime }}=\frac{2m_1}{m_2+m_1}(v_{1^{\prime }}-v_{2^{\prime }})\\ v_{2^{″}}=\frac{2m_1}{m_2+m_1}{v}_{1^{\prime }}+(1-\frac{2m_1}{m_2+m_1})v_{2^{\prime }}\\ v_{2^{″}}=\frac{2m_1}{m_2+m_1}{v}_{1^{\prime }}+\frac{m_2-m_1}{m_2+m_1}{v}_{2^{\prime }}\end{array}$$Um caso especial ocorre se $m_1=m_2,(2)$ então temos simplesmente: $${v_1}^{″}={v_2}^{\prime },\text{e}{v_2}^{″}={v_1}^{\prime }$$Em duas dimensões, podemos reduzir o problema a uma dimensão, considerando que toda a alteração na velocidade devido a colisão entre partículas ocorre apenas na componente paralela a reta que liga o centro das duas esferas. Considerando que a posição de cada partícula é dada por ${\overrightarrow{r}}_i$, então um vetor entre as partículas pode ser escrito como $\overrightarrow{d}={\overrightarrow{r}}_2-{\overrightarrow{r}}_1$. Podemos projetar ambas as velocidades então fazendo: $$u_{i^{\prime }}=\frac{{\overrightarrow{v}}_i\cdot \overrightarrow{d}}{\left|\overrightarrow{d}\right|}={\overrightarrow{v}}_i\cdot \hat{d}$$Obtemos o módulo da velocidade da partícula $i$ na direção $\overrightarrow{d}$ e podemos trabalhar em uma única dimensão para encontrarmos a velocidade de ambas partículas após a colisão nesta dimensão. Ao fim podemos decompor novamente a velocidade final ${u_i}^{″}$ em ambos os eixos mantendo a mesma direção utilizando $\theta =\arctan \left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)$, onde $\mathrm{\Delta }z=z_2-z_1$ é a diferença entre a posição das partículas $2$ e $1$ na componente $z$.

Além disso vale lembrar que há a componente da velocidade perpendicular a $\overrightarrow{d}$, que vamos denotar como ${\overrightarrow{w}}_i={\overrightarrow{v}}_i-u_{i^{\prime }}(\cos \theta ,\sin \theta )$. Esta componente perpendicular permanece inalterada e pode ser visualizada na figura ao lado.

Colisão entre duas partículas mostrando explicitamente os vetores relacionados à partícula 1.

Sendo assim, a velocidade final é dada por: $$\begin{array}{c}{\overrightarrow{v}}_i^{\left(f\right)}=(u_{i^{″}}\cos \theta ,u_{i^{″}}\sin \theta )+{\overrightarrow{w}}_i\\ {\overrightarrow{v}}_i^{\left(f\right)}=u_{i^{″}}(\cos \theta ,\sin \theta )+{\overrightarrow{v}}_i-u_{i^{\prime }}(\cos \theta ,\sin \theta )\\ {\overrightarrow{v}}_i^{\left(f\right)}={\overrightarrow{v}}_i+(u_{i^{″}}-u_{i^{\prime }})(\cos \theta ,\sin \theta )\\ {\overrightarrow{v}}_i^{\left(f\right)}={\overrightarrow{v}}_i^{\left(f\right)}={\overrightarrow{v}}_i+(u_{i^{″}}-u_{i^{\prime }})\overrightarrow{a}\end{array}$$Onde $\overrightarrow{a}$ é um vetor unitário que nos dá a direção entre os centros das partículas. Utilizando ad identidades: $$\cos (\arctan \left(x\right))=\frac{1}{\sqrt{1+x^2}},\sin (\arctan \left(x\right))=\frac{x}{\sqrt{1+x^2}}$$ ficamos então com: $$\begin{array}{c}\overrightarrow{a}=(\cos \theta ,\sin \theta )\\ \overrightarrow{a}=(\cos \arctan \left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right),\sin \arctan \left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right))\\ \overrightarrow{a}=(\frac{1}{\sqrt{1+{\left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)}^2}},\frac{\left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)}{\sqrt{1+{\left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)}^2}})\\ \overrightarrow{a}=\frac{1}{\mathrm{\Delta }x}\frac{1}{\sqrt{1+{\left(\frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}\right)}^2}}(\mathrm{\Delta }x,\mathrm{\Delta }y)\\ \overrightarrow{a}=\frac{1}{\frac{\mathrm{\Delta }x}{\mathrm{\Delta }x}}\frac{1}{\sqrt{\mathrm{\Delta }{x}^2+\mathrm{\Delta }{y}^2}}(\mathrm{\Delta }x,\mathrm{\Delta }y)\end{array}$$Logo: $$\overrightarrow{a}=\frac{(\mathrm{\Delta }x,\mathrm{\Delta }y)}{\sqrt{\mathrm{\Delta }{x}^2+\mathrm{\Delta }{y}^2}}$$ E uma vez que $\overrightarrow{d}={\overrightarrow{r}}_2-{\overrightarrow{r}}_1=(x_2-x_1,y_2-y_1)=(\mathrm{\Delta }x,\mathrm{\Delta }y)$, então $\overrightarrow{a}=\overrightarrow{d}/\left|d\right|=\hat{d}$. Logo: $${\overrightarrow{v}}_i^{\left(f\right)}={\overrightarrow{v}}_i+(u_{i^{″}}-u_{i^{\prime }})\hat{d}$$Ou ainda mais explícito, se fizermos $m_1=m_2$, sendo as partículas $i$ e $j$, onde $i\ne j$, usando (2): $$\begin{array}{cc}{\overrightarrow{v}}_i^{\left(f\right)}& ={\overrightarrow{v}}_i+(u_{j^{\prime }}-u_{i^{\prime }})\hat{d}\\ {\overrightarrow{v}}_i^{\left(f\right)}& ={\overrightarrow{v}}_i+({\overrightarrow{v}}_j\cdot \hat{d}-{\overrightarrow{v}}_i\cdot \hat{d})\hat{d}\end{array}$$ Temos então que: $${\overrightarrow{v}}_i^{\left(f\right)}={\overrightarrow{v}}_i+[({\overrightarrow{v}}_j-{\overrightarrow{v}}_i)\cdot \hat{d}]\hat{d}$$Todo o cálculo exibido foi para uma partícula, para a segunda partícula, o cálculo é análogo.

Distribuição de velocidades para um gás a certa temperatura

A distribuição das velocidads das molécula de um gá é dada pela distribuição de Maxwell-Boltmann

$$f\left(v\right)={\left(\frac{m}{2\pi kT}\right)}^{3/2}4\pi v^2\exp (-\frac{m{v}^2}{kT})$$

Foi realizado uma simulação com 100 partículas por 500.000 de passos em uma grade 100x100 com fronteira toroidal, o resultado encontra-se na figura abaixo assim como o código.

Resultado da simulação com a distribuição de Maxwell-Boltmann para kT=0.7 e m=1.

Código

Aqui chamamos a atenção que a principal diferença do código discutido para autômatos celulares é que agora não temos mais uma grade, mas um modelo contínuo. Também temos a possibilidade de configuração se as fronteiras serão toroidais (extremidades conectadas) ou não. Quanto aos resultados, precisamos considerar que é um modelo simplificado. O Calculo foi realizado para uma colisão entre duas partículas se dentro do mesmo passo deveria ocorrer mais de uma colisão, ou com mais de uma partícula simultaneamente o código falha, para remediar isso, o ideal é utilizar passos de tempo pequenos e/ou uma caixa relativamente grande em relação a quantidade de partículas. Além disso ao optarmos por uma fronteira toroidal ou não, estamos considerando uma caixa fechada ou não, ao considerarmos precisamos considerar a colisão com as fronteira também. Nesse último caso, uma simplificação simples é considerar que esta colisão ocorre primeiro. Novamente o problema surge se houver duas colisões no mesmo passo (agora podendo colidir com a fronteira), logo a recomendação é a mesma.

#Bibliotecas necessárias
from mesa import Agent, Model                                      #Classes Agente e Modelo
from mesa.time import SimultaneousActivation                       #Agendador simultâneo
from mesa.space import ContinuousSpace                             #Malha multigrid
import random                                                      #Número aleatórios
from mesa.datacollection import DataCollector                      #Importamos um coletor de dados
import matplotlib.pyplot as plt                                    #Plotar gráicos             
#from IPython.display import clear_output                           #Configurações de saída do Colab         
import numpy as np    

#AGENTE---------------------------------------------------------------------------------
class Agente(Agent):
    """Classe do agente"""
    def __init__(self, modelo,largura,altura):
        """Bibliotecas necessárias"""
        #modelo     - Modelo que ao qual o agente pertence
        super().__init__(self, modelo)        #Necessário para funcionar o modelo
        self.ppos     = (0,0)                 #Poição que irá se mover
        self.rai= 0.2                         #Raio
        self.lim=(largura,altura)             #Limites da caixa
        v=1                                   #Módulo da velocidade inicial
        ang = 2*random.random()*np.pi
        self.vel      = (v*np.cos(ang),v*np.sin(ang)) #Velocidade decomposta
        self.colisao=False
        self.ac=False
                
    def prox_pos(self):
      dt=0.01
      nx=self.pos[0]+self.vel[0]*0.01
      ny=self.pos[1]+self.vel[1]*0.01
      #Testar os limites
      if (nx < 0):
        nx = - nx; self.vel=(-self.vel[0],self.vel[1])
      if (nx > self.lim[0]):
        nx = 2*self.lim[0] - nx; self.vel=(-self.vel[0],self.vel[1])
      if (ny < 0):
        ny = - ny; self.vel=(self.vel[0], - self.vel[1])
      if (ny > self.lim[1]):
        ny = 2*self.lim[1] - ny;self.vel=(self.vel[0], - self.vel[1])

      return((nx,ny))
    
    def prox_vel(self):
        """Movimentação dos animais"""
        #Checamos se há colisão com outro animal:
        agentes = self.model.continuous.get_neighbors(pos=self.ppos, radius = 3*self.rai, include_center=True)
        for a in agentes:                                         #Percorremos a lista
          if (a.unique_id!=self.unique_id): 
              d=np.sqrt((a.ppos[0]-self.ppos[0])**2+(a.ppos[1]-self.ppos[1])**2)
              if (d<=2*self.rai):
                #Direção da colisão
                d=(a.ppos[0]-self.ppos[0],a.ppos[1]-self.ppos[1])/np.sqrt((a.ppos[0]-self.ppos[0])**2+(a.ppos[1]-self.ppos[1])**2)
                #Partícula 1
                u1=(a.vel[0]-self.vel[0])*d[0]+(a.vel[1]-self.vel[1])*d[1]
                nvel1=(self.vel[0]+u1*d[0],self.vel[1]+u1*d[1])
                #partícula 2,pela conservação de momento:
                nvel2=(self.vel[0]+a.vel[0]-nvel1[0],self.vel[1]+a.vel[1]-nvel1[1])   
                #Se nenhuma partícula tinha se colidido
                if (self.colisao==False and a.colisao==False):          
                  #Atualiza as velocidades
                  self.vel=nvel1; a.vel=nvel2
                  self.ac=True;a.ac=True   
                #Registra a colisão
                self.colisao=True; a.colisao=True       
                break
        return ()                                          

    def movimento(self):
      if(self.colisao):
        return (self.pos)
      else:
        return(self.ppos)
    
    def step(self):
        """Método obrigatório que prepara as mudanças"""
        self.ppos = self.prox_pos()
             
    def advance(self):
        """Método obrigatório que aplica as mudanças"""
        self.prox_vel()
        self.model.continuous.move_agent(self, self.movimento())

#MODELO
class Modelo(Model):
    """Modelo geral"""
    def gerar_particulas(self,N,largura,altura):
      """Função para gerar animais iniciais"""    
      #N        - Número de animais
      #largura  - Largura da grade
      #altura   - Altura da grade
      for n in range(N):
          con=True
          while(con):
            con=False
            a = Agente(self,largura,altura)
            X= largura*random.random()
            Y= altura*random.random()
            if (len(self.schedule.agents)==0):
              self.schedule.add(a)
              self.continuous.place_agent(a, (X, Y))
            else:
              agentes = self.continuous.get_neighbors(pos=(X,Y), radius = 2*a.rai, include_center=True)
              for viz in agentes:
                con = (True) if (len(agentes)==0) else (con)  #Precisamos checar se nasceu no mesmo lugar de outro animal
              if(con==False):
                self.schedule.add(a)
                self.continuous.place_agent(a, (X, Y))
      return()

    def __init__(self, modelo,N,seed=None):
        """Função chamada quando o modelo é inicializazdo"""
        # Modelo   - Dicionário com especificações do modelo
        # N        - Número de partículas
        # seed     - Seed dos números aleatórios do modelo do mesa
        
        # Desempacotar o dicionário
        largura = modelo["Largura"];altura=modelo["Altura"]; random.seed(modelo["Seed"])                     
        self.continuous = ContinuousSpace(x_max= largura, y_max=altura, torus= True, x_min=0, y_min= 0) #Configura a grade
        self.schedule   = SimultaneousActivation(self)            #Configura o agendador
        self.running   = True                                     #Condiçao para seguir executando o modelo
        self.gerar_particulas(N,largura,altura)    #Distribuir as ovelhas
           
    def step(self):
        """Avançar um passo do modelo"""
        self.schedule.step()                              #Avançamos os agentes
        #Limpamos as colisões
        for a in self.schedule.agents:                    #Vamos percorrer os agentes
          a.colisao=False
#Parâmetros
MAX =1000000
N=100
modelo  =  {"Largura":100     ,"Altura":100      ,"Seed":1} 
M = Modelo(modelo,N)
for i in range(MAX):
    M.step()
    if ((i+1)%(MAX/100)==0):
      #clear_output() 
      print(str(100*(1+i)/MAX)+"%")

vel=[]
E=0
for a in M.schedule.agents:
    vel.append(np.sqrt(a.vel[0]**2+a.vel[1]**2))
    E+=a.vel[0]**2+a.vel[1]**2
print(E)
with open('velocidades0.2R100K1000.txt', 'w') as f:
    f.write(str(vel))

Para plotar o gráfico, podemos então usar simplesmente:

import numpy as np
import matplotlib.pyplot as plt

a,b,c=plt.hist(vel,8, density=True, facecolor='g', alpha=0.75)
m=1
K=0.7
x=np.arange(0,2.5, 0.00001)
y = (4/np.sqrt(np.pi))*((m/K)**(3/2))*(x**2)*np.exp(-m*(x**2)/K)
plt.plot(x,y)

Citações

Anterior: Por que usar e o que são modelos baseados em indivíduos | Índice: Ecologia | Próximo: MBA: Caminhante aleatório