Equação de Cahn-Hilliard

Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David, Guilherme Hoss

O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição spinodal de uma mistura binária, utilizando o método FTCS (Forward Time Centered Space).

Decomposição Espinodal

Decomposição espinodal é o nome dado ao processo no qual uma pequena perturbação de um sistema faz com que, uma fase homogênea termodinamicamente instável, diminua sua energia e separe-se espontaneamente em duas outras fases coexistentes, esse é um processo que ocorre sem nucleação, ou seja, é instantâneo. Ela é observada, por exemplo, em misturas de metais ou polímeros e pode ser modelada pela equação de Cahn-Hilliard.

A Equação de Cahn-Hilliard

A equação de Cahn-Hilliard descreve o processo de decomposição espinodal de uma mistura binária. Em outras palavras, é uma equação que descreve o processo de separação de fase entre dois componentes de um fluido binário que se separam de maneira espontânea.

Consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos $c_{a}(x,t)$ e $c_{b}(x,t)$, respectivamente.

Além disso, podemos considerar que - para uma mistura binária - $c_{a}(x,t)+c_{b}(x,t)=1$ e portanto podemos simplificar para apenas uma concentração $c(x,t)$:

$c_{a}(x,t)=c(x,t),c_{b}(x,t)=1-c(x,t)$

Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:

$J=-D\nabla c$

juntamente da equação da continuidade:

${\frac {\partial c}{\partial t}}+\nabla \cdot {\vec {J}}=0$

Onde $D$ é o coeficiente de difusão e ${\vec {J}}$ é o fluxo de difusão. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:

${\frac {\partial c}{\partial t}}=D{\nabla }^{2}c$

Por definição, verificou-se que a concentração não poderia ser a razão da difusão, portanto outra força estaria presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico. Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:

$J=-M\nabla \mu$

Onde $M$ é a mobilidade das partículas (análoga à D) e $\mu$ é o potencial químico. Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick:

${\frac {\partial c}{\partial t}}=M{\nabla }^{2}\mu$

Essa equação também é conhecida como equação de Cahn-Hilliard.

Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:

$\mu ={\frac {\partial g}{\partial c}}$

Onde $g$ é a densidade da energia livre de Gibbs e $c$ é a concentração.

Tendo em vista a substituição do termo $\mu$ por um termo que envolva a concentração dos fluidos, utiliza-se uma equação que descreve a densidade de energia desse sistema através da concentração dos mesmos:

$G=\int _{V}^{}f(c)+{\kappa |\nabla c|}^{2}dV$

Nesse caso, $G$ é a energia livre de Gibbs, $f(c)$ é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e ${\kappa |\nabla c|}^{2}$ é a densidade de energia livre devido ao gradiente de concentração na interface (ou energia de interface).

Além disso, a função $f(c)$ tem o formato de um poço de potencial duplo, que pode ser representado pela seguinte equação:

$f(c)={\frac {(c^{2}-1)^{2}}{4}}$

Levando essas informações em conta e - utilizando um parâmetro $\gamma$ análogo à largura da interface - que é descrito por $\kappa =\gamma ^{2}$ é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema duplo-fásico:

$g(c)=f(c)+{\gamma }^{2}{|\nabla c|}^{2}$

Com essas igualdades agora se torna possível o cálculo de $\mu$ em função da concentração dos fluidos:

$\mu ={\frac {\partial g}{\partial c}}={\frac {\partial f(c)}{\partial c}}+{\frac {\partial {(\gamma }^{2}{|\nabla c|}^{2})}{\partial c}}={\frac {\partial ({\frac {(c^{2}-1)^{2}}{4}})}{\partial c}}+{\gamma }^{2}{\nabla }^{2}c=c^{3}-c+{\gamma }^{2}{\nabla }^{2}c$

Finalmente - utilizando a última expressão encontrada - torna-se possível reescrever o potencial químico em função da mobilidade de suas partículas ($M$) e a concetração do fluido:

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

Essa equação final é chamada de equação de Cahn-Hilliard.

Método FTCS (Forward Time Centered Space)

