%% Antes de processar este arquivo LaTeX (LaTeX2e) deve
%% verificar que o arquivo TEMA.cls estah no mesmo
%% diretorio. O arquivo TEMA.cls pode ser obtido do
%% endereco www.sbmac.org.br/tema.

\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}

\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
    {Modelo e Simula\c{c}\~ao Num\'erica de Intera\c{c}\~oes Envolvendo Bolhas e Gotas}

\author
    {R. MANICA\thanks{manicar@ihpc.a-star.edu.sg; Author para correspond\^encia}\,, E. KLASEBOER\thanks{evert@ihpc.a-star.edu.sg}\,,   Institute of High Performance Computing, 1 Fusionopolis Way, \#16-16 Connexis, Singapore 138632 \\
      \\ D.Y.C. CHAN\thanks{D.Chan@unimelb.edu.au}, Department of Mathematics and Statistics, The University of Melbourne, VIC 3010, Australia. Faculty of Life and Social Sciences, Swinburne University of Technology, VIC 3122, Australia}

\criartitulo

\runningheads {Manica, Klaseboer e Chan}{Intera\c{c}\~oes entre bolhas e gotas}

\begin{abstract}
{\bf Resumo}. Formulamos um modelo te\'orico para o estudo de escoamento em filmes finos envolvendo superf\'icies deform\'aveis (por exemplo bolhas e gotas) que est\~ao interagindo a baixas velocidades. Assumimos que a condi\c{c}\~ao de contorno na interface entre os fluidos \'e de n\~ao-deslizamento tangencial. As equa\c{c}\~oes evolutivas resultantes constituem um sistema alg\'ebrico-diferencial na qual a posi\c{c}\~ao da fronteira avan\c{c}a e deforma ao mesmo tempo e depende da solu\c{c}\~ao global do sistema. O foco principal do trabalho \'e na deriva\c{c}\~ao do modelo e nos detalhes da implementa\c{c}\~ao num\'erica. As equa\c{c}\~oes s\~ao resolvidas usando uma rotina do Matlab e os resultados num\'ericos s\~ao comparados a dados experimentais da literatura que foram produzidos por pesquisadores em differentes laborat\'orios e usando diversas t\'ecnicas, comprovando que o modelo \'e adequado para uma variedade de sistemas.

{\bf Palavras-chave}. Sistema alg\'ebrico-diferencial, equa\c{c}\~oes de Stokes-Reynolds e Young-Laplace, escoamento em filmes finos, coalesc\^encia, gotas e bolhas.
\end{abstract}

\newsec{Introdu\c{c}\~ao}

Intera\c{c}\~oes envolvendo bolhas e gotas est\~ao presentes em tarefas do nosso dia a dia como lavar os pratos ou tomar banho, mas tamb\'em podem ser encontradas em diversos setores industrias relevantes para a economia brasileira como por exemplo extra\c{c}\~ao de \'oleo e g\'as, estabilidade de emuls\~oes, ind\'ustria farmac\^eutica e de cosm\'eticos. Entender a intera\c{c}\~ao entre duas bolhas ou gotas em sistemas multif\'asicos desempenha um papel fundamental na determina\c{c}\~ao das caracter\'isticas de fluxo de todo o sistema. No entanto, a grande diferen\c{c}a de escalas de comprimento caracter\'isticas de tais sistemas implicam em desafios significativos na interpreta\c{c}\~ao dos resultados experimentais e na formula\c{c}\~ao de modelos te\'oricos. Por exemplo, para bolhas e gotas com tamanho em torno de 100 $\mu$m, as colis\~oes que podem levar \`a coalesc\^encia s\~ao afetadas em parte por for\c{c}as de superf\'icie (por exemplo Van der Waals-Lifshitz e el\'etrica) que operam em escalas de nan\^ometros. Por isso \'e importante descrever com precis\~ao as propriedades de fluxo do sistema evolutivo com o mesmo n\'ivel de resolu\c{c}\~ao espacial. Al\'em disso, deforma\c{c}\~oes da superf\'icie que s\~ao da ordem de nan\^ometros em bolhas e gotas s\~ao tamb\'em associadas \`as condi\c{c}\~oes de fluxo e \`a magnitude das for\c{c}as de superf\'icie. Fatores como a natureza da condi\c{c}\~ao de contorno hidrodin\^amica na superf\'icie das bolhas tamb\'em contribuem para a sua intera\c{c}\~ao. Finalmente, \'e importante ter em conta a natureza din\^amica das intera\c{c}\~oes, tais como a depend\^encia das for\c{c}as de colis\~ao entre bolhas e gotas e determinar se as colis\~oes s\~ao est\'aveis ou levam \`a coalesc\^encia.

O uso de din\^amica de fluidos computacional na qual o sistema global \'e estudado considerando o movimento das bolhas em um l\'iquido apresenta uma s\'erie de dificuldades de identifica\c{c}\~ao da interface e tamb\'em da necessidade de malhas extremamente finas para capturar detalhes quando as interfaces est\~ao muito pr\'oximas. O uso da teoria da lubrifica\c{c}\~ao, na qual somente o filme fino \'e considerado, permite estudar a fase final do processo de intera\c{c}\~ao e determinar a estabilidade do sistema. A dificuldade na formula\c{c}\~ao do modelo est\'a no fato de termos um dom\'inio limitado na qual a fronteira est\'a situada em uma posi\c{c}\~ao da bolha fora da zona de intera\c{c}\~ao, mas menor do que seu raio. Essa fronteira avan\c{c}a e deforma simultaneamente e a deforma\c{c}\~ao depende da solu\c{c}\~ao de todo o sistema. Detalhes dessa formula\c{c}\~ao te\'orica foram revisados recentemente por Chan et al. \cite{CKM11ACIS} bem como compara\c{c}\~oes com dados experimentais \cite{CKM11SM}. Solu\c{c}\~oes anal\'iticas aproximadas para esse sistema s\'o s\~ao conhecidas para pequenas deforma\c{c}\~oes \cite{CKM09}.

\newsec{Equa\c{c}\~oes governantes}
Estamos interessados na intera\c{c}\~ao envolvendo superf\'icies deform\'aveis como mostradas na Figura 1. No primeiro exemplo (Figura 1a) vemos uma bolha ancorada a uma seringa, sendo esta aproximada de uma superf\'icie plana atrav\'es do movimento da seringa ou da superf\'cie. Outra possibilidade \'e injetar ar para que a bolha cres\c{c}a e interaja com a superf\'icie. Um experimento como esse foi feito por Fisher et al \cite{FMHRW91,FHMRW92, HFRF93} e a compara\c{c}\~ao com a teoria desenvolvida nesse trabalho vai ser mostrada na se\c{c}\~ao de resultados. J\'a na Figura 2b, temos o caso de duas bolhas ou gotas sendo aproximadas uma da outra devido a uma for\c{c}a exterior, por exemplo a gravidade ou fluxo de cisalhamento. Elas tamb\'em podem ser aproximadas de maneira mec\^anica na qual as gotas est\~ao ancoradas em algum objeto s\'olido, como por exemplo uma superf\'icie ou uma uma viga como no caso do microsc\'opio de for\c{c}a at\^omica \cite{DMC06, VMT10}. 

\begin{figure}[h] %
\centering
\subfigure{\scalebox{1}{\includegraphics[width=7.8cm]{Fig1a}} }
\hspace{1cm}  %
\subfigure{\scalebox{1}{\includegraphics[width=3.cm]{Fig1b}} }
\caption{a) Esquema da intera\c{c}\~ao entre uma bolha ancorada numa seringa interagindo com uma superf\'icie plana, onde $z(r,t)$ \'e a superf\'icie da bolha, enquanto  $h(r,0)$ e $h(r,T)$ representam a separa\c{c}\~ao no tempo inicial 0 e final $T$.  b) Intera\c{c}\~ao entre duas bolhas ou gotas com raios $R_1$ e $R_2$ sendo aproximadas por for\c{c}a exterior.} 
\label{fig1}
\end{figure}

Assumindo que a fase cont\'inua \'e constitu\'ida de um fluido incompress\'ivel Newtoniano, o fluxo \'e  governado pelas equa\c{c}\~oes de Navier-Stokes e pela continuidade escritas na forma \cite{Batchelor67, MB04} 
\begin{eqnarray} \label{NS}
\rho \left( \frac{\partial  \textbf{u}}{\partial t}+ \textbf{u}  \cdot \nabla \textbf{u} \right) & = & - \nabla p + \mu \triangle \textbf{u}\\
\nabla \cdot  \textbf{u} & = &0 \label{NS1}
\end{eqnarray}
onde $\rho$ \'e a massa espec\'ifica, $\mu$ a viscosidade da fase cont\'inua (assumida constante mesmo para separa\c{c}\~oes muito pequenas na ordem de nan\^ometros), $p$ \'e a press\~ao no filme fino e $\textbf{u} = (u, v, w)^T$ \'e a velocidade em fun\c{c}\~ao do tempo $t$ e de tr\^es coordenadas espaciais ($x, y, z$).

Quando a bolha est\'a pr\'oxima da superf\'icie plana ela se deforma de acordo com a Figura 1a.  Nesse caso, usando coordenadas polares, as varia\c{c}\~oes em $z$ s\~ao muito menores que em $r$. Tomando coordenadas cartesianas $x$ e $y$, isso implica
\[\frac{\partial h}{\partial x}\ll 1; \quad \frac{\partial h}{\partial y} \ll 1\]
que \'e uma condi\c{c}\~ao necess\'aria na teoria da lubrifica\c{c}\~ao \cite{Reynolds86} e no estudo de filmes finos em geral \cite{Leal92}. A validade dessa condi\c{c}\~ao foi discutida nos trabalhos de Chesters \cite{AC94, Chesters91, SGC95}. Assumimos tamb\'em que as velocidades de intera\c{c}\~ao s\~ao baixas e consideramos que ocorre fluxo de Stokes na qual o n\'umero de Reynolds \'e pequeno, ou seja, $Re=2R_o \rho  V/\mu  \ll 1$, onde $R_o$ \'e o raio da bolha n\~ao deformada e $V$ a velocidade de aproxima\c{c}\~ao da bolha. Adotamos as escalas sugeridas por Klaseboer et al. \cite{KCGM00} que ser\~ao definidas no cap\'itulo sobre implementa\c{c}\~ao num\'erica. Definimos  o n\'umero capilar  $Ca = \mu V / \sigma$ como  a rela\c{c}\~ao entre for\c{c}as viscosas e a tens\~ao superficial $\sigma$. Inserindo essas escalas nas equa\c{c}\~oes de Navier-Stokes e continuidade (\ref{NS}) e (\ref{NS1}), assumindo que $Ca \ll 1$ e desprezando termos de segunda ordem, resulta nas equa\c{c}\~oes de lubrifica\c{c}\~ao escritas na forma
\begin{eqnarray}
\frac{\partial p}{\partial x}&=&\mu \frac{\partial^2 u}{\partial z^2} \label{lub1}\\
\frac{\partial p}{\partial y}&=&\mu \frac{\partial^2 v}{\partial z^2} \label{lub2}\\
\frac{\partial p}{\partial z}&=&0 \label{lub3}\\
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}&=&0 \label{con_mass}
\end{eqnarray}
Considerando a equa\c{c}\~ao (\ref{lub3}) vemos que $p$ n\~ao depende de $z$. Isso permite obter $u$ e $v$ ao integrar as equa\c{c}\~oes (\ref{lub1}) and (\ref{lub2}) com rela\c{c}\~ao a $z$ onde assumimos que $u=0$ e $v=0$ nas fronteiras $z=0$ e $z=h$. Note que essa condi\c{c}\~ao de contorno assume n\~ao-deslizamento tangencial na interface entre os fluidos, que \'e contr\'aria \`a condi\c{c}\~ao usual, mas que est\'a em acordo com os dados experimentais \cite{CKM11SM}. A  bolha est\'a se aproximando da superf\'icie a uma velocidade $w=-dh/ dt$, onde $h$ \'e a dist\^ancia da bolha \`a superf\'icie. Substituindo $u$, $v$ e $w$ na equa\c{c}\~ao (\ref{con_mass}) obtemos a equa\c{c}\~ao de Stokes-Reynolds

