http://fiscomp.if.ufrgs.br/api.php?action=feedcontributions&user=Arthurdorn&feedformat=atomFísica Computacional - Contribuições do usuário [pt-br]2024-03-29T07:04:02ZContribuições do usuárioMediaWiki 1.39.4http://fiscomp.if.ufrgs.br/index.php?title=Equa%C3%A7%C3%A3o_de_Cahn-Hilliard&diff=4013Equação de Cahn-Hilliard2021-04-05T22:18:59Z<p>Arthurdorn: /* Discussão de Resultados */</p>
<hr />
<div>'''Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David e Guilherme Hoss'''<br />
<br />
O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição espinodal de uma mistura binária, e analisar como é seu comportamento com diferentes coeficientes de difusão, utilizando o método FTCS (''Forward Time Centered Space''). O trabalho foi inspirado no artigo de Sibbing[1].<br />
<br />
== Decomposição Espinodal ==<br />
<br />
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 [2].<br />
<br />
==A Equação de Cahn-Hilliard==<br />
<br />
A equação de Cahn-Hilliard é 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.<br />
Com o intuito de deduzirmos essa equação, consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos <math> c_a(x,t) </math> e <math> c_b(x,t) </math>, respectivamente. [1]<br />
<br />
Além disso, podemos considerar que - para uma mistura binária - <math> c_a(x,t) + c_b(x,t) = 1</math> e portanto podemos simplificar para apenas uma concentração <math> c(x,t) </math>:<br />
:<math><br />
c_a(x,t) = c(x,t), c_b(x,t) = 1 - c(x,t).<br />
</math><br />
<br />
Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:<br />
:<math><br />
J = -D\nabla c<br />
</math><br />
<br />
juntamente da equação da continuidade:<br />
:<math><br />
\frac{\partial c}{\partial t} + \nabla \cdot \vec J = 0.<br />
</math><br />
<br />
Onde <math>D</math> é o coeficiente de difusão e <math>\vec J</math> é o fluxo de difusão de concentração da mistura. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 c.<br />
</math><br />
<br />
A partir dessa equação - como não há a existência de um gradiente de concentração espacial - pode-se esperar que não ocorra mudança na concentração da mistura. No entanto, observa-se que quando a separação de fases ocorre, a difusão demonstra ser contrária ao gradiente de concentração, o que não condiz com a equação anterior.<br />
Tendo isso em vista, conclui-se que a concentração não pode ser a razão da difusão, portanto outra força deve estar presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico (de acordo com Cahn e Hilliard, 1958). Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:<br />
<br />
:<math><br />
J = -M \nabla \mu.<br />
</math><br />
<br />
Onde <math>M</math> é a mobilidade das partículas (análoga à D) e <math>\mu</math> é o potencial químico.<br />
Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick [3]:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 \mu.<br />
</math><br />
<br />
Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c}.<br />
</math><br />
<br />
Onde <math>g</math> é a densidade da energia livre de Gibbs e <math>c</math> é a concentração (de acordo com Schroeder, 1999).<br />
<br />
Tendo em vista a substituição do termo <math>\mu</math> 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 (derivado em [3]):<br />
<br />
:<math><br />
G = \int_{V}^{} f(c) + {\kappa |\nabla c|}^2 dV.<br />
</math><br />
<br />
Nesse caso, <math>G</math> é a energia livre de Gibbs, <math>f(c)</math> é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e <math>{\kappa|\nabla c|}^2</math> é a densidade de energia livre devido ao gradiente de concentração.<br />
<br />
Além disso, a função <math>f(c)</math> - de acordo com [8] - possui o potencial de um poço de potencial duplo. Neste poço, <math>c</math> representa a concentração em escala e está relacionada à temperatura da mistura, que decide se a separação de fases irá - ou não - ocorrer. Esta função pode ser representada pela seguinte equação:<br />
<br />
:<math><br />
f(c) = \frac{(c^2 - 1)^2}{4}.<br />
</math><br />
<br />
Levando essas informações em conta e - utilizando um parâmetro <math>\gamma</math> análogo à largura da interface - que é descrito por <math>\kappa = \gamma ^2</math> é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema com duas fases.<br />
:<math><br />
g(c) = f(c) + {\gamma}^2{|\nabla c |}^2.<br />
</math><br />
<br />
Com essas igualdades agora se torna possível o cálculo de <math>\mu</math> em função da concentração dos fluidos:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c} = \frac{\partial f(c)}{\partial c} + \frac{\partial }{\partial c} (\gamma^2{|\nabla c |}^2) = \frac{\partial }{\partial c} \frac{(c^2 - 1)^2}{4} + {\gamma}^2 {\nabla}^2 c = c^3 - c + {\gamma}^2 {\nabla}^2 c.<br />
</math><br />
<br />
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 (<math>M</math>) e a concentração do fluido:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
Essa equação final é chamada de equação de Cahn-Hilliard. <br />
A equação dependente da difusão é análoga e também funcional:<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
== Método FTCS (Forward Time Centered Space) ==<br />
<br />
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, centrado no espaço". Uma das formas que o método pode ser utilizado é a forma explícita que está descritas abaixo.<br />
<br />
:<math>n\to\Delta t</math><br />
<br />
:<math>j\to\Delta x</math><br />
<br />
===FTCS Explicito===<br />
<br />
:<math>\frac{\partial f}{\partial t}\to \frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}</math><br />
<br />
:<math>\frac{\partial ^2 f}{\partial x^2}\to \frac{f_{j-1}^{n}-2 f_{j}^{n} + f_{j+1}^n}{\Delta x^2}</math><br />
<br />
Para difusão:<br />
<br />
:<math>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})</math><br />
<br />
== Implementação da equação de Cahn-Hilliard 1D pelo método FTCS explicito ==<br />
<br />
<br />
:<math>\displaystyle \frac{\partial c}{\partial t} = D\nabla^{2}(c^{3}-c-\gamma^2\nabla^{2}c)</math> <br />
<br />
:<math>\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})</math><br />
<br />
:<math>\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}</math><br />
<br />
<br />
:<math>\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-2 c_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)</math><br />
<br />
:<math> 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-2 c_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</math><br />
<br />
:<math>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+2 c_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 </math><br />
<br />
'''Condição de Estabilidade '''<br />
<br />
A estabilidade dessa equação mostra-se muito mais complicada de se estipular por ela ser uma equação diferencial de quarta ordem se comparada a equação de difusão, Portanto só iremos analisar a seção 3.3 (Experimental and theoretical stability conditions) do artigo Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions [1]. <br />
<br />
Após a linearização e aplicação do teorema de Gershgorin temos que a condição para estabilidade da equação linear para <math>D = 1</math> é:<br />
<br />
<br />
:<center><math>\Delta t < \displaystyle\frac{(\Delta x)^2}{4+\frac{8\gamma^2}{\Delta x^2}}</math> <br />
<br />
Importante atentar que essa é a condição de estabilidade somente para a equação de Cahn-Hilliard linearizada, não para a original. Tanto que a literatura sobre a equação propõem que <math>\Delta t</math> <math>\varpropto</math> <math>(\Delta x)^4 </math>, que é o que acontece na condição estabilidade linear quando <math>\gamma >> \Delta x</math>.<br />
<br />
O artigo compara os dados experimentais de estabilidade com a estabilidade da equação linearizada relacionado na seção 3.3.4 e conclui que para valores de <math>\frac{\gamma}{\Delta x} \in [0.25,8]</math> a condição teórica encontrada a partir da linearização é uma boa aproximação.<br />
<br />
== Resultados ==<br />
<br />
Com o intuito de testar como o fator de difusão D afeta a evolução da equação de Cahn-Hilliard, comparamos os resultados para os coeficientes de difusão 1, 0.1, 0.01 e 0.001 e analisamos seus gráficos. Para esse estudo foram utilizadas como condições iniciais a concentração -1 para a primeira metade do espaço e 1 para a segunda metade, e o sistema é livre de condições de contorno. <br />
<br />
<br />
[[Arquivo:coeficientes.png|1000px|thumb|center| Gráfico 1: Resultados da simulação variando os coeficientes de difusão (D) para tempos máximos diferentes, <math>\gamma=3.4/128,\Delta t = 1/500000</math>]]<br />
<br />
<br />
Nos gráficos, é possível observar que quanto maior o coeficiente de difusão maior é a velocidade em que a mistura atinge a estabilidade. Além disso, vemos que valores baixos de t produzem soluções mais íngremes que valores altos de t.<br />
<br />
== Discussão de Resultados ==<br />
<br />
Com o objetivo de resolvermos computacionalmente a equação de Cahn-Hilliard de uma dimensão, era de se esperar que ao longo das iterações a inclinação das concentrações fossem menos acentuadas. Podemos observar em todas as situações que conforme evoluímos temporalmente, a derivada primeira da curva vai diminuindo próxima ao centro do gráfico.<br />
<br />
Uma propriedade observada no gráfico 1 é a de que os valores das soluções obtidas utilizando o método FTCS excedem os valores máximos e mínimos permitidos ( C=1 e C=-1), se estivesse modelando uma situação real isso iria contra a lei de conservação de massa, o que pode ocasionar erros nos resultados que exigem uma grande precisão. <br />
<br />
O método FTCS explícito limita-se por causa da condição de estabilidade, por isso esse método não é recomendado para modelos com alto <math>\gamma</math>. Para tais modelos, o método FTCS implícito é mais recomendado por ser incondicionalmente estável.<br />
<br />
== Implementção ==<br />
<br />
<source lang="python"><br />
<br />
import matplotlib.pyplot as plt <br />
<br />
def vector_declaration(L,dx):<br />
c = [[],[]] # vetor concentração<br />
espaco = []<br />
# Condições iniciais<br />
for i in range(int(L/dx)+4):##+4 pois usaremos dois valores antes e depois do ultimo elemento do vetor c<br />
if (i<1/2*L/dx+2):<br />
c[0].append(-1)<br />
c[1].append(-1)<br />
else:<br />
c[0].append(1)<br />
c[1].append(1)<br />
espaco.append(round(i/150,3))<br />
return c, espaco<br />
<br />
def CH_equation(gamma, D, dx, dt, L, TEMPO_MAX): # resolução numérica da equação<br />
<br />
c, espaco = vector_declaration(L, dx)<br />
i = 0<br />
for time in [t*dt/TEMPO_MAX for t in range(int(TEMPO_MAX/dt))]:<br />
for l in range(2,len(c[1][2:-2])):<br />
c1 = c[i][l-1]**3 - 2*c[i][l]**3 + c[i][l+1]**3<br />
c2 = -c[i][l-1] + 2*c[i][l] - c[i][l+1]<br />
c3 = c[i][l-2] -4*c[i][l-1] + 6*c[i][l] - 4*c[i][l+1] + c[i][l+2]<br />
c[1-i][l] = D*dt/(dx**2)*(c1+c2-gamma**2*c3/(dx**2)) + c[i][l]<br />
i = 1-i<br />
<br />
return(espaco[2:-2],c[1-i][2:-2]) ##retirando os elementos a mais do vetor<br />
<br />
tempo=1<br />
tamanho=1/3.5<br />
Difuse= 1<br />
gamma=3.4*1/128<br />
dt=1/500000<br />
dx = 1/128<br />
<br />
figure, axis = plt.subplots(2, 2)<br />
plt.figure(dpi=500)<br />
<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
plt.show()<br />
plt.savefig('graficosDIFUSAO.png', dpi=500)<br />
<br />
</source><br />
<br />
== Referências ==<br />
<br />
[1] SIBBING, Zimo. '''Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions''', tese de bacharelado, 2015.<br />
<br />
[2] Spinodal Decomposition, disponível em: https://en.wikipedia.org/wiki/Spinodal_decomposition<br />
<br />
[3] MARKUS, Wilczek. '''The Cahn-Hilliard Equation''', 2015.<br />
<br />
[4] CAHN, John W.; HILLIARD, John E. '''Free Energy of a Nonuniform System. I. Interfacial Free Energy'''. The Journal of Chemical Physics, 1958.<br />
<br />
[5] CAHN, John W.; HILLIARD, John E. '''Spinodal decomposition: A reprise'''Acta Metallurgica, Volume 19, Issue 2, 1971<br />
<br />
[6] Lei de Fick, disponível em: https://pt.wikipedia.org/wiki/Lei_de_Fick<br />
<br />
[7] Cahn-Hilliard Equation, disponível em: https://pt.qaz.wiki/wiki/Cahn%E2%80%93Hilliard_equation<br />
<br />
[8] Daniel V. Schroeder. An Introduction to Thermal Physics. Addison-Wesley, 1999.<br />
<br />
[9] Dongsun Lee, Joo-Youl Huh, Darae Jeong, Jaemin Shin, Ana Yun, and Junseok Kim. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Computational Materials Science, 2014.<br />
<br />
[10] Neumann Boundary Condition, disponível em: https://en.wikipedia.org/wiki/Neumann_boundary_condition</div>Arthurdornhttp://fiscomp.if.ufrgs.br/index.php?title=Equa%C3%A7%C3%A3o_de_Cahn-Hilliard&diff=4012Equação de Cahn-Hilliard2021-04-05T22:00:51Z<p>Arthurdorn: /* Referências */</p>
<hr />
<div>'''Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David e Guilherme Hoss'''<br />
<br />
O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição espinodal de uma mistura binária, e analisar como é seu comportamento com diferentes coeficientes de difusão, utilizando o método FTCS (''Forward Time Centered Space''). O trabalho foi inspirado no artigo de Sibbing[1].<br />
<br />
== Decomposição Espinodal ==<br />
<br />
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 [2].<br />
<br />
==A Equação de Cahn-Hilliard==<br />
<br />
A equação de Cahn-Hilliard é 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.<br />
Com o intuito de deduzirmos essa equação, consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos <math> c_a(x,t) </math> e <math> c_b(x,t) </math>, respectivamente. [1]<br />
<br />
Além disso, podemos considerar que - para uma mistura binária - <math> c_a(x,t) + c_b(x,t) = 1</math> e portanto podemos simplificar para apenas uma concentração <math> c(x,t) </math>:<br />
:<math><br />
c_a(x,t) = c(x,t), c_b(x,t) = 1 - c(x,t).<br />
</math><br />
<br />
Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:<br />
:<math><br />
J = -D\nabla c<br />
</math><br />
<br />
juntamente da equação da continuidade:<br />
:<math><br />
\frac{\partial c}{\partial t} + \nabla \cdot \vec J = 0.<br />
</math><br />
<br />
Onde <math>D</math> é o coeficiente de difusão e <math>\vec J</math> é o fluxo de difusão de concentração da mistura. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 c.<br />
</math><br />
<br />
A partir dessa equação - como não há a existência de um gradiente de concentração espacial - pode-se esperar que não ocorra mudança na concentração da mistura. No entanto, observa-se que quando a separação de fases ocorre, a difusão demonstra ser contrária ao gradiente de concentração, o que não condiz com a equação anterior.<br />
Tendo isso em vista, conclui-se que a concentração não pode ser a razão da difusão, portanto outra força deve estar presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico (de acordo com Cahn e Hilliard, 1958). Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:<br />
<br />
:<math><br />
J = -M \nabla \mu.<br />
</math><br />
<br />
Onde <math>M</math> é a mobilidade das partículas (análoga à D) e <math>\mu</math> é o potencial químico.<br />
Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick [3]:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 \mu.<br />
</math><br />
<br />
Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c}.<br />
</math><br />
<br />
Onde <math>g</math> é a densidade da energia livre de Gibbs e <math>c</math> é a concentração (de acordo com Schroeder, 1999).<br />
<br />
Tendo em vista a substituição do termo <math>\mu</math> 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 (derivado em [3]):<br />
<br />
:<math><br />
G = \int_{V}^{} f(c) + {\kappa |\nabla c|}^2 dV.<br />
</math><br />
<br />
Nesse caso, <math>G</math> é a energia livre de Gibbs, <math>f(c)</math> é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e <math>{\kappa|\nabla c|}^2</math> é a densidade de energia livre devido ao gradiente de concentração.<br />
<br />
Além disso, a função <math>f(c)</math> - de acordo com [8] - possui o potencial de um poço de potencial duplo. Neste poço, <math>c</math> representa a concentração em escala e está relacionada à temperatura da mistura, que decide se a separação de fases irá - ou não - ocorrer. Esta função pode ser representada pela seguinte equação:<br />
<br />
:<math><br />
f(c) = \frac{(c^2 - 1)^2}{4}.<br />
</math><br />
<br />
Levando essas informações em conta e - utilizando um parâmetro <math>\gamma</math> análogo à largura da interface - que é descrito por <math>\kappa = \gamma ^2</math> é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema com duas fases.<br />
:<math><br />
g(c) = f(c) + {\gamma}^2{|\nabla c |}^2.<br />
</math><br />
<br />
Com essas igualdades agora se torna possível o cálculo de <math>\mu</math> em função da concentração dos fluidos:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c} = \frac{\partial f(c)}{\partial c} + \frac{\partial }{\partial c} (\gamma^2{|\nabla c |}^2) = \frac{\partial }{\partial c} \frac{(c^2 - 1)^2}{4} + {\gamma}^2 {\nabla}^2 c = c^3 - c + {\gamma}^2 {\nabla}^2 c.<br />
</math><br />
<br />
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 (<math>M</math>) e a concentração do fluido:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
Essa equação final é chamada de equação de Cahn-Hilliard. <br />
A equação dependente da difusão é análoga e também funcional:<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
== Método FTCS (Forward Time Centered Space) ==<br />
<br />
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, centrado no espaço". Uma das formas que o método pode ser utilizado é a forma explícita que está descritas abaixo.<br />
<br />
:<math>n\to\Delta t</math><br />
<br />
:<math>j\to\Delta x</math><br />
<br />
===FTCS Explicito===<br />
<br />
:<math>\frac{\partial f}{\partial t}\to \frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}</math><br />
<br />
:<math>\frac{\partial ^2 f}{\partial x^2}\to \frac{f_{j-1}^{n}-2 f_{j}^{n} + f_{j+1}^n}{\Delta x^2}</math><br />
<br />
Para difusão:<br />
<br />
:<math>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})</math><br />
<br />
== Implementação da equação de Cahn-Hilliard 1D pelo método FTCS explicito ==<br />
<br />
<br />
:<math>\displaystyle \frac{\partial c}{\partial t} = D\nabla^{2}(c^{3}-c-\gamma^2\nabla^{2}c)</math> <br />
<br />
:<math>\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})</math><br />
<br />
:<math>\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}</math><br />
<br />
<br />
:<math>\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-2 c_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)</math><br />
<br />
:<math> 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-2 c_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</math><br />
<br />
:<math>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+2 c_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 </math><br />
<br />
'''Condição de Estabilidade '''<br />
<br />
A estabilidade dessa equação mostra-se muito mais complicada de se estipular por ela ser uma equação diferencial de quarta ordem se comparada a equação de difusão, Portanto só iremos analisar a seção 3.3 (Experimental and theoretical stability conditions) do artigo Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions [1]. <br />
<br />
Após a linearização e aplicação do teorema de Gershgorin temos que a condição para estabilidade da equação linear para <math>D = 1</math> é:<br />
<br />
<br />
:<center><math>\Delta t < \displaystyle\frac{(\Delta x)^2}{4+\frac{8\gamma^2}{\Delta x^2}}</math> <br />
<br />
Importante atentar que essa é a condição de estabilidade somente para a equação de Cahn-Hilliard linearizada, não para a original. Tanto que a literatura sobre a equação propõem que <math>\Delta t</math> <math>\varpropto</math> <math>(\Delta x)^4 </math>, que é o que acontece na condição estabilidade linear quando <math>\gamma >> \Delta x</math>.<br />
<br />
O artigo compara os dados experimentais de estabilidade com a estabilidade da equação linearizada relacionado na seção 3.3.4 e conclui que para valores de <math>\frac{\gamma}{\Delta x} \in [0.25,8]</math> a condição teórica encontrada a partir da linearização é uma boa aproximação.<br />
<br />
== Resultados ==<br />
<br />
Com o intuito de testar como o fator de difusão D afeta a evolução da equação de Cahn-Hilliard, comparamos os resultados para os coeficientes de difusão 1, 0.1, 0.01 e 0.001 e analisamos seus gráficos. Para esse estudo foram utilizadas como condições iniciais a concentração -1 para a primeira metade do espaço e 1 para a segunda metade, e o sistema é livre de condições de contorno. <br />
<br />
<br />
[[Arquivo:coeficientes.png|1000px|thumb|center| Gráfico 1: Resultados da simulação variando os coeficientes de difusão (D) para tempos máximos diferentes, <math>\gamma=3.4/128,\Delta t = 1/500000</math>]]<br />
<br />
<br />
Nos gráficos, é possível observar que quanto maior o coeficiente de difusão maior é a velocidade em que a mistura atinge a estabilidade. Além disso, vemos que valores baixos de t produzem soluções mais íngremes que valores altos de t.<br />
<br />
== Discussão de Resultados ==<br />
<br />
Uma propriedade observada no gráfico 1 é a de que os valores das soluções obtidas utilizando o método FTCS excedem os valores máximos e mínimos permitidos ( C=1 e C=-1), se estivesse modelando uma situação real isso iria contra a lei de conservação de massa, o que pode ocasionar erros nos resultados que exigem uma grande precisão. <br />
<br />
O método FTCS explícito limita-se por causa da condição de estabilidade, por isso esse método não é recomendado para modelos com alto <math>\gamma</math>. Para tais modelos, o método FTCS implícito é mais recomendado por ser incondicionalmente estável.<br />
<br />
== Implementção ==<br />
<br />
<source lang="python"><br />
<br />
import matplotlib.pyplot as plt <br />
<br />
def vector_declaration(L,dx):<br />
c = [[],[]] # vetor concentração<br />
espaco = []<br />
# Condições iniciais<br />
for i in range(int(L/dx)+4):##+4 pois usaremos dois valores antes e depois do ultimo elemento do vetor c<br />
if (i<1/2*L/dx+2):<br />
c[0].append(-1)<br />
c[1].append(-1)<br />
else:<br />
c[0].append(1)<br />
c[1].append(1)<br />
espaco.append(round(i/150,3))<br />
return c, espaco<br />
<br />
def CH_equation(gamma, D, dx, dt, L, TEMPO_MAX): # resolução numérica da equação<br />
<br />
c, espaco = vector_declaration(L, dx)<br />
i = 0<br />
for time in [t*dt/TEMPO_MAX for t in range(int(TEMPO_MAX/dt))]:<br />
for l in range(2,len(c[1][2:-2])):<br />
c1 = c[i][l-1]**3 - 2*c[i][l]**3 + c[i][l+1]**3<br />
c2 = -c[i][l-1] + 2*c[i][l] - c[i][l+1]<br />
c3 = c[i][l-2] -4*c[i][l-1] + 6*c[i][l] - 4*c[i][l+1] + c[i][l+2]<br />
c[1-i][l] = D*dt/(dx**2)*(c1+c2-gamma**2*c3/(dx**2)) + c[i][l]<br />
i = 1-i<br />
<br />
return(espaco[2:-2],c[1-i][2:-2]) ##retirando os elementos a mais do vetor<br />
<br />
tempo=1<br />
tamanho=1/3.5<br />
Difuse= 1<br />
gamma=3.4*1/128<br />
dt=1/500000<br />
dx = 1/128<br />
<br />
figure, axis = plt.subplots(2, 2)<br />
plt.figure(dpi=500)<br />
<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
plt.show()<br />
plt.savefig('graficosDIFUSAO.png', dpi=500)<br />
<br />
</source><br />
<br />
== Referências ==<br />
<br />
[1] SIBBING, Zimo. '''Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions''', tese de bacharelado, 2015.<br />
<br />
[2] Spinodal Decomposition, disponível em: https://en.wikipedia.org/wiki/Spinodal_decomposition<br />
<br />
[3] MARKUS, Wilczek. '''The Cahn-Hilliard Equation''', 2015.<br />
<br />
[4] CAHN, John W.; HILLIARD, John E. '''Free Energy of a Nonuniform System. I. Interfacial Free Energy'''. The Journal of Chemical Physics, 1958.<br />
<br />
[5] CAHN, John W.; HILLIARD, John E. '''Spinodal decomposition: A reprise'''Acta Metallurgica, Volume 19, Issue 2, 1971<br />
<br />
[6] Lei de Fick, disponível em: https://pt.wikipedia.org/wiki/Lei_de_Fick<br />
<br />
[7] Cahn-Hilliard Equation, disponível em: https://pt.qaz.wiki/wiki/Cahn%E2%80%93Hilliard_equation<br />
<br />
[8] Daniel V. Schroeder. An Introduction to Thermal Physics. Addison-Wesley, 1999.<br />
<br />
[9] Dongsun Lee, Joo-Youl Huh, Darae Jeong, Jaemin Shin, Ana Yun, and Junseok Kim. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Computational Materials Science, 2014.<br />
<br />
[10] Neumann Boundary Condition, disponível em: https://en.wikipedia.org/wiki/Neumann_boundary_condition</div>Arthurdornhttp://fiscomp.if.ufrgs.br/index.php?title=Equa%C3%A7%C3%A3o_de_Cahn-Hilliard&diff=4011Equação de Cahn-Hilliard2021-04-05T21:59:37Z<p>Arthurdorn: /* Referências */</p>
<hr />
<div>'''Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David e Guilherme Hoss'''<br />
<br />
O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição espinodal de uma mistura binária, e analisar como é seu comportamento com diferentes coeficientes de difusão, utilizando o método FTCS (''Forward Time Centered Space''). O trabalho foi inspirado no artigo de Sibbing[1].<br />
<br />
== Decomposição Espinodal ==<br />
<br />
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 [2].<br />
<br />
==A Equação de Cahn-Hilliard==<br />
<br />
A equação de Cahn-Hilliard é 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.<br />
Com o intuito de deduzirmos essa equação, consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos <math> c_a(x,t) </math> e <math> c_b(x,t) </math>, respectivamente. [1]<br />
<br />
Além disso, podemos considerar que - para uma mistura binária - <math> c_a(x,t) + c_b(x,t) = 1</math> e portanto podemos simplificar para apenas uma concentração <math> c(x,t) </math>:<br />
:<math><br />
c_a(x,t) = c(x,t), c_b(x,t) = 1 - c(x,t).<br />
</math><br />
<br />
Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:<br />
:<math><br />
J = -D\nabla c<br />
</math><br />
<br />
juntamente da equação da continuidade:<br />
:<math><br />
\frac{\partial c}{\partial t} + \nabla \cdot \vec J = 0.<br />
</math><br />
<br />
Onde <math>D</math> é o coeficiente de difusão e <math>\vec J</math> é o fluxo de difusão de concentração da mistura. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 c.<br />
</math><br />
<br />
A partir dessa equação - como não há a existência de um gradiente de concentração espacial - pode-se esperar que não ocorra mudança na concentração da mistura. No entanto, observa-se que quando a separação de fases ocorre, a difusão demonstra ser contrária ao gradiente de concentração, o que não condiz com a equação anterior.<br />
Tendo isso em vista, conclui-se que a concentração não pode ser a razão da difusão, portanto outra força deve estar presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico (de acordo com Cahn e Hilliard, 1958). Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:<br />
<br />
:<math><br />
J = -M \nabla \mu.<br />
</math><br />
<br />
Onde <math>M</math> é a mobilidade das partículas (análoga à D) e <math>\mu</math> é o potencial químico.<br />
Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick [3]:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 \mu.<br />
</math><br />
<br />
Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c}.<br />
</math><br />
<br />
Onde <math>g</math> é a densidade da energia livre de Gibbs e <math>c</math> é a concentração (de acordo com Schroeder, 1999).<br />
<br />
Tendo em vista a substituição do termo <math>\mu</math> 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 (derivado em [3]):<br />
<br />
:<math><br />
G = \int_{V}^{} f(c) + {\kappa |\nabla c|}^2 dV.<br />
</math><br />
<br />
Nesse caso, <math>G</math> é a energia livre de Gibbs, <math>f(c)</math> é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e <math>{\kappa|\nabla c|}^2</math> é a densidade de energia livre devido ao gradiente de concentração.<br />
<br />
Além disso, a função <math>f(c)</math> - de acordo com [8] - possui o potencial de um poço de potencial duplo. Neste poço, <math>c</math> representa a concentração em escala e está relacionada à temperatura da mistura, que decide se a separação de fases irá - ou não - ocorrer. Esta função pode ser representada pela seguinte equação:<br />
<br />
:<math><br />
f(c) = \frac{(c^2 - 1)^2}{4}.<br />
</math><br />
<br />
Levando essas informações em conta e - utilizando um parâmetro <math>\gamma</math> análogo à largura da interface - que é descrito por <math>\kappa = \gamma ^2</math> é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema com duas fases.<br />
:<math><br />
g(c) = f(c) + {\gamma}^2{|\nabla c |}^2.<br />
</math><br />
<br />
Com essas igualdades agora se torna possível o cálculo de <math>\mu</math> em função da concentração dos fluidos:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c} = \frac{\partial f(c)}{\partial c} + \frac{\partial }{\partial c} (\gamma^2{|\nabla c |}^2) = \frac{\partial }{\partial c} \frac{(c^2 - 1)^2}{4} + {\gamma}^2 {\nabla}^2 c = c^3 - c + {\gamma}^2 {\nabla}^2 c.<br />
</math><br />
<br />
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 (<math>M</math>) e a concentração do fluido:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
Essa equação final é chamada de equação de Cahn-Hilliard. <br />
A equação dependente da difusão é análoga e também funcional:<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
== Método FTCS (Forward Time Centered Space) ==<br />
<br />
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, centrado no espaço". Uma das formas que o método pode ser utilizado é a forma explícita que está descritas abaixo.<br />
<br />
:<math>n\to\Delta t</math><br />
<br />
:<math>j\to\Delta x</math><br />
<br />
===FTCS Explicito===<br />
<br />
:<math>\frac{\partial f}{\partial t}\to \frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}</math><br />
<br />
:<math>\frac{\partial ^2 f}{\partial x^2}\to \frac{f_{j-1}^{n}-2 f_{j}^{n} + f_{j+1}^n}{\Delta x^2}</math><br />
<br />
Para difusão:<br />
<br />
:<math>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})</math><br />
<br />
== Implementação da equação de Cahn-Hilliard 1D pelo método FTCS explicito ==<br />
<br />
<br />
:<math>\displaystyle \frac{\partial c}{\partial t} = D\nabla^{2}(c^{3}-c-\gamma^2\nabla^{2}c)</math> <br />
<br />
:<math>\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})</math><br />
<br />
:<math>\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}</math><br />
<br />
<br />
:<math>\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-2 c_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)</math><br />
<br />
:<math> 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-2 c_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</math><br />
<br />
:<math>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+2 c_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 </math><br />
<br />
'''Condição de Estabilidade '''<br />
<br />
A estabilidade dessa equação mostra-se muito mais complicada de se estipular por ela ser uma equação diferencial de quarta ordem se comparada a equação de difusão, Portanto só iremos analisar a seção 3.3 (Experimental and theoretical stability conditions) do artigo Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions [1]. <br />
<br />
Após a linearização e aplicação do teorema de Gershgorin temos que a condição para estabilidade da equação linear para <math>D = 1</math> é:<br />
<br />
<br />
:<center><math>\Delta t < \displaystyle\frac{(\Delta x)^2}{4+\frac{8\gamma^2}{\Delta x^2}}</math> <br />
<br />
Importante atentar que essa é a condição de estabilidade somente para a equação de Cahn-Hilliard linearizada, não para a original. Tanto que a literatura sobre a equação propõem que <math>\Delta t</math> <math>\varpropto</math> <math>(\Delta x)^4 </math>, que é o que acontece na condição estabilidade linear quando <math>\gamma >> \Delta x</math>.<br />
<br />
O artigo compara os dados experimentais de estabilidade com a estabilidade da equação linearizada relacionado na seção 3.3.4 e conclui que para valores de <math>\frac{\gamma}{\Delta x} \in [0.25,8]</math> a condição teórica encontrada a partir da linearização é uma boa aproximação.<br />
<br />
== Resultados ==<br />
<br />
Com o intuito de testar como o fator de difusão D afeta a evolução da equação de Cahn-Hilliard, comparamos os resultados para os coeficientes de difusão 1, 0.1, 0.01 e 0.001 e analisamos seus gráficos. Para esse estudo foram utilizadas como condições iniciais a concentração -1 para a primeira metade do espaço e 1 para a segunda metade, e o sistema é livre de condições de contorno. <br />
<br />
<br />
[[Arquivo:coeficientes.png|1000px|thumb|center| Gráfico 1: Resultados da simulação variando os coeficientes de difusão (D) para tempos máximos diferentes, <math>\gamma=3.4/128,\Delta t = 1/500000</math>]]<br />
<br />
<br />
Nos gráficos, é possível observar que quanto maior o coeficiente de difusão maior é a velocidade em que a mistura atinge a estabilidade. Além disso, vemos que valores baixos de t produzem soluções mais íngremes que valores altos de t.<br />
<br />
== Discussão de Resultados ==<br />
<br />
Uma propriedade observada no gráfico 1 é a de que os valores das soluções obtidas utilizando o método FTCS excedem os valores máximos e mínimos permitidos ( C=1 e C=-1), se estivesse modelando uma situação real isso iria contra a lei de conservação de massa, o que pode ocasionar erros nos resultados que exigem uma grande precisão. <br />
<br />
O método FTCS explícito limita-se por causa da condição de estabilidade, por isso esse método não é recomendado para modelos com alto <math>\gamma</math>. Para tais modelos, o método FTCS implícito é mais recomendado por ser incondicionalmente estável.<br />
<br />
== Implementção ==<br />
<br />
<source lang="python"><br />
<br />
import matplotlib.pyplot as plt <br />
<br />
def vector_declaration(L,dx):<br />
c = [[],[]] # vetor concentração<br />
espaco = []<br />
# Condições iniciais<br />
for i in range(int(L/dx)+4):##+4 pois usaremos dois valores antes e depois do ultimo elemento do vetor c<br />
if (i<1/2*L/dx+2):<br />
c[0].append(-1)<br />
c[1].append(-1)<br />
else:<br />
c[0].append(1)<br />
c[1].append(1)<br />
espaco.append(round(i/150,3))<br />
return c, espaco<br />
<br />
def CH_equation(gamma, D, dx, dt, L, TEMPO_MAX): # resolução numérica da equação<br />
<br />
c, espaco = vector_declaration(L, dx)<br />
i = 0<br />
for time in [t*dt/TEMPO_MAX for t in range(int(TEMPO_MAX/dt))]:<br />
for l in range(2,len(c[1][2:-2])):<br />
c1 = c[i][l-1]**3 - 2*c[i][l]**3 + c[i][l+1]**3<br />
c2 = -c[i][l-1] + 2*c[i][l] - c[i][l+1]<br />
c3 = c[i][l-2] -4*c[i][l-1] + 6*c[i][l] - 4*c[i][l+1] + c[i][l+2]<br />
c[1-i][l] = D*dt/(dx**2)*(c1+c2-gamma**2*c3/(dx**2)) + c[i][l]<br />
i = 1-i<br />
<br />
return(espaco[2:-2],c[1-i][2:-2]) ##retirando os elementos a mais do vetor<br />
<br />
tempo=1<br />
tamanho=1/3.5<br />
Difuse= 1<br />
gamma=3.4*1/128<br />
dt=1/500000<br />
dx = 1/128<br />
<br />
figure, axis = plt.subplots(2, 2)<br />
plt.figure(dpi=500)<br />
<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
plt.show()<br />
plt.savefig('graficosDIFUSAO.png', dpi=500)<br />
<br />
</source><br />
<br />
== Referências ==<br />
<br />
[1] SIBBING, Zimo. '''Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions''', tese de bacharelado, 2015.<br />
<br />
[2] https://en.wikipedia.org/wiki/Spinodal_decomposition<br />
<br />
[3] MARKUS, Wilczek. '''The Cahn-Hilliard Equation''', 2015.<br />
<br />
[4] CAHN, John W.; HILLIARD, John E. '''Free Energy of a Nonuniform System. I. Interfacial Free Energy'''. The Journal of Chemical Physics, 1958.<br />
<br />
[5] CAHN, John W.; HILLIARD, John E. '''Spinodal decomposition: A reprise'''Acta Metallurgica, Volume 19, Issue 2, 1971<br />
<br />
[6] https://pt.wikipedia.org/wiki/Lei_de_Fick<br />
<br />
[7] https://pt.qaz.wiki/wiki/Cahn%E2%80%93Hilliard_equation<br />
<br />
[8] Daniel V. Schroeder. An Introduction to Thermal Physics. Addison-Wesley, 1999.<br />
<br />
[9] Dongsun Lee, Joo-Youl Huh, Darae Jeong, Jaemin Shin, Ana Yun, and Junseok Kim. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Computational Materials Science, 2014.<br />
<br />
[10] https://en.wikipedia.org/wiki/Neumann_boundary_condition</div>Arthurdornhttp://fiscomp.if.ufrgs.br/index.php?title=Equa%C3%A7%C3%A3o_de_Cahn-Hilliard&diff=4007Equação de Cahn-Hilliard2021-04-05T21:21:16Z<p>Arthurdorn: /* Decomposição Espinodal */</p>
<hr />
<div>'''Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David e Guilherme Hoss'''<br />
<br />
O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição espinodal de uma mistura binária, e analisar como é seu comportamento com diferentes coeficientes de difusão, utilizando o método FTCS (''Forward Time Centered Space''). O trabalho foi inspirado no artigo de Sibbing[1].<br />
<br />
== Decomposição Espinodal ==<br />
<br />
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 [2].<br />
<br />
==A Equação de Cahn-Hilliard==<br />
<br />
A equação de Cahn-Hilliard é 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.<br />
Com o intuito de deduzirmos essa equação, consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos <math> c_a(x,t) </math> e <math> c_b(x,t) </math>, respectivamente. [1]<br />
<br />
Além disso, podemos considerar que - para uma mistura binária - <math> c_a(x,t) + c_b(x,t) = 1</math> e portanto podemos simplificar para apenas uma concentração <math> c(x,t) </math>:<br />
:<math><br />
c_a(x,t) = c(x,t), c_b(x,t) = 1 - c(x,t).<br />
</math><br />
<br />
Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:<br />
:<math><br />
J = -D\nabla c.<br />
</math><br />
<br />
juntamente da equação da continuidade:<br />
:<math><br />
\frac{\partial c}{\partial t} + \nabla \cdot \vec J = 0.<br />
</math><br />
<br />
Onde <math>D</math> é o coeficiente de difusão e <math>\vec J</math> é a difusão da mistura. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 c.<br />
</math><br />
<br />
A partir dessa equação - como não há a existência de um gradiente de concentração espacial - pode-se esperar que não ocorra mudança na concentração da mistura. No entanto, observa-se que quando a separação de fases ocorre, a difusão demonstra ser contrária ao gradiente de concentração, o que não condiz com a equação anterior.<br />
Tendo isso em vista, conclui-se que a concentração não pode ser a razão da difusão, portanto outra força deve estar presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico (de acordo com Cahn e Hilliard, 1958). Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:<br />
<br />
:<math><br />
J = -M \nabla \mu.<br />
</math><br />
<br />
Onde <math>M</math> é a mobilidade das partículas (análoga à D) e <math>\mu</math> é o potencial químico.<br />
Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick [3]:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 \mu.<br />
</math><br />
<br />
Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c}.<br />
</math><br />
<br />
Onde <math>g</math> é a densidade da energia livre de Gibbs e <math>c</math> é a concentração (de acordo com Schroeder, 1999).<br />
<br />
Tendo em vista a substituição do termo <math>\mu</math> 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 (derivado em [3]):<br />
<br />
:<math><br />
G = \int_{V}^{} f(c) + {\kappa |\nabla c|}^2 dV.<br />
</math><br />
<br />
Nesse caso, <math>G</math> é a energia livre de Gibbs, <math>f(c)</math> é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e <math>{\kappa|\nabla c|}^2</math> é a densidade de energia livre devido ao gradiente de concentração.<br />
<br />
Além disso, a função <math>f(c)</math> - de acordo com [8] - possui o potencial de um poço de potencial duplo. Neste poço, <math>c</math> representa a concentração em escala e está relacionada à temperatura da mistura, que decide se a separação de fases irá - ou não - ocorrer. Esta função pode ser representada pela seguinte equação:<br />
<br />
:<math><br />
f(c) = \frac{(c^2 - 1)^2}{4}.<br />
</math><br />
<br />
Levando essas informações em conta e - utilizando um parâmetro <math>\gamma</math> análogo à largura da interface - que é descrito por <math>\kappa = \gamma ^2</math> é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema com duas fases.<br />
:<math><br />
g(c) = f(c) + {\gamma}^2{|\nabla c |}^2.<br />
</math><br />
<br />
Com essas igualdades agora se torna possível o cálculo de <math>\mu</math> em função da concentração dos fluidos:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c} = \frac{\partial f(c)}{\partial c} + \frac{\partial }{\partial c} (\gamma^2{|\nabla c |}^2) = \frac{\partial }{\partial c} \frac{(c^2 - 1)^2}{4} + {\gamma}^2 {\nabla}^2 c = c^3 - c + {\gamma}^2 {\nabla}^2 c.<br />
</math><br />
<br />
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 (<math>M</math>) e a concentração do fluido:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
Essa equação final é chamada de equação de Cahn-Hilliard. <br />
A equação dependente da difusão é análoga e também funcional:<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
== Método FTCS (Forward Time Centered Space) ==<br />
<br />
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, centrado no espaço". Uma das formas que o método pode ser utilizado é a forma explícita que está descritas abaixo.<br />
<br />
:<math>n\to\Delta t</math><br />
<br />
:<math>j\to\Delta x</math><br />
<br />
===FTCS Explicito===<br />
<br />
:<math>\frac{\partial f}{\partial t}\to \frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}</math><br />
<br />
:<math>\frac{\partial ^2 f}{\partial x^2}\to \frac{f_{j-1}^{n}-2 f_{j}^{n} + f_{j+1}^n}{\Delta x^2}</math><br />
<br />
Para difusão:<br />
<br />
:<math>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})</math><br />
<br />
== Implementação da equação de Cahn-Hilliard 1D pelo método FTCS explicito ==<br />
<br />
<br />
:<math>\displaystyle \frac{\partial c}{\partial t} = D\nabla^{2}(c^{3}-c-\gamma^2\nabla^{2}c)</math> <br />
<br />
:<math>\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})</math><br />
<br />
:<math>\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}</math><br />
<br />
<br />
:<math>\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-2 c_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)</math><br />
<br />
:<math> 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-2 c_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</math><br />
<br />
:<math>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+2 c_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 </math><br />
<br />
'''Condição de Estabilidade '''<br />
<br />
A estabilidade dessa equação mostra-se muito mais complicada de se estipular por ela ser uma equação diferencial de quarta ordem se comparada a equação de difusão, Portanto só iremos analisar a seção 3.3 (Experimental and theoretical stability conditions) do artigo Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions [1]. <br />
<br />
Após a linearização e aplicação do teorema de Gershgorin temos que a condição para estabilidade da equação linear para <math>D = 1</math> é:<br />
<br />
<br />
:<center><math>\Delta t < \displaystyle\frac{(\Delta x)^2}{4+\frac{8\gamma^2}{\Delta x^2}}</math> <br />
<br />
Importante atentar que essa é a condição de estabilidade somente para a equação de Cahn-Hilliard linearizada, não para a original. Tanto que a literatura sobre a equação propõem que <math>\Delta t</math> <math>\varpropto</math> <math>(\Delta x)^4 </math>, que é o que acontece na condição estabilidade linear quando <math>\gamma >> \Delta x</math>.<br />
<br />
O artigo compara os dados experimentais de estabilidade com a estabilidade da equação linearizada relacionado na seção 3.3.4 e conclui que para valores de <math>\frac{\gamma}{\Delta x} \in [0.25,8]</math> a condição teórica encontrada a partir da linearização é uma boa aproximação.<br />
<br />
== Resultados ==<br />
<br />
Com o intuito de testar como o fator de difusão D afeta a evolução da equação de Cahn-Hilliard, comparamos os resultados para os coeficientes de difusão 1, 0.1, 0.01 e 0.001 e analisamos seus gráficos. Para esse estudo foram utilizadas como condições iniciais a concentração -1 para a primeira metade do espaço e 1 para a segunda metade, e o sistema é livre de condições de contorno. <br />
<br />
<br />
[[Arquivo:coeficientes.png|1000px|thumb|center| Gráfico 1: Resultados da simulação variando os coeficientes de difusão (D) para tempos máximos diferentes, <math>\gamma=3.4/128,\Delta t = 1/500000</math>]]<br />
<br />
<br />
Nos gráficos, é possível observar que quanto maior o coeficiente de difusão maior é a velocidade em que a mistura atinge a estabilidade. Além disso, vemos que valores baixos de t produzem soluções mais íngremes que valores altos de t.<br />
<br />
== Discussão de Resultados ==<br />
<br />
Uma propriedade observada no gráfico 1 é a de que os valores das soluções obtidas utilizando o método FTCS excedem os valores máximos e mínimos permitidos ( C=1 e C=-1), se estivesse modelando uma situação real isso iria contra a lei de conservação de massa, o que pode ocasionar erros nos resultados que exigem uma grande precisão. <br />
<br />
O método FTCS explícito limita-se por causa da condição de estabilidade, por isso esse método não é recomendado para modelos com alto <math>\gamma</math>. Para tais modelos, o método FTCS implícito é mais recomendado por ser incondicionalmente estável.<br />
<br />
== Implementção ==<br />
<br />
<source lang="python"><br />
<br />
import matplotlib.pyplot as plt <br />
<br />
def vector_declaration(L,dx):<br />
c = [[],[]] # vetor concentração<br />
espaco = []<br />
# Condições iniciais<br />
for i in range(int(L/dx)+4):##+4 pois usaremos dois valores antes e depois do ultimo elemento do vetor c<br />
if (i<1/2*L/dx+2):<br />
c[0].append(-1)<br />
c[1].append(-1)<br />
else:<br />
c[0].append(1)<br />
c[1].append(1)<br />
espaco.append(round(i/150,3))<br />
return c, espaco<br />
<br />
def CH_equation(gamma, D, dx, dt, L, TEMPO_MAX): # resolução numérica da equação<br />
<br />
c, espaco = vector_declaration(L, dx)<br />
i = 0<br />
for time in [t*dt/TEMPO_MAX for t in range(int(TEMPO_MAX/dt))]:<br />
for l in range(2,len(c[1][2:-2])):<br />
c1 = c[i][l-1]**3 - 2*c[i][l]**3 + c[i][l+1]**3<br />
c2 = -c[i][l-1] + 2*c[i][l] - c[i][l+1]<br />
c3 = c[i][l-2] -4*c[i][l-1] + 6*c[i][l] - 4*c[i][l+1] + c[i][l+2]<br />
c[1-i][l] = D*dt/(dx**2)*(c1+c2-gamma**2*c3/(dx**2)) + c[i][l]<br />
i = 1-i<br />
<br />
return(espaco[2:-2],c[1-i][2:-2]) ##retirando os elementos a mais do vetor<br />
<br />
tempo=1<br />
tamanho=1/3.5<br />
Difuse= 1<br />
gamma=3.4*1/128<br />
dt=1/500000<br />
dx = 1/128<br />
<br />
figure, axis = plt.subplots(2, 2)<br />
plt.figure(dpi=500)<br />
<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
plt.show()<br />
plt.savefig('graficosDIFUSAO.png', dpi=500)<br />
<br />
</source><br />
<br />
== Referências ==<br />
<br />
[1] SIBBING, Zimo. '''Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions''', tese de bacharelado, 2015.<br />
<br />
[2] https://en.wikipedia.org/wiki/Spinodal_decomposition<br />
<br />
[3] MARKUS, Wilczek. '''The Cahn-Hilliard Equation''', 2015.<br />
<br />
[4] CAHN, John W.; HILLIARD, John E. '''Free Energy of a Nonuniform System. I. Interfacial Free Energy'''. The Journal of Chemical Physics, 1958.<br />
<br />
[5] CAHN, John W.; HILLIARD, John E. '''Spinodal decomposition: A reprise'''Acta Metallurgica, Volume 19, Issue 2, 1971<br />
<br />
[6] https://pt.wikipedia.org/wiki/Lei_de_Fick<br />
<br />
[7] https://pt.qaz.wiki/wiki/Cahn%E2%80%93Hilliard_equation<br />
<br />
[8] Daniel V. Schroeder. An Introduction to Thermal Physics. Addison-Wesley, 1999.<br />
<br />
[9] Dongsun Lee, Joo-Youl Huh, Darae Jeong, Jaemin Shin, Ana Yun, and Junseok Kim. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Computational Materials Science, 2014.</div>Arthurdornhttp://fiscomp.if.ufrgs.br/index.php?title=Equa%C3%A7%C3%A3o_de_Cahn-Hilliard&diff=4006Equação de Cahn-Hilliard2021-04-05T21:21:07Z<p>Arthurdorn: /* Referências */</p>
<hr />
<div>'''Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David e Guilherme Hoss'''<br />
<br />
O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição espinodal de uma mistura binária, e analisar como é seu comportamento com diferentes coeficientes de difusão, utilizando o método FTCS (''Forward Time Centered Space''). O trabalho foi inspirado no artigo de Sibbing[1].<br />
<br />
== Decomposição Espinodal ==<br />
<br />
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 [6].<br />
<br />
==A Equação de Cahn-Hilliard==<br />
<br />
A equação de Cahn-Hilliard é 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.<br />
Com o intuito de deduzirmos essa equação, consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos <math> c_a(x,t) </math> e <math> c_b(x,t) </math>, respectivamente. [1]<br />
<br />
Além disso, podemos considerar que - para uma mistura binária - <math> c_a(x,t) + c_b(x,t) = 1</math> e portanto podemos simplificar para apenas uma concentração <math> c(x,t) </math>:<br />
:<math><br />
c_a(x,t) = c(x,t), c_b(x,t) = 1 - c(x,t).<br />
</math><br />
<br />
Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:<br />
:<math><br />
J = -D\nabla c.<br />
</math><br />
<br />
juntamente da equação da continuidade:<br />
:<math><br />
\frac{\partial c}{\partial t} + \nabla \cdot \vec J = 0.<br />
</math><br />
<br />
Onde <math>D</math> é o coeficiente de difusão e <math>\vec J</math> é a difusão da mistura. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 c.<br />
</math><br />
<br />
A partir dessa equação - como não há a existência de um gradiente de concentração espacial - pode-se esperar que não ocorra mudança na concentração da mistura. No entanto, observa-se que quando a separação de fases ocorre, a difusão demonstra ser contrária ao gradiente de concentração, o que não condiz com a equação anterior.<br />
Tendo isso em vista, conclui-se que a concentração não pode ser a razão da difusão, portanto outra força deve estar presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico (de acordo com Cahn e Hilliard, 1958). Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:<br />
<br />
:<math><br />
J = -M \nabla \mu.<br />
</math><br />
<br />
Onde <math>M</math> é a mobilidade das partículas (análoga à D) e <math>\mu</math> é o potencial químico.<br />
Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick [3]:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 \mu.<br />
</math><br />
<br />
Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c}.<br />
</math><br />
<br />
Onde <math>g</math> é a densidade da energia livre de Gibbs e <math>c</math> é a concentração (de acordo com Schroeder, 1999).<br />
<br />
Tendo em vista a substituição do termo <math>\mu</math> 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 (derivado em [3]):<br />
<br />
:<math><br />
G = \int_{V}^{} f(c) + {\kappa |\nabla c|}^2 dV.<br />
</math><br />
<br />
Nesse caso, <math>G</math> é a energia livre de Gibbs, <math>f(c)</math> é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e <math>{\kappa|\nabla c|}^2</math> é a densidade de energia livre devido ao gradiente de concentração.<br />
<br />
Além disso, a função <math>f(c)</math> - de acordo com [8] - possui o potencial de um poço de potencial duplo. Neste poço, <math>c</math> representa a concentração em escala e está relacionada à temperatura da mistura, que decide se a separação de fases irá - ou não - ocorrer. Esta função pode ser representada pela seguinte equação:<br />
<br />
:<math><br />
f(c) = \frac{(c^2 - 1)^2}{4}.<br />
</math><br />
<br />
Levando essas informações em conta e - utilizando um parâmetro <math>\gamma</math> análogo à largura da interface - que é descrito por <math>\kappa = \gamma ^2</math> é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema com duas fases.<br />
:<math><br />
g(c) = f(c) + {\gamma}^2{|\nabla c |}^2.<br />
</math><br />
<br />
Com essas igualdades agora se torna possível o cálculo de <math>\mu</math> em função da concentração dos fluidos:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c} = \frac{\partial f(c)}{\partial c} + \frac{\partial }{\partial c} (\gamma^2{|\nabla c |}^2) = \frac{\partial }{\partial c} \frac{(c^2 - 1)^2}{4} + {\gamma}^2 {\nabla}^2 c = c^3 - c + {\gamma}^2 {\nabla}^2 c.<br />
</math><br />
<br />
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 (<math>M</math>) e a concentração do fluido:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
Essa equação final é chamada de equação de Cahn-Hilliard. <br />
A equação dependente da difusão é análoga e também funcional:<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
== Método FTCS (Forward Time Centered Space) ==<br />
<br />
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, centrado no espaço". Uma das formas que o método pode ser utilizado é a forma explícita que está descritas abaixo.<br />
<br />
:<math>n\to\Delta t</math><br />
<br />
:<math>j\to\Delta x</math><br />
<br />
===FTCS Explicito===<br />
<br />
:<math>\frac{\partial f}{\partial t}\to \frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}</math><br />
<br />
:<math>\frac{\partial ^2 f}{\partial x^2}\to \frac{f_{j-1}^{n}-2 f_{j}^{n} + f_{j+1}^n}{\Delta x^2}</math><br />
<br />
Para difusão:<br />
<br />
:<math>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})</math><br />
<br />
== Implementação da equação de Cahn-Hilliard 1D pelo método FTCS explicito ==<br />
<br />
<br />
:<math>\displaystyle \frac{\partial c}{\partial t} = D\nabla^{2}(c^{3}-c-\gamma^2\nabla^{2}c)</math> <br />
<br />
:<math>\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})</math><br />
<br />
:<math>\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}</math><br />
<br />
<br />
:<math>\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-2 c_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)</math><br />
<br />
:<math> 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-2 c_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</math><br />
<br />
:<math>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+2 c_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 </math><br />
<br />
'''Condição de Estabilidade '''<br />
<br />
A estabilidade dessa equação mostra-se muito mais complicada de se estipular por ela ser uma equação diferencial de quarta ordem se comparada a equação de difusão, Portanto só iremos analisar a seção 3.3 (Experimental and theoretical stability conditions) do artigo Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions [1]. <br />
<br />
Após a linearização e aplicação do teorema de Gershgorin temos que a condição para estabilidade da equação linear para <math>D = 1</math> é:<br />
<br />
<br />
:<center><math>\Delta t < \displaystyle\frac{(\Delta x)^2}{4+\frac{8\gamma^2}{\Delta x^2}}</math> <br />
<br />
Importante atentar que essa é a condição de estabilidade somente para a equação de Cahn-Hilliard linearizada, não para a original. Tanto que a literatura sobre a equação propõem que <math>\Delta t</math> <math>\varpropto</math> <math>(\Delta x)^4 </math>, que é o que acontece na condição estabilidade linear quando <math>\gamma >> \Delta x</math>.<br />
<br />
O artigo compara os dados experimentais de estabilidade com a estabilidade da equação linearizada relacionado na seção 3.3.4 e conclui que para valores de <math>\frac{\gamma}{\Delta x} \in [0.25,8]</math> a condição teórica encontrada a partir da linearização é uma boa aproximação.<br />
<br />
== Resultados ==<br />
<br />
Com o intuito de testar como o fator de difusão D afeta a evolução da equação de Cahn-Hilliard, comparamos os resultados para os coeficientes de difusão 1, 0.1, 0.01 e 0.001 e analisamos seus gráficos. Para esse estudo foram utilizadas como condições iniciais a concentração -1 para a primeira metade do espaço e 1 para a segunda metade, e o sistema é livre de condições de contorno. <br />
<br />
<br />
[[Arquivo:coeficientes.png|1000px|thumb|center| Gráfico 1: Resultados da simulação variando os coeficientes de difusão (D) para tempos máximos diferentes, <math>\gamma=3.4/128,\Delta t = 1/500000</math>]]<br />
<br />
<br />
Nos gráficos, é possível observar que quanto maior o coeficiente de difusão maior é a velocidade em que a mistura atinge a estabilidade. Além disso, vemos que valores baixos de t produzem soluções mais íngremes que valores altos de t.<br />
<br />
== Discussão de Resultados ==<br />
<br />
Uma propriedade observada no gráfico 1 é a de que os valores das soluções obtidas utilizando o método FTCS excedem os valores máximos e mínimos permitidos ( C=1 e C=-1), se estivesse modelando uma situação real isso iria contra a lei de conservação de massa, o que pode ocasionar erros nos resultados que exigem uma grande precisão. <br />
<br />
O método FTCS explícito limita-se por causa da condição de estabilidade, por isso esse método não é recomendado para modelos com alto <math>\gamma</math>. Para tais modelos, o método FTCS implícito é mais recomendado por ser incondicionalmente estável.<br />
<br />
== Implementção ==<br />
<br />
<source lang="python"><br />
<br />
import matplotlib.pyplot as plt <br />
<br />
def vector_declaration(L,dx):<br />
c = [[],[]] # vetor concentração<br />
espaco = []<br />
# Condições iniciais<br />
for i in range(int(L/dx)+4):##+4 pois usaremos dois valores antes e depois do ultimo elemento do vetor c<br />
if (i<1/2*L/dx+2):<br />
c[0].append(-1)<br />
c[1].append(-1)<br />
else:<br />
c[0].append(1)<br />
c[1].append(1)<br />
espaco.append(round(i/150,3))<br />
return c, espaco<br />
<br />
def CH_equation(gamma, D, dx, dt, L, TEMPO_MAX): # resolução numérica da equação<br />
<br />
c, espaco = vector_declaration(L, dx)<br />
i = 0<br />
for time in [t*dt/TEMPO_MAX for t in range(int(TEMPO_MAX/dt))]:<br />
for l in range(2,len(c[1][2:-2])):<br />
c1 = c[i][l-1]**3 - 2*c[i][l]**3 + c[i][l+1]**3<br />
c2 = -c[i][l-1] + 2*c[i][l] - c[i][l+1]<br />
c3 = c[i][l-2] -4*c[i][l-1] + 6*c[i][l] - 4*c[i][l+1] + c[i][l+2]<br />
c[1-i][l] = D*dt/(dx**2)*(c1+c2-gamma**2*c3/(dx**2)) + c[i][l]<br />
i = 1-i<br />
<br />
return(espaco[2:-2],c[1-i][2:-2]) ##retirando os elementos a mais do vetor<br />
<br />
tempo=1<br />
tamanho=1/3.5<br />
Difuse= 1<br />
gamma=3.4*1/128<br />
dt=1/500000<br />
dx = 1/128<br />
<br />
figure, axis = plt.subplots(2, 2)<br />
plt.figure(dpi=500)<br />
<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
plt.show()<br />
plt.savefig('graficosDIFUSAO.png', dpi=500)<br />
<br />
</source><br />
<br />
== Referências ==<br />
<br />
[1] SIBBING, Zimo. '''Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions''', tese de bacharelado, 2015.<br />
<br />
[2] https://en.wikipedia.org/wiki/Spinodal_decomposition<br />
<br />
[3] MARKUS, Wilczek. '''The Cahn-Hilliard Equation''', 2015.<br />
<br />
[4] CAHN, John W.; HILLIARD, John E. '''Free Energy of a Nonuniform System. I. Interfacial Free Energy'''. The Journal of Chemical Physics, 1958.<br />
<br />
[5] CAHN, John W.; HILLIARD, John E. '''Spinodal decomposition: A reprise'''Acta Metallurgica, Volume 19, Issue 2, 1971<br />
<br />
[6] https://pt.wikipedia.org/wiki/Lei_de_Fick<br />
<br />
[7] https://pt.qaz.wiki/wiki/Cahn%E2%80%93Hilliard_equation<br />
<br />
[8] Daniel V. Schroeder. An Introduction to Thermal Physics. Addison-Wesley, 1999.<br />
<br />
[9] Dongsun Lee, Joo-Youl Huh, Darae Jeong, Jaemin Shin, Ana Yun, and Junseok Kim. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Computational Materials Science, 2014.</div>Arthurdornhttp://fiscomp.if.ufrgs.br/index.php?title=Equa%C3%A7%C3%A3o_de_Cahn-Hilliard&diff=4001Equação de Cahn-Hilliard2021-04-05T21:09:53Z<p>Arthurdorn: </p>
<hr />
<div>'''Grupo: Arthur Dornelles, Bruno Zanette, Gabriel De David e Guilherme Hoss'''<br />
<br />
O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, que descreve o processo de decomposição espinodal de uma mistura binária, e analisar como é seu comportamento com diferentes coeficientes de difusão, utilizando o método FTCS (''Forward Time Centered Space''). O trabalho foi inspirado no artigo de Sibbing[1].<br />
<br />
== Decomposição Espinodal ==<br />
<br />
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 [5].<br />
<br />
==A Equação de Cahn-Hilliard==<br />
<br />
Como dito previamente, a equação de Cahn-Hilliard descreve o processo de decomposição espinodal de uma mistura binária.<br />
Com o intuito de deduzirmos essa equação, consideraremos - de início - uma mistura binária de dois componentes A e B descritas pelas concentrações dos fluidos <math> c_a(x,t) </math> e <math> c_b(x,t) </math>, respectivamente. [1]<br />
<br />
Além disso, podemos considerar que - para uma mistura binária - <math> c_a(x,t) + c_b(x,t) = 1</math> e portanto podemos simplificar para apenas uma concentração <math> c(x,t) </math>:<br />
:<math><br />
c_a(x,t) = c(x,t), c_b(x,t) = 1 - c(x,t).<br />
</math><br />
<br />
Tendo isso em vista, podemos agora utilizar a primeira lei de Fick da difusão:<br />
:<math><br />
J = -D\nabla c.<br />
</math><br />
<br />
juntamente da equação da continuidade:<br />
:<math><br />
\frac{\partial c}{\partial t} + \nabla \cdot \vec J = 0.<br />
</math><br />
<br />
Onde <math>D</math> é o coeficiente de difusão e <math>\vec J</math> é a difusão da mistura. Em seguida, ao combinarmos ambas as equações anteriores o resultado gera a segunda lei de Fick da difusão:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 c.<br />
</math><br />
<br />
A partir dessa equação - como não há a existência de um gradiente de concentração espacial - pode-se esperar que não ocorra mudança na concentração da mistura. No entanto, observa-se que quando a separação de fases ocorre, a difusão demonstra ser contrária ao gradiente de concentração, o que não condiz com a equação anterior.<br />
Tendo isso em vista, conclui-se que a concentração não pode ser a razão da difusão, portanto outra força deve estar presente. E, nesse caso, encontrou-se que a principal força responsável pela difusão negativa é o potencial químico (de acordo com Cahn e Hilliard, 1958). Portanto, outra equação pode ser derivada para generalizar a primeira lei de Fick:<br />
<br />
:<math><br />
J = -M \nabla \mu.<br />
</math><br />
<br />
Onde <math>M</math> é a mobilidade das partículas (análoga à D) e <math>\mu</math> é o potencial químico.<br />
Com essa nova equação podemos agora também deduzir uma nova equação para a segunda lei de Fick [3]:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 \mu.<br />
</math><br />
<br />
Nessa equação, podemos usar a definição do potencial químico através da densidade da energia livre de Gibbs como:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c}.<br />
</math><br />
<br />
Onde <math>g</math> é a densidade da energia livre de Gibbs e <math>c</math> é a concentração (de acordo com Schroeder, 1999).<br />
<br />
Tendo em vista a substituição do termo <math>\mu</math> 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 (derivado em [3]):<br />
<br />
:<math><br />
G = \int_{V}^{} f(c) + {\kappa |\nabla c|}^2 dV.<br />
</math><br />
<br />
Nesse caso, <math>G</math> é a energia livre de Gibbs, <math>f(c)</math> é a densidade de energia livre devido à contribuições de ambas as fases homogêneas e <math>{\kappa|\nabla c|}^2</math> é a densidade de energia livre devido ao gradiente de concentração.<br />
<br />
Além disso, a função <math>f(c)</math> - de acordo com [8] - possui o potencial de um poço de potencial duplo. Neste poço, <math>c</math> representa a concentração em escala e está relacionada à temperatura da mistura, que decide se a separação de fases irá - ou não - ocorrer. Esta função pode ser representada pela seguinte equação:<br />
<br />
:<math><br />
f(c) = \frac{(c^2 - 1)^2}{4}.<br />
</math><br />
<br />
Levando essas informações em conta e - utilizando um parâmetro <math>\gamma</math> análogo à largura da interface - que é descrito por <math>\kappa = \gamma ^2</math> é possível encontrar uma equação que descreve a densidade de energia livre de Gibbs para um sistema com duas fases.<br />
:<math><br />
g(c) = f(c) + {\gamma}^2{|\nabla c |}^2.<br />
</math><br />
<br />
Com essas igualdades agora se torna possível o cálculo de <math>\mu</math> em função da concentração dos fluidos:<br />
<br />
:<math><br />
\mu = \frac{\partial g}{\partial c} = \frac{\partial f(c)}{\partial c} + \frac{\partial }{\partial c} (\gamma^2{|\nabla c |}^2) = \frac{\partial }{\partial c} \frac{(c^2 - 1)^2}{4} + {\gamma}^2 {\nabla}^2 c = c^3 - c + {\gamma}^2 {\nabla}^2 c.<br />
</math><br />
<br />
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 (<math>M</math>) e a concentração do fluido:<br />
<br />
:<math><br />
\frac{\partial c}{\partial t} = M {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
Essa equação final é chamada de equação de Cahn-Hilliard. <br />
A equação dependente da difusão é análoga e também funcional:<br />
:<math><br />
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c).<br />
</math><br />
<br />
== Método FTCS (Forward Time Centered Space) ==<br />
<br />
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, centrado no espaço". Uma das formas que o método pode ser utilizado é a forma explícita que está descritas abaixo.<br />
<br />
:<math>n\to\Delta t</math><br />
<br />
:<math>j\to\Delta x</math><br />
<br />
===FTCS Explicito===<br />
<br />
:<math>\frac{\partial f}{\partial t}\to \frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}</math><br />
<br />
:<math>\frac{\partial ^2 f}{\partial x^2}\to \frac{f_{j-1}^{n}-2 f_{j}^{n} + f_{j+1}^n}{\Delta x^2}</math><br />
<br />
Para difusão:<br />
<br />
:<math>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})</math><br />
<br />
== Implementação da equação de Cahn-Hilliard 1D pelo método FTCS explicito ==<br />
<br />
<br />
:<math>\displaystyle \frac{\partial c}{\partial t} = D\nabla^{2}(c^{3}-c-\gamma^2\nabla^{2}c)</math> <br />
<br />
:<math>\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})</math><br />
<br />
:<math>\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}</math><br />
<br />
<br />
:<math>\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-2 c_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)</math><br />
<br />
:<math> 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-2 c_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</math><br />
<br />
:<math>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+2 c_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 </math><br />
<br />
'''Condição de Estabilidade '''<br />
<br />
A estabilidade dessa equação mostra-se muito mais complicada de se estipular por ela ser uma equação diferencial de quarta ordem se comparada a equação de difusão, Portanto só iremos analisar a seção 3.3 (Experimental and theoretical stability conditions) do artigo Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions [1]. <br />
<br />
Após a linearização e aplicação do teorema de Gershgorin temos que a condição para estabilidade da equação linear para <math>D = 1</math> é:<br />
<br />
<br />
:<center><math>\Delta t < \displaystyle\frac{(\Delta x)^2}{4+\frac{8\gamma^2}{\Delta x^2}}</math> <br />
<br />
Importante atentar que essa é a condição de estabilidade somente para a equação de Cahn-Hilliard linearizada, não para a original. Tanto que a literatura sobre a equação propõem que <math>\Delta t</math> <math>\varpropto</math> <math>(\Delta x)^4 </math>, que é o que acontece na condição estabilidade linear quando <math>\gamma >> \Delta x</math>.<br />
<br />
O artigo compara os dados experimentais de estabilidade com a estabilidade da equação linearizada relacionado na seção 3.3.4 e conclui que para valores de <math>\frac{\gamma}{\Delta x} \in [0.25,8]</math> a condição teórica encontrada a partir da linearização é uma boa aproximação.<br />
<br />
== Resultados ==<br />
<br />
Com o intuito de testar como o fator de difusão D afeta a evolução da equação de Cahn-Hilliard, comparamos os resultados para os coeficientes de difusão 1, 0.1, 0.01 e 0.001 e analisamos seus gráficos. Para esse estudo foram utilizadas como condições iniciais a concentração -1 para a primeira metade do espaço e 1 para a segunda metade, e o sistema é livre de condições de contorno. <br />
<br />
<br />
[[Arquivo:coeficientes.png|1000px|thumb|center| Gráfico 1: Resultados da simulação variando os coeficientes de difusão (D) para tempos máximos diferentes, <math>\gamma=3.4/128,\Delta t = 1/500000</math>]]<br />
<br />
<br />
Nos gráficos, é possível observar que quanto maior o coeficiente de difusão maior é a velocidade em que a mistura atinge a estabilidade. Além disso, vemos que valores baixos de t produzem soluções mais íngremes que valores altos de t.<br />
<br />
== Discussão de Resultados ==<br />
<br />
Uma propriedade observada no gráfico 1 é a de que os valores das soluções obtidas utilizando o método FTCS excedem os valores máximos e mínimos permitidos ( C=1 e C=-1), se estivesse modelando uma situação real isso iria contra a lei de conservação de massa, o que pode ocasionar erros nos resultados que exigem uma grande precisão. <br />
<br />
O método FTCS explícito limita-se por causa da condição de estabilidade, por isso esse método não é recomendado para modelos com alto <math>\gamma</math>. Para tais modelos, o método FTCS implícito é mais recomendado por ser incondicionalmente estável.<br />
<br />
== Implementção ==<br />
<br />
<source lang="python"><br />
<br />
import matplotlib.pyplot as plt <br />
<br />
def vector_declaration(L,dx):<br />
c = [[],[]] # vetor concentração<br />
espaco = []<br />
# Condições iniciais<br />
for i in range(int(L/dx)+4):##+4 pois usaremos dois valores antes e depois do ultimo elemento do vetor c<br />
if (i<1/2*L/dx+2):<br />
c[0].append(-1)<br />
c[1].append(-1)<br />
else:<br />
c[0].append(1)<br />
c[1].append(1)<br />
espaco.append(round(i/150,3))<br />
return c, espaco<br />
<br />
def CH_equation(gamma, D, dx, dt, L, TEMPO_MAX): # resolução numérica da equação<br />
<br />
c, espaco = vector_declaration(L, dx)<br />
i = 0<br />
for time in [t*dt/TEMPO_MAX for t in range(int(TEMPO_MAX/dt))]:<br />
for l in range(2,len(c[1][2:-2])):<br />
c1 = c[i][l-1]**3 - 2*c[i][l]**3 + c[i][l+1]**3<br />
c2 = -c[i][l-1] + 2*c[i][l] - c[i][l+1]<br />
c3 = c[i][l-2] -4*c[i][l-1] + 6*c[i][l] - 4*c[i][l+1] + c[i][l+2]<br />
c[1-i][l] = D*dt/(dx**2)*(c1+c2-gamma**2*c3/(dx**2)) + c[i][l]<br />
i = 1-i<br />
<br />
return(espaco[2:-2],c[1-i][2:-2]) ##retirando os elementos a mais do vetor<br />
<br />
tempo=1<br />
tamanho=1/3.5<br />
Difuse= 1<br />
gamma=3.4*1/128<br />
dt=1/500000<br />
dx = 1/128<br />
<br />
figure, axis = plt.subplots(2, 2)<br />
plt.figure(dpi=500)<br />
<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 0].plot(*CH_equation(gamma, Difuse, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[0, 1].plot(*CH_equation(gamma, Difuse/10, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[0, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 0].plot(*CH_equation(gamma, Difuse/100, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 0].legend(loc="upper left", prop={'size': 6})<br />
<br />
<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100000), label = "tempo: " + str(tempo/100000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10000), label = "tempo: " + str(tempo/10000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/1000), label = "tempo: " + str(tempo/1000))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/100), label = "tempo: " + str(tempo/100))<br />
axis[1, 1].plot(*CH_equation(gamma, Difuse/1000, dx, dt, tamanho, tempo/10), label = "tempo: " + str(tempo/10))<br />
axis[1, 1].legend(loc="upper left", prop={'size': 6})<br />
<br />
plt.show()<br />
plt.savefig('graficosDIFUSAO.png', dpi=500)<br />
<br />
</source><br />
<br />
== Referências ==<br />
<br />
[1] SIBBING, Zimo. '''Numerical methods for the implentation of the Cahn-Hilliard equation in one dimension and a dynamic boundary condition in two dimensions''', tese de bacharelado, 2015.<br />
<br />
[2] MARKUS, Wilczek. '''The Cahn-Hilliard Equation''', 2015.<br />
<br />
[3] CAHN, John W.; HILLIARD, John E. '''Free Energy of a Nonuniform System. I. Interfacial Free Energy'''. The Journal of Chemical Physics, 1958.<br />
<br />
[4] CAHN, John W.; HILLIARD, John E. '''Spinodal decomposition: A reprise'''Acta Metallurgica, Volume 19, Issue 2, 1971<br />
<br />
[5] https://pt.wikipedia.org/wiki/Lei_de_Fick<br />
<br />
[6] https://en.wikipedia.org/wiki/Spinodal_decomposition<br />
<br />
[7] https://pt.qaz.wiki/wiki/Cahn%E2%80%93Hilliard_equation<br />
<br />
[8] Daniel V. Schroeder. An Introduction to Thermal Physics. Addison-Wesley, 1999.<br />
<br />
[9] Dongsun Lee, Joo-Youl Huh, Darae Jeong, Jaemin Shin, Ana Yun, and Junseok Kim. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Computational Materials Science, 2014.</div>Arthurdorn