Acessibilidade / Reportar erro

Simulação tridimensional adaptativa da separação das fases de uma mistura bifásica usando a equação de Cahn-Hilliard

Resumos

Simulamos a separação dos componentes de uma mistura bifásica com a equação de Cahn-Hilliard. Esta equação contém intrincados termos não lineares e derivadas de alta ordem. Além disso, a delgada região de transição entre os componentes da mistura requer muita resolução. Assim, determinar a solução numérica da equação de Cahn-Hilliard não é uma tarefa fácil, principalmente em três dimensões. Conseguimos a resolução exigida no tempo usando uma discretização semi-implícita de segunda ordem. No espaço, obtemos a precisão requerida utilizando malhas refinadas localmente com a estratégia AMR. Essas malhas se adaptam dinamicamente para recobrir a região de transição. O sistema linear proveniente da discretização é solucionado por intermédio de técnicas multinível-multigrid.

Equação biharmônica; malhas adaptativas com refinamento local; métodos semi-implícitos; modelos de campo de fase; multigrid multinível


We simulate 3D spinodal decomposition modeled by the Cahn-Hilliard equation. This equation has intricate nonlinear terms and high order derivatives. Moreover, the thin transition region between the components of the mix requires high resolution. Thus, the computation of the Cahn-Hilliard equation is not an easy task, specially in three dimensions. We get the required resolution in time using a semi-implicit second order discretization scheme. In space, we obtain high accuracy utilizing meshes locally refined with the AMR strategy. These meshes adapt dynamically to cover the transition region. The linear system obtained from the discretization is solved by multilevel- multigrid techniques.

Biharmonic equation; adaptive mesh refinements; semi-implicit methods; phase-field models; multi-level multigrid


Simulação tridimensional adaptativa da separação das fases de uma mistura bifásica usando a equação de Cahn-Hilliard1 1 Parte desta pesquisa foi financiada pela National Science Foundation (projeto DMS 0609996), pela FAPESP (projetos 04/13781-1 e 06/57099-5) e pelo CNPq.

R.L. NósI; H.D. CenicerosII; A.M. RomaIII

IDepartamento Acadêmico de Matemática, UTFPR, 80230-901 Curitiba, PR, Brasil. rudimarnos@utfpr.edu.br

IIDepartment of Mathematics, UCSB, 93106 Santa Barbara, CA, USA. hdc@math.ucsb.edu

IIIDepartamento de Matemática Aplicada, IME, USP, 05311-970 São Paulo, SP, Brasil. roma@ime.usp.br

RESUMO

Simulamos a separação dos componentes de uma mistura bifásica com a equação de Cahn-Hilliard. Esta equação contém intrincados termos não lineares e derivadas de alta ordem. Além disso, a delgada região de transição entre os componentes da mistura requer muita resolução. Assim, determinar a solução numérica da equação de Cahn-Hilliard não é uma tarefa fácil, principalmente em três dimensões. Conseguimos a resolução exigida no tempo usando uma discretização semi-implícita de segunda ordem. No espaço, obtemos a precisão requerida utilizando malhas refinadas localmente com a estratégia AMR. Essas malhas se adaptam dinamicamente para recobrir a região de transição. O sistema linear proveniente da discretização é solucionado por intermédio de técnicas multinível-multigrid.

Palavras-chave: Equação biharmônica, malhas adaptativas com refinamento local, métodos semi-implícitos, modelos de campo de fase, multigrid multinível.

ABSTRACT

We simulate 3D spinodal decomposition modeled by the Cahn-Hilliard equation. This equation has intricate nonlinear terms and high order derivatives. Moreover, the thin transition region between the components of the mix requires high resolution. Thus, the computation of the Cahn-Hilliard equation is not an easy task, specially in three dimensions. We get the required resolution in time using a semi-implicit second order discretization scheme. In space, we obtain high accuracy utilizing meshes locally refined with the AMR strategy. These meshes adapt dynamically to cover the transition region. The linear system obtained from the discretization is solved by multilevel- multigrid techniques.

Keywords: Biharmonic equation, adaptive mesh refinements, semi-implicit methods, phase-field models, multi-level multigrid.

