Equação de Kuramoto–Sivashinsky

Nesta página apresentamos uma simulação de Dinâmica Molecular utilizando o potencial de Lennard-Jones como potencial de interação entre as partículas.

Equação

A equação de Kuramoto-Sivashinsky (KS) foi inicialmente proposta por Yoshiki Kuramoto e Gregory Sivashinsky para a modelagem das instabilidades da frente de uma chama laminar.

$u_{t}+u_{xx}+u_{xxxx}+{\frac {1}{2}}u_{x}^{2}=0$

No entanto, a forma diferencial $v=u_{x}$ da equação foi descoberta independentemente na área de dinâmica de fluidos, onde ficou mais conhecida para a modelagem de filmes finos.

$v_{t}+v_{xx}+v_{xxxx}+vv_{x}=0$

A equação apresenta comportamentos caóticos ao ser estudada em geometrias anelares, ou seja, condições de contorno periódicas; por isso esse será o caso estudado aqui.

Métodos de Solução

Diferenças Finitas (FTCS)

O método mais simples foi aplicado a fim de entender o comportamento da equação. Tentativas de linearização resultaram em uma equação de comportamento muito distinto do descrito na literatura; e portanto, não é de nosso interesse estudá-la. É importante relembrar que a aplicação do método FTCS não é recomendado no caso de não linearidades pois o método tem baixa precisão e pode acumular erros.

Como esperado, a solução numérica teve sua estabilidade atrelada à suas condições iniciais; ou seja, a solução é instável (exceto em casos bem específicos, onde ela se comporta de forma análoga à equação da difusão).

Solução FTCS de condição inicial $\sin {\left({\frac {\pi }{20}}x\right)}^{2}$

Solução FTCS de condição inicial $\sin {\left({\frac {\pi }{10}}x\right)}^{2}$

Método Espectral (FFT)

Utilizando o método fft em conjunto com as propriedades da transformada de fourier, pode-se resolver a EDP utilizando a forma de derivada espectral. Esse método possui alta precisão e é capaz de incorporar o comportamento ondulatório de soluções, algo característico da equação de Kuramoto-Sivashinsky.

A derivada espectral é caracterizada pela transformada de fourier de uma derivada, o que facilita a resolução de equações diferenciais parciais:

${\mathcal {F}}(v_{x})=ik{\mathcal {F}}(v)\longrightarrow v_{x}={\mathcal {F^{-1}}}\left(ik{\mathcal {F}}(v)\right)$.

Essa propriedade permite que a resolução da parte espacial da equação de Kuramoto-Sivashinsky diferencial seja dada por

$v_{t}={\mathcal {F^{-1}}}(k^{2}{\hat {v}})-{\mathcal {F^{-1}}}(k^{4}{\hat {v}})-{\frac {v{\mathcal {F^{-1}}}(ik{\hat {v}})}{2}}$,

onde ${\hat {v}}\equiv {\mathcal {F}}(v)$. A parte temporal também pode ser resolvida da mesma forma, mas o método de Runge-Kutta (RK45) foi aplicado em seu lugar por sua maior simplicidade e velocidade de processamento.

Soluções

Através do método espectral foi possível estudar o comportamento da Equação de Kuramoto-Sivashinsky em sua forma diferencial. Para uma melhor interpretação das soluções, os dados são dispostos em dois modelos de gráfico: um gif em que a passagem de tempo é representada pela mudança de frames, e um gráfico com escala de cores; onde os eixos representam o tempo e o espaço.

Solução espectral de condição inicial $\sin {\left({\frac {2\pi }{3}}x\right)}^{2}$

Mapa temporal de condição inicial $\sin {\left({\frac {2\pi }{3}}x\right)}^{2}$

Ao partir de uma condição inicial senoidal, a equação se mantém estável por longos períodos de tempo, o que atrasa a chegada de seu comportamento mais característico: a formação de faixas. Esse comportamento é mais perceptível ao se analisar o mapa temporal da solução, deixando clara a tendência da equação de concentrar mínimos e máximos próximos, gerando picos e vales que se relacionam de forma complexa. No entanto, é fácil perceber que o comportamento simulado não é caótico e sim quase periódico. Isso está atrelado a um fator muito importante para o estudo da equação, o tamanho da periodicidade simulada, ou seja, seu domínio (L).

Caos gerado pelo Domínio (L)

É esperado que a solução para a equação dependa das condições iniciais escolhidas, mas uma relação entre a solução e o tamanho do domínio estudado é algo bem mais surpreendente. A equação KS apresenta esse comportamento inesperado e por isso é objeto de estudo até mesmo em teorias de caos. Sua forma simples apresenta três elementos importantes:

Um termo crescente, com funcionamento inverso à equação do calor:

$-v_{xx}$;

um termo dissipador, responsável pelo amortecimento de picos e vales estreitos:

