Equação de Cahn-Hilliard em 2D

Leonardo Dasso Migottto WORK IN PROGRESS

O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, utilizando a Transformada Rápida de Fourier em uma e (principalmente) em duas dimensões. Será explorado as variações em concentração inicial e seus respectivos padrões formados, dados coeficientes de difusão e largura da superfície fixos.

Esta equação já foi tratada em detalhes por colegas anteriores a mim , e a leitura do trabalho por eles desenvolvido é recomendada para maior entendimento da equação. O foco deste trabalho é explorar a solução numérica para a equação quando tratada em duas dimensões, onde a formação de padrões apresenta resultados mais interessantes. No entanto, a fim de facilitar a implementação e entendimento em duas dimensões, também será exibido uma implementação em uma dimensão.

Equação de Cahn-Hiliiard utilizando Transformada de Fourier

Para encontrar a equação que implementaremos com o uso da Transformada Rápida de Fourier, precisamos encontrar a nossa equação representada no espaço de Fourier. Seguirei a literatura de S. Bulent Biner , onde há um capítulo dedicado a resolver equações de difusão com métodos que utilizam esta transformada. Primeiro, resolveremos em uma dimensão a equação, que segue abaixo:

${\frac {\partial c}{\partial t}}=D{\nabla }^{2}(c^{3}-c-\gamma {\nabla }^{2}c)$

Em uma dimensão, os laplacianos podem ser substituídos pela derivada segunda em relação a $x$, resultando na seguinte equação:

${\frac {\partial c}{\partial t}}=D{\Bigl (}{\frac {\partial ^{2}(c^{3}-c)}{\partial x^{2}}}-\gamma {\frac {\partial ^{4}c}{\partial x^{4}}}{\Bigr )}$

Para solucioná-la numericamente, aplicaremos a Transformada de Fourier à frente em ambos os lados, da maneira descrita abaixo, onde k é o respectivo coeficiente de Fourier):

${\frac {\partial \{c\}_{k}}{\partial t}}=D{\Bigl (}{\Bigl \{}{\frac {\partial ^{2}c^{3}-c}{\partial x^{2}}}{\Bigr \}}_{k}-\gamma {\Bigl \{}{\frac {\partial ^{4}c}{\partial x^{4}}}{\Bigr \}}_{k}{\Bigr )}$

Em seguida, substituimos as derivadas espaciais pela sua equivalente no espaço de Fourier:

${\Bigl \{}{\frac {\partial ^{n}c}{\partial x^{n}}}{\Bigr \}}_{k}\to (ik)^{n}\{c\}_{k}$

Assim, obtemos a seguinte equação:

${\frac {\partial \{c\}_{k}}{\partial t}}=D{\bigl (}-k^{2}(\{c^{3}\}_{k}-\{c\}_{k})-\gamma k^{4}\{c\}_{k}{\bigr )}$

O próximo passo é fazer a derivada à direita quanto ao tempo da seguinte maneira:

${\frac {\partial \{c\}_{k}}{\partial t}}\to {\frac {\{c\}_{k}^{n+1}-\{c\}_{k}^{n}}{\Delta t}}$

Substituindo na equação e reescrevendo-a a fim de isolar $\{c\}_{k}^{n+1}$, obtemos a equação final:

$\{c\}_{k}^{n+1}=\{c\}_{k}^{n}+D\Delta t{\bigl (}-k^{2}(\{c^{3}\}_{k}^{n}-\{c\}_{k}^{n})-\gamma k^{4}\{c\}_{k}^{n}{\bigr )}$

Dado que conhecemos a forma da equação em uma dimensão, podemos encontrar sua equivalente bidimensional com maior facilidade. A única diferença entre as duas equações está no laplaciano, que resultará na derivada no eixo $y$ aparecer também. No entanto, a notação da transformada permanece a mesma, e $k$ representará um vetor com coordenada $(k_{1},k_{2})$ com módulo ${\sqrt {k_{1}^{2}+k_{2}^{2}}}$, onde $k_{1}$ e $k_{2}$ são os coeficientes em $x$ e $y$ respectivamente.

Código em uma dimensão

O código completo estará disponível no final desta página.

Resultados em uma dimensão e discussão

Como já há um trabalho que trata em detalhes a implementação unidimensional e seus resultados, irei comparar aqui ambas implementações. Abaixo, vemos uma animação comparando ambos métodos a partir de uma condição inicial aleatória, utillizando condições de contorno periódicas, com o maior valor de erro destacado no topo. O valor das constantes relevantes são: $\Delta t={\frac {10^{-6}}{22}},\Delta x={\frac {1}{128}},D=1,\gamma ={\frac {3.4}{128}}$.

Como podemos ver, a diferença dos valores entre os resultados obtidos pelo método FTCS e o método das transformadas é minúscula (após poucos instantes o maior módulo da diferença é menor que 0,01). É interessante citar que a escolha de $dx$ como um valor tão específico não é em vão: a função RFFT utilizada no código apresenta maior rapidez de execução quando o vetor utilizado é de tamanho $2^{n},n=1,2,3,...$. Com este valor de $dx$, nosso vetor é composto de 128 elementos, ou $2^{7}$ elementos, fazendo uso dessa vantagem. Dado a escala do sistema na qual ambos métodos foram utilizados, os tempos de processamento são mínimos, e não há merito em compará-los. No entanto, para sistemas maiores, onde o número de operações necessárias aumenta drasticamente, o método das transformadas, sob as condições corretas, será extremamente rápido. E isso será observado com clareza quando o sistema for bidimensional.

Resultados em duas dimensões

Para todos resultados abaixo, foram utilizados as seguintes constantes relevantes: $\Delta t={\frac {1.8}{10^{7}}},\Delta x={\frac {1}{128}},D=1,\gamma =0.01$. Também foram utilizadas condições de contorno periódicas. Como observado anteriormente, o valor de $dx$ escolhido serve para que nosso array de valores possua $10^{8}$ elementos, também fazendo uso da eficiência completa da função RFFT2 (equivalente bidimensional da função RFFT). Segue abaixo uma IMAGEM ? da difusão sob estes parâmetros com concentração média 0.

INSERIR IMAGEM AQUI (NÃO PODE GIF, MUITO PESADO)

Aqui há uma animação de um gráfico com concentração média inicial 0.5 (como se houvesse 75% em uma fase e 25% em outra).

Discussão dos resultados em duas dimensões

Podemos notar que há a formação de "listras" conforme passa o tempo, de modo que estas listras tornam-se mais largas. Também nota-se que, quanto mais tempo se passa, menor é o movimento das concentrações, por elas estarem se aproximando da estabilidade. É intuitivo que, neste caso, o estado de equilíbrio do sistema seja quando há somente dois "blocos", um de cada cor, de modo que as "listras" (dois grandes blocos, neste caso) não tenham como se tornar mais largas.