# Phase-field model for the Rayleigh–Taylor instability of immiscible fluids

###### Abstract

The Rayleigh–Taylor instability of two immiscible fluids in the limit of small Atwood numbers is studied by means of a phase-field description. In this method the sharp fluid interface is replaced by a thin, yet finite, transition layer where the interfacial forces vary smoothly. This is achieved by introducing an order parameter (the phase field) whose variation is continuous across the interfacial layers and is uniform in the bulk region. The phase field model obeys a Cahn–Hilliard equation and is two-way coupled to the standard Navier–Stokes equations. Starting from this system of equations we have first performed a linear analysis from which we have analytically rederived the known gravity-capillary dispersion relation in the limit of vanishing mixing energy density and capillary width. We have performed numerical simulations and identified a region of parameters in which the known properties of the linear phase (both stable and unstable) are reproduced in a very accurate way. This has been done both in the case of negligible viscosity and in the case of nonzero viscosity. In the latter situation only upper and lower bounds for the perturbation growth-rate are known. Finally, we have also investigated the weakly-nonlinear stage of the perturbation evolution and identified a regime characterized by a constant terminal velocity of bubbles/spikes. The measured value of the terminal velocity is in perfect agreement with available theoretical prediction. The phase-field approach thus appears to be a valuable tecnhique for the dynamical description of the stages where hydrodynamic turbulence and wave-turbulence enter into play.

A. Celani, A. Mazzino, P. Muratore-Ginanneschi and L. Vozella] ANTONIOCELANI, ANDREAMAZZINO, PAOLOMURATORE-GINANNESCHIand LARAVOZELLA

## 1 Introduction

The Rayleigh-Taylor (RT) instability is a fluid-mixing
mechanism occurring when a heavy, denser, fluid
is pushed into a lighter one. For a fluid in a gravitational field,
such a mechanism was first discovered by Lord Rayleigh in the 1880s
Rayleigh (1883)
and later applied to all accelerated fluids by
Sir Geoffrey Taylor in 1950 Taylor (1950).
The relevance of this mixing mechanism embraces many different
phenomena occurring in completely different contexts.
We just mention, among the many,
astrophysical supernova explosions and geophysical formations
like salt domes and volcanic islands Di Prima & Swinney (1981); Dimonte & Schneider (2000), continental magmatism caused by lithospheric gravitational instability
Lee, Rudnick & Brimhall Jr. (2001); Ducea & Saleeby (1998), inertial confinement fusion Cook & Zhou (2002) and cloud formation in atmospheric sciences Schultz et al. (2006).

Back to classical fluids applications, RT instability is the first step
eventually leading to a fully developed turbulent regime.
A deeper understanding of the mechanism of flows driven by RT instability
thus would shed light on the many processes that underpin fully developed
turbulence.

The difficulty inherent in sustaining an unstable density stratification
has challenged experimentalists for over half a century. Several innovative
approaches have been recently developed (see e.g., Ramaprabhu & Andrews 2004).

With the advent of supercomputers, high-resolution numerical
simulations of RT at high Reynolds numbers have become a reality. However,
simulations using many different benchmark codes and
experiments disagree
already on apparently innocent observables like,
for instance, the value of the
growth constant, , associated to the spread of the
turbulent mixing zone (see, e.g., Di Prima & Swinney 1981). The differences can be as high as .

Despite the long history of RT turbulence, a consistent phenomenological theory has been presented only very recently by Chertkov (2003) for the miscible case. The theoretical predictions by Chertkov have been verified by Celani, Mazzino & Vozella (2006) exploiting numerical simulations in two spatial dimensions. For the three-dimensional miscible case we refer, e.g, to Young et al. (2001).

In many of the aforementioned situations where the RT instability
has an important role, the two fluids are immiscible owing to a non
negligible surface tension. At level of linear analysis the role played
by a non zero surface tension was addressed by Chandrasekhar (1961). The successive
dynamics falling in a
turbulent regime has been recently analyzed by Chertkov, Kolokolov & Lebedev (2005).
Using a phenomenological approach, the authors suggest
the existence of a Kolmogorov cascade between the integral scale
and a time-dependent scale related to the
typical drop size. Below the latter scale, associated
to an emulsion-like region,
a wave energy cascade takes in. This is mediated by weakly
interacting capillary waves propagating on top of the drop surface.
Eventually, the energy is dissipated by viscous forces.

RT instability and RT turbulence of immiscible fluids
thus appear richer than the corresponding miscible situations.
The existence of two different cascades poses a serious challenge to
numerical investigations of the immiscible RT problems. The emulsion-like
phase indeed occurs at very small scales and the energy transfer
takes place on the interfaces. These are geometrical objects
close to singularities and thus difficult to describe appropriately in
a numerical scheme.
Accuracy and efficiency are thus fundamental requirements
to reproduce the correct statistical features characterizing
immiscible RT turbulence.