O FTCS é um método numérico utilizado para resolver equações diferenciais parciais, tais como a difusão do calor e do transporte de massa, traduzindo, significa "Progressivo no tempo, avançado no espaço". Esse método pode ser utilizado em sua forma implícita ou explícita que estão descritas abaixo.

$n\to \Delta t$

$j\to \Delta x$

FTCS Explicito

${\frac {\partial f}{\partial t}}\to {\frac {f_{j}^{n+1}-f_{j}^{n}}{\Delta t}}$

${\frac {\partial ^{2}f}{\partial x^{2}}}\to {\frac {f_{j-1}^{n}-2f_{j}^{n}+f_{j+1}^{n}}{\Delta x^{2}}}$

Para difusão:

$f_{j}^{n+1}=f_{j}^{n}+{\frac {D\Delta t}{(\Delta x)^{2}}}(f_{j-1}^{n}-2f_{j}^{n}+f_{j+1}^{n})$

FTCS Implicito (BTCS)

${\frac {\partial f}{\partial t}}\to {\frac {f_{j}^{n}-f_{j}^{n-1}}{\Delta t}}$

${\frac {\partial ^{2}f}{\partial x^{2}}}\to {\frac {f_{j-1}^{n}-2f_{j}^{n}+f_{j+1}^{n}}{\Delta x^{2}}}$

Para difusão:

$f_{j}^{n+1}=f_{j}^{n}+{\frac {D\Delta t}{(\Delta x)^{2}}}(f_{j-1}^{n+1}-2f_{j}^{n+1}+f_{j+1}^{n+1})$

Resolução do Cahn-Hilliard Equation para FTCS explicito para x somente:

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

$\displaystyle {\frac {c_{j}^{n+1}-c_{j}^{n}}{\Delta t}}=D\displaystyle {\frac {\partial ^{2}}{\partial x^{2}}}(c^{3}-c-\gamma ^{2}\displaystyle {\frac {\partial ^{2}c}{\partial x^{2}}})$

$\displaystyle {\frac {c_{j}^{n+1}-c_{j}^{n}}{\Delta t}}=D{\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{(\Delta x)^{2}}}$

$\displaystyle {\frac {c_{j}^{n+1}-c_{j}^{n}}{\Delta t}}=D\left({\frac {(c_{j-1}^{n})^{3}-2(c_{j}^{n})^{3}+(c_{j+1}^{n})^{3}}{(\Delta x)^{2}}}-{\frac {c_{j-1}^{n}-2c_{i}^{n}+c_{j+1}^{n}}{(\Delta x)^{2}}}-\gamma ^{2}{\frac {c_{j-2}^{n}-4c_{j-1}^{n}+6c_{j}^{n}-4c_{j+1}^{n}+c_{j+2}^{n}}{(\Delta x)^{4}}}\right)$

$c_{j}^{n+1}=D\Delta t\left({\frac {(c_{j-1}^{n})^{3}-2(c_{i}^{n})^{3}+(c_{j+1}^{n})^{3}}{(\Delta x)^{2}}}-{\frac {c_{j-1}^{n}-2c_{j}^{n}+c_{j+1}^{n}}{(\Delta x)^{2}}}-\gamma ^{2}{\frac {c_{j-2}^{n}-4c_{j-1}^{n}+6c_{j}^{n}-4c_{j+1}^{n}+c_{j+2}^{n}}{(\Delta x)^{4}}}\right)+c_{j}^{n}$

$c_{j}^{n+1}={\frac {D\Delta t}{(\Delta x)^{2}}}\left((c_{j-1}^{n})^{3}-2(c_{i}^{n})^{3}+(c_{j+1}^{n})^{3}-{c_{j-1}^{n}+2c_{j}^{n}-c_{j+1}^{n}}-\gamma ^{2}{\frac {c_{j-2}^{n}-4c_{j-1}^{n}+6c_{j}^{n}-4c_{j+1}^{n}+c_{j+2}^{n}}{(\Delta x)^{2}}}\right)+c_{j}^{n}$

Condição de Estabilidade

$\Delta t<\displaystyle {\frac {(\Delta x)^{2}}{2}}$