\documentclass{TEMA}

\usepackage[brazil]{babel}      % para texto em Português
%\usepackage[english]{babel}    % para texto em Inglês

\usepackage[latin1]{inputenc}   % para acentuação em Português
%\input{P-margin.inf}

\usepackage[dvips]{graphics}
\usepackage{subfigure}
\usepackage{graphicx}
\usepackage{epsfig}

\usepackage{amsmath}
\usepackage{wrapfig}
\usepackage{algorithm}
\usepackage{algorithmic}
\floatname{algorithm}{Algoritmo}
\renewcommand{\listalgorithmname}{Lista de Algoritmos}

\newcommand{\B}{{\tt\symbol{92}}}
\newcommand{\til}{{\tt\symbol{126}}}
\newcommand{\chap}{{\tt\symbol{94}}}
\newcommand{\agud}{{\tt\symbol{13}}}
\newcommand{\crav}{{\tt\symbol{18}}}

\begin{document}

%********************************************************
\title
    {Simulação tridimensional adaptativa da separação das fases de uma mistura bifásica usando a equação de Cahn-Hilliard\thanks{Parte desta pesquisa 
     foi financiada pela {\em National Science Foundation} 
     (projeto DMS 0609996) e pela FAPESP (projetos 04/13781-1 e 06/57099-5).}}

\author
     {R.L. NÓS%
     \thanks{rudimarnos@utfpr.edu.br}\,,
     Departamento Acadêmico de Matemática,
     UTFPR, 80230-901 Curitiba, PR, Brasil \\
     H.D. CENICEROS%
     \thanks{hdc@math.ucsb.edu}\,,
     Department of Mathematics,
     UCSB, 93106 Santa Barbara, CA, USA \\
     A.M. ROMA%
     \thanks{roma@ime.usp.br}\,,
     Departamento de Matemática Aplicada,
     IME,
     USP, 05311-970 São Paulo, SP, Brasil.}

\criartitulo

\runningheads {Ceniceros, Nós e Roma}{Separação dos componentes de uma mistura bifásica}

\vspace*{-0.5cm}