Our aim here is to perform a first step along this direction by focusing on direct numerical simulations of immiscible RT instability. The numerical strategy we exploit here is known as phase-field model Bray (2002); Cahn & Hilliard (1958); Badalassi, Ceniceros & Banerjee (2003); Ding, Spelt & Shu (2007). The main idea of the method is to treat the interface between two immiscible fluids as a thin mixing layer across which physical properties vary steeply but continuously. The evolution of the mixing layer is ruled by an order parameter (the phase field) that obeys a Cahn-Hilliard equation Cahn & Hilliard (1958). The method permits to avoid a direct tracking of the interface and easily produces the correct interfacial tension from the mixing layer free energy.

We present here an accurate numerical study that validates the phase–field approach by testing known results of immiscible RT instability both at level of linear and weakly nonlinear analysis. From our results, it turns out that this strategy is a valuable option for a quantitative treatment of the turbulent regime characterized by the interplay between hydrodynamic and interface degrees of freedom.

The paper is organized as follows.
In Sec. 2 we introduce the Rayleigh–Taylor problem and discuss
the related phase-field approach.
A detailed analysis
of the
energy balance between purely hydrodynamic degrees of freedom and interface
degrees of freedom is presented. Finally, the dispersion relation
for gravity-capillary waves is obtained by analytical calculations
starting from the phase-field equation coupled to the Navier–Stokes equations.

In Sec. 3 the results from the direct numerical simulations are presented and
compared with known results for the linear analysis. We focus both on the case of zero viscosity and on that of negligible viscosity. Both
stable and unstable configurations are considered.
Finally, the weakly nonlinear
regime is considered and the resulting terminal velocity of bubbles/spikes
compared with existing theoretical predictions.

Sec. 4 is devoted to some conclusions and perspectives.

## 2 System configuration and phase-field model

Our system consists of two immiscible, incompressible fluids (labeled by 1 and 2) having different densities, and , with the denser fluid placed, e.g., above the less dense one (see Fig. 1).

In the absence of gravity,
this flow configuration is stable.
In presence of the gravitational force, surface tension may be able
to keep the system in equilibrium, provided the density contrast
is not too large.

Let us start by describe the equilibrium configuration and then
pass to the
evolution (RT instability) that occurs
when a perturbation is imposed to the interface separating the two
fluids.

### 2.1 Equilibrium state

Let us consider an equilibrium state where fluid 1 is placed below fluid 2 and they are separated by a sharp interface. The fact that the interface is sharp (i.e. a discontinuity in the fluid properties) poses a serious challenge to numerical simulations. Indeed, for sharp interfaces, the evolution equations are obtained by following fluid 1 and 2 separately with the appropriate boundary condition at interface (see, for instance, Smolianski, Haario & Luukka 2005; Sethian 1999). Other approaches follow the interface alone. In this latter case, the movement of the interface is naturally amenable to a Lagrangian description, while the bulk flow is conventionally solved in an Eulerian framework. These approaches employ a mesh that has grid points on the interfaces and deforms according to the flow. A major shortcoming of these approaches is in that they cannot handle properly topological changes such as breakup, coalescence and reconnections (see Yue et al. 2004, and references therein). In this respect, the phase–field method is, by far and large, more effective, at the expense of a larger number of grid points required.

The idea of the phase–field method is to replace the sharp interface with a diffuse one in such a way that the numerical computation of interface movement and deformation can be carried out on fixed grids Anderson, McFadden & Wheeler (1998); Jacqmin (1999). More quantitatively, this amounts to assigning to the system a Ginzburg–Landau free energy, , espressed in term of the order parameter as Cahn & Hilliard (1958); Bray (2002); Yue et al. (2004):

(2.0) |

where is the region of space occupied by the system,
is a mixing energy density
and is the capillary
width, representative of the interface thickness.
The order parameter is a field which
serves to identify fluid 1 and 2. We assume
in the region occupied by fluid 1 and
in those where fluid 2 is present.

The equilibrium state is the minimizer of the free energy .
The mechanism which keeps the system in this configuration is
the competition between two effects due to the two addends
in (2.1).
The first term
favours a perfect mixing (i.e. in ,
this term being the interface energy contribution) whereas the second one
one drives the system towards demixing (the associated
term in , the
bulk contribution, has indeed a minimum for ).
The nontrivial final equilibrium
state is just the results of this competition.
More quantitatively, the final state is obtained by minimizing the free-energy
functional with respect to variations of the function , i.e., solving:

(2.0) |

