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

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

[
7 July 2019
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).

Figure 1: Fluids configuration corresponding to a heavier fluid of density placed above a lighter one of density .

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

Figure 2: The square growth-rate (see Eq. (2.0)) for three different values of corresponding to three different values of the critical wave number : (solid circle), (solid triangle) and (solid rhombus). The dashed line is the linear-theory prediction expressed by the relation (2.0).
Figure 3: The square growth-rate for (solid circles), (solid triangles) and (solid rhombus), all smaller than , for six different values of ranging from to . The dashed line corresponds to the linear-theory prediction.

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 .

Figure 4: Time behavior of for (in Fig. 2 corresponding to the solid triangle) and for four values of . (a) , (b) , (c) and (d) . The numerical results (symbols) are compared with the corresponding best fit expressions (see the text for details). In the insets the interface perturbation, , is plotted at revealing a very accurate linear analysis prediction.

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.

Figure 5: Behavior of the dimensionless perturbation growth-rate, , for and corresponding to . Dotted lines correspond to upper and lower bounds for the growth-rate (see Eqs. (A 0) and (A 0)). The arrow selects a value of the viscosity for which the time evolution of is reported in the inset. The continuous line is the best fit slope (see text).

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.

Figure 6: Time behavior of the perturbation maximum, , for and . The critical wave number is . Numerical results (symbols) are compared with the prediction from linear theory (see Eq. (2.0)).

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).

Figure 7: Time evolution of amplitude perturbation . The dots are our numerical results, the dashed line is the prediction by Eq. (B 0).
Figure 8: Two-color snapshots of the phase field. Black (white) corresponds to (). Frames are equally spaced in time in the interval (see also Fig. 7).

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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
76591
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description