1. Introdução

A equação de Cahn-Hilliard é um modelo de campo de fase [3, 7, 9, 10, 11, 12, 14, 15, 16, 18, 20, 24] que permite simular a separação dos componentes de uma mistura bifásica. Os modelos de campo de fase constituem uma classe particular dos modelos de interface difusiva. Nesses modelos, as transições abruptas nas interfaces entre os diferentes fluidos ou materiais são substituídas por camadas delgadas nas quais as forças interfaciais são distribuídas suavemente. A idéia básica é introduzir um parâmetro de ordem ou campo de fase que descreve em cada instante o estado do fluido. Esse parâmetro de ordem varia continuamente sobre as finas camadas interfaciais, tornando-se mais uniforme no interior das fases.

Fundamentados nas pesquisas de Badalassi, Ceniceros e Banerjee [3] e de Ceniceros, Nós e Roma [12, 14, 20], apresentamos neste trabalho uma metodologia para efetuar simulações tridimensionais completamente adaptativas de um modelo de campo de fase descrito pela equação de Cahn-Hilliard. A metodologia numérica aplicada à solução dessa equação consiste no uso de uma discretização temporal linear semi-implícita de segunda ordem e de uma acurada discretização espacial em malhas refinadas localmente que se adaptam dinamicamente para recobrir a interface de transição. Isto garante precisão tanto no tempo quanto no espaço.

Na discretização temporal da equação de Cahn-Hilliard aplicamos o Método de Gear com coeficientes variáveis no tempo [12, 14, 20] e o sistema linear proveniente da discretização é solucionado com técnicas multinível-multigrid [1, 2, 8, 13, 14, 17, 19, 20, 21, 23]. A equação de Cahn-Hilliard é reescrita na forma de um sistema de equações diferenciais parciais para evitar a discretização direta do termo biharmônico e possibilitar o emprego de um multigrid linear. Quanto à discretização do domínio, utilizamos o refinamento local adaptativo (AMR - Adaptive Mesh Refinement) introduzido por Berger [4, 5, 6] na solução numérica de equações diferenciais parciais hiperbólicas. Optamos por essa técnica devido à sua eficiência e baixo custo computacional na remalhagem.

2. Modelo Matemático

2.1. Modelo de campo de fase

A energia livre do sistema permite avaliar o comportamento da separação dos componentes de uma mistura. No modelo de campo de fase que adotamos, essa energia é definida como um funcional de [3, 7, 9, 14, 15, 20], ou seja,

onde Ω é a região do espaço ocupada pelo sistema, é o termo que representa a energia interfacial, c > 0 é uma constante associada à espessura da interface de transição e f (Φ) é a densidade de energia típica. Esta tem dois mínimos que correspondem as duas fases estáveis do fluido.

A energia livre (2.1) está sujeita à restrição , sendo Φm o valor médio do campo de fase Φ e |Ω| o volume do domínio Ω.

A primeira derivada funcional da energia livre (2.1)

é o potencial químico, uma medida da tendência das partículas à difusão. Quando o potencial químico é constante, temos a minimização do funcional F (Φ) com respeito às variações de Φ. Essa minimização fornece o perfil de equilíbrio na interface.

Para a densidade de energia típica, usamos a função double well potential ("poço duplo") [3, 14, 15, 20]

com α = β = 1 (Helmholtz free energy). Assim, f' (Φ) = Φ3 - Φ. Com essa escolha, os valores 1 e -1 representam as duas fases da mistura e a região de transição é identificada pela variação dos valores de Φ de -1 a 1.

O perfil de equilíbrio em Rn é dado pelas soluções da equação

As soluções Φ± = ± 1 são estáveis, uniformes e representam as fases coexistentes. A terceira solução

é não uniforme, unidimensional na direção z e satisfaz as condições de fronteira .

A solução (2.2) [3, 15, 20] descreve o perfil de equilíbrio de uma interface plana normal à direção z, de espessura proporcional a ξ = , a qual separa as duas fases.

A espessura da interface de equilíbrio é definida como sendo a distância de . Dessa forma temos que