where is the so-called chemical potential (see, for instance, Cahn & Hilliard 1958; Bray 2002; Yue et al. 2004). If one considers an one-dimensional interface, varying along the gravitational direction , one easily finds the solution of Eq. (2.1) as Cahn & Hilliard (1958); Bray (2002); Yue et al. (2004):

(2.0) |

This solution exists and is stable in all dimensions although the decay rate of perturbations depends upon the dimensionality Korvola, Kupiainen & Taskinen (2005). From (2.1) one immediately realizes that the sharp-interface limit is obtained for : in this case . Moreover, the surface tension is equal to the integral of the free-energy density along the interface (see, for example, Landau & Lifshitz 2000). For a plane interface, this integral yields Cahn & Hilliard (1958); Bray (2002); Yue et al. (2004):

(2.0) |

It is now easy to verify how the sharp interface limit is obtained: it suffices to take the limits and to zero keeping fixed to the value prescribed by surface tension Liu & Shen (2003).

### 2.2 Perturbation evolution

Let us now suppose to impose a small perturbation on the (finite thickness) interface separating the two fluids. Such perturbation will displace the phase field from the previous equilibrium configuration, which minimized the free-energy , to a new configuration for which in general, . The system will react so as to try to reach again an equilibrium configuration. In formulae:

(2.0) |

being the so-called mobility (see, for instance, Bray 2002; Yue et al. 2004).
Notice the presence of the Laplacian
operator in front of .
Notice that the mass of each fluid
is conserved, as imposed by the physics of
the problem under consideration.

The dynamics of the velocity field
is governed by the usual Boussinesq Navier-Stokes equations
Kundu & Cohen (2001) plus
an additional stress contribution arising at the interface where
the effect of surface tension enters into play Bray (2002); Yue et al. (2004); Berti et al. (2005).
The equations of motion are:

(2.0) | |||||

(2.0) |

In the first equation and is the kinematic viscosity. The quantity is the coupling term that accounts for capillary forces. It is easy to verify that it can be rewritten as plus a gradient term which can be absorbed into the pressure term. Finally, is the buoyancy contribution, being the deviation of the actual density, , from the mean density :

The buoyancy contribution can be rewritten in terms of , and as:

(2.0) | |||||

where is the Atwood number.

### 2.3 Energetics

Let us define the kinetic energy (per unit volume), , and the potential
energy (per unit volume), , for our system ruled by Eqs. (2.2), (2.2) and (2.2).

By definition of potential energy, we have:

(2.0) | |||||

being the total volume occupied by the fluids and brackets,
, denote spatial averages. In Eq. (2.0)
the constant is chosen such to set the potential energy to zero for
vanishing Atwood number.

In a similar way, one can define the kinetic energy per unit volume as:

(2.0) | |||||

From Eqs. (2.2) and (2.2) we immediately realize that such equations are left invariant under the simultaneous transformation , . As a consequence, and the resulting kinetic energy simply reads:

(2.0) |

By defining the total energy of the two-fluid system is

The equation for is obtained by multiplying Eq. (2.2) by and then taking spatial average. We easily get:

(2.0) |

Let us now take Eq. (2.2), multiply it by , and take the average:

(2.0) |

by translational invariance and Leibniz rule. We thus have:

(2.0) |

where we have used the fact that . The free–energy variation is

(2.0) | |||||

i.e.,

(2.0) |

where we have introduced the strain tensor and assumed boundary
conditions suitable to justify integrations by parts.

The energy balance takes then the form:

(2.0) |

The global system in thus intimately dissipative,
even for a vanishing kinetic viscosity.

It is worth emphasizing the cancellation of
by the kinetic and the free–energy contributions, due to
exchanges between the velocity field and the interface.

### 2.4 Dispersion relation for the phase-field model

The aim of this section is to show that the well-known
dispersion relation for gravity-capillary waves Chandrasekhar (1961) can
be easily obtained within the phase-field formalism. To do that,
let us
concentrate our attention on a two-dimensional problem and indicate
by the gravity direction.
Moreover, we will assume heavier fluid to be placed below the lighter
one, in a way to have a stable situation.
For a given perturbation imposed to the interface,
the problem is to determine how the perturbation evolves in time.

Denoting by a small perturbation imposed to a planar interface, we
can rewrite as:

(2.0) |

where can be larger than , yet it has to
be smaller than the scale
of variation of (small amplitudes).

Locally, the interface is in equilibrium, i.e.:

(2.0) |

where . In this limit we have:

(2.0) |

Linearizing Eq. (2.2) for small interface velocity we have, neglecting the viscous term:

(2.0) |

The integration in the vertical direction interpreted in the principle value sense

(2.0) | |||||

(2.0) |

yields:

(2.0) |

having used the relations , and

(2.0) |

The height variation of the interface has to match the vertical fluid velocity, thus giving:

(2.0) |

