'''Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara'''

'''Objetivo:''' O trabalho tem como objetivo descrever a dinâmica molecular de um conjunto de partículas-planares infinitas paralelas. o modelo simplificado para estudar o comportamento de muitos corpos interagentes ao longo de um único eixo espacial é usado para feixes de partículas carregadas e também para sistemas autogravitantes. O estudo desse tipo de interação é importante para ter insights de sistemas de maior dimensionalidades.

= O modelo =

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível [1] e também para sistemas auto-gravitantes [2]. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

== Descrição do Modelo ==

=== Disposição espacial ===

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

[[Arquivo:1paralelas.gif]]

=== Ausência de Efeitos de Colisão ===

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

[[Arquivo:2diagrama.png]]

:<math>
\begin{cases} mv_1 - mv_2 = -mv_1' + mv_2' \\ \\ \frac{v_1' + v_2'}{v_1 + v_2} = 1 \end{cases}
</math>

O sistema de equações resulta em: <math>v_1' = v_2 \ \ \ ; \ \ \ v_2' = v_1</math>

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x.

=== Interação entre as partículas ===

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

[[Arquivo:3corpolivre.png]]

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial <math>\sigma</math>:
:<math>
\oint \overrightarrow{E}. d\overrightarrow{S} = \frac{Q_{int}}{\epsilon_0}
</math>
:<math>
E.A - (-E.A) = \frac{\sigma . A }{\epsilon_0}
</math>
:<math>
E= \frac{\sigma }{2\epsilon_0}
</math>

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

=== Interferência externa do sistema ===

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

[[Arquivo:4graficoforça.png]]

= Comportamento de uma Partícula no Sistema =

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

:<math>
E_i(x) = -\lambda\Theta[-(x-x_i)] + \lambda\Theta(x-x_i)
</math>

[[Arquivo:5perfilcampo.png]]

Em que a constante <math>\lambda</math> pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função <math>\Theta</math> equivale a função degrau de Heaviside:

:<math>
\Theta(x) = \begin{cases} 0, \ \ x<0 \\ 1, \ \ x>0 \end{cases}
</math>

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

[[Arquivo:6forçaresultante.png]]

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a <math>(N-1)</math>.

É possível notar que para uma distribuição de partículas-plano do tipo:

:<math>
n(x) = \sum_{i=1}^N \delta(x_i)
</math>

com o índice <math>\ \ ''i''\ </math> sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

:<math>
F_i = [-(N-1) + 2(i-1)]F
</math>
:<math>
F_i = (2i - N-1)F
</math>

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

:<math>
F^{grav.}_{R,i} = (2i - N-1)F
</math>

:<math>
F^{el.}_{R,i} = (2i - N-1)F - \kappa x_i
</math>

= Conservação da Energia do Sistema =

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

== Energia Cinética ==

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

:<math>
En_{cin} = \frac{m}{2}\sum_{i=1}^N \dot{x_i}^2
</math>

== Energia Potencial Restauradora ==

Observada no sistema eletricamente carregado citado, pode ser descrita como:

:<math>
En_{pot} = \frac{\kappa}{2}\sum_{i=1}^N x_i^2
</math>

