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 dimensões da superfície fixos.

Cahn-Hilliard bidimensional com condições de contorno periódicas.

Esta equação já foi tratada em detalhes por colegas anteriores a mim , e a leitura do trabalho por eles desenvolvido é recomendada para melhor 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.

Todos códigos desenvolvidos para duas dimensões foram feitos a fim de que qualquer pessoa, com Python e as bibliotecas necessárias instaladas, possa executá-los sem nenhum conhecimento prévio, e estes códigos podem ser encontrados no final desta página. As bibliotecas utilizadas que requerem instalação são: Numpy, Scipy e Matplotlib, e o renderizador de vídeo FFMPEG foi utilizado, porém na maioria dos computadores modernos executando Windows, na ausência deste, outro renderizador será utilizado sem necessidade de alterar o código. Porém, caso seja necessário instalar o renderizador, este link possui um tutorial simples em inglês para Windows: [1]

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 desse tipo. Primeiro, resolveremos a equação em uma dimensã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, encontramos sua equivalente bidimensional com 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 equação final é escrita da mesma forma, de modo que $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 está disponível no final desta página. Abaixo há o excerto da funcão, em Python, que calcula o instante seguinte do nosso sistema, utilizando as funções RFFT e IRFFT do pacote Scipy.

def cahnfourier1d(cc, k2, k4):
    cct = rfft(cc)
    cct3 = rfft(cc**3 - cc)
    cct = cct + difd*dt*(-k2*(cct3) - k4*cct)
    cc = irfft(cct)
    return cc

No código, k2 e k4 são os coeficientes elevados às respectivas potências, porém k4 está multiplicado pelo valor de gamma. Calculando estas constantes préviamente é possível evitar cálculos reduntantes durante o código. É importante destacar que ambas funções funcionam utilizando o vetor completo de pontos para fazer os cálculos, eliminando qualquer necessidade de iteração sobre o vetor de valores (uma diferença notável do método FTCS).

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 alguns instantes 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 {1}{22\cdot 10^{6}}},\Delta x={\frac {1}{128}},D=1,\gamma ={\Bigl (}{\frac {3.4}{128}}{\Bigr )}^{2},x_{max}=1$.

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 aproxima-se de 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.

A solução estacionária do sistema bidimensional, que será avaliada a posteriormente, também nos permite inferir a forma da solução estácionaria do sistema unidimensional dado um tempo infinito: $c(x)=sin{\Bigl (}{\frac {2\pi x}{L}}+\phi {\Bigr )}$, onde $\phi$ é uma constante de fase a fim de compensar pelas condições de contorno periódicas.

Código em Duas Dimensões

Todos códigos estão disponíveis na íntegra no final desta página. A fim de poupar tempo, quatro programas foram criados:

- Um programa que gera, a partir de valores aleatórios parametrizados, valores até um certo instante de tempo, salvando-os em um arquivo.

- Um programa que, lê um arquivo de valores já gerados e gera mais valores até um certo instante de tempo, salvando-os neste mesmo arquivo.

- Um programa que lê um arquivo gerado préviamente, criando uma animação dos valores salvos no arquivo.

- Um programa que lê um arquivo gerado préviamente, criando frames dos valores salvos no arquivo.

Um padrão de arquivo de formato .npy é foi criado, que armazena as informações necessárias para gerar mais valores ou plotar os gráficos, atém de um .txt que serve como uma consulta local dos valores pelo usuário (arquivos .npy não podem ser lidos no bloco de notas). Todos arquivos são salvos em pastas individuias cujo nome é um número aleatório único, sendo este número o nome do arquivo de texto, o do arquivo .npy e o arquivo com os valores calculados. Abaixo há o excerto da funcão, em Python, que calcula o instante seguinte do nosso sistema, utilizando as funções RFFT2 e IRFFT2 do pacote Scipy.

def cahnfourier2d(aa, kk2, kk4):
    cct = rfft2(aa)
    cct3 = rfft2(aa**3)
    cct = cct + difd*dt*(-kk2*(cct3 - cct) - kk4*cct)
    ccn = irfft2(cct)
    return ccn

No código, kk2 e kk4 são os coeficientes elevados às respectivas potências, porém kk4 está multiplicado pelo valor de gamma. Vale destacar a semelhança entre este código e o unidimensional, fruto da discretização quase idêntica.

Resultado em Duas Dimensões: Concentração Média 0

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^{2},x_{max}=1$. 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^{14}$ elementos, também fazendo uso da eficiência completa da função RFFT2 (equivalente bidimensional da função RFFT).

