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_{i}=-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 (aqui considerado, $\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$ 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_{i}=-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})^{2}-(x_{j}-x_{j-1})^{2}\right)$,

Sendo que as partículas de índice zero e N estão fixas, 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:

em que $a$ é o vetor das posições projetado num vetor de um seno com a frequência do modo que queremos plotar:

e:

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

Resultados

O sistema foi iniciado com o modo normal de oscilação 1 (seno com frequência $\omega =$), 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.

Segue abaixo a evolução do sistema ao longo do tempo, apenas para gerar o .gif considerando:

${\begin{cases}\alpha =1,2\\k=0,95\\m=1,05\\N=30\\t_{max}=4000\\\end{cases}}$

Simulação FPU para o problema proposto.

Já os resultados das energias, consideramos todas constantes iguais à 1, mais de 100 partículas e o número de iterações está em cada imagem. procedemos de duas formas. Calculamos a porcentagem de cada modo normal de oscilação pela fft, o que representa um comportamento muito similar às energias por modo:

Porcentagem de cada modo vibracional ao longo do tempo. A escala do tempo está reduzida em 20 vezes.

E também calculamos pela equação de energia citada acima, obtivemos os seguintes comportamentos:

• 

  Energias por modo de oscilação. A escala do tempo está aumentada em 2 vezes.
• 

  Energias por modo de oscilação com fft. A escala do tempo está aumentada em 5 vezes.

Discussões

Como o intuito era replicar os resultados e o comportamento das energias através da simulação com dinâmica molecular, obtivemos resultados muito parecidos comparando estudos já realizados sobre este problema .

A questão de simular a energia usando a porcentagens dos senos também faz sentido pensando quanta contribuição o modo normal de oscilação tem no movimento atual, como pudemos ver que o comportamento foi bem similar ao das energias calculadas.

Implementação

Implementamos a iteração do movimento das partículas por Velocity-Verlet.

A questão de variar os parâmetros do problema, como já mencionado previamente nesta wiki, e também apresentado nos gráficos e no gif, vimos que não interferia muito no caráter dos resultados (e o quanto isso poderia influenciar no paradoxo de equipartição de energia), porém o fizemos para gerar imagens melhores e o gif com menos de 2 MB.

### 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)
    calcula_energias(pos)
  
  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()

Link para Códigos

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