== Energia Potencial Intrínseca ==

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições <math>x=L</math> e <math>x= -L</math> (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são <math>V(-L) = V(L) = 0</math>.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, <math>x_0</math> e <math>x_{N+1}</math>, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

:<math>
\Delta V_i = (x_{i} - x_{i - 1})E_i
</math>

com <math>E_i</math> sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

:<math>
E_i = (2i - N - 1)\lambda
</math>

[[Arquivo:7perfilcampoN4.png]]

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

:<math>
En_k = \alpha \sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante <math>\alpha</math> que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

:<math>
En_{int} = \alpha \sum_{k = 1}^{N}\sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

:<math>
E_{total} = En_{cin} + En_{pot} + En_{int} = cte.
</math>

= Modelagem Computacional =

== Criando a Distribuição do Sistema ==

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas <math>(x_i,y_i)</math> existirá uma partícula espelhada
com coordenadas <math>(-x_i,-y_i)</math>.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

<source lang='python'>

# Importando as bibliotecas necessárias
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[] #criamos uma lista vazia

N = 1000 # Determinamos o número de partículas utilizado

for i in range(int(N)):

A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados

A = np.array(A)
X= np.array([]) #Arrays usados para exemplificar
V= np.array([]) #distribuição inicial

for i in range (len(A)): #simetrização da distribuição
X=np.append(X,A[i][0])
V=np.append(V,A[i][1])
V=np.append(V,-A[i][1])
X=np.append(X,-A[i][0])

plt.xlabel('posição') #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

</source>

[[Arquivo:8espaçofaseinicial.png]]

== Função de Indexação ==

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo <math>x</math>.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e <math>2N</math> linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

<source lang='python'>

teste = np.array([2,6,3,9,0,1,4]) #criado um vetor de teste para que seus valores sejam indexados

def indexador(arr): #define-se a função

ordenador = [] #é criada uma lista vazia

for i in range(len(arr)):

ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual

ordenador.sort() # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento

for i in range(len(arr)):

ordenador[i].append(i) # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência

for i in range(len(arr)):

a = ordenador[i][0]
b = ordenador[i][1] # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
ordenador[i][0] = b
ordenador[i][1] = a

ordenador.sort() # reordenar os elementos de acordo com a ordem antiga

lista = []

for i in range(len(arr)):

lista.append(ordenador[i][2]) #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável

indexador = np.array(lista)

return indexador # e essa ordem é retornada como saída da função

print(indexador(teste))

</source>

== Função de movimento e aceleração das partículas ==

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado (<math>y</math>).

Para as coordenadas espaciais x (<math>y[i], \ 0\leq i< N </math>) temos que:

:<math>
\dot{x_i} = y_{[i + N]}
</math>

Para as coordenadas de velocidade v (<math>y[i], \ N\leq i< 2N-1 </math>) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias (<math>m = q = \kappa = 1</math>) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave (<math>F = 1/N</math>). O programa será escrito levando em conta o modelo de cargas elétricas:

:<math>
\dot{v_i} = \ddot{x_i} = \frac{1}{m} \Big[ -\kappa y_{[i - N]} + (2i - N - 1)F \Big]
</math>

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice <math>1</math> possui <math>N</math> partículas à esquerda (x negativo) e <math>(N-1)</math> partículas a direita (x positivo), o que resulta uma força de interação de <math>F</math> para a direita. A de índice 2 terá força de interação <math>3F</math>, a de índice 3 terá <math>5F</math> e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

:<math>
\ddot{x_i} = \frac{1}{m}\Big[ -\kappa y_{[i - N]} + \frac{y_{[i - N]}}{|y_{[i - N]}|}(2i- 1)F \Big]
</math>

<source lang='python'>

def deriv(t, y):

p = len(y)
k = indexador(y[0:int(p/2)])
sinal=abs(y[0:int(p/2)])/y[0:int(p/2)]

dydt = np.zeros(int(p))

dydt[0:int(p/2)] = y[int(p/2):int(p)] #definindo as derivadas temporais de posição

#definindo as derivadas temporais de velocidade
dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*sinal

deriv = dydt[0:int(p)]

return deriv

</source>

== Função de Energia do estado ==

<source lang='python'>

def Energia_total(y): # Para calcular a energia precisamos apenas do vetor de estado y do sistema de partículas utilizado

l = int(len(y)/2) # a variável l representa a quantidade de valores que teremos
pot = np.zeros(l+1) # cria-se um array com l + 1 dimensões para representar as l partículas e os respectivos potenciais

r = abs(np.array(y[0:l])) # cria-se um array r com as coordenadas de posição do estado y
r = np.append(r,4) # adiciona-se uma coordenada para representar o contorno externo do experimento
r.sort() # Faz-se o ordenamento desse array por distância da origem

pot_0 = 0 # O potencial do contorno deve ser nulo

for i in range(l,0,-1): # Loop decrescente que define o potencial de posição para cada partícula

pot[i] = pot_0

pot_0 = pot_0 + (1/2)*(i/l)*abs((r[i] - r[i-1]))

pot[0] = pot_0

v = np.array(y[l:-1]) # Cria-se um array com os valores de velocidade

E_total = 0

for i in range(l-1): #Soma-se as energias de cada partícula

E_total = E_total + pot[i] + 0.5*v[i]**2 + 0.5*r[i]**2

return 2*E_total, r[0:l+1], pot[0:l+1]

</source>

== Evolução da Distribuição no Tempo ==

Para simular a evolução temporal da distribuição inicial escolhida, precisamos de um método de integração numérica no tempo. O método escolhido deverá partir do estado da distribuição e da função de derivada temporal atribuída para estimar como a evolução ocorre.

<source lang='python'>

from scipy.integrate import solve_ivp #integrador numérico escolhido

def nova_fase(y,i): #função que define a função do espaço de fase em cada instante de tempo, adaptada para o integrador numérico escolhido

fase = np.array([])

for j in range(len(y)):

fase = np.append(fase,y[j][i])

return fase

Fase = [] # Cria uma lista para reservar a evolução do espaço de fase no tempo
T = np.array([]) # Cria uma lista que registra os tempos de marcação do integrador numérico
E = np.array([]) # Cria uma lista que registra a energia da distribuição no tempo
E_pot = [] # Cria uma lista que registra o comportamento da energia potencial ao longo de x no tempo
X_pot = [] # O análogo ao anterior para a posição
tempo = 0.0 # o tempo começa em 0
delta_t = 1 # o intervalo de tempo é 0.1

X = np.random.rand(4000)
V = np.random.rand(4000) - np.ones(4000)*0.5
y = np.append(X,2*V) #define-se a fase inicial

for j in range(50):

sol = solve_ivp(deriv, (tempo,tempo + delta_t), y , rtol = 0.00001) # integrador numérico realiza a integração

T= np.append(T, sol.t)

tempo = tempo + delta_t
for i in range(len(sol.t)): #o integrador possui um número de vetores intermediários ao vetor da solução final sol.y[-1]

y = nova_fase(sol.y,i) # aqui são registradas as funções definidas para cada estado de integração
Fase.append(y)

Ei = Energia_total(y)

E = np.append(E, Ei[0])

E_pot.append(Ei[2])
X_pot.append(Ei[1])

y = nova_fase(sol.y, -1) # A distribuição atual é atualizada para nova integração

</source>

=== Animações ===

Foi utilizado o código abaixo para gerar os gifs da evolução temporal

<source lang='python'>

from matplotlib.animation import FuncAnimation
from IPython import display

# código base para montar gifs da evolução temporal
fig, (ax1, ax2, ax3) = plt.subplots(3,1 ,figsize = (6,14)) #cria uma figura de três linhas

lines = [] #objeto padrão da biblioteca FuncAnimation

def animate(frame):

# As entradas são limpas para evitar superposição de imagens nos frames
ax1.cla()
ax2.cla()
ax3.cla()

# Plot do espaço de fases
ax1.scatter(Fase[frame][0:4000], Fase[frame][4000:8000], color = 'k', s =5)
ax1.scatter(-Fase[frame][0:4000], -Fase[frame][4000:8000], color = 'k', s = 5)
ax1.set_xlabel('$posição$')
ax1.set_ylabel('$velocidade$')
ax1.set_xlim(-4,4)
ax1.set_ylim(-2.5,2.5)

# Plot da distribuição de Potencial
ax2.plot(X_pot[frame],E_pot[frame], color = 'grey')
ax2.plot(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'grey')
ax2.scatter(X_pot[frame],E_pot[frame], color = 'b', s =5)
ax2.scatter(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'b', s =5)
ax2.set_ylabel("potencial")
ax2.set_xlabel("posição")

# Plot da energia total no tempo
ax3.plot(T[0:frame],E[0:frame], color = 'k')
ax3.set_ylabel("Energia")
ax3.set_xlabel("tempo")
ax3.set_xlim(0, T[-1])
ax3.set_ylim(E[0]- 10, E[0] + 10)

return lines

# Execução do programa e geração do gif inline
anima = FuncAnimation(fig, animate,frames=len(T), interval=20, blit=True)
video = anima.to_html5_video()
html = display.HTML(video)
display.display(html)

plt.close()

</source>

[[Arquivo:Mom-pot-e.gif]]

Arquivo:Mom-pot-e.gif

2021-11-29T16:11:30Z

Renanrs:

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-29T16:03:38Z

2021-11-22T11:26:31Z

Renanrs: /* Função de movimento e aceleração das partículas */

Em construção

'''Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara'''

'''Objetivo:'''

= O modelo =

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível [1] e também para sistemas auto-gravitantes [2]. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

== Descrição do Modelo ==

=== Disposição espacial ===

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

ADD IMAGEM

=== Ausência de Efeitos de Colisão ===

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

ADD IMAGEM

:<math>
\begin{cases} mv_1 - mv_2 = -mv_1' + mv_2' \\ \\ \frac{v_1' + v_2'}{v_1 + v_2} = 1 \end{cases}
</math>

O sistema de equações resulta em: <math>v_1' = v_2 \ \ \ ; \ \ \ v_2' = v_1</math>

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x.

=== Interação entre as partículas ===

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

ADD IMAGEM

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial <math>\sigma</math>:
:<math>
\oint \overrightarrow{E}. d\overrightarrow{S} = \frac{Q_{int}}{\epsilon_0}
</math>
:<math>
E.A - (-E.A) = \frac{\sigma . A }{\epsilon_0}
</math>
:<math>
E= \frac{\sigma }{2\epsilon_0}
</math>

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

=== Interferência externa do sistema ===

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

ADD IMAGEM

= Comportamento de uma Partícula no Sistema =

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

:<math>
E_i(x) = -\lambda\Theta[-(x-x_i)] + \lambda\Theta(x-x_i)
</math>

ADD IMAGEM

Em que a constante <math>\lambda</math> pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função <math>\Theta</math> equivale a função degrau de Heaviside:

:<math>
\Theta(x) = \begin{cases} 0, \ \ x<0 \\ 1, \ \ x>0 \end{cases}
</math>

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

ADD IMAGEM

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a <math>(N-1)</math>.

É possível notar que para uma distribuição de partículas-plano do tipo:

:<math>
n(x) = \sum_{i=1}^N \delta(x_i)
</math>

com o índice <math>\ \ ''i''\ </math> sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

:<math>
F_i = [-(N-1) + 2(i-1)]F
</math>
:<math>
F_i = (2i - N-1)F
</math>

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

:<math>
F^{grav.}_{R,i} = (2i - N-1)F
</math>

:<math>
F^{el.}_{R,i} = (2i - N-1)F - \kappa x_i
</math>

= Conservação da Energia do Sistema =

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

== Energia Cinética ==

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

:<math>
En_{cin} = \frac{m}{2}\sum_{i=1}^N \dot{x_i}^2
</math>

== Energia Potencial Restauradora ==

Observada no sistema eletricamente carregado citado, pode ser descrita como:

:<math>
En_{pot} = \frac{\kappa}{2}\sum_{i=1}^N x_i^2
</math>

== Energia Potencial Intrínseca ==

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições <math>x=L</math> e <math>x= -L</math> (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são <math>V(-L) = V(L) = 0</math>.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, <math>x_0</math> e <math>x_{N+1}</math>, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

:<math>
\Delta V_i = (x_{i} - x_{i - 1})E_i
</math>

com <math>E_i</math> sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

:<math>
E_i = (2i - N - 1)\lambda
</math>

ADD IMAGEM

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

:<math>
En_k = \alpha \sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante <math>\alpha</math> que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

:<math>
En_{int} = \alpha \sum_{k = 1}^{N}\sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

:<math>
E_{total} = En_{cin} + En_{pot} + En_{int} = cte.
</math>

= Modelagem Computacional =

== Criando a Distribuição do Sistema ==

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas <math>(x_i,y_i)</math> existirá uma partícula espelhada
com coordenadas <math>(-x_i,-y_i)</math>.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

<source lang='python'>

# Importando as bibliotecas necessárias
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[] #criamos uma lista vazia

N = 1000 # Determinamos o número de partículas utilizado

for i in range(int(N)):

A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados

A = np.array(A)
X= np.array([]) #Arrays usados para exemplificar
V= np.array([]) #distribuição inicial

for i in range (len(A)): #simetrização da distribuição
X=np.append(X,A[i][0])
V=np.append(V,A[i][1])
V=np.append(V,-A[i][1])
X=np.append(X,-A[i][0])

plt.xlabel('posição') #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

</source>

ADD IMAGEM - <matplotlib.collections.PathCollection at 0x7f4ab8366ad0>

== Função de Indexação ==

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo <math>x</math>.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e <math>2N</math> linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

<source lang='python'>

teste = np.array([2,6,3,9,0,1,4]) #criado um vetor de teste para que seus valores sejam indexados

def indexador(arr): #define-se a função

ordenador = [] #é criada uma lista vazia

for i in range(len(arr)):

ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual

ordenador.sort() # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento

for i in range(len(arr)):

ordenador[i].append(i) # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência

for i in range(len(arr)):

a = ordenador[i][0]
b = ordenador[i][1] # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
ordenador[i][0] = b
ordenador[i][1] = a

ordenador.sort() # reordenar os elementos de acordo com a ordem antiga

lista = []

for i in range(len(arr)):

lista.append(ordenador[i][2]) #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável

indexador = np.array(lista)

return indexador # e essa ordem é retornada como saída da função

print(indexador(teste))

</source>

== Função de movimento e aceleração das partículas ==

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado (<math>y</math>).

Para as coordenadas espaciais x (<math>y[i], \ 0\leq i< N </math>) temos que:

:<math>
\dot{x_i} = y_{[i + N]}
</math>

Para as coordenadas de velocidade v (<math>y[i], \ N\leq i< 2N-1 </math>) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias (<math>m = q = \kappa = 1</math>) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave (<math>F = 1/N</math>). O programa será escrito levando em conta o modelo de cargas elétricas:

:<math>
\dot{v_i} = \ddot{x_i} = \frac{1}{m} \Big[ -\kappa y_{[i - N]} + (2i - N - 1)F \Big]
</math>

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice <math>1</math> possui <math>N</math> partículas à esquerda (x negativo) e <math>(N-1)</math> partículas a direita (x positivo), o que resulta uma força de interação de <math>F</math> para a direita. A de índice 2 terá força de interação <math>3F</math>, a de índice 3 terá <math>5F</math> e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

:<math>
\ddot{x_i} = \frac{1}{m}\Big[ -\kappa y_{[i - N]} + \frac{y_{[i - N]}}{|y_{[i - N]}|}(2i- 1)F \Big]
</math>

<source lang='python'>

def deriv(t, y):

p = len(y)
k = indexador(y[0:int(p/2)])
sinal=abs(y[0:int(p/2)])/y[0:int(p/2)]

dydt = np.zeros(int(p))

dydt[0:int(p/2)] = y[int(p/2):int(p)] #definindo as derivadas temporais de posição

#definindo as derivadas temporais de velocidade
dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*sinal

deriv = dydt[0:int(p)]

return deriv

</source>

== Função de Energia do estado ==

<source lang='python'>

def Energia_total(y): # Para calcular a energia precisamos apenas do vetor de estado y do sistema de partículas utilizado

l = int(len(y)/2) # a variável l representa a quantidade de valores que teremos
pot = np.zeros(l+1) # cria-se um array com l + 1 dimensões para representar as l partículas e os respectivos potenciais

r = abs(np.array(y[0:l])) # cria-se um array r com as coordenadas de posição do estado y
r = np.append(r,4) # adiciona-se uma coordenada para representar o contorno externo do experimento
r.sort() # Faz-se o ordenamento desse array por distância da origem

pot_0 = 0 # O potencial do contorno deve ser nulo

for i in range(l,0,-1): # Loop decrescente que define o potencial de posição para cada partícula

pot[i] = pot_0

pot_0 = pot_0 + (1/2)*(i/l)*abs((r[i] - r[i-1]))

pot[0] = pot_0

v = np.array(y[l:-1]) # Cria-se um array com os valores de velocidade

E_total = 0

for i in range(l-1): #Soma-se as energias de cada partícula

E_total = E_total + pot[i] + 0.5*v[i]**2 + 0.5*r[i]**2

return 2*E_total, r[0:l+1], pot[0:l+1]

</source>

== Evolução da Distribuição no Tempo ==

Para simular a evolução temporal da distribuição inicial escolhida, precisamos de um método de integração numérica no tempo. O método escolhido deverá partir do estado da distribuição e da função de derivada temporal atribuída para estimar como a evolução ocorre.

<source lang='python'>

from scipy.integrate import solve_ivp #integrador numérico escolhido

def nova_fase(y,i): #função que define a função do espaço de fase em cada instante de tempo, adaptada para o integrador numérico escolhido

fase = np.array([])

for j in range(len(y)):

fase = np.append(fase,y[j][i])

return fase

Fase = [] # Cria uma lista para reservar a evolução do espaço de fase no tempo
T = np.array([]) # Cria uma lista que registra os tempos de marcação do integrador numérico
E = np.array([]) # Cria uma lista que registra a energia da distribuição no tempo
E_pot = [] # Cria uma lista que registra o comportamento da energia potencial ao longo de x no tempo
X_pot = [] # O análogo ao anterior para a posição
tempo = 0.0 # o tempo começa em 0
delta_t = 1 # o intervalo de tempo é 0.1

X = np.random.rand(4000)
V = np.random.rand(4000) - np.ones(4000)*0.5
y = np.append(X,2*V) #define-se a fase inicial

for j in range(50):

sol = solve_ivp(deriv, (tempo,tempo + delta_t), y , rtol = 0.00001) # integrador numérico realiza a integração

T= np.append(T, sol.t)

tempo = tempo + delta_t
for i in range(len(sol.t)): #o integrador possui um número de vetores intermediários ao vetor da solução final sol.y[-1]

y = nova_fase(sol.y,i) # aqui são registradas as funções definidas para cada estado de integração
Fase.append(y)

Ei = Energia_total(y)

E = np.append(E, Ei[0])

E_pot.append(Ei[2])
X_pot.append(Ei[1])

y = nova_fase(sol.y, -1) # A distribuição atual é atualizada para nova integração

</source>

=== Animações ===

Foi utilizado o código abaixo para gerar os gifs da evolução temporal

<source lang='python'>

from matplotlib.animation import FuncAnimation
from IPython import display

# código base para montar gifs da evolução temporal
fig, (ax1, ax2, ax3) = plt.subplots(3,1 ,figsize = (6,14)) #cria uma figura de três linhas

lines = [] #objeto padrão da biblioteca FuncAnimation

def animate(frame):

# As entradas são limpas para evitar superposição de imagens nos frames
ax1.cla()
ax2.cla()
ax3.cla()

# Plot do espaço de fases
ax1.scatter(Fase[frame][0:4000], Fase[frame][4000:8000], color = 'k', s =5)
ax1.scatter(-Fase[frame][0:4000], -Fase[frame][4000:8000], color = 'k', s = 5)
ax1.set_xlabel('$posição$')
ax1.set_ylabel('$velocidade$')
ax1.set_xlim(-4,4)
ax1.set_ylim(-2.5,2.5)

# Plot da distribuição de Potencial
ax2.plot(X_pot[frame],E_pot[frame], color = 'grey')
ax2.plot(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'grey')
ax2.scatter(X_pot[frame],E_pot[frame], color = 'b', s =5)
ax2.scatter(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'b', s =5)
ax2.set_ylabel("potencial")
ax2.set_xlabel("posição")

# Plot da energia total no tempo
ax3.plot(T[0:frame],E[0:frame], color = 'k')
ax3.set_ylabel("Energia")
ax3.set_xlabel("tempo")
ax3.set_xlim(0, T[-1])
ax3.set_ylim(E[0]- 10, E[0] + 10)

return lines

# Execução do programa e geração do gif inline
anima = FuncAnimation(fig, animate,frames=len(T), interval=20, blit=True)
video = anima.to_html5_video()
html = display.HTML(video)
display.display(html)

plt.close()

</source>

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-22T11:23:39Z

Renanrs: /* Função de movimento e aceleração das partículas */

Em construção

'''Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara'''

'''Objetivo:'''

= O modelo =

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível [1] e também para sistemas auto-gravitantes [2]. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

== Descrição do Modelo ==

=== Disposição espacial ===

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

ADD IMAGEM

=== Ausência de Efeitos de Colisão ===

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

ADD IMAGEM

:<math>
\begin{cases} mv_1 - mv_2 = -mv_1' + mv_2' \\ \\ \frac{v_1' + v_2'}{v_1 + v_2} = 1 \end{cases}
</math>

O sistema de equações resulta em: <math>v_1' = v_2 \ \ \ ; \ \ \ v_2' = v_1</math>

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x.

=== Interação entre as partículas ===

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

ADD IMAGEM

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial <math>\sigma</math>:
:<math>
\oint \overrightarrow{E}. d\overrightarrow{S} = \frac{Q_{int}}{\epsilon_0}
</math>
:<math>
E.A - (-E.A) = \frac{\sigma . A }{\epsilon_0}
</math>
:<math>
E= \frac{\sigma }{2\epsilon_0}
</math>

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

=== Interferência externa do sistema ===

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

ADD IMAGEM

= Comportamento de uma Partícula no Sistema =

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

:<math>
E_i(x) = -\lambda\Theta[-(x-x_i)] + \lambda\Theta(x-x_i)
</math>

ADD IMAGEM

Em que a constante <math>\lambda</math> pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função <math>\Theta</math> equivale a função degrau de Heaviside:

:<math>
\Theta(x) = \begin{cases} 0, \ \ x<0 \\ 1, \ \ x>0 \end{cases}
</math>

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

ADD IMAGEM

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a <math>(N-1)</math>.

É possível notar que para uma distribuição de partículas-plano do tipo:

:<math>
n(x) = \sum_{i=1}^N \delta(x_i)
</math>

com o índice <math>\ \ ''i''\ </math> sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

:<math>
F_i = [-(N-1) + 2(i-1)]F
</math>
:<math>
F_i = (2i - N-1)F
</math>

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

:<math>
F^{grav.}_{R,i} = (2i - N-1)F
</math>

:<math>
F^{el.}_{R,i} = (2i - N-1)F - \kappa x_i
</math>

= Conservação da Energia do Sistema =

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

== Energia Cinética ==

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

:<math>
En_{cin} = \frac{m}{2}\sum_{i=1}^N \dot{x_i}^2
</math>

== Energia Potencial Restauradora ==

Observada no sistema eletricamente carregado citado, pode ser descrita como:

:<math>
En_{pot} = \frac{\kappa}{2}\sum_{i=1}^N x_i^2
</math>

== Energia Potencial Intrínseca ==

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições <math>x=L</math> e <math>x= -L</math> (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são <math>V(-L) = V(L) = 0</math>.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, <math>x_0</math> e <math>x_{N+1}</math>, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

:<math>
\Delta V_i = (x_{i} - x_{i - 1})E_i
</math>

com <math>E_i</math> sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

:<math>
E_i = (2i - N - 1)\lambda
</math>

ADD IMAGEM

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

:<math>
En_k = \alpha \sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante <math>\alpha</math> que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

:<math>
En_{int} = \alpha \sum_{k = 1}^{N}\sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

:<math>
E_{total} = En_{cin} + En_{pot} + En_{int} = cte.
</math>

= Modelagem Computacional =

== Criando a Distribuição do Sistema ==

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas <math>(x_i,y_i)</math> existirá uma partícula espelhada
com coordenadas <math>(-x_i,-y_i)</math>.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

<source lang='python'>

# Importando as bibliotecas necessárias
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[] #criamos uma lista vazia

N = 1000 # Determinamos o número de partículas utilizado

for i in range(int(N)):

A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados

A = np.array(A)
X= np.array([]) #Arrays usados para exemplificar
V= np.array([]) #distribuição inicial

for i in range (len(A)): #simetrização da distribuição
X=np.append(X,A[i][0])
V=np.append(V,A[i][1])
V=np.append(V,-A[i][1])
X=np.append(X,-A[i][0])

plt.xlabel('posição') #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

</source>

ADD IMAGEM - <matplotlib.collections.PathCollection at 0x7f4ab8366ad0>

== Função de Indexação ==

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo <math>x</math>.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e <math>2N</math> linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

<source lang='python'>

teste = np.array([2,6,3,9,0,1,4]) #criado um vetor de teste para que seus valores sejam indexados

def indexador(arr): #define-se a função

ordenador = [] #é criada uma lista vazia

for i in range(len(arr)):

ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual

ordenador.sort() # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento

for i in range(len(arr)):

ordenador[i].append(i) # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência

for i in range(len(arr)):

a = ordenador[i][0]
b = ordenador[i][1] # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
ordenador[i][0] = b
ordenador[i][1] = a

ordenador.sort() # reordenar os elementos de acordo com a ordem antiga

lista = []

for i in range(len(arr)):

lista.append(ordenador[i][2]) #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável

indexador = np.array(lista)

return indexador # e essa ordem é retornada como saída da função

print(indexador(teste))

</source>

== Função de movimento e aceleração das partículas ==

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado (<math>y</math>).

Para as coordenadas espaciais x (<math>y[i], \ 0\leq i< N </math>) temos que:

:<math>
\dot{x_i} = y_{[i + N]}
</math>

Para as coordenadas de velocidade v (<math>y[i], \ N\leq i< 2N-1 </math>) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias (<math>m = q = \kappa = 1</math>) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave (<math>F = 1/N</math>). O programa será escrito levando em conta o modelo de cargas elétricas:

:<math>
\dot{v_i} = \ddot{x_i} = \frac{1}{m}[ -\kappa y_{[i - N]} + (2i - N - 1)F ]
</math>

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice <math>1</math> possui <math>N</math> partículas à esquerda (x negativo) e <math>(N-1)</math> partículas a direita (x positivo), o que resulta uma força de interação de <math>F</math> para a direita. A de índice 2 terá força de interação <math>3F</math>, a de índice 3 terá <math>5F</math> e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

:<math>
\ddot{x_i} = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + \frac{y_{[i - N]}}{|y_{[i - N]}|}(2i- 1)F \Big{]}
</math>

<source lang='python'>

def deriv(t, y):

p = len(y)
k = indexador(y[0:int(p/2)])
sinal=abs(y[0:int(p/2)])/y[0:int(p/2)]

dydt = np.zeros(int(p))

dydt[0:int(p/2)] = y[int(p/2):int(p)] #definindo as derivadas temporais de posição

#definindo as derivadas temporais de velocidade
dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*sinal

deriv = dydt[0:int(p)]

return deriv

</source>

== Função de Energia do estado ==

<source lang='python'>

def Energia_total(y): # Para calcular a energia precisamos apenas do vetor de estado y do sistema de partículas utilizado

l = int(len(y)/2) # a variável l representa a quantidade de valores que teremos
pot = np.zeros(l+1) # cria-se um array com l + 1 dimensões para representar as l partículas e os respectivos potenciais

r = abs(np.array(y[0:l])) # cria-se um array r com as coordenadas de posição do estado y
r = np.append(r,4) # adiciona-se uma coordenada para representar o contorno externo do experimento
r.sort() # Faz-se o ordenamento desse array por distância da origem

pot_0 = 0 # O potencial do contorno deve ser nulo

for i in range(l,0,-1): # Loop decrescente que define o potencial de posição para cada partícula

pot[i] = pot_0

pot_0 = pot_0 + (1/2)*(i/l)*abs((r[i] - r[i-1]))

pot[0] = pot_0

v = np.array(y[l:-1]) # Cria-se um array com os valores de velocidade

E_total = 0

for i in range(l-1): #Soma-se as energias de cada partícula

E_total = E_total + pot[i] + 0.5*v[i]**2 + 0.5*r[i]**2

return 2*E_total, r[0:l+1], pot[0:l+1]

</source>

== Evolução da Distribuição no Tempo ==

Para simular a evolução temporal da distribuição inicial escolhida, precisamos de um método de integração numérica no tempo. O método escolhido deverá partir do estado da distribuição e da função de derivada temporal atribuída para estimar como a evolução ocorre.

<source lang='python'>

from scipy.integrate import solve_ivp #integrador numérico escolhido

def nova_fase(y,i): #função que define a função do espaço de fase em cada instante de tempo, adaptada para o integrador numérico escolhido

fase = np.array([])

for j in range(len(y)):

fase = np.append(fase,y[j][i])

return fase

Fase = [] # Cria uma lista para reservar a evolução do espaço de fase no tempo
T = np.array([]) # Cria uma lista que registra os tempos de marcação do integrador numérico
E = np.array([]) # Cria uma lista que registra a energia da distribuição no tempo
E_pot = [] # Cria uma lista que registra o comportamento da energia potencial ao longo de x no tempo
X_pot = [] # O análogo ao anterior para a posição
tempo = 0.0 # o tempo começa em 0
delta_t = 1 # o intervalo de tempo é 0.1

X = np.random.rand(4000)
V = np.random.rand(4000) - np.ones(4000)*0.5
y = np.append(X,2*V) #define-se a fase inicial

for j in range(50):

sol = solve_ivp(deriv, (tempo,tempo + delta_t), y , rtol = 0.00001) # integrador numérico realiza a integração

T= np.append(T, sol.t)

tempo = tempo + delta_t
for i in range(len(sol.t)): #o integrador possui um número de vetores intermediários ao vetor da solução final sol.y[-1]

y = nova_fase(sol.y,i) # aqui são registradas as funções definidas para cada estado de integração
Fase.append(y)

Ei = Energia_total(y)

E = np.append(E, Ei[0])

E_pot.append(Ei[2])
X_pot.append(Ei[1])

y = nova_fase(sol.y, -1) # A distribuição atual é atualizada para nova integração

</source>

=== Animações ===

Foi utilizado o código abaixo para gerar os gifs da evolução temporal

<source lang='python'>

from matplotlib.animation import FuncAnimation
from IPython import display

# código base para montar gifs da evolução temporal
fig, (ax1, ax2, ax3) = plt.subplots(3,1 ,figsize = (6,14)) #cria uma figura de três linhas

lines = [] #objeto padrão da biblioteca FuncAnimation

def animate(frame):

# As entradas são limpas para evitar superposição de imagens nos frames
ax1.cla()
ax2.cla()
ax3.cla()

# Plot do espaço de fases
ax1.scatter(Fase[frame][0:4000], Fase[frame][4000:8000], color = 'k', s =5)
ax1.scatter(-Fase[frame][0:4000], -Fase[frame][4000:8000], color = 'k', s = 5)
ax1.set_xlabel('$posição$')
ax1.set_ylabel('$velocidade$')
ax1.set_xlim(-4,4)
ax1.set_ylim(-2.5,2.5)

# Plot da distribuição de Potencial
ax2.plot(X_pot[frame],E_pot[frame], color = 'grey')
ax2.plot(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'grey')
ax2.scatter(X_pot[frame],E_pot[frame], color = 'b', s =5)
ax2.scatter(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'b', s =5)
ax2.set_ylabel("potencial")
ax2.set_xlabel("posição")

# Plot da energia total no tempo
ax3.plot(T[0:frame],E[0:frame], color = 'k')
ax3.set_ylabel("Energia")
ax3.set_xlabel("tempo")
ax3.set_xlim(0, T[-1])
ax3.set_ylim(E[0]- 10, E[0] + 10)

return lines

# Execução do programa e geração do gif inline
anima = FuncAnimation(fig, animate,frames=len(T), interval=20, blit=True)
video = anima.to_html5_video()
html = display.HTML(video)
display.display(html)

plt.close()

</source>

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-22T11:22:17Z

Renanrs:

Em construção

'''Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara'''

'''Objetivo:'''

= O modelo =

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível [1] e também para sistemas auto-gravitantes [2]. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

== Descrição do Modelo ==

=== Disposição espacial ===

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

ADD IMAGEM

=== Ausência de Efeitos de Colisão ===

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

ADD IMAGEM

:<math>
\begin{cases} mv_1 - mv_2 = -mv_1' + mv_2' \\ \\ \frac{v_1' + v_2'}{v_1 + v_2} = 1 \end{cases}
</math>

O sistema de equações resulta em: <math>v_1' = v_2 \ \ \ ; \ \ \ v_2' = v_1</math>

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x.

=== Interação entre as partículas ===

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

ADD IMAGEM

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial <math>\sigma</math>:
:<math>
\oint \overrightarrow{E}. d\overrightarrow{S} = \frac{Q_{int}}{\epsilon_0}
</math>
:<math>
E.A - (-E.A) = \frac{\sigma . A }{\epsilon_0}
</math>
:<math>
E= \frac{\sigma }{2\epsilon_0}
</math>

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

=== Interferência externa do sistema ===

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

ADD IMAGEM

= Comportamento de uma Partícula no Sistema =

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

:<math>
E_i(x) = -\lambda\Theta[-(x-x_i)] + \lambda\Theta(x-x_i)
</math>

ADD IMAGEM

Em que a constante <math>\lambda</math> pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função <math>\Theta</math> equivale a função degrau de Heaviside:

:<math>
\Theta(x) = \begin{cases} 0, \ \ x<0 \\ 1, \ \ x>0 \end{cases}
</math>

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

ADD IMAGEM

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a <math>(N-1)</math>.

É possível notar que para uma distribuição de partículas-plano do tipo:

:<math>
n(x) = \sum_{i=1}^N \delta(x_i)
</math>

com o índice <math>\ \ ''i''\ </math> sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

:<math>
F_i = [-(N-1) + 2(i-1)]F
</math>
:<math>
F_i = (2i - N-1)F
</math>

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

:<math>
F^{grav.}_{R,i} = (2i - N-1)F
</math>

:<math>
F^{el.}_{R,i} = (2i - N-1)F - \kappa x_i
</math>

= Conservação da Energia do Sistema =

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

== Energia Cinética ==

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

:<math>
En_{cin} = \frac{m}{2}\sum_{i=1}^N \dot{x_i}^2
</math>

== Energia Potencial Restauradora ==

Observada no sistema eletricamente carregado citado, pode ser descrita como:

:<math>
En_{pot} = \frac{\kappa}{2}\sum_{i=1}^N x_i^2
</math>

== Energia Potencial Intrínseca ==

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições <math>x=L</math> e <math>x= -L</math> (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são <math>V(-L) = V(L) = 0</math>.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, <math>x_0</math> e <math>x_{N+1}</math>, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

:<math>
\Delta V_i = (x_{i} - x_{i - 1})E_i
</math>

com <math>E_i</math> sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

:<math>
E_i = (2i - N - 1)\lambda
</math>

ADD IMAGEM

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

:<math>
En_k = \alpha \sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante <math>\alpha</math> que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

:<math>
En_{int} = \alpha \sum_{k = 1}^{N}\sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

:<math>
E_{total} = En_{cin} + En_{pot} + En_{int} = cte.
</math>

= Modelagem Computacional =

== Criando a Distribuição do Sistema ==

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas <math>(x_i,y_i)</math> existirá uma partícula espelhada
com coordenadas <math>(-x_i,-y_i)</math>.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

<source lang='python'>

# Importando as bibliotecas necessárias
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[] #criamos uma lista vazia

N = 1000 # Determinamos o número de partículas utilizado

for i in range(int(N)):

A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados

A = np.array(A)
X= np.array([]) #Arrays usados para exemplificar
V= np.array([]) #distribuição inicial

for i in range (len(A)): #simetrização da distribuição
X=np.append(X,A[i][0])
V=np.append(V,A[i][1])
V=np.append(V,-A[i][1])
X=np.append(X,-A[i][0])

plt.xlabel('posição') #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

</source>

ADD IMAGEM - <matplotlib.collections.PathCollection at 0x7f4ab8366ad0>

== Função de Indexação ==

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo <math>x</math>.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e <math>2N</math> linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

<source lang='python'>

teste = np.array([2,6,3,9,0,1,4]) #criado um vetor de teste para que seus valores sejam indexados

def indexador(arr): #define-se a função

ordenador = [] #é criada uma lista vazia

for i in range(len(arr)):

ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual

ordenador.sort() # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento

for i in range(len(arr)):

ordenador[i].append(i) # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência

for i in range(len(arr)):

a = ordenador[i][0]
b = ordenador[i][1] # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
ordenador[i][0] = b
ordenador[i][1] = a

ordenador.sort() # reordenar os elementos de acordo com a ordem antiga

lista = []

for i in range(len(arr)):

lista.append(ordenador[i][2]) #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável

indexador = np.array(lista)

return indexador # e essa ordem é retornada como saída da função

print(indexador(teste))

</source>

== Função de movimento e aceleração das partículas ==

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado (<math>y</math>).

Para as coordenadas espaciais x (<math>y[i], \ 0\leq i< N </math>) temos que:

:<math>
\dot{x_i} = y_{[i + N]}
</math>

Para as coordenadas de velocidade v (<math>y[i], \ N\leq i< 2N-1 </math>) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias (<math>m = q = \kappa = 1</math>) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave (<math>F = 1/N</math>). O programa será escrito levando em conta o modelo de cargas elétricas:

:<math>
\dot{v_i} = \ddot{x_i} = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + (2i - N - 1)F \Big{]}
</math>

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice <math>1</math> possui <math>N</math> partículas à esquerda (x negativo) e <math>(N-1)</math> partículas a direita (x positivo), o que resulta uma força de interação de <math>F</math> para a direita. A de índice 2 terá força de interação <math>3F</math>, a de índice 3 terá <math>5F</math> e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

:<math>
\ddot{x_i} = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + \frac{y_{[i - N]}}{|y_{[i - N]}|}(2i- 1)F \Big{]}
</math>

<source lang='python'>

def deriv(t, y):

p = len(y)
k = indexador(y[0:int(p/2)])
sinal=abs(y[0:int(p/2)])/y[0:int(p/2)]

dydt = np.zeros(int(p))

dydt[0:int(p/2)] = y[int(p/2):int(p)] #definindo as derivadas temporais de posição

#definindo as derivadas temporais de velocidade
dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*sinal

deriv = dydt[0:int(p)]

return deriv

</source>

== Função de Energia do estado ==

<source lang='python'>

def Energia_total(y): # Para calcular a energia precisamos apenas do vetor de estado y do sistema de partículas utilizado

l = int(len(y)/2) # a variável l representa a quantidade de valores que teremos
pot = np.zeros(l+1) # cria-se um array com l + 1 dimensões para representar as l partículas e os respectivos potenciais

r = abs(np.array(y[0:l])) # cria-se um array r com as coordenadas de posição do estado y
r = np.append(r,4) # adiciona-se uma coordenada para representar o contorno externo do experimento
r.sort() # Faz-se o ordenamento desse array por distância da origem

pot_0 = 0 # O potencial do contorno deve ser nulo

for i in range(l,0,-1): # Loop decrescente que define o potencial de posição para cada partícula

pot[i] = pot_0

pot_0 = pot_0 + (1/2)*(i/l)*abs((r[i] - r[i-1]))

pot[0] = pot_0

v = np.array(y[l:-1]) # Cria-se um array com os valores de velocidade

E_total = 0

for i in range(l-1): #Soma-se as energias de cada partícula

E_total = E_total + pot[i] + 0.5*v[i]**2 + 0.5*r[i]**2

return 2*E_total, r[0:l+1], pot[0:l+1]

</source>

== Evolução da Distribuição no Tempo ==

Para simular a evolução temporal da distribuição inicial escolhida, precisamos de um método de integração numérica no tempo. O método escolhido deverá partir do estado da distribuição e da função de derivada temporal atribuída para estimar como a evolução ocorre.

<source lang='python'>

from scipy.integrate import solve_ivp #integrador numérico escolhido

def nova_fase(y,i): #função que define a função do espaço de fase em cada instante de tempo, adaptada para o integrador numérico escolhido

fase = np.array([])

for j in range(len(y)):

fase = np.append(fase,y[j][i])

return fase

Fase = [] # Cria uma lista para reservar a evolução do espaço de fase no tempo
T = np.array([]) # Cria uma lista que registra os tempos de marcação do integrador numérico
E = np.array([]) # Cria uma lista que registra a energia da distribuição no tempo
E_pot = [] # Cria uma lista que registra o comportamento da energia potencial ao longo de x no tempo
X_pot = [] # O análogo ao anterior para a posição
tempo = 0.0 # o tempo começa em 0
delta_t = 1 # o intervalo de tempo é 0.1

X = np.random.rand(4000)
V = np.random.rand(4000) - np.ones(4000)*0.5
y = np.append(X,2*V) #define-se a fase inicial

for j in range(50):

sol = solve_ivp(deriv, (tempo,tempo + delta_t), y , rtol = 0.00001) # integrador numérico realiza a integração

T= np.append(T, sol.t)

tempo = tempo + delta_t
for i in range(len(sol.t)): #o integrador possui um número de vetores intermediários ao vetor da solução final sol.y[-1]

y = nova_fase(sol.y,i) # aqui são registradas as funções definidas para cada estado de integração
Fase.append(y)

Ei = Energia_total(y)

E = np.append(E, Ei[0])

E_pot.append(Ei[2])
X_pot.append(Ei[1])

y = nova_fase(sol.y, -1) # A distribuição atual é atualizada para nova integração

</source>

=== Animações ===

Foi utilizado o código abaixo para gerar os gifs da evolução temporal

<source lang='python'>

from matplotlib.animation import FuncAnimation
from IPython import display

# código base para montar gifs da evolução temporal
fig, (ax1, ax2, ax3) = plt.subplots(3,1 ,figsize = (6,14)) #cria uma figura de três linhas

lines = [] #objeto padrão da biblioteca FuncAnimation

def animate(frame):

# As entradas são limpas para evitar superposição de imagens nos frames
ax1.cla()
ax2.cla()
ax3.cla()

# Plot do espaço de fases
ax1.scatter(Fase[frame][0:4000], Fase[frame][4000:8000], color = 'k', s =5)
ax1.scatter(-Fase[frame][0:4000], -Fase[frame][4000:8000], color = 'k', s = 5)
ax1.set_xlabel('$posição$')
ax1.set_ylabel('$velocidade$')
ax1.set_xlim(-4,4)
ax1.set_ylim(-2.5,2.5)

# Plot da distribuição de Potencial
ax2.plot(X_pot[frame],E_pot[frame], color = 'grey')
ax2.plot(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'grey')
ax2.scatter(X_pot[frame],E_pot[frame], color = 'b', s =5)
ax2.scatter(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'b', s =5)
ax2.set_ylabel("potencial")
ax2.set_xlabel("posição")

# Plot da energia total no tempo
ax3.plot(T[0:frame],E[0:frame], color = 'k')
ax3.set_ylabel("Energia")
ax3.set_xlabel("tempo")
ax3.set_xlim(0, T[-1])
ax3.set_ylim(E[0]- 10, E[0] + 10)

return lines

# Execução do programa e geração do gif inline
anima = FuncAnimation(fig, animate,frames=len(T), interval=20, blit=True)
video = anima.to_html5_video()
html = display.HTML(video)
display.display(html)

plt.close()

</source>

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-21T15:06:40Z

Renanrs:

Em construção

'''Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara'''

'''Objetivo:'''

= O modelo =

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível [1] e também para sistemas auto-gravitantes [2]. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

== Descrição do Modelo ==

=== Disposição espacial ===

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

ADD IMAGEM

=== Ausência de Efeitos de Colisão ===

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

ADD IMAGEM

:<math>
\begin{cases} mv_1 - mv_2 = -mv_1' + mv_2' \\ \\ \frac{v_1' + v_2'}{v_1 + v_2} = 1 \end{cases}
</math>

O sistema de equações resulta em: <math>v_1' = v_2 \ \ \ ; \ \ \ v_2' = v_1</math>

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x.

=== Interação entre as partículas ===

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

ADD IMAGEM

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial <math>\sigma</math>:
:<math>
\oint \overrightarrow{E}. d\overrightarrow{S} = \frac{Q_{int}}{\epsilon_0}
</math>
:<math>
E.A - (-E.A) = \frac{\sigma . A }{\epsilon_0}
</math>
:<math>
E= \frac{\sigma }{2\epsilon_0}
</math>

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

=== Interferência externa do sistema ===

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

ADD IMAGEM

= Comportamento de uma Partícula no Sistema =

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

:<math>
E_i(x) = -\lambda\Theta[-(x-x_i)] + \lambda\Theta(x-x_i)
</math>

ADD IMAGEM

Em que a constante <math>\lambda</math> pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função <math>\Theta</math> equivale a função degrau de Heaviside:

:<math>
\Theta(x) = \begin{cases} 0, \ \ x<0 \\ 1, \ \ x>0 \end{cases}
</math>

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

ADD IMAGEM

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a <math>(N-1)</math>.

É possível notar que para uma distribuição de partículas-plano do tipo:

:<math>
n(x) = \sum_{i=1}^N \delta(x_i)
</math>

com o índice <math>\ \ ''i''\ </math> sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

:<math>
F_i = [-(N-1) + 2(i-1)]F
</math>
:<math>
F_i = (2i - N-1)F
</math>

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

:<math>
F^{grav.}_{R,i} = (2i - N-1)F
</math>

:<math>
F^{el.}_{R,i} = (2i - N-1)F - \kappa x_i
</math>

= Conservação da Energia do Sistema =

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

== Energia Cinética ==

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

:<math>
En_{cin} = \frac{m}{2}\sum_{i=1}^N \dot{x_i}^2
</math>

== Energia Potencial Restauradora ==

Observada no sistema eletricamente carregado citado, pode ser descrita como:

:<math>
En_{pot} = \frac{\kappa}{2}\sum_{i=1}^N x_i^2
</math>

== Energia Potencial Intrínseca ==

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições <math>x=L</math> e <math>x= -L</math> (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são <math>V(-L) = V(L) = 0</math>.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, <math>x_0</math> e <math>x_{N+1}</math>, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

:<math>
\Delta V_i = (x_{i} - x_{i - 1})E_i
</math>

com <math>E_i</math> sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

:<math>
E_i = (2i - N - 1)\lambda
</math>

ADD IMAGEM

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

:<math>
En_k = \alpha \sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante <math>\alpha</math> que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

:<math>
En_{int} = \alpha \sum_{k = 1}^{N}\sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

:<math>
E_{total} = En_{cin} + En_{pot} + En_{int} = cte.
</math>

= Modelagem Computacional =

== Criando a Distribuição do Sistema ==

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas <math>(x_i,y_i)</math> existirá uma partícula espelhada
com coordenadas <math>(-x_i,-y_i)</math>.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

<source lang='python'>

# Importando as bibliotecas necessárias
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[] #criamos uma lista vazia

N = 1000 # Determinamos o número de partículas utilizado

for i in range(int(N)):

A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados

A = np.array(A)
X= np.array([]) #Arrays usados para exemplificar
V= np.array([]) #distribuição inicial

for i in range (len(A)): #simetrização da distribuição
X=np.append(X,A[i][0])
V=np.append(V,A[i][1])
V=np.append(V,-A[i][1])
X=np.append(X,-A[i][0])

plt.xlabel('posição') #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

</source>

ADD IMAGEM - <matplotlib.collections.PathCollection at 0x7f4ab8366ad0>

== Função de Indexação ==

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo <math>x</math>.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e <math>2N</math> linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

<source lang='python'>

teste = np.array([2,6,3,9,0,1,4]) #criado um vetor de teste para que seus valores sejam indexados

def indexador(arr): #define-se a função

ordenador = [] #é criada uma lista vazia

for i in range(len(arr)):

ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual

ordenador.sort() # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento

for i in range(len(arr)):

ordenador[i].append(i) # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência

for i in range(len(arr)):

a = ordenador[i][0]
b = ordenador[i][1] # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
ordenador[i][0] = b
ordenador[i][1] = a

ordenador.sort() # reordenar os elementos de acordo com a ordem antiga

lista = []

for i in range(len(arr)):

lista.append(ordenador[i][2]) #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável

indexador = np.array(lista)

return indexador # e essa ordem é retornada como saída da função

print(indexador(teste))

</source>

== Função de movimento e aceleração das partículas ==

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado (<math>y</math>).

Para as coordenadas espaciais x (<math>y[i], \ 0\leq i< N </math>) temos que:

:<math>
\dot{x_i} = y_{[i + N]}
</math>

Para as coordenadas de velocidade v (<math>y[i], \ N\leq i< 2N-1 </math>) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias (<math>m = q = \kappa = 1</math>) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave (<math>F = 1/N</math>). O programa será escrito levando em conta o modelo de cargas elétricas:

:<math>
\dot{v}_i = \ddot{x}_i = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + (2i - N - 1)F \Big{]}
</math>

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice <math>1</math> possui <math>N</math> partículas à esquerda (x negativo) e <math>(N-1)</math> partículas a direita (x positivo), o que resulta uma força de interação de <math>F</math> para a direita. A de índice 2 terá força de interação <math>3F</math>, a de índice 3 terá <math>5F</math> e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

:<math>
\ddot{x}_i = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + \frac{y_{[i - N]}}{|y_{[i - N]}|}(2i- 1)F \Big{]}
</math>