O motivo de não explorarmos valores diferentes de $\gamma$ nem de $D$ está na sensibilidade da estabilidade da equação: pequenas mudanças na ordem de grandeza de qualquer uma das constantes pode resultar num sistema instável.

Segue abaixo mídias da difusão sob estes parâmetros com concentração média 0. A "baixa qualidade" das imagens se dá pelo tamanho de array utilizado (128x128). É possível executar o programa com um $\Delta x$ menor, logo um array maior, porém o tamanho escolhido já permite analisarmos os dados, além de ter um tempo de processamento menor (lembre-se que estamos utilizando um $\Delta t$ minúsculo, além de que, utilizar metade do $\Delta x$ resultará em um array com 4 vezes o número de células).

Instante inicial com concentrações aleatórias.

Tempo máximo da animação: 0.012. Repare na breve homogenização do sistema, seguido da rápida formação das fases.

Instante expandido a fim de demonstrar as condições de contorno periódicas. Esta imagem é composta por uma grade 3x3 das de tamanho menor.

Tempo máximo da animação: 0.3. Este vídeo está acelerado em relação ao anterior. Repare na difusão das "bolhas" que não entraram em contato com outras concentrações.

Tempo máximo da animação: 3. Este vídeo está acelerado em relação ao anterior. Repare no arredondamento da concentração se formou.

Tempo máximo da animação: 6. Este vídeo está acelerado em relação ao anterior. Utilizando da visualização expandida, podemos ver que, ao atingir a estabilidade, uma das fases formou uma concentração circular, enquanto a outra ocupou o espaço remanescente.

Último instante da simulação, no qual as fases atingiram a estabilidade.

Resultado em Duas Dimensões: Concentração Média 0.5

Abaixo seguem mídias da simulação de condições iniciais idênticas às anteriores, porém com concentração média inicial 0.5 (equivalente a 75% uma fase e 25% em outra).

Instante inicial com concentrações aleatórias. Repare nas cores serem majoritariamente verdes ou azuis, relativas à concentrações de valor positivo.

Tempo máximo da animação: 0.012. Repare na breve homogenização do sistema, seguido da rápida formação das fases pequenas isoladas.

Instante expandido a fim de demonstrar as condições de contorno periódicas. Esta imagem é composta por uma grade 3x3 das de tamanho menor.

Tempo máximo da animação: 0.3. Este vídeo está acelerado em relação ao anterior. Repare na difusão das "bolhas" que não entraram em contato com outras concentrações.

Tempo máximo da animação: 3. Este vídeo está acelerado em relação ao anterior. Repare na formação de 3 concentrações já circulares de forma, além do início da difusão de uma delas.

Tempo máximo da animação: 7. Este vídeo está acelerado em relação ao anterior. Utilizando da visualização expandida, podemos ver novamente que, ao atingir a estabilidade, uma das fases formou uma concentração circular, enquanto a outra ocupou o espaço remanescente.

Último instante da simulação, no qual as fases atingiram a estabilidade.

Discussão dos resultados em duas dimensões

Para tempos curtos, vemos o fenômeno da homogenização do sistema, seguida da formação das regiões das diferentes fases, conforme as concentrações aleatórias difundem. Como temos uma condição inicial de valores aleatórios, ela é equivalente ao cenário real de uma mistura homogênea com perturbações. Logo a homogenização que percebemos não é parte da realidade, mas sim uma breve adequação indireta da simulação.

Após poucos instantes, já percebemos a separação completa das fases, com exceção das fronteiras, que apresentam uma mistura de ambas. É importante destacar que, conforme maiores as concentrações de fases se tornam, menor a velocidade que elas se "movem". Ao lembrar que a equação de Cahn-Hilliard descreve a separação de fases em um sistema homogêneo termodiamicamente instável, que reduz sua energia conforme a separação, a desaçeleração da separação ocorre pela aproximação da sua energia atual à energia mínima do sistema.

No entanto, vemos que a tendência é a redução dessas fronteiras, de modo que as cocentrações maiores tendem a encostar nos menores (ou os menores difundem no ambiente, caso não haja uma concentração maior próximo o suficiente). Mesmo sem saber que o estado de equilíbrio, que neste caso, é um círculo, poderíamos intuir que seria ele: o círculo é a forma que possui maior área por perímetro, e vimos que as concentrações morfam de modo a diminuir esse perímetro. No caso de concentrações diferentes entre as fases, a de menor volume (ou neste caso, área) total, forma o círculo, caso contrário, é "aleatório" (depende da posição das concentrações iniciais).