$-v_{xxxx}$;

e um termo não linear, provavelmente responsável pelo comportamento ondulatório da equação, transferindo energia entre os estados crescentes e amortecidos:

$-vv_{x}$.

A soma desses comportamentos faz com que em domínios periódicos pequenos, a equação consiga se manter estável, atingindo resultados constantes ou ondas planas. Isso pode ser identificado nas soluções abaixo, que apesar de possuírem a mesma condição inicial ( $\sin {\left({\frac {2\pi }{L}}x\right)}^{2}$ ), são mapeadas em domínios de tamanhos diferentes.

Solução espectral para $L=1$

Mapa temporal para $L=1$

Solução espectral para $L=2$

Mapa temporal para $L=2$

Solução espectral para $L=5$

Mapa temporal para $L=5$

É notável a diferença entre as soluções encontradas, mas alguns pontos merecem destaque. Primeiramente existe uma forte relação entre a mudança de fase da equação e seu domínio, ou seja, domínios pequenos são atraídos mais rapidamente para a solução final; enquanto domínios maiores tendem a se manter estáveis por mais tempo. Por outro lado, a solução final tem forma mais simples para domínios pequenos, enquanto domínios maiores tendem a apresentar padrões mais intrincados.

A fim de estudar o caos gerado pela equação, é necessário que pequenas variações (ruído) seja inserido. Para tanto, a condição inicial se manteve a mesma já estudada, mas somada à valores gerados aleatoriamente a partir do intervalo $[-0.01,0.01]$ para cada ponto do domínio.

Mapa temporal para $L=1$ sem ruído

Mapa temporal para $L=2$ sem ruído

Mapa temporal para $L=5$ sem ruído

Mapa temporal para $L=1$ com ruído

Mapa temporal para $L=2$ com ruído

Mapa temporal para $L=5$ com ruído

Ao analisar as soluções com ruído é possível perceber a sensibilidade de cada solução às condições iniciais. A solução de domínio $L=1$ se manteve quase inalterada, atingindo um resultado muito similar ao da solução original. Já o domínio $L=2$ apresentou uma sensibilidade muito maior às pequenas variações do ruído; inicialmente a solução aparenta seguir o mesmo caminho, mas aos poucos deformidades começam a surgir e a solução perde suas simetrias e periodicidades. Por fim, a solução de domínio $L=5$ teve seu comportamento quase completamente alterado; sua maior estabilidade deixa de estar presente e seus padrões não apresentam qualquer simetria.

Implementação do Código em Python

#Importação de bibliotecas
from scipy.fft import fft, ifft, fftshift
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import glob

#Parâmetros Gerais
L = 3
N = L*100
t_max = 60000

dx = L/N

X = np.arange(0,L,dx)
K = 4*np.pi/L * np.arange(-L/2,L/2,dx)
K = fftshift(K)

#Forma diferencial da Equação KS
def KS(t,F,K=K):
  fourier = fft(F)
  return ifft(K**2 * fourier) - ifft(K**4 * fourier) - F/2 * ifft(K * 1j* fourier)

#Forma integral da Equação KS
def IKS(t,F,K=K):
  fourier = fft(F)
  return ifft(K**2 * fourier) - ifft(K**4 * fourier) - (ifft(K * 1j* fourier))**2

#Declaração das condições iniciais
F0 = np.sin(2*np.pi*X/L)**2 + 0j
t = 0
tempo = 0
plot = 1

V = []
T = []


#Solução por RK45
sol=integrate.RK45(KS,0,F0,t_max)


while t < t_max:
  passo = sol.step_size
  if passo == None:
    passo = 0
  tempo += passo
  T.append(tempo)
  V.append(sol.y)

  if t%1000==0:
    fig = plt.figure(figsize=(5,5))
    plt.plot(X,sol.y.real,label='tempo:%.2f'% (tempo))
    plt.xlabel('X')
    plt.ylabel('v(x,t)')
    plt.legend(loc='upper right')
    plt.ylim(-5,5)
    fig.savefig('./grafs/%c.PNG'% (plot))
    plt.close()
    plot+=1

  sol.step()
  t +=1


# Criação do GIF
frames = (Image.open(f) for f in sorted(glob.glob("./grafs/*.PNG")))
frame_one = next(frames)
frame_one.save("animado.gif", format="GIF", append_images=frames, save_all=True, duration=150, loop=0)


#Criação mapa de cores real
fig = plt.figure(figsize=(5,5))
graf = plt.pcolormesh(X,T,np.array(V).real,cmap="inferno",vmin=-5,vmax=5)
plt.xlabel('x')
plt.ylabel('t')
cbar = fig.colorbar(graf)
cbar.ax.set_title('v(x,t)')
fig.savefig("mapa_temporal.png")
plt.close()