<source lang='python'>

def deriv(t, y):

p = len(y)
k = indexador(y[0:int(p/2)])
sinal=abs(y[0:int(p/2)])/y[0:int(p/2)]

dydt = np.zeros(int(p))

dydt[0:int(p/2)] = y[int(p/2):int(p)] #definindo as derivadas temporais de posição

#definindo as derivadas temporais de velocidade
dydt[int(p/2):int(p)] = -y[0:int(p/2)] + ((np.ones(int(p/2)) + 2*k )/(p))*sinal

deriv = dydt[0:int(p)]

return deriv

</source>

== Função de Energia do estado ==

<source lang='python'>

def Energia_total(y): # Para calcular a energia precisamos apenas do vetor de estado y do sistema de partículas utilizado

l = int(len(y)/2) # a variável l representa a quantidade de valores que teremos
pot = np.zeros(l+1) # cria-se um array com l + 1 dimensões para representar as l partículas e os respectivos potenciais

r = abs(np.array(y[0:l])) # cria-se um array r com as coordenadas de posição do estado y
r = np.append(r,4) # adiciona-se uma coordenada para representar o contorno externo do experimento
r.sort() # Faz-se o ordenamento desse array por distância da origem

pot_0 = 0 # O potencial do contorno deve ser nulo

for i in range(l,0,-1): # Loop decrescente que define o potencial de posição para cada partícula

