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)$

Lembrando que, conforme explicado no trabalho anterior a este, $c$ é a soma da concentração de duas fases, e seu valor varia de -1 (presença de só uma das fases) a 1 (presença somente da outra fase). 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. As cores variam de vermelho (-1) a violeta (1), de modo que as cores variam de acordo com o espectro visível, onde a cor amarela corresponde a 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 delas, a desaçeleração da separação ocorre pela aproximação da sua energia atual à energia mínima do sistema.

Percebemos que o estado de menor energia é uma concentração circular, enquanto o restante do meio é ocupado pela outra. Mesmo sem saber que este é o estado estacionário, poderíamos intuir que seria ele ao ver os instantes iniciais da simulação: o círculo é a forma que possui maior área por perímetro, e vimos que as fases diferentes morfam de modo a diminuir as fronteiras entre elas. No caso de uma quantidade diferente entre as fases, a de menor volume (ou neste caso, área) total, forma o círculo, caso contrário, é "aleatório" (depende dos locais das perturbações).

Códigos

Todos códigos referentes a duas dimensões possuem, nas primeiras linhas, comentários que destacam quais valores o usuário é recomendado a alterar. Caso você não possua experiência com programação, é recomendado não alterar as linhas de código abaixo dos locais indicados.

Mídia:cahn2dfourier.py Programa utilizado para gerar valores sem partir de nenhum arquivo anterior. Os dados são identificados por um código numérico único, e serão salvos em uma pasta com nome de tal código, onde ali se encontrarão:

- O arquivo código.npy principal, que consiste nos arrays.

- O arquivo valcódigo.npy, que consiste em um arquivo que é lido pelos outros programas contendo informações cruciais.

- O arquivo código.txt, que contém tais informações cruciais em um formato legível ao usuário.

Utilizado em conjunto com os programas "framer" e "animador", você pode gerar imagens como as vistas neste trabalho.

Mídia:cahn2dfouriercontinue.py Programa utilizado para gerar valores partindo de outros já antes calculados. O arquivo com os dados novos (concatenados aos antigos) se encontrará na pasta original, substituindo o arquivo antigo. Recomendado para gerar valores por tempos longos sem precisar deixar o programa executando por longas horas sem intervalo.

Mídia:cahn2dframer.py Programa utilizado para gerar as imagens estáticas vistas neste trabalho.

Mídia:cahn2danimador.py programa utilizado para gerar as animações vistas neste trabalho. Atente que, para computadores sem placa de vídeo dedicada, a renderização dos gifs pode levar muito tempo, em casos como este recomendo gerar as animações em menor resolução ou com intervalos maiores entre as imagens.

Código utilizado para gerar o frame de comparação em uma dimensão. Não recomendo acessá-lo, por não possuir um resultado interessante, além de não estar otimizado nem comentado. Mídia:cahn1dcompare.py