2.2. Equação de Cahn-Hilliard

No contexto de problemas transientes e assumindo que os fluxos difusivos interfaciais são aproximadamente proporcionais ao potencial químico [9], a dinâmica do campo de fase é descrita pela equação de Cahn-Hilliard

onde Φ = Φ (t,x,y,z), com -1 ≤ Φ (t,x,y,z) ≤ 1, é o campo de fase, µ(Φ) é o potencial químico e M (Φ) > 0 é a mobilidade ou coeficiente Onsager. Para este usamos, assim como Badalassi et al [3], a forma M (Φ) = Mc (1 - ϒ Φ2), onde Mc e ϒ são constantes não negativas. Como simulamos a separação das fases com mobilidade constante, adotamos Mc = 1 e ϒ = 0.

3. Metodologia Numérica

3.1. Estratégia semi-implícita

Seguimos a estratégia de Badalassi et al [3] para tratar a não linearidade proveniente do potencial químico na equação (2.3).

Lembrando que , podemos reescrever (2.3) como

sendo

Obtemos a equação (3.1) somando e subtraindo o termo

à equação (2.3). No termo (3.3), f' (Φ) é linearizado por Φ e a mobilidade M é majorada por M.

Na equação (3.1), o termo (3.3) será implicitado e o termo (3.2) explicitado, caracterizando assim o método semi-implícito [3, 14, 20].

3.2. Discretização temporal

A equação (3.1) contém o termo biharmônico ∇4Φ . Para evitar a discretização direta desse termo, consideramos as igualdades , as quais permitem reescrever a equação (3.1) na forma do sistema de equações

com .

Na solução das equações (3.4) e (3.5) adotamos condições de contorno de Neumann homogêneas ou de fluxo nulo, isto é, n . ∇Φ1 = 0 e n . ∇Φ2 = 0, onde n é o vetor unitário normal à fronteira do domínio de solução, ou condições de contorno periódicas.

Na discretização temporal das equações (3.4) e (3.5) utilizamos o Método de Gear com coeficientes variáveis no tempo [12, 14, 20], um método de segunda ordem. Assim,

onde

Como o Método de Gear é um método de dois passos no tempo, precisamos de Φ1 e Φ2 no instante de tempo t1. Determinamos Φ11 e Φ21 por intermédio do Método de Euler Implícito, o qual pode ser obtido das equações (3.6) e (3.7) considerando-se α2 = 1, α1 = -1, α0= 0, β1 = 1 e β0 = 0.

Usamos uma adaptação do ciclo W [13, 20] para solucionar o sistema formado pelas equações (3.6) e (3.7). Nessa técnica multigrid, empregamos Gauss-Seidel red-black [23] para relaxar a correção uma vez na descida e uma vez na subida. Na malha mais grossa, ao invés de solucionarmos exatamente a equação residual, relaxamos a correção algumas vezes. No ciclo W adaptado, utilizamos uma média simples na restrição e interpolações trilineares no prolongamento.

3.3. Discretização espacial

O sistema (3.6)-(3.7) é solucionado em um domínio computacional Ω = [A1, B1] x [A2, B2] x [A3, B3] refinado localmente, ou seja, apenas a região de transição na separação dos componentes da mistura é recoberta por malhas espaciais com resolução mais elevada. A Figura (1) mostra um exemplo com dois níveis de refinamento.


A malha com refinamento local (malha composta) ilustrada na Figura 1-(b) é denotada como 64 x 64 x 64 L2 ou 643 L2, ou seja, uma malha com dois níveis de refinamento sendo que as malhas mais grossas estão contidas em uma malha 64 x 64 x 64 e as malhas mais finas em uma malha 128 x 128 x 128. Nessa malha refinada localmente, calculamos o campo de fase Φ no centro da célula computacional, um paralelepípedo de arestas ∆x, ∆y e ∆z, e discretizamos os operadores gradiente, divergente e Laplaciano empregando diferenças finitas progressivas, regressivas e centradas de segunda ordem [20].

3.4. Estabilidade