pot[i] = pot_0

pot_0 = pot_0 + (1/2)*(i/l)*abs((r[i] - r[i-1]))

pot[0] = pot_0

v = np.array(y[l:-1]) # Cria-se um array com os valores de velocidade

E_total = 0

for i in range(l-1): #Soma-se as energias de cada partícula

E_total = E_total + pot[i] + 0.5*v[i]**2 + 0.5*r[i]**2

return 2*E_total, r[0:l+1], pot[0:l+1]

</source>

== Evolução da Distribuição no Tempo ==

Para simular a evolução temporal da distribuição inicial escolhida, precisamos de um método de integração numérica no tempo. O método escolhido deverá partir do estado da distribuição e da função de derivada temporal atribuída para estimar como a evolução ocorre.

<source lang='python'>

from scipy.integrate import solve_ivp #integrador numérico escolhido

def nova_fase(y,i): #função que define a função do espaço de fase em cada instante de tempo, adaptada para o integrador numérico escolhido

fase = np.array([])

for j in range(len(y)):

fase = np.append(fase,y[j][i])

return fase

Fase = [] # Cria uma lista para reservar a evolução do espaço de fase no tempo
T = np.array([]) # Cria uma lista que registra os tempos de marcação do integrador numérico
E = np.array([]) # Cria uma lista que registra a energia da distribuição no tempo
E_pot = [] # Cria uma lista que registra o comportamento da energia potencial ao longo de x no tempo
X_pot = [] # O análogo ao anterior para a posição
tempo = 0.0 # o tempo começa em 0
delta_t = 1 # o intervalo de tempo é 0.1