\begin{abstract}
{\bf 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.

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


%********************************************************
\newsec{Introdução}

\hspace*{0.4cm} A equação de Cahn-Hilliard é um modelo de campo de fase \cite{ceniceros,bray,cahn,ceniceros2,cenroma3,chella,gurtin,kim2} 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 {\em parâmetro de ordem} ou {\em campo de fase} $\displaystyle{\phi}$ 
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 \cite{ceniceros} e 
de Ceniceros, Nós e Roma \cite{ceniceros2,cenroma3,nos1}, 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 \cite{ceniceros2,cenroma3,nos1} e o sistema linear prove\-niente da discretização é 
solucionado com técnicas {\it multinível-multigrid} \cite{altas,altas2,briggs,nos,cenroma3,kim,mccormick,nos1,trottenberg,yavneh}. 
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 
{\it multigrid} linear. Quanto à discretização do domínio, utilizamos o refinamento local 
adaptativo (AMR - {\em Adaptive Mesh Refinement}) introduzido por Berger \cite{berger1,berger2,berger5} 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.

%********************************************************

\newsec{Modelo matemático}

\subsection{Modelo de campo de fase}

\hspace*{0.4cm} 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 $\displaystyle{\phi}$ \cite{ceniceros,bray,cahn,cenroma3,chella,nos1}, ou seja,
\begin{equation}
\displaystyle{F\left[\phi\right] = \int_{\Omega} \left(f\left(\phi \left(\mathbf{x}\right)\right)
    + \frac{1}{2}\epsilon^2 \left|\nabla \phi\left(\mathbf{x}\right)\right|^2\right)d\mathbf{x}}, 
\label{free}
\end{equation}
onde $\displaystyle{\Omega}$ é a região do espaço ocupada pelo sistema, 
$\displaystyle{\frac{1}{2}\epsilon^2 \left|\nabla \phi\right|^2}$ é o termo que representa a energia 
interfacial, $\displaystyle{\epsilon >0}$ é uma constante associada à espessura da interface de 
transição e $\displaystyle{f\left(\phi \right)}$ é a densidade de energia típica. Esta tem dois 
mínimos que correspondem as duas fases estáveis do fluido.

A energia livre (\ref{free}) está sujeita à restrição 
$\displaystyle{\int_{\Omega} \phi \left(t,\mathbf{x}\right)d\mathbf{x} = \phi_m |\Omega |}$,
sendo $\displaystyle{\phi_m}$ o valor médio do campo de fase $\displaystyle{\phi}$ e 
$\displaystyle{|\Omega |}$ o volume do domínio $\displaystyle{\Omega}$.

A primeira derivada funcional da energia livre (\ref{free})
\[
\displaystyle{\mu\left(\phi\right)=\frac{\delta F\left[\phi\right]}{\delta
           \phi\left(\mathbf{x}\right)}=
           f^{\prime}\left(\phi\left(\mathbf{x}\right) \right)-\epsilon^2 \nabla^2\phi\left(
           \mathbf{x}\right)}  
\]
é o {\em 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 $\displaystyle{F\left[\phi\right]}$ com 
respeito às variações de $\displaystyle{\phi}$. Essa minimização fornece o perfil de equilíbrio na 
interface.

Para a densidade de energia típica, usamos a função {\em double well potential} (``poço duplo'') 
\cite{ceniceros,cenroma3,chella,nos1}
\[
\displaystyle{f\left(\phi\right) = \frac{\alpha}{4}\left(\frac{\beta}{\alpha} - \phi^2\right)^2} 
\]
com $\displaystyle{\alpha=\beta=1}$ ({\em Helmholtz free energy}). Assim, $\displaystyle{f^{\prime}\left(\phi\right) =  
\phi^{3} - \phi}$. Com essa escolha, os valores $\displaystyle{1}$ e $\displaystyle{-1}$ 
representam as duas fases da mistura e a região de transição é identificada pela variação dos 
valores de $\displaystyle{\phi}$ de $\displaystyle{-1}$ a $\displaystyle{1}$.

O perfil de equilíbrio em $R^n$ é dado pelas soluções da equação
\[
\displaystyle{\mu\left(\phi\right) = \phi^3 - \phi - \epsilon^2 \nabla^2\phi = 0}. 
\]

As soluções $\displaystyle{\phi_{\pm} = \pm 1}$ são estáveis, uniformes e representam as fases 
coexistentes. A terceira solução 
\begin{equation}
\displaystyle{\phi_0\left(z\right) =
  \tanh\left(\frac{\sqrt{2}}{2\xi}z\right)}\label{onesol} 
\end{equation}
é não uniforme, unidimensional na direção $\displaystyle{z}$ e satisfaz as condições de fronteira 
$\displaystyle{\lim_{z\rightarrow \pm \infty}\phi_0\left(z\right) = \pm 1}$.

A solução (\ref{onesol}) \cite{ceniceros,chella,nos1} descreve o perfil de equilíbrio de uma interface plana normal 
à direção $\displaystyle{z}$, de espessura proporcional a $\displaystyle{\xi = \epsilon}$, a qual separa 
as duas fases.

A espessura da interface de equilíbrio $\displaystyle{\delta}$ é definida como sendo a 
distância de $\displaystyle{\phi_0\left(z_1\right) = }$ $\displaystyle{0,9\phi_{-}}$ a 
$\displaystyle{\phi_0\left(z_2\right) = 0,9\phi_{+}}$. Dessa forma temos que  
\[
\displaystyle{\delta = |z_2 - z_1| = 2\sqrt{2}\xi \tanh^{-1}(0,9) \approx{4,164\xi}}.
\]

\subsection{Equação de Cahn-Hilliard}

\hspace*{0.4cm} No contexto de problemas transientes e assumindo que os fluxos difusivos interfaciais são aproximadamente proporcionais ao potencial químico \cite{cahn}, a dinâmica do campo de fase é descrita pela equação de Cahn-Hilliard
\begin{equation}
\displaystyle{\frac{\partial \phi}{\partial t} = \nabla \cdot
  \left[M\left(\phi\right)\nabla \mu\left(\phi\right)\right]}, 
\label{chnc} 
\end{equation}
onde $\displaystyle{\phi = \phi \left(t,x,y,z\right)}$, com $\displaystyle{-1 \leq 
\phi(t,x,y,z)\leq 1}$, é o campo de fase, $\displaystyle{\mu\left(\phi\right)}$ é o potencial 
químico e $\displaystyle{M\left(\phi\right)>0}$ é a {\em mobilidade} ou 
{\em coeficiente Onsager}. Para este usamos, assim como Badalassi et al \cite{ceniceros}, 
a forma $\displaystyle{M\left(\phi\right) = M_c\left(1 - \gamma \phi^2\right)}$, onde 
$\displaystyle{M_c}$ e $\displaystyle{\gamma}$ são constantes não negativas. Como simulamos a separação das fases com mobilidade constante, adotamos $M_c=1$ e $\gamma=0$.

%********************************************************

\newsec{Metodologia numérica}

\subsection{Estratégia semi-implícita}

\hspace{0.4cm} Seguimos a estratégia de Badalassi et al \cite{ceniceros} para tratar a não 
linearidade proveniente do potencial químico na equação (\ref{chnc}). 

Lembrando que $\displaystyle{\nabla f^{\prime}\left(\phi\right) = f^{\prime
    \prime}\left(\phi\right) \nabla \phi}$, podemos reescrever (\ref{chnc}) como
\begin{equation}
\displaystyle{\frac{\partial \phi}{\partial t} = \tilde{M}\nabla^2
  \left(\tau \phi - \epsilon^2 \nabla ^2 \phi\right) +
  g\left(\phi\right)},\label{chncmod}   
\end{equation}
sendo
\begin{equation} 
\displaystyle{g\left(\phi\right) = \nabla \cdot \left[M \left(\phi\right) \nabla \mu 
\left(\phi\right)\right] - \tilde{M}\nabla^2 \left(\tau \phi - \epsilon^2 \nabla
  ^2 \phi\right)}, \label{chncmodg} \\
\end{equation}
$\displaystyle{\tau \geq c_1 \max_{-1 \leq \phi \leq 1}\left\{ f^{\prime \prime}
  \left(\phi\right)\right\},\:\:c_1 \in ]0,\infty[\:\:\:e\:\:\: \tilde{M} = c_2 
\max_{-1 \leq \phi \leq 1} \left\{M \left(\phi\right) \right\},\:\:c_2 \in ]0,\infty[}$. 
\vspace*{0.3cm}

Obtemos a equação (\ref{chncmod}) somando e subtraindo o termo 
\begin{equation}
\displaystyle{\tilde{M} \nabla \cdot 
  \left[\nabla \left(\tau \phi - \epsilon^2 \nabla ^2 \phi\right)\right]}
  \label{chncm3}      
\end{equation}
à equação (\ref{chnc}). No termo (\ref{chncm3}), $\displaystyle{f^{\prime}\left(\phi\right)}$ é 
linearizado por $\displaystyle{\tau \phi}$ e a mobilidade $\displaystyle{M}$ é majorada por 
$\displaystyle{\tilde{M}}$.

Na equação (\ref{chncmod}), o termo (\ref{chncm3}) será implicitado e o termo 
(\ref{chncmodg}) explicitado, caracterizando assim o método semi-implícito 
\cite{ceniceros,cenroma3,nos1}.  

\subsection{Discretização temporal}

\hspace{0.4cm} A equação (\ref{chncmod}) contém o termo biharmônico 
$\displaystyle{\nabla^4 \phi}$ . Para evitar a discretização direta desse termo, 
consideramos as igualdades $\displaystyle{\phi_1 = \phi} \label{equalphi1}$ e 
$\displaystyle{\phi_2= \tau \phi_1 - \epsilon^2 \nabla^2 \phi_1}$, as quais permitem reescrever 
a equação (\ref{chncmod}) na forma do sistema de equações 
\begin{align}
\displaystyle{\frac{\partial \phi_1}{\partial t}} & \displaystyle{ = \tilde{M}\nabla^2
  \phi_2 + g_1\left(\phi_1,\phi_2\right)}, \label{chphi1} \\    
\displaystyle{\phi_2} & \displaystyle{ = \tau \phi_1 - \epsilon^2 \nabla^2 \phi_1},
  \label{chphi2}  
\end{align}
com $\displaystyle{g_1\left(\phi_1,\phi_2\right) = \nabla \cdot \left[M \left(\phi_1\right) 
\nabla \mu \left(\phi_1\right)\right] - \tilde{M}\nabla^2 \phi_2}$. 
\vspace*{0.3cm}

Na solução das equações (\ref{chphi1}) e (\ref{chphi2}) adotamos condições de contorno de Neumann homogêneas ou de fluxo nulo, isto é, $\displaystyle{\mathbf{n} \cdot \nabla \phi_1 = 0}$ e 
$\displaystyle{\mathbf{n} \cdot\nabla \phi_2 = 0}$, onde $\displaystyle{\mathbf{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 (\ref{chphi1}) e (\ref{chphi2}) utilizamos o Método de 
Gear com coeficientes variáveis no tempo \cite{ceniceros2,cenroma3,nos1}, um método de segunda ordem. Assim,
\begin{align} 
\displaystyle{\frac{\alpha_2\phi_1^{n+1}+\alpha_1\phi_1^n+\alpha_0\phi_1^{n-1}}
    {\Delta t}} & \displaystyle{= \tilde{M} \nabla ^2 \phi_2^{n+1} + 
    \beta_1g_1\left(\phi_1^{n},\phi_2^{n}\right) + \beta_0 g_1\left(\phi_1^{n-1},
    \phi_2^{n-1}\right)}, \label{chphi1g} \\
\displaystyle{\phi_2^{n+1}} & \displaystyle{= \tau \phi_1^{n+1} - 
    \epsilon ^2\nabla ^2 \phi_1^{n+1}}, \label{chphi2g}
\end{align} 
onde $\displaystyle{g_1\left(\phi_1^m,\phi_2^m\right) = \nabla \cdot \left[M\left(
\phi_1^m\right)\nabla \mu\left(\phi_1^m\right)\right] - \tilde{M} \nabla^2 \left(
\phi_2^m\right)}$. 
\vspace*{0.3cm}

Como o Método de Gear é um método de dois passos no tempo, precisamos de 
$\displaystyle{\phi_1}$ e $\displaystyle{\phi_2}$ no instante de tempo $\displaystyle{t^1}$. 
Determinamos $\displaystyle{\phi_1^1}$ e $\displaystyle{\phi_2^1}$ por intermédio do Método 
de Euler Implícito, o qual pode ser obtido das equações (\ref{chphi1g}) e (\ref{chphi2g}) 
considerando-se
$\displaystyle{\alpha_2 = 1}$, $\displaystyle{\alpha_1 = -1}$, 
$\displaystyle{\alpha_0 = 0}$, $\displaystyle{\beta_1 = 1}$ e
$\displaystyle{\beta_0 = 0}$. 

Usamos uma adaptação do ciclo W \cite{nos,nos1} para solucionar o sistema formado pelas 
equações (\ref{chphi1g}) e (\ref{chphi2g}). Nessa técnica {\it multigrid}, empregamos 
Gauss-Seidel {\it red-black} \cite{yavneh} 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.

\subsection{Discretização espacial}

\hspace{0.4cm} O sistema (\ref{chphi1g})-(\ref{chphi2g}) é solucionado em um domínio 
computacional $\displaystyle{\Omega = \left[A_1,B_1\right] \times}$ $\displaystyle{
\left[A_2,B_2\right] \times \left[A_3,B_3\right]}$ 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 (\ref{disdomain}) mostra um exemplo com dois níveis de 
refinamento. 

\begin{figure}[htbp!]
\centering{
\epsfig{file=spag16.ps, height=4.5cm,width=4.5cm}
\epsfig{file=spagmeshes16.ps, height=4.5cm,width=4.5cm}
}
\centerline{
{\bf (a)} \qquad \qquad \qquad \qquad \qquad \quad {\bf (b)} 
}
\vspace{-0.8cm}
\caption{ {\small  Separação dos componentes de uma mistura bifásica com dois níveis 
de refinamento local: (a) componentes da mistura e região de transição no domínio 
$\displaystyle{\Omega = [0,1] \times [0,1] \times [0,1]}$; (b) refinamento local. O nível 
mais grosso corresponde a uma malha $64 \times 64 \times 64$.}}
\label{disdomain}
\end{figure}

A malha com refinamento local (malha composta) ilustrada na Figura \ref{disdomain}-(b) é denotada 
como $64 \times 64 \times 64$ L$2$ ou $64^3$ L$2$, ou seja, uma malha com dois níveis de refinamento 
sendo que as malhas mais grossas estão contidas em uma malha $64 \times 64 \times 64$ e as malhas 
mais finas em uma malha $128 \times 128 \times 128$. Nessa malha refinada localmente, calculamos 
o campo de fase $\displaystyle{\phi}$ no centro da célula computacional, um paralelepípedo de 
arestas $\displaystyle{\Delta x}$, $\displaystyle{\Delta y}$ e $\displaystyle{\Delta z}$, e 
discretizamos os operadores gradiente, divergente e Laplaciano empregando diferenças finitas 
\cite{nos1}.  

\subsection{Estabilidade}

\hspace{0.4cm} 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
\begin{equation}
\displaystyle{\Delta t = \min\left\{\Delta t_1,\Delta t_2\right\}}, 
\label{estabchnc}
\end{equation}
onde $\displaystyle{\Delta t_1 = t_{final}-t}$, $\displaystyle{\Delta t_2 = \min\left\{\Delta x,
\Delta y,\Delta z\right\}}$, $t_{final}$ é o tempo necessário para que o estado estacionário seja alcançado e $\displaystyle{\Delta x}$, $\displaystyle{\Delta y}$ e 
$\displaystyle{\Delta 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 (\ref{estabchnc}) temos que $\displaystyle{\Delta t}$ é da ordem 
de $\displaystyle{\Delta x}$ (ou $\displaystyle{\Delta y}$ ou $\displaystyle{\Delta z}$), o que torna 
os erros no tempo e no espaço proporcionais. Alteramos o passo temporal (\ref{estabchnc}) quando 
queremos capturar a separação das fases a partir de uma condição inicial estabelecida por uma 
perturbação randômica ({\it spinodal decomposition}). Neste caso, usamos a estratégia adaptativa para 
o passo temporal descrita no Algoritmo \ref{spdecomp}, a qual possibilita capturar com precisão também 
a parte inicial e rápida da separação.

\begin{algorithm}[h!]
\caption{\small Seleção do passo temporal de integração para a equação de Cahn-Hilliard.}
\label{spdecomp}
\vspace{0.5cm}
se $0 \le t \le c\epsilon,\:\:c \in ]0,\infty[$ então                                          \\
   $\displaystyle{\:\:\:\:\:\Delta t_1 = t_{final}-t}$                                    \\
   $\displaystyle{\:\:\:\:\:\Delta t_2 = \epsilon^2}$                                     \\
   $\displaystyle{\:\:\:\:\:\Delta t = \min\left\{\Delta t_1,\Delta t_2\right\}}$         \\
senão                                                                                     \\ 
   $\displaystyle{\:\:\:\:\:\Delta t_1 = t_{final}-t}$                                    \\
   $\displaystyle{\:\:\:\:\:\Delta t_2 = \min\left\{\Delta x,\Delta y,\Delta z\right\}}$  \\
   $\displaystyle{\:\:\:\:\:\Delta t = \min\left\{\Delta t_1,\Delta t_2\right\}}$         \\
fim se                                                                                    \\
\end{algorithm}

%********************************************************
\vspace{-0.5cm}
\newsec{Refinamento adaptativo}

\hspace{0.4cm} 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) \cite{berger5}. Dessa maneira, uma célula computacional no nível imediatamente abaixo do nível mais 
fino é marcada se o campo de fase $\displaystyle{\phi}$ calculado no seu centro satisfaz a propriedade
\begin{equation}
\displaystyle{\phi \in \left[\phi_{trans} - p\:\phi_{max}, \phi_{trans} + p\:\phi_{max}\right]},
\label{phiref}
\end{equation}
sendo $\displaystyle{\phi_{trans}=0}$, $\displaystyle{\phi_{max} = 1}$  o valor máximo do campo 
de fase e $\displaystyle{p}$ um porcentual de $\displaystyle{\phi_{max}}$, como por exemplo, 
$\displaystyle{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 (\ref{phiref}) a cada 
$\displaystyle{n\:\Delta t}$ passos no tempo, com $\displaystyle{n=2}$ ou $\displaystyle{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.

%********************************************************

\newsec{Verificação da ordem do esquema numérico}

\hspace{0.4cm} 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 
(\ref{chnc}). Na malha com refinamento local, utilizamos dois níveis de refinamento - Figura 
\ref{compositgrid} - 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. 

\begin{figure}[htbp!]
\centering{
\epsfig{file=meshesx.ps, height=4.5cm,width=4.5cm}
\epsfig{file=meshesy.ps, height=4.5cm,width=4.5cm}
}
\centerline{
{\bf (a)} \qquad \qquad \qquad \qquad \qquad \quad {\bf (b)} 
}
\vspace{-0.8cm}
\caption{ {\small Malha composta com dois níveis de refinamento (o nível mais grosso corresponde a uma malha $32 \times 32 \times 32$)}: (a) corte em $x=0,1$; (b) corte em $y=0,2$.}
\label{compositgrid}
\end{figure}  

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 \cite{nos1}. Isto nos permite usar a estimativa
\begin{equation}
\displaystyle{\log_2 R \rightarrow n} 
\label{converg}
\end{equation}
para determinar a ordem de convergência do campo de fase $\displaystyle{\phi}$. Em (\ref{converg}), 
$\displaystyle{n}$ é a ordem de convergência e $\displaystyle{R}$ é a razão entre os erros absolutos 
cometidos nas malhas com espaçamento $\displaystyle{h}$ e $\displaystyle{\frac{h}{2}}$, respectivamente.

Verificamos a ordem de discretização do esquema numérico proposto para a equação de Cahn-Hilliard 
(\ref{chnc}) testando-o na solução da equação 
\begin{equation}
\displaystyle{\frac{\partial \phi}{\partial t}\left(t,\mathbf{x}\right) = 
\nabla \cdot \left[M\left(\phi\left(t,\mathbf{x}\right)\right)\nabla \mu\left(\phi
\left(t,\mathbf{x}\right)\right)\right] + F_{ch}\left(t,\mathbf{x}\right)}, 
\label{chforcing}      
\end{equation}
onde $\displaystyle{F_{ch}\left(t,\mathbf{x}\right)}$ é um termo forçante.

Construimos um problema teste suave escolhendo a função
\[
\displaystyle{\phi_e \left(t,\mathbf{x}\right) = \sin^{3}\left(2\pi x + 2\pi y + 2\pi z 
 + t\right)},
\]
com $\displaystyle{0 \leq t \leq 10}$ e $\displaystyle{\left(x,y,z\right) = 
\mathbf{x} \in [0,1] \times [0,1] \times [0,1]}$, como solução exata de (\ref{chforcing}) 
(solução manufaturada), de tal maneira que $\displaystyle{-1 \leq \phi_e \leq 1\:\:\:\forall 
\left(t,\mathbf{x}\right)}$. A partir dessa escolha, o termo forçante é dado por
\begin{displaymath}
\displaystyle{F_{ch}\left(t,\mathbf{x}\right) = \frac{\partial\phi_e}{\partial t}
\left(t,\mathbf{x}\right) - \nabla \cdot \left[M\left(\phi_e\left(t,\mathbf{x}
\right)\right)\nabla\mu\left(\phi_e\left(t,\mathbf{x}\right)\right)\right]}  
\end{displaymath}
e a condição inicial por $\displaystyle{\phi\left(0,\mathbf{x}\right) = \phi_e\left(0,\mathbf{x}\right)}$.

Adotamos mobilidade constante, isto é, $\displaystyle{M\left(\phi_e\left(t,\mathbf{x}\right)\right) = 1}$, 
condições de contorno periódicas e o passo temporal (\ref{estabchnc}). Para as constantes presentes nas 
equações discretas (\ref{chphi1g}) e (\ref{chphi2g}), atribuimos os valores $\displaystyle{\tilde{M} = 1}$, 
$\displaystyle{\tau = 2}$ \cite{xu} e $\displaystyle{\epsilon^2 = 0,01}$.

A Tabela \ref{chtable1} 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.

\begin{table}[htbp!]
\begin{center}
\begin{tabular}{|c|c|c||c|c|c|} 
\hline
Malhas & $\|\phi-\phi_e\|_2$ & R & Malhas & $\|\phi-\phi_e\|_2$ & R \\ 
\hline
\hline
$32^3$ L$1$ & $7,475653 \times 10^{-2}$ &  & $32^3$ L$2$ & $8,085192 \times 10^{-2}$ & \\
\hline
$64^3$ L$1$ & $1,843576 \times 10^{-2}$ & $4,05$ & $64^3$ L$2$ & $2,376079 \times 10^{-2}$ & $3,40$ \\
\hline
$128^3$ L$1$ & $4,584894 \times 10^{-3}$ & $4,02$ & $128^3$ L$2$ & $5,963587 \times 10^{-3}$ & $3,98$ \\
\hline
\end{tabular}
\centerline{
{\bf (a)} \qquad \qquad \qquad \qquad \qquad \qquad \quad {\bf (b)} 
}
\vspace{-0.8cm}
\caption{ {\small Razão de convergência para a equação de Cahn-Hilliard com condições de contorno periódicas e mobilidade constante: (a) malha uniforme; (b) malha composta.}}
\label{chtable1}  
\end{center}
\end{table}

%********************************************************
\vspace{-1.0cm}
\newsec{Simulação computacional}

\hspace{0.4cm} Na simulação da separação das fases usamos como condição inicial para o 
campo de fase $\displaystyle{\phi}$ uma distribuição randômica fornecida por
\[
\displaystyle{\phi \left(0,x,y,z\right) = \phi_m + \epsilon\:r(x,y,z)},
\label{randomdist}
\]
onde $\displaystyle{\phi_m = 0}$ é a concentração constante da mistura uniforme, 
$\displaystyle{\epsilon}$ é a constante relacionada à espessura da interface de 
transição e $\displaystyle{r(x,y,z) \in [-1,1]}$ com média igual a zero.

As equações (\ref{chphi1}) e (\ref{chphi2}) são solucionadas no domínio 
$\displaystyle{\Omega = [0,1] \times [0,1] \times [0,1]}$, com condições de 
contorno de Neumann homogêneas, ou seja, \\ \centerline{$\displaystyle{\mathbf{n} \cdot \nabla \phi_1 = 0}$ e 
$\displaystyle{\mathbf{n} \cdot \nabla \phi_2 = 0}$.}

O passo temporal empregado é aquele mencionado no Algoritmo \ref{spdecomp}, 
sendo $\displaystyle{c=10}$, e as equações discretas (\ref{chphi1g}) e (\ref{chphi2g}) são 
integradas durante um tempo $\displaystyle{t=23}$. Este tempo é suficiente para que ocorra 
a separação completa dos dois componentes da mistura.

Utilizamos uma malha refinada localmente $\displaystyle{32^3}$ L$\displaystyle{3}$,  
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 $\displaystyle{128 \times 128 \times 128}$. Como a fase 
inicial da separação $\displaystyle{\left(0 \leq t \leq 10\epsilon \right)}$ ocorre rapidamente, 
ativamos nessa etapa o refinamento localizado a cada $\displaystyle{2}$ passos no tempo. Após essa 
fase, a cada $\displaystyle{4}$ passos. Usamos a estratégia (\ref{phiref}) com 
$\displaystyle{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 (\ref{chphi1g}) e (\ref{chphi2g}), 
adotamos para elas os valores \\ \centerline{$\displaystyle{\tau = 2}$, $\displaystyle{\tilde{M} = 1}$, 
$\displaystyle{M_c=1\:e\:\gamma = 0 \rightarrow M = 1}$ e 
$\displaystyle{\epsilon^2 = 0,00095 \approx \left(\frac{4,0}{128}\right)^2}$.}

A Figura \ref{sdgraphs} 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 $\displaystyle{\phi_m}$ é preservado pelo menos até a segunda casa decimal 
\cite{ceniceros}. Quanto ao número de células computacionais, percebemos a redução do custo 
computacional no decorrer do tempo. 

\begin{figure}[htbp!]
\centering{
\epsfig{file=free_energy.ps, height=4.0cm,width=4.0cm}
\epsfig{file=phiaverage.ps, height=4.0cm,width=4.0cm}
\epsfig{file=cellsnumber.ps, height=4.0cm,width=4.0cm}
}
\centerline{
{\bf (a)} \qquad \qquad \qquad \qquad \qquad {\bf (b)} 
\qquad \qquad \qquad \qquad \qquad {\bf (c)}
}
\vspace{-0.8cm}
\caption{ {\small Simulação da separação dos componentes de uma mistura 
bifásica empregando a equação de Cahn-Hilliard: (a) energia livre do sistema; 
(b) valor médio $\displaystyle{\phi_m}$ do campo de fase $\displaystyle{\phi}$
; (c) número de células computacionais.}}
\label{sdgraphs}
\end{figure}

A Figura \ref{sp1} 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 \ref{sp2}, a malha composta em um instante do estado estacionário.  
\vspace{1.0cm}

\begin{figure}[htbp!]
\centering{
\epsfig{file=spag1.ps, height=4.0cm,width=4.0cm}
\epsfig{file=spag3.ps, height=4.0cm,width=4.0cm}
\epsfig{file=spag4.ps, height=4.0cm,width=4.0cm}
}
\centerline{
{\bf (a)} \qquad \qquad \qquad \qquad \qquad {\bf (b)} 
\qquad \qquad \qquad \qquad \qquad {\bf (c)}
}
\centering{
\epsfig{file=spag5.ps, height=4.0cm,width=4.0cm}
\epsfig{file=spag8.ps, height=4.0cm,width=4.0cm}
\epsfig{file=spag12.ps, height=4.0cm,width=4.0cm}
}
\centerline{
{\bf (d)} \qquad \qquad \qquad \qquad \qquad {\bf (e)} 
\qquad \qquad \qquad \qquad \qquad {\bf (f)}
}
\centering{
\epsfig{file=spag15.ps, height=4.0cm,width=4.0cm}
\epsfig{file=spag17.ps, height=4.0cm,width=4.0cm}
\epsfig{file=spag20.ps, height=4.0cm,width=4.0cm}
}
\centerline{
{\bf (g)} \qquad \qquad \qquad \qquad \qquad {\bf (h)} 
\qquad \qquad \qquad \qquad \qquad {\bf (i)}
}
\vspace{-0.8cm}
\caption{ {\small Simulação da separação dos componentes de uma mistura 
bifásica empregando a equação de Cahn-Hilliard: estado da separação nos tempos 
(a) $t=0$, (b) $t \approx 3,8$x$10^{-2}$, (c) $t \approx 5,3$x$10^{-2}$, 
(d) $t \approx 7,6$x$10^{-2}$, (e) $t \approx 0,16$, (f) $t \approx 0,46$,
(g) $t \approx 1,89$, (h) $t \approx 5,58$ e (i) $t \approx 23,00$.}}
\label{sp1}
\end{figure}

\begin{figure}[htbp!]
\centering{
\epsfig{file=spaghetti_contourmeshes_3232.ps, height=4.2cm,width=4.2cm}
}
\vspace{-0.5cm}
\caption{ {\small Malha composta no instante $t \approx 23,00$.}}
\label{sp2}
\end{figure}

%********************************************************

\vspace*{-1.1cm}
\newsec{Conclusão}

\hspace{0.4cm} 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 {\it 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.

%********************************************************

\newsec{Ferramentas computacionais}

\hspace{0.4cm} Desenvolvemos o código computacional em linguagem Fortran 90 e executamos a 
simulação em uma {\it 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 
{\it absoft} para Fortran 90 e sistema operacional {\it Linux} (ydl).  Visualizamos os 
resultados com o Tecplot 360.

%********************************************************

\vspace*{.3cm}

\begin{abstract}
{\bf 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 require 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.
\end{abstract}

\begin{thebibliography}{8}

\bibitem{altas}
I. Altas, J. Dym, M.M. Gupta, R.P. Manohar, Multigrid solution of automatically generated high-order discretizations for the biharmonic equation, {\em SIAM J. Sci. Comput.}, {\bf 19(5)} (1998), 1575-1585.

\bibitem{altas2}
I. Altas, J. Erhel, M.M. Gupta, High accuracy solution of three-dimensional biharmonic equations, {\em Num. Algorith.}, {\bf 29} (2002), 1-19.

\bibitem{ceniceros} 
V.E. Badalassi, H.D. Ceniceros, S. Banerjee, Computation of multiphase systems with phase field models, {\em J. Comput. Phys.}, {\bf 190} (2003), 371-397.

\bibitem{berger1}
J. Bell, M.J. Berger, J. Saltzman, M. Welcome, Three-dimensional adaptive mesh refinement for hyperbolic conservation laws, {\em SIAM J. Sci. Comput.}, {\bf 15(1)} (1994), 127-138.

\bibitem{berger2}
M.J. Berger, Data structures for adaptive grid generation, {\em SIAM J. Sci. Stat. Comput.}, {\bf 7(3)} (1986), 904-916.

\bibitem{berger5} 
M.J. Berger, I. Rigoutsos, An algorithm for point clustering and grid generation, {\em IEEE Trans. Syst. Man Cybernet.}, {\bf 21(5)} (1991), 1278-1286.

\bibitem{bray}
A.J. Bray, Theory of phase-ordering kinetics, {\em Adv. Phys.}, {\bf 43(3)} (1994), 357-459.

\bibitem{briggs}
W.L. Briggs, ``A multigrid tutorial'', Society for Industrial and Applied Mathematics, Philadelphia, 1987.

\bibitem{cahn}
J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system - I. Interfacial free energy, {\em J. Chem. Phys.}, {\bf 28(2)} (1958), 258-267.

\bibitem{ceniceros2}
H.D. Ceniceros, R.L. Nós, A.M. Roma, Three-dimensioal, fully adaptive simulations of phase-field fluid models, {\em J. Comput. Phys.}, {\bf 229(17)} (2010), 6135-6155.

\bibitem{nos} 
H.D. Ceniceros, R.L. Nós, A.M. Roma, Solução de equações diferenciais par\-ciais 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.

\bibitem{cenroma3} 
H.D. Ceniceros, A.M. Roma, A nonstiff, adaptive mesh refinement-based method for the Cahn-Hilliard equation, {\em J. Comput. Phys.}, {\bf 225(2)} (2007), 1849-1862.

\bibitem{chella}
R. Chella, J. Viñals, Mixing of a two-phase fluid by a cavity flow, {\em Phys. Rev. E}, {\bf 53(4)} (1996), 3832-3840.

\bibitem{gurtin}
M.E. Gurtin, D. Polignone, J. Viñals, Two-phase binary fluids and immiscible fluids described by an order parameter, {\em Math. Models Methods Appl. Sci.}, {\bf 6(6)} (1996), 815-831.

\bibitem{kim}
J. Kim, K. Kang, J.S. Lowengrub, Conservative multigrid methods for Cahn-Hilliard fluids, {\em J. Comput. Phys.}, {\bf 193} (2004), 511-543.

\bibitem{kim2}
J. Kim, A numerical method for the Cahn-Hilliard equation with a variable mobility, {\em Commun. Nonlinear Sci. Numer. Simul.}, {\bf 12} (2007), 1560-1571.

\bibitem{mccormick}
S.F. Mccormick, Multigrid methods, em ``Frontiers in Applied Mathematics'', Vol. 3, SIAM Books, Philadelphia, 1987.

\bibitem{nos1} 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.

\bibitem{trottenberg}
U. Trottenberg, C. Oosterlee, A. Schüller, ``Multigrid'', Academic Press, London, 2001.

\bibitem{xu}
C. Xu, T. Tang, Stability analysis of a large time-stepping methods for epitaxial growth models, {\em SIAM J. Numer. Anal.}, {\bf 44(4)} (2006), 1759-1779.

\bibitem{yavneh}
I. Yavneh, Multigrid smoothing factors for red-black Gauss-Seidel relaxation applied to a class of elliptic operators, {\em SIAM J. Numer. Anal.}, {\bf 32(4)} (1995), 1126-1138.

\end{thebibliography}


\end{document}
\newpage
$ \  \  $  \thispagestyle{myheadings}  \markboth{      }{   }