\begin{equation}\label{SRbolha}
\frac{\partial h}{\partial t}=\frac{1}{12 \mu} \frac{\partial}{\partial x}\left(h^{3} \frac{\partial p}{\partial x}\right)+\frac{1}{12 \mu} \frac{\partial}{\partial y}\left(h^{3} \frac{\partial p}{\partial y}\right)
\end{equation}

Note que a equa\c{c}\~ao (\ref{SRbolha}) \'e v\'alida tanto para uma bolha interagindo com uma superf\'icie plana como para duas bolhas interagindo entre si. Para a press\~ao usamos a equa\c{c}\~ao de Young-Laplace (tamb\'em conhecida como balan\c{c}o de tens\~ao normal) escrita na forma
\begin{equation}\label{YL999}
p+\Pi=\frac{2 \sigma}{R} -\frac{\sigma}{n} \left( \frac{\partial^2h}{\partial x^2}+\frac{\partial^2h}{\partial y^2}\right)
\end{equation}
onde $\Pi (h(r,t))$ \'e uma for\c{c}a de superf\'icie, por exemplo van der Waals ou el\'etrica que depende somente da dist\^ancia entre as interfaces e $2 \sigma/R$ \'e a press\~ao de Laplace. Quando a deforma\c{c}\~ao \'e pequena, $R \sim R_o$ e $2 \sigma/R_o$ \'e a diferen\c{c}a de press\~ao entre o interior e exterior da bolha. O valor constante $n$ \'e 1 para uma bolha ou 2 para duas bolhas. Nesse segundo caso, $R_o$ \'e dado pela m\'edia harm\^onica entre $R_1$ e $R_2$ definido como
\[ \frac{1}{R_o}=\frac{1}{2}\left(\frac{1}{R_1}+\frac{1}{R_2}\right)\]