The last step is to relate the velocity at the interface with the integral . This is done by restricting to potential flows:

(2.0) |

For , denoting with “” the Fourier Transform, we have:

(2.0) |

(2.0) |

(2.0) |

(2.0) |

Therefore:

(2.0) |

so that in space we have:

(2.0) |

From these two equations we immediately get:

(2.0) |

with:

(2.0) |

that is the expected dispersion relation Chandrasekhar (1961).
For the stable configuration we have, for all values of :
,
i.e. any initially imposed perturbation will not grow indefinitely.

From Eq. (2.0) and the initial condition:

(2.0) |

we immediately have:

(2.0) |

and the velocity at the interface reads:

(2.0) |

Assuming an initial perturbation of the form , from Eqs. (2.0) and (2.0) we obtain:

(2.0) |

and the velocity components, for , read:

(2.0) | |||||

(2.0) |

where we used the relation .

For , in a similar way we obtain the velocity field components:

(2.0) | |||||

(2.0) |

When in the initial configuration the heavier fluid placed above the lighter one, the dispersion relation (2.0) trasforms in:

(2.0) |

which is readly obtained by flipping the sign of . For surface tension is not able to contrast gravity-induced vertical motion with the final result that amplitude perturbations grows exponentially: the flow is unstable. More precisely, from relation (2.0) and for we have:

(2.0) |

and Eqs. (2.0) - (2.0) transform in:

(2.0) | |||||

(2.0) |

for , and:

(2.0) | |||||

(2.0) |

for .

## 3 Numerical investigation

In this section we report results we have obtained exploiting
direct numerical simulations (DNS) of the phase-field
model for the Rayleigh–Taylor problem
described in the preceeding sections.
Our attention will be focused
both on the linear phase of the perturbation evolution and on the
weakly nonlinear regime governed by plumes, for .

In the present study we will
consider initial perturbations imposed to the interface varying along
one of the horizontal directions, say the -axis, and invariant along the
other horizontal direction, say the -axis.
The perturbation is thus intimately two-dimensional
a fact that allows us to solve the original Navier–Stokes equations coupled
to the phase field in two dimensions. This clearly permits to obtain
high accuracy and thus to properly test the phase-field approach against
known results for both the linear and the nonlinear evolution stage.

For a two-dimensional flow it is
convenient to introduce the vorticity field
[]
and study
the equations

(3.0) |

(3.0) |

In order to efficiently and accurately solve those equations we exploit a
pseudospectral method C88 (). Accordingly, periodic boundary conditions
have to be assumed along the two directions. For the horizontal direction
it is a natural choice (see e.g. Cabot & Cook (2006); Liu & Shen (2003)) while along the vertical one
this choice deserves some comments.
As initial condition we started from the
hyperbolic-tangent profile, Eq. (2.1),
for with
the interface placed in the middle of the
domain.
The fact that we have periodic boundary conditions along
simply means that far from the middle of the domain the
hyperbolic-tangent profile
has to be distorted in order to satisfy periodic boundary conditions.
However, both in the linear and in the weakly nonlinear
regimes the amplitude of the interface perturbation is always much
smaller than
the vertical size of the box, so that the actual choice of boundary conditions at the top and bottom can be safely neglected.

Such a strategy has been
already exploited for the miscible case by Celani, Mazzino & Vozella (2006).

The box has a horizontal to vertical
aspect ratio for the linear analysis stage and
for the weakly nonlinear evolution.
In the latter case we take a smaller aspect ratio owing to the fact that
the perturbation can reach a higher amplitude (with respect to case of the
linear analysis).

In both cases the resolution is collocation
points.
We need such a high resolution (despite the fact that we focus on
a linear and weakly nonlinear study)
in order to have a well described interface separating the two phases.
In our simulations the mixing width () is 6 mesh points.

The time evolution is implemented by a standard second-order
Runge–Kutta scheme.

The physically relevant parameters in the present problem
are the kinematic viscosity
, the buoyancy intensity and the surface tension .
Both and will be
varied in our study, while will be kept fixed to a fixed value (see
below).
The surface tension is related to the ratio with
(and thus )
sufficiently small in order to have a finite value for the
surface tension and, at the same time,
to reproduce the correct sharp-interface limit.
Finally, the parameter appearing in the relaxation term in
Eq. (3.0) must satisfy the requirement that be
small, so as to enforce
‘istantaneous’ local equilibrium between flow and interface. Here we used
the value (model units) .

All simulations presented here start from an initial condition corresponding to an equilibrium configuration: velocity identically zero and hyperbolic tangent profile for the phase field , expressed by the relation of the form: with

For a given we choose the initial amplitude in a way that (where ) is sufficiently small to fall in the linear phase (i.e. ) and is sufficiently large for the wave disturbance to see an almost infinitesimal mixing width (i.e. ). Specific numerical values are reported in the next sections.