X = np.random.rand(4000)
V = np.random.rand(4000) - np.ones(4000)*0.5
y = np.append(X,2*V) #define-se a fase inicial

for j in range(50):

sol = solve_ivp(deriv, (tempo,tempo + delta_t), y , rtol = 0.00001) # integrador numérico realiza a integração

T= np.append(T, sol.t)

tempo = tempo + delta_t
for i in range(len(sol.t)): #o integrador possui um número de vetores intermediários ao vetor da solução final sol.y[-1]

y = nova_fase(sol.y,i) # aqui são registradas as funções definidas para cada estado de integração
Fase.append(y)

Ei = Energia_total(y)

E = np.append(E, Ei[0])

E_pot.append(Ei[2])
X_pot.append(Ei[1])

y = nova_fase(sol.y, -1) # A distribuição atual é atualizada para nova integração

</source>

=== Animações ===

Foi utilizado o código abaixo para gerar os gifs da evolução temporal

<source lang='python'>

from matplotlib.animation import FuncAnimation
from IPython import display

# código base para montar gifs da evolução temporal
fig, (ax1, ax2, ax3) = plt.subplots(3,1 ,figsize = (6,14)) #cria uma figura de três linhas

lines = [] #objeto padrão da biblioteca FuncAnimation

def animate(frame):

