Numerical modeling of Rayleigh-Taylor instabilty vith interFoam

The Rayleigh-Taylor instability is an instability of the interface separating two immiscible fluids ; the densest of which rests on the lightest. Following a disturbance of the interface, the dense fluid enters the light fluid.

The interface growths with a characteristic “mushroom” shape. For example, Rayleigh-Taylor instability manifests itself in several areas such as astrophysics with the Crab Nebula, thermonuclear fusion with laser implosion of deuterium-tritium targets, or meteorology with cloud formation.

References :

[1] : D.H. SHARP, An overview of Rayleigh-Taylor Instability. 0167-2789/84/$03.00 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

[2] : Stéphane Popinet. Stabilite et formation de jets dans les bulles cavitantes : developpement d’une methode de chaine de marqueurs adaptee au traitement numerique des equations de navier-stokes avec surfaces libres. http://gfs.sourceforge.net/examples/examples/node7.html

[3] : Rusche, H. Computational fluid dynamics of dispersed two-phase flows at high phase fractions Imperial College of Science, Technology & Medicine, Department of Mechanical Engineering, 2002.

The linear theory describes from Navier-Stokes equations the space and time evolution of the initial interface perturbation. The case is assumed two dimensional and the parameters are shown in the figure below (ref[1]).

The initial perturbation is defined by the following expression :

\(y_{s} = \eta(t) \cos(kx)\)         (1)

It’s possible to show using various arguments from Bernoulli equations or potential theory that the perturbation evolution is described by the following differential equation:

\(\frac{d^2\eta}{dt^2}(t)= \alpha(k)^2\eta(t)\)         (2)

With:

\( \alpha(k)^2\ = g(\frac{\rho_{h}-\rho_{l}}{\rho_{h}+\rho_{l}})k-(\frac{\sigma}{\rho_{h}+\rho_{l}})k^3\)         (3)

\(\sigma\) being the surface tension (N/m) and\(\rho_{h},\rho_{l}\) the heavier and lighter density.

The differential equation (2) has for solution :

\(\eta(t)=\eta(0)cosh(\alpha t)\)         (4)

The calculation assumptions are those conventionally used (ref [2]) for modeling Rayleigh-Taylor instability.

• The domain size is rectangular 1 m(width)x 4 m(height).
• The mesh is structured (generated with blockMesh) with a cells number of 128×512.
• Lighter fluid properties: density 1.255 kg/m^3 ; dynamic viscosity 3.13 10^-3 Pa.s.
• Heavier fluid properties: density 0.1694  kg/m^3 ; dynamic viscosity 3.13 10^-3 Pa.s.
• Surface tension coefficient : 0.01 N/m.
• Initial perturbation shape (t = 0 sec) : cosinus with 0.05 m of amplitude and 1m of wave lenght.
• Laminar regime.
• Lateral patchs of slip wall type. Upper boundary of “atmosphere” type.

Setting up the case requires the initialization of the interface between the heavy fluid and the light fluid with a cosine-like shape. There are several possibilities to perform this operation:

• Use swak4Foam: This requires compiling the utility which is not necessarily easy depending on the version of OpenFOAM® you use.
• Use setFields. setFields is an OpenFOAM® utility for initializing a given field, including the possibility to use a stl format surface (generated for example with Salome).
• Use codeStream. codeStream offers the possibility of coding directly in the dictionaries of boundary conditions its instructions.

In the context of this post we will limit ourselves to the use of the second method. The use of the setFields command requires the completion of a dictionary located in the /system directory and detailed below:

The following figures show initial interface after setFields :

The discretization schemes are second-order except for the temporal derivative scheme (1st implicit order Euler). The term div(phi, alpha) corresponds to the convection term of the transport equation of the volume fraction. The term div(phirb, alpha) corresponds to the second convection term called numerical compression term of interface. This term [3] makes it possible to artificially compress the interface. interFoam solver uses the following numerical equation for the volume fraction:

\(\frac{\partial\alpha}{\partial t} + \nabla (\alpha \overrightarrow{U}) + \nabla (\alpha (1 – \alpha) \overrightarrow{U_r}) \)         (5)

The compression effect can be change by adjusting the cAlpha value (which influence the value of the relative speed \(\overrightarrow{U_r}\)) in fvSolution.

• A zero value do not compress the interface \(\overrightarrow{U_r} = \overrightarrow{0}\)(this will lead to a diffusive interface).
• A value of 1 corresponds to a conservative value. In most cases this value is recomended.
• A value greater than one lead to an enhancement of the compression.

It should be noted that this term compression tends to deform the interface.

The pressure / velocity coupling of interFoam is handled through the PISO loop (ie a PIMPLE loop with only one external corrector (nOuterCorrector = 1).

The following animated picture show the time evolution of the phase fraction and velocity fields calculated by interFoam. The results show that the perturbation develops downwards generating the “mushroom” characteristic shape.

The comparison of interFoam results with linear theory gives reasonable agreement.