\documentclass{article}
\usepackage[utf8]{inputenc}
\title{Lattice-Boltzmann modeling of multi-component systems: An \author{okuksen }
\date{June 2017}
\usepackage{natbib}
\usepackage{graphicx}
\begin{document}
\maketitle
\begin{abstract}
\end{abstract}
\maketitle
\section {LBM for multiphase fluids }
\subsection{Introduction}
Flows of multicomponent fluids mixtures occurs in a variety of natural as well as
technologically relevant processes, from ink-jet printing and processing of multi-component polymer blends to multiphase flows of oil-water mixtures in porous medium during enhanced oil recovery processes, hence modeling of such flows is of interest for numerous applications; we will provide a few examples by the end
of this sections. The major challenge in capturing the dynamics of multiphase fluids is
tracking or capturing the dynamics of the interface between the fluid components. The interface can be represented either as a sharp interface, or as so called diffuse interfaces, where the interface between the phases is relatively wide (diffuse) and is described
through the effective phase field. Hence all the available methods can be divided into sharp-interface methods and diffuse-interface methods. A number of approaches such as as boundary integral and boundary element methods can be used to track a sharp moving interface; in such methods, a grid undergoes deformation as the interface is deformed and re-meshing of the interface is often required.
Keeping track of moving interface can be computationally expensive and might not be feasible for the cases where morphological transitions are of interest (such as phase separation between the components). The multiphase LBM approach belongs to the class
of the diffusive interface methods \cite{Anderson},\cite{diffusive2004}. The important advantage of these class of methods is that the interface does not need to be tracked but the interfacial flows including dynamics of phase separation are captured through the
interactions between the different components.
There are three most well known LBM approaches that can be used to model
multiphase flows and capture the phase separation between the components. In
this chapter we will focus in details on the free-energy lattice Boltzmann approach
proposed in 1996 by Swift at el {\cite{swift}}and on the practical application of this
approach. The advantage of this approach is that the equilibrium distribution
functions are defined based on the system's free energy that also includes a gradient
term defining an interacial tension as we derive below. This allows one to define and to vary in interfacial tension
more easily that in any other multiphase, multicomponent LBM. We note that two other well known multiphase LBM approaches include color gradient method
proposed by Gunstensen et al. and a pseudo-potential model by Shan and Chen. The
color gradient approach was in fact the first multiphase LBM approach; in this method,
instead of a single distribution function as for example defined for the single
component fluid in the section above, two particle distribution functions were introduced for the first time,
red and blue distribution functions for two different immiscible red and blue
fluid phases. Local equilibrium distribution functions were defined by the local macroscopic parameters for each component and so-called "re-coloring" step updated the equilibrium distribution functions based on the color gradients; the
phase separation in this approach was driven by the repulsive interactions based on the color
gradient and color momentum. In the Shan-Chen pseudo-potential model, the non-
local interaction were introduced that controls the equation of state and result in the spontaneous phase separation between the components when the equation of state is appropriately chosen. This model produces more stable and accurate results
then the model by Gunstensen et al [Hou,Bao] and currently one of the most
commonly used multiphase LBM approaches.
An excellent review that compares all these approaches for multiphase flows is given
by Chen and Doolen. A number of more recent reviews on all three of these LBM
approaches exist that are given with respect to the specific questions such as flow in
the porous medium {YangBoek2013} or the multi-component multiphase flows
with high density ratios {BaoShaffler2013}. ADD MORe
\{subsection}{Governing continuum equations}
To describe the dynamics of the binary mixtures using the free-energy lattice
Boltzmann approach, we first define an order parameter in the binary system, which
is the difference between the densities of the phases A and B and can be written as
$\phi=\rho^A({\bf r})-\rho^B({\bf r})$, where
$\rho^A({\bf r})$ and $\rho^B({\bf r})$ are the densities of each of the
components. Correspondingly, the total density of system is defined as
$\phi=\rho^A({\bf r})+\rho^B({\bf r})$ The interface between the phases
varies smoothly across a number of lattice cites, so that in equilibrium across the AB
interface the order parameter changes gradually typically from the value set to the "-1" within the A phase to the "+1" within the B as we show below. It is worth
noting that the width of the interface between the fluids phases in LBM as well as in
others diffuse interface models is typically significantly larger than the interface width in
real experimental systems when one relates dimensionless values used in simulations and corresponding experimental values. Hence the connection between the length scales and time
scales should be made carefully based on the interfacial properties of the system (see below).
We first specify the full set of the continuum equations that is effectively numerically integrated
by using multiphase LBM approach with respective continuum boundary conditions.
Various CFD techniques could also be used to solve this system of continuum
equations, and the major advantage of using LBM are the relative ease of
implementation including implementations in the cases of complicated boundary conditions encompassing
chemically and topographically patterned walls as well as low computational cost of
simulations and a straightforward ability to parallelize the code.
The starting point of the LBM free energy approach is defining a
suitable free energy potential; the specific form of this
potential depends on the system of interest. One of the most common common choices is a Landau free energy potential:
\begin{equation}\label{F}
F = \int dV\Big[\Psi(\phi, \rho,T) + \kappa/2
\left | \bf \nabla \phi)^2\right |\Big] + \int dS\Psi_{s}(\phi_{s}).
\end{equation}
The first term, $\Psi(\phi, \rho,T)$, represents the bulk free energy and depends only on the densities and order parameter at a given moment in time. The integral here is taken over the entire volume of the fluid. The particular choice of this functional depends on the system of interests, often this functional is taken as \begin{equation}\label{psilandau}
\Psi(\phi, \rho,T)= \frac{c^2}{3}\rho ln{\rho}-\frac{a}{2}\phi^2+\frac{b}{4}\phi^4
\end{equation}
In general, the coefficient $a$ in the above equation depends on the temperature within the system \cite{Bray}. If the temperature is below the corresponding critical temperature, the fluids are immisible and undergo phase-separation when quenched from the homogeneous mixture; in this case, the coefficients in eq. (\ref{psilandau}) are chosen to be positive and the potential above represents so-called double-well potential\cite{Bray}. Typically the depth of the minima corresponding to the pure A and pure B phases are chosen to be equal, moreover, one typically chooses $b=a$, so that the equilibrium values of the phase-separating binary fluids are set at "-1" and "+1" for the A and B phases, respectively as we have already mentioned above.
We note that one could increase the temperature above the critical value which in turn would result in a change in sign of the parameter $a$ and instead of the double-well potential, the free energy will only have a single minimum corresponding to the uniformly intermixed phase. The first term in eq. (\ref{psilandau}) that depends on the total density of the fluid, is chosen in this particular functional form to improve numerical accuracy of LBM implementation, where $c=\Delta x/ \Delta t$ is the speed of lattice, with $\Delta x $ and $\Delta t$ being defined as the spacing between the nearest neighboring cites and the time step, respectively.
Notably, one could also add additional interactions to this system in a straightforward manner. For example, if the nanoparticles with preferential wetting interactions are dispersed within the system, corresponding potential energy favoring one of phases surrounding the nanoparticles \cite{AnnaNP} could be added as an additional term in eq. (\ref{F}). ADD Another common form of the free energy potential that is used in the free-energy multiphase LBM explicitly includes the repulsive energy term between the two components, $\lambda\rho^A\rho^B$,
where $ \lambda $ is the strength of the repulsion, reads:
\begin{equation}\label{psi0}
\psi(\phi, \rho,T)= \frac{{\lambda\rho}}{4}
(1-\frac{\phi^2}{\rho^2}) - T \rho +\frac{T}{2}
(\rho+\phi)
ln (\frac{\rho + \phi}{2}) + \frac{T}{2}
(\rho - \phi) ln (\frac{\rho - \phi}{2}).
\end{equation}
Here again the free energy potential takes a double-well potential form similar to that in (\ref{psilandau}), but the entropic and enthalpic contributions are accounted for explicitly so that the energy.
The second term in eq. (\ref{F}) describes the cost of forming the interface between the two phases; the interfacial tension between the two phases is proportional to the coefficient $\kappa$ as we show below. Finally, the last term in (\ref{F}) accounts for the interaction between the fluids and bounding surfaces; this is through this term one could account for the specific wetting interactions in these systems. This energy could be written as a power series in $\phi_{s}$, where $\phi_{s}$ is the value of the order parameter on the surface[18,19- Yongting], however, typically only linear term is taken into account [13,20,21,22]
\begin{equation}\label{psis}
\Psi_{s}(\phi_{s})=-h\phi_{s}
\end{equation}
Under the isothermal conditions, the dynamics of the binary fluid is described by the system of equations encompassing the continuity equation, the convection-diffusion equation for the order parameter and Navier-Stokes equation
and {\cite{Landau}}. This system of equation ensues conservation of total density, order parameter and the fluid momentum, respectively, and and can be written as
\begin{equation}\label{econt0}
\partial_t \rho + \bf \nabla \cdot (\rho \bf u) =0.
\end{equation}
\begin{equation}\label{NS}
\partial_t (\rho \bf u) + \bf \nabla \cdot (\rho (\bf u \bf u)
=
- \nabla \cdot \bf P + \eta \bf \nabla^2\bf u ,
\end{equation}
\begin{equation} \label{diff}
\partial_t (\phi) + \bf \nabla \cdot (\phi \bf u)=
M \bf \nabla^2\mu
\end {equation}
where $\bf u$ is the fluid velocity, $\eta$ is the viscosity and $M$ is the mobility of the order parameter. The chemical potential $ \mu$ and pressure tensor ${\bf P}$
are determined based on the chosen expression for the free energy. Taking the functional derivative of the free energy given in \ref{F}, one derives both, the chemical potential and the boundary conditions on the surface on the box. For the bulk free energy expression given in \ref{psilandau}, the chemical potential is written as \cite{Julia,Chris,Briant}:
\begin{equation} \label{muphi}
\mu = \delta F / \delta \phi = -a \phi+ b\phi^3-k \partial_\alpha \partial_\alpha \phi
\end{equation}
Finally, the pressure tensor for the same free energy can be written as \cite{}:
\begin{equation}\label{Pab}
P_{\alpha \beta} =
p_0\delta_{\alpha \beta} +
k (\partial_ \alpha \phi \partial_\beta
\phi -
1/2 \partial_ \gamma \phi \partial_
\gamma
\phi \delta_{\alpha \beta}
-\phi \partial_ \gamma \partial_ \gamma
\phi \delta_{\alpha \beta})
\end{equation}
where
\begin{equation} \label{po}
p_0 =
-\psi(\phi,\rho)+
\rho \frac{\partial \Psi(\phi,\rho)}
{\partial \rho} +
\phi \frac{\partial \psi(\phi,\rho)}
{\partial \phi}=\frac{c^2\rho}{3}-\frac{a}{2}\phi^2+\frac{3b}{4}\phi^4
\end{equation}
Notably, while the LBM studies focusing on the hydrodynamics of phase separation are relatively recent, the same evolution equation as in \ref{diff}, also refereed to as a Cahn-Hilliad equation, with the free energy given by the eq. (\ref{psilandau}) and the total density of the system set to the constant value, was used to model the dynamics of phase separation with conserved order parameter in a number of earlier studies focusing on the phase-separation in the system in the absence of the hydrodynamic effects.
The above systems of dynamic equations should be supplemented by the respective boundary conditions; these boundary conditions are formulated based on the problem in hand. Multiphase LBM is successfully used to model a class of problems that focuses on the phase separation in bulk for the uniform mixtures quenched into the phase-separating region; in this case, periodic boundary conditions are used. In this class of problems, it is important to identify possible finite size effects on the simulation results. The conclusions from the simulations, for example, a domain growth exponent calculated from the simulations, should not depend on the size of the simulation box. To make sure the results of the simulations indeed are not affected by the finite size effect one needs to ensure that the domain size remains significantly smaller than the size of the simulation box and that the growth dynamics remains unchanged within the considered set of parameters when one increases the box size. In many cases, binary fluid is bounded by the surrounding surfaces in one or more dimensions. Figure 1 gives the reader a few examples of the implementaion of LBM approaches to solve various problems in the multiphase fluids. The images in the left panel show in those cases, the last term in eq. \ref{F} allows one to specify the boundary conditions for the order parameter on these surfaces.
If the surface potential is taken in the form \ref{psis}, the boundary conditions for the order parameter can be written as \begin{equation}\label{bcphi}\bf n \cdot \bf \nabla \phi = -h/\kappa
\end{equation},where $\bf n$ is the unit vector normal to the surface bounding the fluid. The parameter $h$ is a model parameter that allows one to set the static contact angle, $\theta$, at the surface using the following relationship
\begin{equation} \label{h}
h = \sqrt{2kb} \cdot sgn(\pi/2-\theta)\Big[cos(\alpha/3)(1-cos\alpha/3)\Big]^(1/2)
\end{equation},where $\alpha=arccos(sin^2\theta)$ and sgn(y) gives the sign of y. The details of the derivation of this equation as well as it's implementation for the case of liquid-vapor interface in a contact with the solid surface are given in \cite{BriantContactPTR}; in essence, by means of setting the value of $h$, one can define the boundary conditions that are capable of accurately reproducing Young's equation in equilibrium. Herein, we remind the reader that the contact angle is the angle that the liqid/vapor or the interface between the two fluids forms with solid surface; the case of $\theta=0$ corresponds to the complete, and the case of $<0\theta=0<\pi$ corresponds to the partial wetting. It was shown that the accuracy of setting the actual contact angle through the eq. \ref{bcphi} and implementing this boundary condition as provided in \cite{BriantContactPTR} is within $3\%$ for the large range of values of $30^0<\theta<150^0$. For the velocity field, one often chooses no-slip boundary conditions (zero velocity at the walls). The density of the binary fluid we are focusing in herein remains close to the constant value (chosen to be equal to unity); for the cases of the liquid-vapor binary systems, one would need to impose the respective densities at the bounding surfaces. It is important to emphasize that in equilibrium, the chemical potential in the system $\mu$ is equal to zero.
It is important to consider an equilibrium solution across the fluid interface, and then compare the numerical solution one can achieve with LBM algorithm with this analytical solution. One can find an exact analytical solution of the order parameter distribution across the fluid interface by analytically solving the equation $\mu=0$ in one dimension across the interface. This solution reads
\begin{equation} \label{phieq}
\phi_{eq}=\sqrt{\frac{a}{b}} tanh \frac{(x-x_0)}{\sqrt{2}\xi}
\end{equation}where $x_0$ marks the center of the interface, and the value of $\xi =\kappa/a$ defines the interface thickness. From the eq. \ref{phieq}, the order parameter varies from $-\sqrt(a/b)$ for the pure A phase to $+\sqrt(a/b)$ for the pure B phase. For simplicity, in majority of the simulations studies one typically chooses the value of $b$ to be equal to $a$ so that the order parameter varies from "-1" to "+1" in the majority of simulations at a constant temperature. The interfacial tension of the system per unit length of the interface is defined as an integral of the total energy over the entire width of the interface (perpendicular to the interface) when the equilibrium is reached. Hence, to calculate the interfacial tension in our system, we integrate over the interface assuming that the coordinate $x$ is taken in the direction normal to the fluid interface and setting the origin of the coordinate system to the middle of the interface between the two phases so that $x_0=0$. Integrating across the interface taking into account the above equilibrium solutions ref{phieq}, we find the interfacial tension as
\begin{equation}\label{sigma}
\sigma=
\int_{-\infty}^{+\infty} dx\left [ -\frac{a}{2} \phi_{eq}^2(x)+\frac{b}{4} \phi_{eq}^4(x) -\frac{\kappa}{2} (\partial_x\phi_{eq}(x)))^2\right ]=\sqrt{}\frac{8\kappa a^3}{9b^2}
\end {equation}
One can easily show that the two contributions to the total energy in eq.(\ref{sigma}), the contribution from the local, bulk free energy and the contribution from (first two terms) and the gradient terms are equal. The
LMB formulation effectively allows us to find the numerical solution of the above system of equations with the given boundary conditions.
In multiphase LBM, the physical variables (density, order parameter and velocity) are described through the two sets of distribution functions that are discrete in space and time as shown below, which is different from the single set of distribution functions we discussed at the beginning of this chapter. We will mainly focus on the cases of binary fluids with the same densities and viscosities and on the single relaxation time approach for these systems. We will however comment on the extension of this approach to the MRT scheme by the end of this section.
Finally, we note that while herein we focus only on the dynamics of the binary fluid, this formulation can be extended in the straightforward manner to model the behavior of the ternary fluid; in this case, the total free energy of the system needs to be extended to the triple-well potential, where each of the potential minima corresponds to one of the three fluids and the interfacial energies need to be defined between all the components within the system. \cite{ternary}.
\section{Lattice Boltzmann Algorithm}
\label{C2}Similar to the single component LBM approach introduced above, LBM simulations for multicomponent fluids consist of two steps, first is collision of fluid "particles" of two types at each lattice cite of the regular lattice and the second is propagation of these particle to the neighboring sites based on the collision rules. Correspondingly, two distribution functions are now defined, $f_i (\bf{r})$ and $g_i
(\bf{r})$ on each lattice site $\bf{r}$. We will provide here an example of a nine-
velocity model defined in two dimensions with velocity vectors
${\bf e}_i=(\pm 1,0),(0,\pm 1), (\pm 1, \pm
1), (0,0)$. The conserved physical variables, density $\rho (\bf x,t)$, order parameter $\phi (\bf x,t)$ and momentum $\bf j (\bf x,t)=\rho (\bf x,t) \bf u (\bf x,t)$ are calculated through these distribution functions as:
\begin{equation}
\rho=\sum_i f_i, \qquad \rho u_\alpha =
\sum_i f_i e_{i\alpha},
\qquad \phi = \sum_i {g_i}.
\label{eq1}
\end{equation}
The time evolution equation for these two particle distribution functions are calculated again in two steps, collision and a streaming step. During the collision step, the values of the distribution functions are updated based on the collision operator (here, we use the BGK approximation)
\begin{equation}\label{fi}
f_i^{'}({\bf x},t)=f_i(\bf x,t)-\frac{f_i-f_i^{eq}}{\tau_{\rho}}
\end{equation}
and
\begin{equation}\label{gi}
g_i^{'}({\bg x},t)=f_i(\bf x,t)-\frac{g_i-g_i^{eq}}{\tau_{\phi}}
\end{equation}
The streaming step that follows the collision step moves the particles along the corresponding lattice velocity directions, $\bf e_i$as
\begin{equation}
f_i({\bf x},t)+{\bf e} \bf \Delta x,t+\Delta t)=f_i^{'}({\bf x},t)
\end{equation}
and
\begin{equation}
g_i({\bf x},t)+{\bf e} \bf \Delta x,t+\Delta t)=g_i^{'}({\bf x},t)
\end{equation}
The functions $f_i^{eq}$ and $g_i^{eq}$ in the above equations denote the corresponding equilibrium distribution functions. These functions are defined as a power series in the local velocity and have to satisfy a number of constraints, First, the conservation of number density for each component and the conservation of momentum for the bulk are imposed as following:
\begin{equation}\label{constr1}
\sum_i f_i^{eq} = \rho,\qquad \sum_i f_i^{eq}
e_{i \alpha} = \rho
u_{\alpha}, \qquad \sum_i g_i^{eq} = \phi.
\end{equation}
In addition, the following constraints are also imposed \cite{Swift},
\cite{Colin}:
\begin{equation}\label{constr2}
\sum_i f_i^{eq} e_{i\alpha}e_{i\beta} =
P_{\alpha\beta}+\rho
u_\alpha u_\beta,
\label{eq6}
\end{equation}
\begin{equation}\label{constr4}
\sum_i g_i^{eq} e_{i \alpha} = \phi
u_{\alpha}, \quad \sum_i g_i^{eq} e_{i\alpha}
e_{i\beta} =
\Gamma \mu \delta_{\alpha\beta}+\phi
u_\alpha u_\beta.
\end{equation}
With the above constraints, expansion of eqs. ({\ref{eqf}})-
({\ref{eqg}}) recovers the above system of equations with the following values of the kinematic viscosity
\begin{equation}
\nu=\Delta t \frac{c^2}{3}(\tau_{\rho}-\frac{1}{2})
\end{equation}
and the mobility of the order parameter
\begin{equation}
M=\Delta t \Gamma (\tau_{\phi}-\frac{1}{2})
\end{equation}
The additional term (last term on the r.h.s. of
the above equation) is common to these Lattice
Boltzmann schemes (\cite{Swift}), but it
is small compared to the other terms.
The equilibrium distribution functions can be
written as a polynomial expansions to satisfy all constraints listed above. An example of the set of the distribution function is given in the reference
\subsection{Spurious velocities}
One of the well-known problems in simulating fluids using multiphase LBM techniques is an existence of unphysical flows near the interfaces when the system reached an equilibrium due to the discretization errors. These currents are referred to as spurious currents or spurious velocities. Below, we review the spurious velocities appearing in the free-energy LBM only, for comparison of spurious velocities in different multiphase LBM approaches we refer the reader to the recent review \cite{Connington} and references therein. An example of spurious velocity field around a liquid droplet in equilibrium with its vapor phase is shown in Fig. SpuriosA. The simulations in this image were conducted using the standard free energy lattice-gas LBM approach, the density field within the droplet is higher and corresponds to the liquid phase (more specifically, this density is set at NN in the given example) and the density of the outer vapor phase corresponds to $NN$. Herein, the size of the arrow corresponds to the magnitude of the velocity so that one can clearly see that the maximum spurious velocities appear at the droplet interface, the flow field forms eight vortexes and the maximum magnitude of these spurious velocities are on the order of $10^{-4}$ dimensionless LB units. What one ideally expects is that when the system reaches equilibrium the system will be at rest in the absence of thermal fluctuations, hence the velocities should approach zero values.
Similar spurious velocities are observed for the binary fluid. An example of the distribution of the spurious velocity field for the droplet of fluid A within the fluid B is given in \cite{WagnerSpurious}, the density of the fluid changes across a few lattice sites across the interfacial region according to the eq. \ref{phieq} from the value close to “-1.” within the center of the droplet to the value close to “+1” in the outer region. The maximum spurious vectors again could reach the values as high as $ 10^{-4}$ dimensionless lattice units, and the spurious velocities again form eight vortexes around the droplet \cite{WagnerSpurious}. We note that the eq. \ref{muphi} is also satisfied only approximately due to the discretization errors, so that after reaching the equilibrium instead reproducing $\mu=0$ solution the LBM approach gives the non-zero, not-uniform distribution of $\mu$ with additional error terms due to the discretization errors‚ the density plot showing $mu$ at equilibrium is also given in \cite{WagnerSpurious}.
The magnitude of the spurious velocities depends on the specifics of the LBM implementations, on the chosen equilibrium distribution functions, surface tension, viscosity as well as on the curvature of the interface. Notably, for the purely diffusive systems, where one essentially needs to integrate only the Chan-Hilliard equation and the hydrodynamics is not taken into account, the LBM approach can be re-formulated with the spurious velocities reduced to the machine precision \cite{WagnerSpurious}. Furthermore, it was also shown that the spurious velocities will reappear in the same system when the hydrodynamics is reintroduced but can again be reduced almost to the machine precision by modifying the forcing term in the Navier-Sockes equation; however, this reduction came at a cost: the momentum was no loner conserved and the method was no longer stable unless additional corrections in viscosity was applied.
A significant reduction in the magnitude of the spurious velocities was achieved for both liquid-gas [] and binary liquid [] systems in the recent work by the optimal choice of stencils for calculating gradients and Laplacian in a combination with the corresponding optimization of the equilibrium distribution functions. These choices resulted in developing numerically stable schemes that satisfied momentum conservation. While the spurious velocities were not reduced to the machine precision, they were reduced significantly for the system in Fig. 1 with these improvements, see Figure 1B. The justification of the optimal choice of the stencils for the gradient and Laplacian are given in ref
\begin{equation}
\label {modified gradients}
\partial_x=\frac{1}{12\Delta x}
\begin{matrix}
-1 & 0 & 1
-4 & 0 & 4
-1 & 0 & 1
\end{matrix},
\Delta^2 =\frac{1}{6\Delta x^2}
\begin{bmatrix}
1 & 4 & 1\\
4 & -20 & 4\\
1 & 4 & 1
\end{bmatrix}
\end{equation}
Such choices improve the isotropy of calculating the corresponding operators and hence lead to the reduction in spurious velocities.
In addition to the spurious velocities appearing at the interfaces,
the spurious velocities would also appear at the surfaces, for example while dealing with the problems simulating the contact line motion along the surface. It was however shown that for the binary fluid with the viscosities of both phases been equal, the standard choice of the equilibrium distribution functions and with the motion of the contact line can be simulated with the high accuracy. In many fundamentally and technologically relevant problems, however, one needs to focus on the binary fluids with the different viscosities.
We note that with the single relaxation time approach, only the binary fluid with the same viscosities are considered; to account for the difference in the fluid viscosities one needs to implement MRT approach.
\begin{thebibliography}{99}
\bibitem{Knight} J. Knight, New Scientist
{\bf 159}, 28, (1998).
\bibitem{study1} S. Dietrich, J. Phys.: Condens.
Matter {\bf 8}, 9127, (1996).
%\bibitem{DeGennes} P. G. de Gennes: Rev. Mod.
Phys. {\bf l57}, N3.
\bibitem{Bauer} C. Bauer and S. Dietrich, Phys.
Rev. E {\bf 60}, 6119, (1999).
\bibitem{study2} S. Puri and S. L. Frisch. J.
Phys.: Condens. Matter {\bf 9}, 2109, (1997).
\bibitem{study3} S. Bastea, S. Puri and J. L.
Lebowitz, Phys. Rev. E {\bf 63}, 041513, (2001).
\bibitem{SC92} M. K. Chaudhury, G. M. Whitesides,
Science {\bf 256}, 1539, (1992).
\bibitem{PatternedSubstrates} A. Karim et al,
Phys. Rev. E; {\bf 56}, R6273, (1998).
%\\K.Fukunaga, H.Elbs, G.Krausch, Langmuir,
{\bf 16}, 3474, (2000)
\bibitem{Rowlinson} J.S.Rowlinson, B. Widom,
{\it Molecular Theory of Capillarity},
Clarendon Press: Oxford, (1982).
\bibitem{Kataouka} D. Kataouka; S. Troian.
Nature, {\bf 402}, 794, (1999).
\bibitem{Papra} A. Papra, et al, Langmuir
{\bf 17}, 4090, (2001).
\bibitem{Bain} C. Bain; G. D. Burnetthall, R.
R. Montgomerie, Nature, {\bf 372}, 414,(1994).
\bibitem{Lee} S. W. Lee, P. E. Laibinis, JACS,
{\bf 122}, 5395, (2000).
\bibitem{Gau} H. Gau, S. Herminghaus, P. Lenz,
R. Lipowsky, Science {\bf 283}, 46, (1999).
\bibitem{Gallardo} B. S. Gallardo et al, Science
1999, {\bf 283}, 57, (1999).
\bibitem{Rascon} C. Rascon, A. O. Parry, Nature
{\bf 407}, 986, (2000).
\bibitem{Chen} S. Chen and G. D. Doolen, Annual
Rev. Fluid Mech. {\bf 30}, 329, (1998).
\bibitem{Martys} N. S. Martys, H. Chen, Phys.
Rev. E {\bf 53}, 743, (1996).
\bibitem{emulsion} We note that for the
pharmaceutical and cosmetics industries,
microfluidic assays commonly require the
formation of stable emulsions of immiscible
fluids, such as oil and water.
\bibitem{Wilbur} J. L. Wilbur et al.,
Nanotechnology {\bf 7}, 452 (1996).
\bibitem{beebe} B. Zhao, J. Moore, D. J.
Beebe, Science, {\bf 291}, 1023, (2001),
and references therein.
\bibitem{Landau} L. D. Landau and E. M.
Lifshitz, {\it Fluid Mechanics, 2nd Ed.},
Pergamon Press, Oxford, (1987).
\bibitem{Swift} M. R. Swift, E. Orlandini,
W. R. Osborn, J. M. Yeomans, Phys. Rev. E
{\bf 54}, 5041, (1996).
\bibitem{Anna} A. C. Balazs et al., J. Phys.
Chem. B {\bf 104}, 3411, (2000).
\bibitem{OurLangm} O. Kuksenok, J. M. Yeomans,
A. Balazs, Langmuir, in press.
\bibitem{Onuki} A. Onuki, J. Phys.: Condens.
Matter {\bf 9}, 6119, (1997).
\bibitem{Colin} C. Denniston, E. Orlandini,
J. M. Yeomans, Europhys. Letts. {\bf 52}, 481, (2000).
\bibitem{Colin} A. Bray, Advances in Physics, 2002, Vol. 51, No. 2, 481, (2002)
\bibitem{D. M. Anderson, Annu. Rev. Fluid Mech. 1998. 30:139â65
\bibitem{Diffusive2004} P. Yue at al, J. Fluid Mech. (2004), vol. 515, pp. 293â317.
\end{thebibliography}
\newpage
\begin{figure}[h!]
\centering
\includegraphics[scale=1.7]{universe.jpg}
\caption{The Universe}
\label{fig:univerise}
\end{figure}
\section{Conclusion}
``I always thought something was fundamentally wrong with the universe'' \citep{adams1995hitchhiker}
\end{document}