# As entradas são limpas para evitar superposição de imagens nos frames
ax1.cla()
ax2.cla()
ax3.cla()

# Plot do espaço de fases
ax1.scatter(Fase[frame][0:4000], Fase[frame][4000:8000], color = 'k', s =5)
ax1.scatter(-Fase[frame][0:4000], -Fase[frame][4000:8000], color = 'k', s = 5)
ax1.set_xlabel('$posição$')
ax1.set_ylabel('$velocidade$')
ax1.set_xlim(-4,4)
ax1.set_ylim(-2.5,2.5)

# Plot da distribuição de Potencial
ax2.plot(X_pot[frame],E_pot[frame], color = 'grey')
ax2.plot(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'grey')
ax2.scatter(X_pot[frame],E_pot[frame], color = 'b', s =5)
ax2.scatter(-np.flip(X_pot[frame]),np.flip(E_pot[frame]), color = 'b', s =5)
ax2.set_ylabel("potencial")
ax2.set_xlabel("posição")

# Plot da energia total no tempo
ax3.plot(T[0:frame],E[0:frame], color = 'k')
ax3.set_ylabel("Energia")
ax3.set_xlabel("tempo")
ax3.set_xlim(0, T[-1])
ax3.set_ylim(E[0]- 10, E[0] + 10)