### 3.1 Linear instability for negligible viscosity

The aim of this section is to verify the growth-rate (2.0)
which holds in the linear phase when the viscosity is negligible.

In order to do so, we take a small value of ( in the model units)
and vary (up to ,
the critical wave-number separating unstable from stable wave-modes) and
and take a fixed value of . The ratio while
ranges from to in the range of
considered.

The behavior of the square growth-rate
is shown in dimensionless form in Fig. 2
as a function of for three different values
of (obtained by varying )
and in Fig. 3 by varying for three different values of
. In both figures, symbols refer to the numerical results and
the dashed line is the theoretical expectation given by (2.0).

The numerical data in Figs. 2 and 3 have been obtained via best-fit

of , the spatial average of as a function of time. The latter average is computed over a horizontal strip containing the interface (placed in the middle of the computational domain) and having an extension of above and below the interface. This has been done to avoid spurious contaminations coming from the upper and lower domain regions affected by the boundary conditions. In formulae:

(3.0) |

where we used the expression (2.0) and (2.0)
for and , respectively.

The best fit has been done with as unique free parameter
and its high accuracy can be verified in Fig. 4
where we show the time evolution of for
(solid triangles in Fig. 2) and for four values of smaller than .
At nonlinear
effects
start to enter into play giving rise to corrections to the linear analysis
(see Sec. 3.4). Up to that time, linear theory is very accurate as
one can also realize by looking at the insets of Fig. 2
where the sinusoidal form of is reported for .

### 3.2 Linear instability for finite viscosity

The aim of this section is to investigate numerically
how the growth-rate, , is modified by
viscosity. As discussed in Appendix A, both an upper and a lower bound
for the perturbation growth-rate are known (see Eqs. (A 0) and (A 0)) and we
want to assess how the actual growth-rates compare with those.

For such purpose,
we choose a surface tension, , and
in such a way to obtain instability for
few (unstable) wavenumbers. Our choice was
(see Sec. 3.1)
thus corresponding to 5 unstable wavenumbers.

As far as the initial perturbation is concerned, we report here the case corresponding
to . Initial perturbations with a larger wavenumber simply
need an initial smaller amplitude (and eventually a larger numerical resolution)
in order to satisfy
and . Here, we have
and . Such ratios turned out to be sufficiently
‘asymptotic’ to produce accurate results.
The effect of viscosity is studied
by considering twelve values of viscosity in the range (model units).

The results of our simulations are summarized in Fig. 5 where the behavior of the square
perturbation growth-rate, ,
is shown as a function of viscosity. The numerical predictions have been compared with the
available theoretical bounds (dashed lines).

Note that the numerical points are always in between the two bounds and also how the relative differences
between the upper bound and the numerical values are . This latter fact is compatible, for example,
with the results of Menikoff et al. (1977).

The value of the growth-rates have been obtained via best of (see Eq. (3.0)).
Unlike what we did in previous section, here we perform the fit within the exponential region.
The reason is that the non-asymptotic form of the perturbation time-evolution is unknown in the present case.

The fit accuracy can be appreciated in the inset of Fig. 5
where the temporal evolution of the pertubation for (model units) is shown together with the best fit slope
(dashed line) from which is determined.
Error bars, estimated by looking at the fit sensitivity by varying the length of the fit interval,
are of the order of the symbol sizes.

### 3.3 Stable configuration: gravity-capillary waves

The performance of the
phase-field approach in the unstable regime predicted
by linear theory both in the presence
and in the absence of viscosity proved to be very good.
As discussed in Sec. 2.4, for sufficiently
large surface tensions and/or
sufficiently small differences between fluids density,
a perturbation initially imposed
to the fluid interface may maintain
its initial amplitude giving rise to the
dispersion relation (2.0).
The waves resulting from the balance between gravity and surface tension
are known as gravity-capillary waves. Our aim here is
to verify their dispersion relation.

To do that, we have fixed the parameters to obtain a
critical wavenumber of order one.
For (model units) and the same as
in the unstable case, one has . The
first accessible wavenumber is thus stable and should evolve
in time according to (2.0).
However, the geometrical/computational configuration used
in the unstable case did not produce
sufficiently accurate results. In particular, using the same
domain aspect ratio
and the same ratio between
perturbation
amplitude and perturbation wave-length we found a dynamics too dissipative with respect to what is
expected. In the absence of viscosity,
dissipation arises in the phase field formulation due to the sole contribution proportional to
in Eq. (2.0). The latter parameter has been taken sufficiently small to ensure a negligible effects
inside a period of oscillation. The specific value was .
To avoid spurious dissipation, as that induced by nonlinear effects, we reduced the amplitude of the
initial perturbation with respect to the unstable case. Also, we increased the size of the periodicity box
along the gravitational direction in a way to reduce possible spurious contribution arising
from the upper/lower part of the computational domain where instabilities, not present in the unstable case,
might now develop.
The above choice on the amplitude of the inital perturbation
implies a consequent reduction
of . The following set
of parameters have been used: , and a resolution
of . For an initial perturbation on , its initial amplitude has been chosen
to have and .
The behavior of the maximum, , of the initial perturbation is shown as a function of time in
Fig. 6.

