Problema de Fermi-Pasta-Ulam

Grupo: Augusto M Giani e Henrique Padovani

O objetivo deste trabalho é replicar os resultados do problema proposto por Fermi-Pasta-Ulam em 1953 sobre sistemas dinâmicos não lineares. As análises serão sobre a solução dos modos de vibração comparados à solução analítica para poucas massas e também sobre a energia do sistema para os modos de oscilação, enquanto o sistema evolui no tempo.

O Problema

O Problema proposto constitui-se de simulações em uma rede de partículas ligadas entre si através de molas que obedecem a Lei de Hooke com uma correção não-linear quadrática ou cúbica

• 

  Figura 1. Problema de Fermi-Pasta-Ulam, molas acopladas.

A lei de forças que rege o comportamento deste sistema é:

$F=-k\Delta x-\alpha k\Delta x^{2}-\beta k\Delta x^{3}$.

Onde $\Delta x$ e a deformação a cada 2 massas acopladas ($x_{i+1}-x_{i}$), $k$ é a constante elástica da mola, $\alpha$ é um parâmetro de deformação arbitrário que controla a correção não linear quadrática e $\beta$ é o parâmetro que controla a correção cúbica. Importante ressaltar que se $\alpha$ é possuir assumir um valor não nulo, real, $\beta$ é igual a zero no nosso sistema, ou vice-versa. Não estamos analisando correções quadráticas somadas com correções cúbicas neste trabalho.

Motivação: O paradoxo Enrico Fermi, John R. Pasta, Stanislaw M. Ulam e Mary Tsingou

A premissa inicial do paradoxo de FPUT consiste no Teorema da Equipartição de Energia. O sistema consistia em uma corrente de partículas, com as extremidades fixas, que interagiam entre seus vizinhos somente com forças elásticas (as forças teriam um termo linear como a Força de Hooke e mais um termo não-linear, podendo ser quadrático ou cúbico). Era esperado que a energia total fosse distribuída igualmente entre as partículas. No caso em questão, a distribuição de energia entre as partículas pode ser descrita através dos seus modos normais de vibração.

A análise do problema gerou um paradoxo que começaria a ser respondido somente 10 anos depois, o que ajudou no desenvolvimento das teorias de sólitons e do caos. Pretendia-se observar a distribuição uniforme de energia entre os diversos modos normais de vibração com o passar do tempo (ao longo das iterações da simulação computacional). Isso significaria que o sistema alcançou um equilíbrio térmico e seria uma exemplificação computacional do Teorema de Equipartição de Energia. Caso as forças entre as partículas fossem estritamente lineares, isso não ocorreria, pois a energia alocada em cada modo não conseguiria acessar outros modos. Imaginava-se que uma componente não-linear na força tornaria acessível qualquer modo de vibração, porém não foi o observado.

A princípio, foi observada a tendência do sistema de distribuir a energia. O primeiro modo de vibração antes estimulado, perdeu energia ao longo do tempo, a qual começou a se alocar nos modos de energia mais baixos. Entretanto, por um descuido, deixaram a simulação decorrer por um tempo maior do que era planejado. Ao retornar ao laboratório para corrigir tal erro, se depararam com um resultado inesperado. A energia, que supostamente deveria estar igualmente partilhada entre os modos de vibração, estava quase completamente alocada no primeiro modo de vibração. De fato, somente 3% da energia não estava presente no primeiro modo. Devido a esta observação, deixaram a simulação correr por ainda mais tempo. Notaram, então, que existia um ciclo, no qual a energia saía do primeiro modo de vibração, começava a se distribuir nos modos mais baixos, para, por fim, voltar quase que inteiramente para o primeiro modo de vibração. Contudo, em 2015, relatou-se que o sistema FPUT poderia atingir equipartição de energia pelo menos entre modos normais livres (interação entre três fônons).

Discretização

A discretização deste problema gira em torno de abrir a equação das forças, e com o termo de aceleração, iterar o movimento das partículas a partir disso . Partimos do problema com correção quadrática, ou seja, $\beta =0$. Partindo de:

$F=-k\Delta x-\alpha k\Delta x^{2}$,

substituímos pelas variáveis discretas:

$m{\ddot {x_{j}}}=-k\left((x_{j+1}-x_{j})-(x_{j}-x_{j-1})\right)-\alpha k\left((x_{j+1}-x_{j})-(x_{j}-x_{j-1})\right)$,

Chegamos em:

$m{\ddot {x_{j}}}=k\left(x_{j+1}-2x_{j}+x_{j-1}\right)\left[1+\alpha \left(x_{j+1}-x_{j-1}\right)\right]$,

Em que ${\ddot {x_{j}}}$ é a aceleração da j-ésima partícula, com ela conseguimos integrar o movimento das partículas.

A Energia do sistema é calculada para cada ciclo de oscilação, porém para obtermos os resultados dos modos de oscilação e compararmos com os estudos atuais e o original de de FPUT, calculamos a energia dos primeiros modos de vibração da corda para demonstrar o comportamento visivelmente periódico destas energias. É possível calcular a energia dos modos de vibração através da equação:

$E_{total}=\sum _{i=2}^{N}{\frac {\omega _{i}^{2}A_{i}^{2}+{\dot {A_{i}^{2}}}}{2}}$,

onde $A_{i}$ corresponde à:

$A_{i}=C\sin {\left(i\theta \right)}$

Conseguimos resolver $C$ com a equação geral de $x(t)$:

$x_{i}(t)=\sum _{m=1}^{N}C_{m}\sin {\left({\frac {im\pi }{N+1}}\right)}\cos {\left(wmt\right)}$,

e as freqûencias $\omega _{i}$ são dadas por: $\omega _{m}=2\omega _{0}\sin \left({\frac {m\pi }{2(N+1)}}\right)$.

$N$ é o número de partículas e $\omega _{0}={\frac {k}{massa}}$

Resultados

O sistema foi iniciado com o modo normal de oscilação 1 (nenhum nodo), com velocidades iniciais igual a zero de cada partícula. O movimento começa apenas pelas forças já presentes entre cada partículas, pelas equações apresentadas acima.

legenda.

NOSSO GIF DEU MAIS DE 7 MEGA -----------------> VER COMO REDUZIR O GIF JÁ NO CÓDIGO!!!! ELE TA BEM RUIM ASSIM

legenda.

TEM QUE MUDAR AS FIGURAS DE BAIXO:

• 

  legenda.
• 

  legenda.

CASO PRECISE: ${\begin{cases}X(x=0,y)=X{0}=.....\\Y(x=L,y)=Y(x,y=0)=Y(x,y=L)=.....\\\end{cases}}$

Discussões

Como o intuito era replicar os resultados através da simulação com dinâmica molecular obtivemos resultados muito parecidos comparando estudos já realizados sobre este problema (Citar estes artigos)

Implementação

Usamos XX partículas, com modo de oscilação YY

### Exemplo da iteração do movimento utilizando forcacom correção quadrática ###
### código em python

def aceleracao(pos,alpha,k,massa):
    size = len(pos)
    acel = [0.0 for i in range(size)]
    for i in range(1,size-1):
        acel[i] = (k/massa) * ((pos[i+1] + pos[i-1] - 2*pos[i]) * ( 1.0 + alpha*(pos[i+1]-pos[i-1]) ) )
    return acel

def velocidade(velo, acel, dt):
    size = len(velo)
    new_velo = [0.0 for i in range(size)]
    for i in range(size):
        new_velo[i] = velo[i] + 0.5*acel[i]*dt
    return new_velo

def posicao(pos, velo, dt):
    size = len(posY)
    new_posY = [0.0 for i in range(size)]
    for i in range(size):
        new_posY[i] = posY[i] + new_veloY[i]*dt
    return new_posY

N = número de partículas
dt = 0.2
x = np.linspace(0, x_final, dt)
pos = np.sin( 2*x*pi / (N*dt))

while t < tmax: # Loop temporal
     
    plt.scatter(x,pos) # plotagem dos gráficos
    acel = aceleracao(pos_old,alpha,k,massa)
    velo = velocidade(veloY, new_aceY,dt)
    pos  = posicao(posY_old, veloY, dt)
  
  pos_old = pos.copy()  #sem ".copy()" o python usa o mesmo endereço de memória para 2 variáveis
  t = t + td

gera_gif()

Falar um pouco mais da implementação?

Link para Códigos

Fizemos no ambiente Colab em .ipynb, segue link do github:[1]