return lines

# Execução do programa e geração do gif inline
anima = FuncAnimation(fig, animate,frames=len(T), interval=20, blit=True)
video = anima.to_html5_video()
html = display.HTML(video)
display.display(html)

plt.close()

</source>

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-21T15:00:23Z

Renanrs:

Em construção

'''Grupo: Antônio Carlan, Gabriela Pereira, Renan Soares e Victor Gandara'''

'''Objetivo:'''

= O modelo =

É possível criar um modelo simplificado de interação de partículas imaginando que todas elas são compostas por planos de mesma massa e sem espessura, paralelos entre si e com liberdade de movimento relevante no eixo perpendicular às superfícies das mesmas.

Tal modelo pode ser utilizado para estudar a dinâmica de um número elevado de partículas e conta com algumas particularidades simplificadoras de complexidade numérica. Alguns estudos utilizam esse modelo para estudar o comportamento de feixes carregados com a menor dimensionalidade possível [1] e também para sistemas auto-gravitantes [2]. As duas aplicações citadas serão as demonstradas a seguir mas o modelo é versátil para diferentes condições iniciais e de contorno.

== Descrição do Modelo ==

=== Disposição espacial ===

Um número "N" de partículas-plano de mesma massa é disposto ao longo do eixo x com uma função de distribuição n(x) e com todas as partículas paralelas ao plano yz. No caso do modelo auto-gravitante, movimentos ao longo do eixo z não têm relevância prática mas para o modelo de um feixe de planos é definido que todos os planos se movem na direção z com mesma velocidade.

A consideração de movimento ao longo do eixo z para o modelo de feixe é importante para a análise do comportamento do espaço de fases ao longo do comprimento do acelerador de partículas. Aqui, podemos imaginar estudar casos em que o feixe gerado inicialmente no espaço de fase terá um comportamento constante no tempo para cada secção transversal em xy para uma coordenada em z pré-estabelecida.

ADD IMAGEM

=== Ausência de Efeitos de Colisão ===

O modelo leva em consideração que todas os planos possuem a mesma massa e que as colisões não levam a perda de energia no sistema.

Dessa forma podemos prever o que ocorre quando a posição de dois planos coincide pois há conservação no momento linear em x e o coeficiente de restituição (razão entre as velocidades de afastamento e aproximação dos corpos) é igual a 1.

ADD IMAGEM

:<math>
\begin{cases} mv_1 - mv_2 = -mv_1' + mv_2' \\ \\ \frac{v_1' + v_2'}{v_1 + v_2} = 1 \end{cases}
</math>

O sistema de equações resulta em: <math>v_1' = v_2 \ \ \ ; \ \ \ v_2' = v_1</math>

Dessa forma, é irrelevante ponderar se houve de fato uma colisão com troca de velocidades ou se as partículas simplesmente cruzam sem interferir de forma relevante. A última consideração é a levada em conta no ponto de vista computacional pois exclui a necessidade de atualizar a velocidade das partículas que passaram pela mesma posição no eixo x.

=== Interação entre as partículas ===

Para analisar a interação entre as partículas é preciso entender como uma única partícula afeta seu entorno.

No caso das partículas-plano carregadas eletricamente podemos considerar os efeitos de campo elétrico de um isolante plano infinito:

ADD IMAGEM

Pela lei de Gauss, para uma dada distribuição de carga com densidade superficial <math>\sigma</math>:
:<math>
\oint \overrightarrow{E}. d\overrightarrow{S} = \frac{Q_{int}}{\epsilon_0}
</math>
:<math>
E.A - (-E.A) = \frac{\sigma . A }{\epsilon_0}
</math>
:<math>
E= \frac{\sigma }{2\epsilon_0}
</math>

O mesmo pode ser feito para um plano com densidade de massa distribuída uniformemente, para o caso do sistema auto-gravitante, com mudança na constante física multiplicativa e no sinal.

O fato importante a ser destacado aqui é a total independência do campo gerado pela partícula e a distância. Dessa forma, a interação entre as partículas em qualquer sistema que seguir esse modelo não terá dependência da distância entre elas e sim da posição relativa entre as partículas consideradas.

=== Interferência externa do sistema ===

Para o caso dos planos no sistema auto-gravitante não é necessária ação externa alguma para manter o equilíbrio em um sistema dinâmico estável. Sempre é levado em conta que todas as partículas têm velocidade em x inferior a velocidade de escape do restante do sistema e, por conta disso e das forças sempre atrativas de interação, o sistema não expande indefinidamente.

Já no caso de planos carregados eletricamente as forças de interação entre os planos serão sempre repulsivas e a energia total do sistema é positiva. Sem uma interação externa, mesmo um sistema composto por um número limitado de partículas de fraca interação expandiria indefinidamente. Por causa disso, nesse caso é considerado que há uma força colimadora no feixe de planos proporcional a distância da origem. A origem dessa força externa vem de modelos de maior dimensionalidade de aceleradores de partícula que contam com um campo magnético externo capaz de colimar um feixe de cargas radialmente, por não existir um paralelo direto para o modelo unidimensional consideramos apenas uma força restauradora que impede a expansão indefinida do sistema como um todo. Também é possível ajustar condições de contorno simulando as paredes externas do acelerador como placas condutoras infinitas e com potencial definido. Esse tipo de consideração pode ser feita para sistemas que não são simétricos e impede que as cargas se afastem demasiadamente da origem.

ADD IMAGEM

= Comportamento de uma Partícula no Sistema =

Graças a caracterização das partículas temos que o campo gerado por uma partícula-plano pode ser caracterizado (para ambos os sistemas citados) da seguinte forma:

:<math>
E_i(x) = -\lambda\Theta[-(x-x_i)] + \lambda\Theta(x-x_i)
</math>

ADD IMAGEM

Em que a constante <math>\lambda</math> pode ser tanto positiva (cargas positivas) quanto negativa (cargas negativas ou interação gravitacional) e a função <math>\Theta</math> equivale a função degrau de Heaviside:

:<math>
\Theta(x) = \begin{cases} 0, \ \ x<0 \\ 1, \ \ x>0 \end{cases}
</math>

Dessa forma, para exemplificação, um sistema de 5 partículas-plano teria uma interação da seguinte forma:

ADD IMAGEM

Aqui, F seria a força de interação entre duas partículas no caso de planos carregados. É possível notar que a força de interação varia de duas unidades dessa constante para cada partícula contada em uma direção e que a força máxima sentida (pelas partículas-plano localizadas nas extremidades) é proporcional a <math>(N-1)</math>.

É possível notar que para uma distribuição de partículas-plano do tipo:

:<math>
n(x) = \sum_{i=1}^N \delta(x_i)
</math>

com o índice <math>\ \ ''i''\ </math> sendo a ordem da partícula contada do sentido negativo para o positivo do eixo x, teremos uma força resultante de interação para uma partícula que pode ser definida por

:<math>
F_i = [-(N-1) + 2(i-1)]F
</math>
:<math>
F_i = (2i - N-1)F
</math>

Essa força é suficiente para manter o sistema auto-gravitante coeso mas para as partículas carregadas devemos adicionar a força restauradora descrita anteriormente:

:<math>
F^{grav.}_{R,i} = (2i - N-1)F
</math>

:<math>
F^{el.}_{R,i} = (2i - N-1)F - \kappa x_i
</math>

= Conservação da Energia do Sistema =

Ao realizar uma simulação numérica com um grande número de entes que representam um sistema físico é importante que seja feita uma verificação sobre o quão apropriados os resultados da simulação de fato são. Uma das formas mais comuns de verificar se uma simulação numérica está correta na dinâmica molecular de vários corpos e efetuar periodicamente a soma da energia do sistema.

No caso de sistemas como os descritos, não há nenhum incremento de energia proveniente de uma fonte externa. Mesmo para o sistema eletricamente carregado, a força restauradora é potencial e não muda a energia total do sistema dada uma distribuição inicial de posição e momento.

Para todos os efeitos, a energia dos sistemas estudados pode ser dividida em três partes:

== Energia Cinética ==

Todas as partículas da distribuição escolhida terão massa e velocidade, de forma que a energia cinética do sistema será dada por:

:<math>
En_{cin} = \frac{m}{2}\sum_{i=1}^N \dot{x_i}^2
</math>