The continuous line is relative to a sinusoidal with pulsation obtained from (2.0). The agreement between theory and numerics is satisfactory both for the amplitude and for the pulsation. Note the small reduction of , in one oscillation period: only 1 grid box over 4096.

### 3.4 Weakly non-linear stage

In this section we investigate the early stages of the nonlinear dynamics.
We focus on the rising/falling velocity of plumes
in the limit of small Atwood numbers when spikes and bubbles are known to coincide.
The theoretical prediction for the terminal velocity is reported in
Appendix B. Our aim here is both to verify the existence of a regime
characterized by a costant ‘terminal’ velocity and, secondly, to compare the prediction (B 0)
for such terminal velocity with our numerical data.

The physical parameters are chosen to magnify the effect of the surface
tension on the terminal velocity. This happens when
the wavenumber
of the initial pertubation (still supposed unimodal) is slightly below .
Here we choose and such that and thus
look at the dynamics associated to the
wavenumber . The initial perturbation has an amplitude ;
the initial dynamics is thus linear.
Although we are interested to investigate the case of zero viscosity, in order to prevent
numerical instabilities we add a small viscosity
(model units).
In Fig. 7 the perturbation amplitude is shown as a function of
time: symbols correspond to our numerical data and the dashed line is the prediction (B 0).
A good agreement is found between numerics and theory in the range
.
At larger times, neighboring plumes start to interact
and the arguments leading to (B 0) do not apply any longer.
In Fig. 8
we show some snapshots of the evolution of the two fluids. Figures
are equally spaced in time in the interval .
Black corresponds to ; white to .
Their shape is similar to that experimentally observed.
Note the aforementioned spike/bubble symmetry
corresponding to the up-down symmetry of our original evolution equations.
by Waddell, Nieserhaus & Jacobs (2001).

## 4 Conclusions and perspectives

In this paper we showed that the phase–field model provides
a valuable numerical instrument for the study of immiscible, convective hydrodynamics. As a testground for this model, we have considered the Rayleigh–Taylor
instability. Numerical results compare very well with known analytical results
both for the linearly stable and unstable case,
and for the weakly nonlinear stages of the latter.

All these results are very encouraging in view of the
next important step that is the
the numerical simulation of immiscible RT turbulence. There, the interplay
of all the fundamental mechanisms that we have illustrated here
(instabilities and wave propagation)
is expected to give rise
to a small-scale emulsion-like phase dominated by gravity-capillary waves
and by a large-scale
hydrodynamic range of scales where classical Kolmogorov
turbulence should appear. This theoretical suggestion still awaits numerical
confirmation, and the phase–field model provides the appropriate method
to pursue this goal.

We acknowledge useful discussions with Hekki Haario. AM and LV have been partially supported by PRIN 2005 project n. 2005027808 and by CINFAI consortium (AM). LV acknowledges support from From Discrete to Continuous models for Multiphase Flows TEKES project n. 40289/05.

## A Bounds for the perturbation growth-rate in the presence of viscosity

The effect of viscosity is to reduce the perturbation growth-rate. However it does not remove the instabilities. Analytically, it is more difficult to consider the effect of viscosity with respect to surface tension (see Eq. (115) at page 443 of Chandrasekhar 1961). Nonetheless, it is possible to determine a lower and an upper bound to the growth-rate . These bounds are the solutions to the following equations Menikoff et al. (1977):

(A 0) | |||||

(A 0) |

where is the growth-rate in the inviscid case (see Eq. (2.0)). The solution of Eq. (A 0) is:

(A 0) |

while only a numerical solution is available for Eq. (A 0).

The goodness of those upper and lower bounds are numerically investigated
in Sec. 3.2 by means of the phase-field model.

## B Models for the terminal bubbles/spike velocities in the weekly nonlinear regime

Substantial deviations from the linear theory are observed
when the perturbation amplitude reaches a size of the
order of 0.1 - 0.4 Sharp (1984).