As equa\c{c}\~oes (\ref{SRbolha}) e (\ref{YL999}) s\~ao v\'alidas na regi\~ao de intera\c{c}\~ao quando existe a presen\c{c}a de um filme fino. Precisamos de condi\c{c}\~oes de contorno numa posi\c{c}\~ao $x$ e $y$ fora do filme mas menor que o raio da bolha. Devido a complexidade do sistema bidimensional em encontrar condi\c{c}\~oes de contorno, pois estas variam dependendo do sistema, vamos mostrar resultados assumindo que o sistema \'e axisim\'etrico, visto que a maioria dos resultados experimentais apresentam simetria axial.

\newsec{Equa\c{c}\~oes para o sistema axisim\'etrico}
Assumindo que o sistema mant\'em simetria axial nos permite resolver o problema em uma \'unica dimens\~ao espacial. Usando o mesmo procedimento do caso bidimensional podemos facilmente chegar \`as equa\c{c}\~oes abaixo \cite{CKM11SM}

\begin{equation}\label{SR}
\frac{\partial h}{\partial t}=\frac{1}{12 \mu r} \frac{\partial}{\partial r}\left(r h^{3} \frac{\partial p}{\partial r}\right)
\end{equation}

\begin{equation}\label{YL}
p+\Pi=\frac{2 \sigma}{R} -\frac{\sigma}{n r} \frac{\partial}{\partial r}\left(r  \frac{\partial h}{\partial r}\right)
\end{equation}
onde $n=1$ para uma bolha e $n=2$ para duas bolhas. Note que se a bolha tem formato esf\'erico implica $p=0$ e no caso de o filme ser plano temos $p= 2\sigma/R$ se assumirmos que a press\~ao de superf\'icie $\Pi$ \'e nula. A for\c{c}a de intera\c{c}\~ao do sistema \'e dada por
\begin{equation}\label{FG}
F=2\pi\int_0^{\infty}r(p+\Pi) dr
\end{equation}
onde $p+\Pi \to 0 $ quando $r \to \infty$. Para completar o formalismo, precisamos de uma condi\c{c}\~ao inicial e quatro condi\c{c}\~oes de contorno adequadas para o sistema. A condi\c{c}\~ao inicial \'e dada por
\begin{equation}
\label{init}
h(r,0) = h_0+\frac{n r^2}{2R_o}
\end{equation}
onde $h_0$ \'e a separa\c{c}\~ao inicial e assumimos que a superf\'icie mant\'em a sua forma esf\'erica original. As condi\c{c}\~oes de contorno no centro de simetria s\~ao dadas por
\begin{eqnarray}
\label{BCaxi}
\frac{\partial h}{\partial r}  & = & 0 \quad \mbox{  em} \quad r = 0 \\
\frac{\partial p}{\partial r}  & = & 0 \quad  \mbox{  em} \quad r = 0
\end{eqnarray}

Escolhemos uma posi\c{c}\~ao $r_{max}$ (Figura 1a) onde aplicamos as outras duas condi\c{c}\~oes de contorno. Sabemos que fora da \'area de intera\c{c}\~ao a press\~ao decai numa fun\c{c}\~ao $r^{-4}$ \cite{YD90} que pode ser implementada como
\begin{equation}\label{dpdr}
 \frac{\partial p}{\partial r} + \frac{4}{r}p   =   0  \quad \mbox{ em} \quad r=r_{max}
\end{equation}

Outra possibilidade \'e simplesmente assumir que a press\~ao \'e nula em $r_{max}$, mas resultados num\'ericos mostraram que a condi\c{c}\~ao (\ref{dpdr}) fornece resultados mais acurados. Precisamos de uma condi\c{c}\~ao de contorno que forne\c{c}a o movimento da superf\'icie $\partial h/\partial t$ na posi\c{c}\~ao $r=r_{max}$. Uma possibilidade \'e assumir que n\~ao existe deforma\c{c}\~ao e que a superf\'icie deform\'avel se aproxima com a velocidade da base s\'olida \cite{AC94, KCGM00}
\begin{equation}\label{wrongbc}
\frac{\partial h}{\partial t} =  \mp V \quad \mbox{ em} \quad r=r_{max}
\end{equation}
onde usamos $-V$ se as superf\'icies est\~ao se aproximando e $V$ se est\~ao se afastando. Todavia, resultados num\'ericos \cite{CCLMD05} mostraram que a condi\c{c}\~ao (\ref{wrongbc}) fornece resultados que dependem da escolha de $r_{max}$, ou seja a solu\c{c}\~ao muda com o tamanho do dom\'inio. Isso se deve ao fato da superf\'icie deform\'avel ter comportamento logar\'itmico quando uma for\c{c}a exterior \'e aplicada. A forma assint\'otica da solu\c{c}\~ao da equa\c{c}\~ao (\ref{YL}), assumindo volume constante, \'e dada por \cite{CKM11SM} 
\begin{equation}
h(r,t) \sim h(0,t)+\frac{F}{2 \pi \sigma}\left[\ln\left(\frac{F}{8 \pi \sigma R_o}\right)+B(\theta)\right]
\end{equation}
onde $B(\theta)$ \'e uma constante que depende do \^angulo $\theta$ que a bolha est\'a ancorada \`a superf\'icie. Esse resultado indica que mesmo uma for\c{c}a pontual tem efeito de deformar toda a bolha, que adecua sua superf\'icie de forma a minimizar energia. 

Durante o doutorado do primeiro autor \cite{Manica07}, desenvolvemos uma condi\c{c}\~ao de contorno que \'e independente da escolha de $r_{max}$ dada por \cite{CCLMD05, CKM11ACIS}

\begin{equation}\label{hdotBC}
\frac{\partial h}{\partial t} +  \frac{n}{4 \pi \sigma} \frac{dF}{dt} \left[ 2 +  \ln\left(\frac{r_{max}^2}{4R_0^2}\right)
  +  \ln\left(\frac{1+ \cos \theta}{1-\cos \theta}\right)\right]  = \mp V  \quad \mbox{ em} \quad r=r_{max}
\end{equation}
Assim, completamos a formula\c{c}\~ao te\'orica do problema. Esse sistema pode ser adaptado para o caso de bolhas ou gotas interagindo sob o efeito da gravidade na qual a for\c{c}a \'e constante, o que gera um sistema de equa\c{c}\~oes adimensionais totalmente independente de qualquer vari\'avel e temos uma solu\c{c}\~ao universal para o problema. Mostraremos a seguir a implementa\c{c}\~ao num\'erica do sistema evolutivo de equa\c{c}\~oes.

\newsec{Implemeta\c{c}\~ao num\'erica}
As equa\c{c}\~oes evolutivas s\~ao implementadas na forma adimensional onde usamos as escalas sugeridas por Klaseboer et al. \cite{KCGM00} que foram baseadas nas escalas propostas por Chesters \cite{AC94}. A escala do tempo \'e dada pela velocidade $t_c = h_c/V$, a escala da press\~ao \'e dada pela press\~ao de Laplace, a escala radial depende da regi\~ao onde ocorre deforma\c{c}\~ao  $r_c^2 = R_0 h_c$ e a separa\c{c}\~ao $h_c$ \'e escolhida para nondimesionalizar as equa\c{c}\~oes.  Com essas escolhas as quantidades adimensionas s\~ao:
\begin{eqnarray*}
h_c & = &   R_0 Ca^{1/2} \\
r_c & = &   R_0  Ca^{1/4}\\
p_c & = &   \sigma/R_0 \\
t_c & = &   \mu  Ca^{-1/2} /p_c
\end{eqnarray*}
onde $Ca = \mu V / \sigma$ \'e o n\'umero capilar --  a rela\c{c}\~ao entre for\c{c}as viscosas e a tens\~ao superficial. A forma adimensional das equa\c{c}\~oes
 (\ref{SR}), (\ref{YL}) e (\ref{FG}) fica (usando os mesmos s\'imbolos para as vari\'aveis adimensionais)
\begin{eqnarray}\label{nthinfilm}
\frac{\partial h}{\partial t} & = & \frac{1}{12  r}\frac{\partial 
}{\partial r}\left(rh^3 \frac{\partial p}{\partial r} \right) \\
p + \Pi & = &2 - \frac{1}{nr} \frac{\partial }{\partial r}\left(r \frac{\partial 
h }{\partial r}\right)\label{npressure} \\
\frac{F}{2\pi\sigma} =G & = & \int_0^{\infty} r [p(r,t)+ \Pi(h(r,t))] dr
\end{eqnarray}
onde $n=1$ para uma bolha e $n=2$ para duas bolhas. Note que definimos a quantidade $G$ relacionada com a for\c{c}a de intera\c{c}\~ao. As condi\c{c}\~oes iniciais e de contorno ficam
\begin{eqnarray}
h(r,0) & = & h_0 + \frac{nr^2}{2} \\
\frac{\partial h}{\partial r}  & = & 0 \quad \mbox{  em \quad $r = 0$} \\
\frac{\partial p}{\partial r}  & = & 0 \quad \mbox{  em \quad $r = 0$}\\
\frac{\partial p}{\partial r} + \frac{4}{r}p  &  = &  0 \quad  \mbox{  em \quad $r=r_{max}$}\\
\frac{\partial h}{\partial t}  + \frac{n\alpha}{2} \frac{d G}{d t}   &=&  \mp 1  \quad \mbox{  em \quad $r=r_{max}$}\label{newBC}
\end{eqnarray}
onde
\begin{equation} \label{BCst}
\alpha=\left[ 2 + \ln\left(r_{max}^2 Ca^{1/2} \over 4\right)
  +  \ln\left(\frac{1+ \cos \theta}{1-\cos \theta}\right)\right]
\end{equation}
para linha de contato fixa em $r=r_1$. Para o caso na qual o \^angulo $\theta=\theta_p$ \'e fixo temos  \cite{CCLMD05}
\begin{equation}
 \alpha=\left[ 2 + \ln\left(r_{max}^2 Ca^{1/2} \over 4\right)
  +  \ln\left(\frac{1+ \cos \theta_p}{1-\cos \theta_p}\right) -\frac{2}{2+\cos \theta_p}\right]
\end{equation}

Note que o uso da condi\c{c}\~ao de contorno (\ref{newBC}) ao inv\'es de (\ref{wrongbc}) altera a natureza do problema, pois o sistema se torna alg\'ebrico-diferencial na qual a condi\c{c}\~ao de contorno depende da solu\c{c}\~ao do sistema e precisa ser resolvida simultaneamente.

Discretizamos o sistema na coordenada espacial $r$ usando o m\'etodo da linhas. Usamos diferen\c{c}as finitas centrais em $r$ nas equa\c{c}\~oes  (\ref{nthinfilm}) e (\ref{npressure}) e obtemos um sistema de equa\c{c}\~oes diferenciais ordin\'arias para $h_j(t) \equiv h(j \Delta r, t), j = 0 \cdots N$ onde $N  =  r_{max}/\Delta r$. Usamos uma malha uniforme $r = [0, r_{max}]$ com $\Delta r = 0.05 $ e $ r_{max} = 10$ produzindo um sistema de 200 equa\c{c}\~oes. As condi\c{c}\~oes de contorno em $r=0$ s\~ao usadas pra obter $\dot{h}_0$ e a  equa\c{c}\~ao (\ref{newBC}) fornece a condi\c{c}\~ao para $\dot{h}_N$. Isso requer $G$ como um vari\'avel extra a ser resolvida. Ao calcular $G$, dividimos a regi\~ao de integra\c{c}\~ao em duas partes:
\begin{equation}
G=\int_0^{r_{max}} r[ p + \Pi]  \ dr +  \int_{r_{max}}^{\infty} rp \ dr
\end{equation}
onde a primeira integral \'e calculada numericamente usando a regra de Simpson. Desta forma,  $G$ est\'a relacionada com todas as outras vari\'aveis $h_j$ como uma restri\c{c}\~ao alg\'ebrica. A segunda integral, onde $\Pi$ pode ser desprezado, pode ser calculada analiticamente baseado na forma $r^{-4}$ da press\~ao \cite{YD90}.

Resumindo, o sistema final de equa\c{c}\~oes tem a forma
\begin{equation}
\label{system}
\left(
\begin{array}{ccccc}
1 & 0 &  \cdots  & 0 & 0 \\
0 & 1 &  0  & \cdots & 0\\
   &     &   \vdots &    &  \\
0 & 0 &  \cdots  & 1 & \alpha \\
0 & 0 &  \cdots  & 0 & 0 \\
\end{array}
\right)
\left(
\begin{array}{c}
\dot{h}_0 \\
\dot{h}_1\\
      \vdots  \\
     \dot{h}_N \\
  \dot{G}
\end{array}
\right) =
\left(
\begin{array}{c}
f_0 \\
f_1\\
    \vdots  \\
   \mp 1 \\
  G - \sum_j \beta_j g(h_j)
\end{array}
\right)
\end{equation}
onde $\beta_j$ s\~ao os coeficientes da regra de Simpson e os valores $f_j$ representam as contribui\c{c}\~oes da equa\c{c}\~ao de evolu\c{c}\~ao do filme fino. Esse sistema tem matriz massa singular e \'e alg\'ebico-diferencial de \'indice 1 e pode ser resolvido usando algum software padr\~ao para EDOs, no nosso caso usando a rotina  \verb+ode15s+ no Matlab.

\section{Exemplos de aplica\c{c}\~oes}
Os resultados das simula\c{c}\~oes num\'ericas s\~ao comparados com dados experimentais da literatura. No primeiro exemplo considerou-se a intera\c{c}\~ao entre uma bolha e uma superf\'icie s\'olida plana de quartzo na qual dados foram publicados por Fisher et al \cite{FMHRW91,FHMRW92, HFRF93}. Detalhes adicionais sobre esse experimento bem como a teoria implementada podem ser encontrados na literatura recente \cite{CKM11SM, MC11}. No segundo exemplo consideramos a intera\c{c}\~ao entre duas gotas (fase discreta) imersas num segundo l\'iquido (fase cont\'inua) na qual os resultados experimentais foram obtidos por Klaseboer et al \cite{KCGM00}. Durante aquele experimento v\'arios materias foram usados incluindo \'agua, diferentes tipos de \'oleo, glicerina. Informa\c{c}\~oes adicionais relacionados com esse experimento, bem como situa\c{c}\~oes semelhantes de intera\c{c}\~ao entre duas superf\'icies deform\'aveis podem ser encontradas em \cite{CKM11ACIS, MKC08}.

\subsection{Intera\c{c}\~ao entre uma bolha e uma superf\'icie s\'olida plana}
O esquema do experimento foi mostrado na Figura 1a na qual uma bolha de ar \'e injetada de uma seringa que est\'a inicialmente pr\'oxima de uma superf\'icie plana. Os resultados experimentais, bem como a solu\c{c}\~ao num\'erica das equa\c{c}\~oes governantes, s\~ao reproduzidos na Figura 2 para uma caso de \'agua natural. Os dados experimentais foram extra\'idos da literatura \cite{FMHRW91, FHMRW92, HFRF93} onde os autores enfatizaram que nenhum modelo te\'orico existente  era capaz de explicar os resultados obtidos. 

\begin{figure}[h]
\centering %
\includegraphics[width=6.5cm]{Fig2} 
\caption{Intera\c{c}\~ao entre uma bolha e uma superf\'icie plana de quartzo em \'agua natural. Linhas cont\'inuas representam a teoria enquanto os s\'imbolos s\~ao experimentos que foram obtidos para $t$ = 5, 50 e 200 segundos.} %
\label{fig2}
\end{figure}

Uma bolha de raio $R_o \sim$ 1.1 mm est\'a inicialmente a uma dist\^ancia $h_0 \sim 40 \mu m$ da superf\'icie de quartzo. Quando ar \'e injetado a bolha cresce e eventualmente interage com a superf\'icie.  A press\~ao na parte central de intera\c{c}\~ao se torna superior \`a sua press\~ao interna, o que ocasiona uma invers\~ao de curvatura (curvas (a)-(c) na Figura 2). Esse processo inicial ocorre numa escala de tempo r\'apida e n\~ao foi capturado durante o experimento. Durante a fase de escoamento do filme fino (curvas (d)-(f) na Figura 2) o processo \'e lento durando em torno de 200 segundos, onde percebemos uma excelente concord\^ancia entre o modelo te\'orico (linhas s\'olidas) e o experimento (s\'imbolos) para diferentes tempos. Como o sistema tem simetria axial, o experimento foi conduzido para um \'unico raio.  Eventualmente um filme est\'avel aparece devido \`a repuls\~ao gerada pela for\c{c}a el\'etrica na superf\'icie (curva (f)). Esse fen\^omeno tamb\'em foi observado em experimentos envolvendo gotas de merc\'urio interagindo com uma superf\'icie plana \cite{MCCHC07}.

\subsection{Intera\c{c}\~ao entre duas gotas}
Um dos experimentos mais completos que estudam o fluxo e a evolu\c{c}\~ao de filmes finos entre duas gotas iguais (raio $\sim$1500 $\mu$m) de mesmo material em um outro l\'iquido foi realizado por Klaseboer et al \cite{KCGM00} durante seu trabalho de doutorado. Os dados experimentais foram analizados usando as mesmas equa\c{c}\~oes do presente trabalho, mas a condi\c{c}\~ao de contorno adotada naquele trabalho era dada pela equa\c{c}\~ao (\ref{wrongbc}). Embora os autores tenham sido capazes de obter resultados qualitativamente e at\'e quantitativamente precisos, estudos posteriores mostraram que a condi\c{c}\~ao (\ref{wrongbc}) n\~ao \'e adequada para este problema \cite{CCLMD05}. Isso ocorre pelo fato dessa condi\c{c}\~ao de contorno n\~ao levar em considera\c{c}\~ao a deforma\c{c}\~ao das gotas mantendo o volume constante.

\begin{figure}[h]
\centering %
\subfigure{\scalebox{1}{\includegraphics[width=5.cm]{Fig3a}} } %
\hspace{1cm}  %
\subfigure{\scalebox{1}{\includegraphics[width=6.2cm]{Fig3b}} } %
\caption{a) Exemplo de an\'eis de Newton obtidos experimentalmente. b) Evolu\c{c}\~ao da parte central $h_o$ e da separa\c{c}\~ao m\'inima $h_{r}$. As linhas representam a teoria enquanto os c\'irculos s\~ao dados experimentais. A flecha indica o momento que as gotas pararam de ser aproximadas em $t$=27 s.} %
\label{fig3}
\end{figure}

O experimento consiste de duas seringas colocadas uma contra a outra na qual duas gotas foram injetadas nas pontas. Ao aproximar as seringas com velocidade constante $V$, as gotas eventualmente interagem e deformam. O escoamento no filme fino pode ser bem lento se a fase cont\'inua tiver viscosidade alta e demorar muitos minutos e at\'e horas at\'e que ocorra coalesc\^encia. Por outro lado, experimentos em \'agua com gotas  \cite{DMC06} ou bolhas \cite{VMT10} muito menores ($R_o \sim 100 \mu m$) mostraram que a coalesc\^encia pode ocorrer em menos de 1 segundo. 

Os dados experimentais foram obtidos usando interfer\^ometro onde os an\'eis de Newton mostrados na Figura 3a podem ser convertidos em separa\c{c}\~oes atrav\'es da rela\c{c}\~ao de Braggs para um anel de ordem $m$: $h= m \lambda / 2n$ onde $\lambda$ = 632.8 nm \'e o comprimento de onda do laser e $n$ = 1.41 \'e o \'indice de refra\c{c}\~ao do \'oleo de silicone. 

Algumas caracter\'sticas importantes do escoamento no filme s\~ao a evolu\c{c}\~ao da parte central $h_0$ e da posi\c{c}\~ao onde a separa\c{c}\~ao \'e m\'inima $h_{r}$ que s\~ao mostrados na Figura 3b. Essa posi\c{c}\~ao \'e onde eventualmente ocorre coalesc\^encia desse sistema. Isso ocorre quando $h_{r}$ se torna suficientemente pequeno que a intera\c{c}\~ao de Van de Waals provoca a ruptura do filme. Esse \'e um exemplo de sistema inst\'avel enquanto no exemplo apresentado em na se\c{c}\~ao 5.1 o sistema era est\'avel devido a repuls\~ao el\'etrica entre a superf\'icie e a bolha. 

\section{Conclus\~ao e perspectivas}
Nesse trabalho apresentamos um sistema de equa\c{c}\~oes alg\'ebrico-diferencial que foi derivado baseado em problemas experimentais e nas caracter\'isticas f\'isicas da evolu\c{c}\~ao de filmes finos para superf\'icies deform\'aveis. Os resultados num\'ericos foram comparados com dados experimentais da literatura. A excelente concord\^ancia mostra que o modelo te\'orico desenvolvido captura as caracter\'isticas f\'isicas dos problemas estudados. O uso de teoria e recursos computacionais proporciona entender diversos fatores que n\~ao podem ser facilmente visualizados a partir dos dados experimentais.

Sistemas multif\'asicos reais encontrados na natureza e em processos industriais s\~ao geralmente mais complicados que os estudados nesse trabalho. A grande maioria das intera\c{c}\~oes n\~ao apresenta simetria axial, pode ocorrer coalesc\^encia prematura devido a presen\c{c}a de impurezas, efeitos relacionados a n\'umeros de Reynolds altos e turbul\^encia. Esses efeitos constituem muitas possibilidades para pesquisas futuras.  

\vspace*{.3cm}

\noindent {\bf {\large Agradecimentos}}\ 
Gostar\'iamos de agradecer ao Prof Dr AL De Bortoli (UFRGS) por fornecer in\'umeras sugest\~oes que melhoraram consideravelmente o manuscrito. DYC Chan \'e professor visitante da Universidade Nacional de Cingapura (NUS) e cientista em resid\^encia no Instituto do Computa\c{c}\~ao de Alta Performance (IHPC), Cingapura. R Manica gostaria de agradecer a UFRGS por fornecer a base matem\'atica durante gradua\c{c}\~ao e mestrado e a Universidade de Melbourne, Australia pela bolsa para fazer doutorado pleno sob supervis\~ao de DYC Chan. E finalmente agradece a NUS que forneceu ajuda financeira para fazer intercambio em Cingapura, onde conheceu E Klaseboer e mais tarde se tornou pesquisador no IHPC. 

\begin{abstract}
{\bf Abstract}. We developed a theoretical model to study the thin films involving deformable interfaces (for example drops and bubbles) which are interacting at low speeds. We assume the tangentially immobile boundary conditions holds at the fluid-fluid interface. The evolution equations for such system are of differential-algebraic nature in which the position of the boundary advances and deforms at the same time and the deformation depends on the solution. We focus on the model derivation and numerical implementation. The equations are solved using a Matlab routine and the numerical results are compared to experimental data from the literature, which were produced by researchers in different labs and using different techniques and shows the model is adequate to solve a variety of problems.

\end{abstract}

\begin{thebibliography}{8}
 
 \bibitem{AC94}  S. Abid, A.K. Chesters, The drainage and rupture of partially-mobile film between colliding drops at constant approach velocity, {\em Int. J. Multiphase Flow}, {\bf 7}, (1994), 613--629. 
 
  \bibitem{Batchelor67} G.K. Batchelor, ``An Introduction to Fluid Dynamics'', Cambridge University Press, 1967.
 
 \bibitem{CCLMD05}  S.L. Carnie, D.Y.C. Chan, C. Lewis, R. Manica, R.R. Dagastine, Measurement of dynamical forces between deformable drops using the atomic force microscope. {I. Theory}, {\em Langmuir}, {\bf 21}, (2005), 2912--2292. 

\bibitem{CKM09}  D.Y.C. Chan, E. Klaseboer, R. Manica, Dynamic deformations and forces in soft matter, {\em Soft Matter}, {\bf 5}, (2009), 2858--2861. 

\bibitem{CKM11SM}  D.Y.C. Chan, E. Klaseboer, R. Manica, Film drainage and coalescence between deformable drops and bubbles, {\em Soft Matter}, {\bf 7}, (2011), 2235--2264. 

\bibitem{CKM11ACIS}  D.Y.C. Chan, E. Klaseboer, R. Manica, Theory of non-equilibrium force measurements involving deformable drops and bubbles,  {\em Adv. Colloid Interface Sci.}, {\bf 165}, (2011), 70--90. 

\bibitem{Chesters91}  A.K. Chesters, Modelling coalescence processes in fluid-liquid dispersions: A review of current understanding,  {\em Trans. Inst. Chem. Engrs. Part A}, {\bf 69}, (1991), 259-270.

 \bibitem{DMC06} R.R. Dagastine, R. Manica, S.L. Carnie,D.Y.C. Chan, G.W. Stevens, F. Grieser, ``Dynamic forces between two deformable oil droplets in water'',   {\em Science}, {\bf 313}, (2006), 210--213.

\bibitem{FMHRW91} L.R. Fisher, E.E. Mitchell, D. Hewitt, J. Ralston, J. Wolfe, The drainage of a thin aqueous film between a solid surface and an approaching gas bubble, {\em Colloids Surf.}, {\bf 52}, (1991), 163--174. 

\bibitem{FHMRW92} L.R. Fisher, D. Hewitt, E.E. Mitchell, J. Ralston, J. Wolfe, The drainage of an aqueous film between a solid plane and an air bubble, {\em Adv. Colloid Interface Sci.}, {\bf 39}, (1992), 397--416.

\bibitem{HFRF93} D. Hewitt, D. Fornasiero,  J. Ralston, L.R. Fisher, Aqueous Film Drainage at the Quartz/Water/Air Interface, {\em J. Chem. Soc. Faraday Trans.}, {\bf 89}, (1993), 817--822.

\bibitem{KCGM00}  E. Klaseboer, J.Ph. Chevaillier,  C. Gourdon, O. Masbernat, Film drainage between colliding drops at constant approach velocity: experiments and modeling, {\em J. Colloid Interface Sci.}, {\bf 229}, (2000), 274--285. 

\bibitem{Leal92} L.G. Leal, ``Laminar flow and convective transport processes'', Butterworth
Heinemann, Boston, 1992.

\bibitem{MB04} R. Manica, A.L. De Bortoli, Simulations of sudden expansion flows for power-law fluids, {\em J. Non-Newtonian Fluid Mech.} {\bf 121}, (2004), 35--40.

\bibitem{Manica07} R. Manica, ``Modelling Hydrodynamic Interactions between Deformable Droplets'', PhD Thesis, The University of Melbourne, 2007.

\bibitem{MCCHC07} R. Manica, J.N. Connor, S.L. Carnie, R.G. Horn, D.Y.C. Chan, Dynamics of interactions involving deformable drops: hydrodynamic dimpling under attractive and repulsive electrical double layer interactions, {\em Langmuir} {\bf 23}, (2007), 626--637.

\bibitem{MC11} R. Manica, D.Y.C. Chan, Drainage of the air-water-quartz film: experiments and theory, {\em Phys. Chem. Chem. Phys.}, {\bf 13}, (2011), 1434--1439.

\bibitem{MKC08} R. Manica, E. Klaseboer, D.Y.C. Chan, Dynamic interactions between drops--a critical assessment, {\em Soft Matter}, {\bf 4}, (2008), 1613--1616. 

\bibitem{Reynolds86} O. Reynolds, On the theory of lubrication and its application, {\em Philos. Trans. R. Soc. London} {\bf A177}, (1886), 157.

\bibitem{SGC95} A. Saboni, C. Gourdon, A. K. Chesters, Drainage and rupture of partially mobile films during coalescence in liquid-liquid systems under a constant interaction force, {\em J. Colloid and Interface Sci.} {\bf 175}, (1995), 27--35.

\bibitem{VMT10} I.U. Vakarelskia, R. Manica, X. Tang, S.J. O'Shea, G.W. Stevens, F. Grieser,
R.R. Dagastine, D.Y.C. Chan, ``Dynamic interactions between microbubbles in water'', {\em PNAS}  {\bf 175}, (2010) 11177--11182.

\bibitem{YD90} S.G. Yiantsios, R.H. Davis, On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface, {\em J. Fluid Mech.}, {\bf 217}, (1990), 547--573.

\end{thebibliography}

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