== Energia Potencial Restauradora ==

Observada no sistema eletricamente carregado citado, pode ser descrita como:

:<math>
En_{pot} = \frac{\kappa}{2}\sum_{i=1}^N x_i^2
</math>

== Energia Potencial Intrínseca ==

Energia contida na interação entre as partículas. A energia de uma partícula depende do potencial na sua posição da sua carga (ou massa).

Como potenciais não são absolutos, é estabelecido um Potencial nulo de referência nas condições de contorno do sistema. Como pode ser visto no gif com as partículas-plano carregadas em azul, são colocadas placas condutoras planas nas posições <math>x=L</math> e <math>x= -L</math> (caso haja simetria no espaço de fase do sistema) e é estabelecido que os potenciais nessas posições são <math>V(-L) = V(L) = 0</math>.

Uma maneira mais simples de determinar o potencial das N partículas é considerar que há nas posições de contorno, <math>x_0</math> e <math>x_{N+1}</math>, partículas não interagentes.

A diferença de potencial entre duas posições vizinhas pode ser determinada por:

:<math>
\Delta V_i = (x_{i} - x_{i - 1})E_i
</math>

com <math>E_i</math> sendo o módulo do campo entre a i-ésima e a (i - 1)-ésima partícula, dado por:

:<math>
E_i = (2i - N - 1)\lambda
</math>

ADD IMAGEM

O que resulta na energia de uma partícula específica (k-ésima) descrita por:

:<math>
En_k = \alpha \sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

Aqui foram aglutinadas as possíveis constantes multiplicativas em uma constante <math>\alpha</math> que pode ser escolhida para se adequar ao sistema estudado e a energia intrínseca do sistema é dada pela soma da energia de interação de cada uma das partículas. Como o sistema é dividido em intervalos com variação de potencial discretos, podemos escrever a energia intrínseca da seguinte forma:

:<math>
En_{int} = \alpha \sum_{k = 1}^{N}\sum_{i = 1}^{k}(x_i - x_{i - 1})(2i - N - 1)
</math>

A soma dessas três energias deve se manter constante (ou com pouca variação relativa significativa, pelo menos) durante a simulação:

:<math>
E_{total} = En_{cin} + En_{pot} + En_{int} = cte.
</math>

= Modelagem Computacional =

== Criando a Distribuição do Sistema ==

O primeiro passo para criar a simulação das N partículas é definir as coordenadas fundamentais das partículas em si (posição e velocidade).

A simulação que será utilizada de exemplo é simétrica no espaço de fase, ou seja, para cada partícula com coordenadas <math>(x_i,y_i)</math> existirá uma partícula espelhada
com coordenadas <math>(-x_i,-y_i)</math>.

Esse artifício é utilizado para que o centro de massa mantenha-se na origem e que o momento linear do sistema em x seja sempre nulo. Com isso é garantido que haverá sempre, também, a simetria par entre as energias das partículas espelhadas e a simetria ímpar entre as forças resultantes.

Com isso garantido, basta fazer o cálculo numérico de uma metade da distribuição levando em conta a presença da outra metade, porém não tornando necessário que os cálculos sejam repetidos pois o comportamento de cada partícula espelhada torna-se perfeitamente previsível.

<source lang='python'>

# Importando as bibliotecas necessárias
from math import *
import matplotlib.pyplot as plt
import numpy as np

A=[] #criamos uma lista vazia

N = 1000 # Determinamos o número de partículas utilizado

for i in range(int(N)):

A.append([np.random.rand() ,2*(np.random.rand() - 1/2)]) #valores aleatórios de x entre 0 e 1 e de v entre -1 e 1 são adicionados

A = np.array(A)
X= np.array([]) #Arrays usados para exemplificar
V= np.array([]) #distribuição inicial

for i in range (len(A)): #simetrização da distribuição
X=np.append(X,A[i][0])
V=np.append(V,A[i][1])
V=np.append(V,-A[i][1])
X=np.append(X,-A[i][0])

plt.xlabel('posição') #plot da distribuição inicial
plt.ylabel('momento linear')
plt.title('Espaço de fase inicial')
plt.grid(True)

plt.scatter(X,V)

</source>

ADD IMAGEM - <matplotlib.collections.PathCollection at 0x7f4ab8366ad0>

== Função de Indexação ==

Em outras modelagens como a interação de cargas em sistemas bidimensionais ou de massas em tridimensionais é necessário calcular a força de interação entre partículas em função das distâncias. Na modelagem proposta não é necessário calcular distâncias mas sim saber o ordenamento das partículas ao longo do eixo <math>x</math>.

No entanto, é muito comum em simulações que utilizam métodos de integração utilizar vetores que possuam apenas uma coluna e <math>2N</math> linhas, sendo a primeira metade de linhas correspondentes às posições e a outra metade correspondente às velocidades.

Para não precisar atualizar as linhas de posição e velocidade do vetor 2N-dimensional, utiliza-se uma função de indexação que será responsável por indicar qual é a ordem da partícula correspondente a uma determinada linha em um dado estado.

Assim, a posição da partícula é atualizada no tempo e seu índice também, porém a posição dela no vetor de fase nunca se modifica. Esse fato é importante para estudos que pretendem analisar a trajetória de um único corpo no tempo, como é feito nos diagramas de Poincaré, por exemplo.

<source lang='python'>

teste = np.array([2,6,3,9,0,1,4]) #criado um vetor de teste para que seus valores sejam indexados

def indexador(arr): #define-se a função

ordenador = [] #é criada uma lista vazia

for i in range(len(arr)):

ordenador.append([abs(arr[i]),i]) # para cada valor criamos uma lista com o valor absoluto e a ordem atual

ordenador.sort() # é feito o ordenamento da lista 'ordenador' de acordo com o primeiro elemento

for i in range(len(arr)):

ordenador[i].append(i) # é adicionado um terceiro valor para cada elemento de ordenador com o valor da sua nova ordem na sequência

for i in range(len(arr)):

a = ordenador[i][0]
b = ordenador[i][1] # É feita a troca do segundo com o primeiro valor de cada elemento de ordenador para...
ordenador[i][0] = b
ordenador[i][1] = a

ordenador.sort() # reordenar os elementos de acordo com a ordem antiga

lista = []

for i in range(len(arr)):

lista.append(ordenador[i][2]) #por fim guarda-se a ordem correspondente a cada valor da lista colocada como variável

indexador = np.array(lista)

return indexador # e essa ordem é retornada como saída da função

print(indexador(teste))

</source>

== Função de movimento e aceleração das partículas ==

Para simular como um estado evolui no tempo, é necessário que tenhamos as derivadas temporais de todas as coordenadas do vetor de estado (<math>y</math>).

Para as coordenadas espaciais x (<math>y[i], \ 0\leq i< N </math>) temos que:

:<math>
\dot{x_i} = y_{[i + N]}
</math>

Para as coordenadas de velocidade v (<math>y[i], \ N\leq i< 2N-1 </math>) teremos que usar a relação de forças que foi descrita no modelo. Algumas simplificações comuns serão feitas como considerar constantes físicas unitárias (<math>m = q = \kappa = 1</math>) e relações de proporcionalidade com o número de partículas do sistema para mantê-lo com um comportamento mais suave (<math>F = 1/N</math>). O programa será escrito levando em conta o modelo de cargas elétricas:

:<math>
\dot{v}_i = \ddot{x}_i = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + (2i - N - 1)F \Big{]}
</math>

A relação acima pode ser descrita de forma diferente se levarmos em conta que a simulação é simétrica em x e que todas as posições iniciais são positivas.

A partícula que recebe índice <math>1</math> possui <math>N</math> partículas à esquerda (x negativo) e <math>(N-1)</math> partículas a direita (x positivo), o que resulta uma força de interação de <math>F</math> para a direita. A de índice 2 terá força de interação <math>3F</math>, a de índice 3 terá <math>5F</math> e assim por diante.

Antes de reescrever a equação da aceleração deve-se levar em consideração de uma das partículas distribuídas no espaço de fase assumirem coordenadas negativas (caso tenham proximidade da origem e velocidade negativa). Novamente usando o argumento de simetria, sabe-se que caso uma partícula assuma uma posição negativa, a partícula espelhada assumirá a posição simétrica em relação a origem. Dessa forma basta calcular a força de interação sobre a partícula espelhada e inverter o sinal para a partícula original. Levando tudo isso em conta teremos a seguinte equação da aceleração:

:<math>
\ddot{x}_i = \frac{1}{m}\Big{[} -\kappa y_{[i - N]} + \frac{y_{[i - N]}}{|y_{[i - N]}|}(2i- 1)F \Big{]}
</math>

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-21T14:39:16Z

Renanrs:

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-21T14:36:04Z

Renanrs:

Dinâmica Molecular de Partículas-Planares Interagentes

2021-11-21T14:35:12Z

Renanrs:

2021-10-15T13:39:12Z

Renanrs:

Arquivo:FAE visualization.jpg

2021-10-15T13:38:59Z

Renanrs: trabalho elementos finitos 2021/2 fonte: https://en.wikipedia.org/wiki/Finite_element_method

trabalho elementos finitos 2021/2

fonte: https://en.wikipedia.org/wiki/Finite_element_method

Arquivo:Matrizglobal.png

2021-10-15T13:38:26Z

Renanrs: trabalho elementos finitos 2021/2

trabalho elementos finitos 2021/2

Arquivo:Nos.png

2021-10-15T13:38:16Z

Renanrs: trabalho elementos finitos 2021/2

trabalho elementos finitos 2021/2

Arquivo:Contribuição elemento.png

2021-10-15T13:37:55Z

Renanrs: trabalho elementos finitos 2021/2

trabalho elementos finitos 2021/2