In that case the perturbation evolution is nonlinear.
Then the disturbance grows non-linearly
and the interface starts to deform.
Indeed, at least for finite values of ,
the interface can be divided into spikes corresponding to the regions where
the heavier fluid penetrates into the lighter one, and
bubbles associated to those regions where
lighter fluid rises in the heavier one.
The roll-up of vortices produces a mushroom-type
shape for bubbles and spikes (see, for instance, Waddell, Nieserhaus & Jacobs 2001).
When the fluid densities are similar
(corresponding to our case )
spikes and bubbles coincide and approach a constant and equal velocity.
In both cases, the exponential growth of the velocity perturbation amplitude
characterizing the linear phase of the evolution is
replaced by a linear-in-time behavior Waddell, Nieserhaus & Jacobs (2001).
Two models are available to describe this stage:
the drag-buoyancy model Alon et al. (1995)
and the “Layzer model” Layzer (1955); Goncharov (2003); Young & Ham (2006).
The former model describes bubble and spike motion by balancing the
buoyancy and drag forces and it assumes that this velocities reach a
constant values for sufficiently long times.
The latter model uses an expansion of the perturbation amplitudes and conservation
equations near the tip of bubbles
and spikes.
This
approach has been first
applied to the fluid-vacuum interface ()
Layzer (1955) and then extended to arbitrary Atwood
numbers Goncharov (2003) and to include
the surface tension contribution Young & Ham (2006).
According to the latter study,
in our case (bidimensional flow, immiscible fluids and small Atwood number)
one expects that the terminal bubble
and spike velocity be equal to Young & Ham (2006):

(B 0) |

This expectation is numerically tested, in Sec. 3.4, by exploiting the phase-field method.