A discretização semi-implícita para a equação de Cahn-Hilliard tem um comportamento incondicionalmente estável nos testes numéricos. Nesse caso, para condições iniciais suaves e para garantir a precisão no espaço, empregamos o passo temporal dado por

onde ∆t1 = tfinal - t, ∆t2 = min {∆x, ∆y e ∆z}, tfinal é o tempo necessário para que o estado estacionário seja alcançado e ∆x, ∆y e ∆z são os passos espaciais em cada uma das direções no nível que contém as malhas mais finas. No passo temporal (3.8) temos que ∆t é da ordem de ∆x (ou ∆y ou ∆z), o que torna os erros no tempo e no espaço proporcionais. Alteramos o passo temporal (3.8) quando queremos capturar a separação das fases a partir de uma condição inicial estabelecida por uma perturbação randômica (spinodal decomposition). Neste caso, usamos a estratégia adaptativa para o passo temporal descrita no Algoritmo 1, a qual possibilita capturar com precisão também a parte inicial e rápida da separação.


4. Refinamento Adaptativo

Queremos refinar localmente a região de transição porque esta requer alta resolução. Em nosso modelo, as fases estáveis do fluido bifásico são representadas pelos valores -1 e 1 e a transição pelos valores entre eles. Usamos esta propriedade para marcar as células computacionais que serão refinadas usando o algoritmo de Berger e Rigoutsos (AMR) [6]. Dessa maneira, uma célula computacional no nível imediatamente abaixo do nível mais fino é marcada se o campo de fase Φ calculado no seu centro satisfaz a propriedade

sendo Φtrans = 0, Φmax = 1 o valor máximo do campo de fase e p um porcentual de Φmax, como por exemplo, p = 0,99. Com esta estratégia refinamos p% da região de transição.

A adaptatividade do refinamento é definida pela verificação da propriedade (4.1) a cada nt passos no tempo, com n = 2 ou n = 4. Se pelo menos uma célula computacional do novo conjunto de células marcadas no nível imediatamente abaixo do nível mais fino não pertencer ao conjunto de células marcadas anteriormente, o processo de refinamento é então novamente ativado.

5. Verificação da Ordem do Esquema Numérico

Apresentamos um teste em malhas uniformes e em malhas refinadas localmente para comprovar numericamente a convergência do esquema de segunda ordem proposto para a equação de Cahn-Hilliard (2.3). Na malha com refinamento local, utilizamos dois níveis de refinamento - Figura 2 - que permanecem inalterados ao longo do tempo (malha composta estática). Queremos dessa forma verificar se os erros de truncamento introduzidos pela interpolação e pela discretização nas interfaces entre as malhas grossas e finas são controlados corretamente a fim de evitar a degradação da precisão global.


Para a análise de convergência, assumimos que a solução numérica obtida possui uma expansão assintótica em potências do espaçamento da malha [20]. Isto nos permite usar a estimativa

para determinar a ordem de convergência do campo de fase Φ. Em (5.1), n é a ordem de convergência e R é a razão entre os erros absolutos cometidos nas malhas com espaçamento , respectivamente.

Verificamos a ordem de discretização do esquema numérico proposto para a equação de Cahn-Hilliard (2.3) testando-o na solução da equação

onde Fch (t,x) é um termo forçante.

Construimos um problema teste suave escolhendo a função

com 0 ≤ t ≤ 10 e (x, y, z) = x Є [0,1] x [0,1] x [0,1], como solução exata de (5.2) (solução manufaturada), de tal maneira que . A partir dessa escolha, o termo forçante é dado por

e a condição inicial por Φ (0,x) = Φe (0,x).

Adotamos mobilidade constante, isto é, M (Φe(t,x)) = 1, condições de contorno periódicas e o passo temporal (3.8). Para as constantes presentes nas equações discretas (3.6) e (3.7), atribuimos os valores , = 2 [22] e e2 = 0,01.

A Tabela 1 mostra os resultados obtidos para o teste em uma malha uniforme e na malha com refinamento local descrita anteriormente, respectivamente. Observando-a, percebemos a convergência de segunda ordem do esquema numérico.

6. Simulação Computacional