- Alon et al. (1995) Alon, U., Hecht, J., Ofer, D. & Shvarts, D. 1995 Power laws and similarity of Rayleigh–Taylor and Richtmyer–Meshkov mixing fronts at all density ratio. Phys. Rev. Lett. 74(4), 534–537
- Anderson, McFadden & Wheeler (1998) Anderson, D. M., McFadden, G. B. & Wheeler, A. A. 1998 Diffuse-interface methods in fluid mechanics. Annu. Rev. Fluid Mech. 30, 139–165
- Badalassi, Ceniceros & Banerjee (2003) Badalassi, V. E., Ceniceros & H. D., Banerjee, S. 2003 Computation of multiphase systems with phase field models. J. Comput. Phys. 190, 371–397
- Berti et al. (2005) Berti, S., Boffetta, G., Cencini, M. & Vulpiani, A. 2005 Turbulence and coarsening in active and passive binary mixtures. Phys. Rev. Lett. 95, 224501-1–224501-4
- Bray (2002) Bray, A. J. 2002 Theory of phase-ordering kinetics. Advances in Physics 51(2), 481–587
- Cabot & Cook (2006) Cabot, W. H. & Cook, A. W. 2006 Reynolds number effects on Rayleigh–Taylor instability with possible implications for type-Ia supernovae. nature physics 2, 562–568
- Cahn & Hilliard (1958) Cahn, J. W. & Hilliard, J. E. 1958 Free energy of a non uniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258–267
- (Canuto et al.) Canuto, C., Hussaini, M. Y., Quarteroni, A. & Zang, T. A. 1988 Spectral Methods in Fluid Dynamics Springer Series in Computational Physics. Springer-Verlag
- Celani, Mazzino & Vozella (2006) Celani, A.,Mazzino, A. & Vozella, L. 2006 Rayleigh–Taylor turbulence in two-dimensions. Phys. Rev. Lett. 96, 134504-1–134504-4
- Chandrasekhar (1961) Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. New York: Dover
- Chertkov (2003) Chertkov, M. 2003 Phenomenology of Rayleigh–Taylor turbulence. Phys. Rev. Lett. 91, 115001-1–115001-4
- Chertkov, Kolokolov & Lebedev (2005) Chertkov, M., Kolokolov, I. & Lebedev, V. 2005 Effects of surface tension on immiscible Rayleigh–Taylor turbulence. Phys. Rev. E 71, 055301-1–055301-4
- Cook & Zhou (2002) Cook, A. W. & Zhou, Y. 2002 Energy transfer in Rayleigh–Taylor instability. Phys. Rev. E 66, 026312-1–026312-12
- Dimonte & Schneider (2000) Dimonte, G. & Schneider, M. 2000 Density ratio dependence of Rayleigh–Taylor mixing for sustained and impulsive acceleration histories. Phys. Fluids 12, 304–321
- Ding, Spelt & Shu (2007) Ding, H., Spelt, P. D. M. and Shu, Chang 2007 Diffuse interface model for incompressible two-phase flows with large density ratios. J. Comput. Phys. 226, 2078–2095
- Di Prima & Swinney (1981) Di Prima, R. C. & Swinney, H. L. 1981 Hydrodynamic Instabilities and the Transition to Turbulence. eds. Swinney, H. L. & Gollup, J. P. Springer, Berlin
- Ducea & Saleeby (1998) Ducea, M. & Saleeby, J. 1998 A case for delamination of the deep batholithic crust beneath the Sierra Nevada, California. Int. Geology Rev. 40, 78–93
- Goncharov (2003) Goncharov, V. N. 2003 Analytical model of nonlinear, single-mode, classical Rayleigh–Taylor instability at arbitrary Atwood numbers. Phys. Rev. Lett. 88(13), 134502-2–134502-4
- Korvola, Kupiainen & Taskinen (2005) Korvola, T., Kupiainen, A. & Taskinen, J. 2005 Anomalous scaling for three-dimensional Cahn-Hilliard fronts. Comm. Pure Appl. Math. 58(8), 1077–1115
- Kull (2002) Kull, H. J. 1991 Theory of the Rayleigh–Taylor instability. Phys. Rep. 206(5), 197–325
- Kundu & Cohen (2001) Kundu, P. K. & Cohen, I. M. 2001 Fluids Mechanics - Second Edition Academic Press
- Jacqmin (1999) Jacqmin, D. 1999 Calculation of two-phase Navier–Stokes flows using Phase-Field modeling. J. Comp. Phys. 155, 96–127
- Layzer (1955) Layzer, D 1955 On the instability of superposed fluids in a gravitational field. Astrophys. J. 122, 1–12
- Lee, Rudnick & Brimhall Jr. (2001) Lee , C.-T., Rudnick, R. L. & Brimhall Jr., G. H. 2001 Deep lithospheric dynamics beneath the Sierra Nevada during the Mesozoic and Cenozoic as inferred from xenolith petrology. Geochem. Geophys. Geosys. 2, 2001GC000152
- Liu & Shen (2003) Liu, C. & Shen, J. 2003 A phase field model for the mixture of two incompressible fluids and its approximation by a Fourier-spectral method. Physica D 179, 211–228
- Menikoff et al. (1977) Menikoff, R., Mjolsness, R. C., Sharp, D. H. & Zemach, C. 1977 Unstable normal mode for Rayleigh–Taylor instability in viscous fluids. Phys. Fluids 20(12), 2000–2004
- Landau & Lifshitz (2000) L.D. Landau & E.M. Lifshitz 2000 Fluid Mechanics Volume 6 of Course of Theoretical Physics Second Edition, Revised Butterworth Heinemann
- Ramaprabhu & Andrews (2004) Ramabrabhu, P. & Andrews, M. J. 2004 Experimental investigation of Rayleigh–Taylor mixing at small Atwood numbers. J. Fluid Mech. 502, 233–271
- Rayleigh (1883) Lord Rayleigh 1883 Investigation of the caracter of the equilibrium of an incompressible heavy fluid of variable density. Proc. London Math. Soc. 14, 170.
- Schultz et al. (2006) Schultz, D.M., Kanak, K. M., Straka, J. M., Trapp, R. J., Gordon, B. A., Zrnić, D. S., Bryan, G. H., Durant, A. J., Garrett, T. J., Klein, P. K. & Lilly, D. K. 2006 The mysteries of Mammatus clouds: observations and formation mechanisms. J. Atmos. Sci. 10, 2409–2435
- Sethian (1999) Sethian, A.J. 1999 Level Set methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science. Cambridge University Press: Cambridge
- Sharp (1984) Sharp, D. H. 1984 An overview of Rayleigh–Taylor instability. Physica D 12, 3–18
- Smolianski, Haario & Luukka (2005) Smolianski, A., Haario, H. & Luukka, P. 2005 Vortex shedding behind a rising bubble and two-bubble coalescence: A numerical approach. Appl. Math. Model 29, 615–632
- Taylor (1950) Taylor, G. I. 1950 The instability of liquid surfaces when accelerated in a direction perpendicular to their planes I. Proc. R. Soc. A 201, 192–197
- Waddell, Nieserhaus & Jacobs (2001) Waddell, J. T., Niederhaus, C. E. & Jacobs, J. W. 2001 Experimental study of Rayleigh–Taylor instability: low Atwood number systems with single-mode initial perturbations. Phys. Fluids 13(5), 1263–1273
- Young et al. (2001) Young, Y. N., Tufo, H., Dubey, A. & Rosner, R. 2001 On the miscibile Rayleigh–Taylor instability: two and three dimensions. J. Fluid Mech. 447, 377–408
- Young & Ham (2006) Young, Y. N.& Ham, F. E. 2006 Surface tension in incompressible Rayleigh–Taylor mixing flow. J. Turbul. 71(7), 1–23.
- Yue et al. (2004) Yue, P., Feng, J.J., Liu, C. & Shen, J. 2004 A diffuse-interface method for simulating two-phase flows of complex fluids. J. Fluid Mech. 515, 293–317