Na simulação da separação das fases usamos como condição inicial para o campo de fase Φ uma distribuição randômica fornecida por

onde Φm = 0 (spaghetti case) é a concentração constante da mistura uniforme, e é a constante relacionada à espessura da interface de transição e r (x, y, z) Є [-1, 1]com média igual a zero.

As equações (3.4) e (3.5) são solucionadas no domínio Ω = [0, 1] x [0, 1] x [0, 1], com condições de contorno de Neumann homogêneas, ou seja,

n . ∇ Φ1 = 0 e n . ∇Φ2 = 0.

O passo temporal empregado é aquele mencionado no Algoritmo 1, sendo c = 10, e as equações discretas (3.6) e (3.7) são integradas durante um tempo t = 23. Este tempo é suficiente para que ocorra a separação completa dos dois componentes da mistura.

Utilizamos uma malha refinada localmente 323 L3, sendo que a malha do nível mais fino coincide inicialmente com o domínio. Desta forma, começamos integrando as equações em uma malha 128 x 128 x 128. Como a fase inicial da separação (0 ≤ t ≤ 10e) ocorre rapidamente, ativamos nessa etapa o refinamento localizado a cada 4 passos no tempo. Após essa fase, a cada p = 0,99 passos. Usamos a estratégia (4.1) com p = 0,99 para marcar as células computacionais que serão submetidas ao AMR e definir assim as regiões refinadas.

Quanto às constantes presentes nas equações discretas (3.6) e (3.7), adotamos para elas os valores

A Figura 3 mostra o bom comportamento do código computacional na execução do teste. Observando-a, constatamos que a energia livre descresce monotonicamente em direção a um valor mínimo e que Φm é preservado pelo menos até a segunda casa decimal [3]. Quanto ao número de células computacionais, percebemos a redução do custo computacional no decorrer do tempo.


A Figura 4 ilustra algumas etapas da simulação da separação dos componentes de uma mistura bifásica empregando a equação de Cahn-Hilliard; a Figura 5, a malha composta em um instante do estado estacionário.



7. Conclusão

Apresentamos uma estratégia computacional para a simulação tridimensional completamente adaptativa de um modelo de campo de fase conservativo descrito pela equação de Cahn-Hilliard. Com ela simulamos a separação dos componentes de uma mistura bifásica.

A metodologia numérica empregada combinou malhas refinadas localmente com uma discretização temporal semi-implícita de segunda ordem e técnicas multinível-multigrid lineares para a solução de sistemas. A discretização proposta é isenta de restrições severas de estabilidade e a ordem da discretização foi comprovada numericamente com o emprego de uma solução manufaturada.

8. Ferramentas Computacionais

Desenvolvemos o código computacional em linguagem Fortran 90 e executamos a simulação em uma Power Mac G5 (modelo M9592LL/A) com processador quad (duplo dual) de 2.5GHz, 16GB de memória RAM, 250GB de disco rígido, aritmética de 64 bits, compilador absoft para Fortran 90 e sistema operacional Linux (ydl). Visualizamos os resultados com o Tecplot 360.

  • [1] I. Altas, J. Dym, M.M. Gupta, R.P. Manohar, Multigrid solution of automatically generated high-order discretizations for the biharmonic equation, SIAM J. Sci. Comput., 19 (1998), 1575-1585.
  • [2] I. Altas, J. Erhel, M.M. Gupta, High accuracy solution of three-dimensional biharmonic equations, Num. Algorith., 29 (2002), 1-19.
  • [3] V.E. Badalassi, H.D. Ceniceros, S. Banerjee, Computation of multiphase systems with phase field models, J. Comput. Phys., 190 (2003), 371-397.
  • [4] J. Bell, M.J. Berger, J. Saltzman, M. Welcome, Three-dimensional adaptive mesh refinement for hyperbolic conservation laws, SIAM J. Sci. Comput., 15 (1994), 127-138.
  • [5] M.J. Berger, Data structures for adaptive grid generation, SIAM J. Sci. Stat. Comput., 7 (1986), 904-916.
  • [6] M.J. Berger, I. Rigoutsos, An algorithm for point clustering and grid generation, IEEE Trans. Syst. Man Cybernet., 21 (1991), 1278-1286.
  • [7] A.J. Bray, Theory of phase-ordering kinetics, Adv. Phys., 43 (1994), 357-459.
  • [8] W.L. Briggs, "A Multigrid Tutorial", Society for Industrial and Applied Mathematics, Philadelphia, 1987.
  • [9] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system - I. Interfacial free energy, J. Chem. Phys., 28 (1958), 258-267.
  • [10] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system - II. Thermodynamical basis, J. Chem. Phys., 30 (1959), 1121-1135.
  • [11] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system - III. Nucleation in a two-component incompressible fluid, J. Chem. Phys., 31 (1959), 688-699.
  • [12] H.D. Ceniceros, R.L. Nós, A.M. Roma, Three-dimensioal, fully adaptive simulations of phase-field fluid models, J. Comput. Phys., 229 (2010), 6135-6155.
  • [13] H.D. Ceniceros, R.L. Nós, A.M. Roma, Solução de equações diferenciais parciais elípticas por técnicas multinível-multigrid em malhas tridimensionais bloco-estruturadas com refinamento localizado, em "XXV Congresso Nacional de Matemática Aplicada e Computacional", São Paulo, SP, 2005.
  • [14] H.D. Ceniceros, A.M. Roma, A nonstiff, adaptive mesh refinement-based method for the Cahn-Hilliard equation, J. Comput. Phys., 225 (2007), 1849-1862.
  • [15] R. Chella, J. Viñals, Mixing of a two-phase fluid by a cavity flow, Phys. Rev. E, 53 (1996), 3832-3840.
  • [16] M.E. Gurtin, D. Polignone, J. Viñals, Two-phase binary fluids and immiscible fluids described by an order parameter, Math. Models Methods Appl. Sci., 6 (1996), 815-831.
  • [17] J. Kim, K. Kang, J.S. Lowengrub, Conservative multigrid methods for Cahn-Hilliard fluids, J. Comput. Phys., 193 (2004), 511-543.
  • [18] J. Kim, A numerical method for the Cahn-Hilliard equation with a variable mobility, Commun. Nonlinear Sci. Numer. Simul., 12 (2007), 1560-1571.
  • [19] S.F. Mccormick, Multigrid methods, em "Frontiers in Applied Mathematics", Vol. 3, SIAM Books, Philadelphia, 1987.
  • [20] R.L. Nós, "Simulações de escoamentos tridimensionais bifásicos empregando métodos adaptativos e modelos de campo de fase", Tese de Doutorado, IME, USP, São Paulo, SP, 2007.
  • [21] U. Trottenberg, C. Oosterlee, A. Schüller, "Multigrid", Academic Press, London, 2001.
  • [22] C. Xu, T. Tang, Stability analysis of a large time-stepping methods for epitaxial growth models, SIAM J. Numer. Anal., 44 (2006), 1759-1779.
  • [23] I. Yavneh, Multigrid smoothing factors for red-black Gauss-Seidel relaxation applied to a class of elliptic operators, SIAM J. Numer. Anal., 32 (1995), 1126-1138.
  • [24] O. Wodo, B. Ganapathysubramanian, Computationally efficient solution to the Cahn-Hilliard equation: adaptive implicit time schemes, mesh sensitivity analysis and the 3D isoperimetric problem, J. Comput. Phys., 230 (2011), 6037-6060.
  • 1
    Parte desta pesquisa foi financiada pela National Science Foundation (projeto DMS 0609996), pela FAPESP (projetos 04/13781-1 e 06/57099-5) e pelo CNPq.
  • Datas de Publicação

    • Publicação nesta coleção
      05 Out 2012
    • Data do Fascículo
      Abr 2012
    Sociedade Brasileira de Matemática Aplicada e Computacional Rua Maestro João Seppe, nº. 900, 16º. andar - Sala 163 , 13561-120 São Carlos - SP, Tel. / Fax: (55 16) 3412-9752 - São Carlos - SP - Brazil
    E-mail: sbmac@sbmac.org.br