The Rayleigh-Taylor instability

Tags: category:benchmark feature:2d feature:cartesian feature:community-benchmark feature:analytical-solution

(sec:benchmarks:rayleigh_taylor_instability)=`

The Rayleigh-Taylor instability#

This section was contributed by Cedric Thieulot.

This benchmark is carried out in Deubelbeiss and Kaus [2008], Gerya [2010], Thieulot [2011]) and is based on the analytical solution by Ramberg [Ramberg, 1968], which consists of a gravitationally unstable two-layer system. Free slip are imposed on the sides while no-slip boundary conditions are imposed on the top and the bottom of the box. Fluid 1 \((\rho_1,\eta_1)\) of thickness \(h_1\) overlays fluid 2 \((\rho_2,\eta_2)\) of thickness \(h_2\) (with \(h_1+h_2=L_y\)). An initial sinusoidal disturbance of the interface between these layers is introduced and is characterized by an amplitude \(\Delta\) and a wavelength \(\lambda=L_x/2\) as shown in Fig. 211.

../../../../../_images/setup.png

Fig. 211 Setup of the Rayleigh-Taylor instability benchmark (taken from Thieulot [2011])#

Under this condition, the velocity of the diapiric growth \(v_y\) is given by the relation

\[\frac{v_y}{\Delta} = - K \frac{\rho_1-\rho_2}{2 \eta_2} h_2 g \qquad \qquad \text{with} \qquad \qquad K=\frac{-d_{12}}{c_{11}j_{22}-d_{12}i_{21}}\]

where \(K\) is the dimensionless growth factor and

\[\begin{split}c_{11} &= \frac{\eta_1 2 \phi_1^2}{\eta_2(\cosh 2\phi_1 - 1 - 2\phi_1^2)} - \frac{2\phi_2^2}{\cosh 2\phi_2 - 1 - 2 \phi_2^2}\\ d_{12} &= \frac{\eta_1(\sinh 2\phi_1 -2\phi_1)}{\eta_2(\cosh 2\phi_1 -1 -2\phi_1^2)} + \frac{\sinh 2\phi_2 - 2\phi_2}{\cosh 2\phi_2 -1 -2\phi_2^2} \\ i_{21} &= \frac{\eta_1\phi_2 (\sinh 2 \phi_1 + 2 \phi_1)}{\eta_2(\cosh 2\phi_1 -1 -2\phi_1^2)} + \frac{\phi_2 (\sinh 2\phi_2 + 2\phi_2)}{\cosh 2\phi_2 -1 -2\phi_2^2} \\ j_{22} &= \frac{\eta_1 2 \phi_1^2 \phi_2}{\eta_2(\cosh 2\phi_1 -1-2\phi_1^2)} - \frac{2\phi_2^3}{ \cosh 2\phi_2 -1 -2\phi_2^2}\\ \phi_1&=\frac{2\pi h_1}{\lambda} \\ \phi_2&=\frac{2\pi h_2}{\lambda}\end{split}\]

We set \(L_x=L_y=512 \text{ km}\), \(h_1=h_2=256 \text{ km}\), \(|\boldsymbol{g}|=10\text{m}/\text{s}^2\), \(\Delta=3\text{ km}\), \(\rho_1=3300\text{ kg}/\text{m}^3\), \(\rho_2=3000\text{kg}/\text{m}^3\), \(\eta_1=1e21\text{ Pa s}\). \(\eta_2\) is varied between \(10^{20}\) and \(10^{23}\) and 3 values of \(\lambda\) (64, 128, and 256km) are used. Adaptive mesh refinement based on density is used to capture the interface between the two fluids, as shown in Fig. 212 and Fig. 213. This translates as follows in the input file:

subsection Mesh refinement
  set Initial global refinement = 4
  set Initial adaptive refinement = 6
  set Strategy = density
  set Refinement fraction = 0.6
end

Description of benchmark files

../../../../../_images/grid.png

Fig. 212 Grid with initial global refinement 4 and adaptive refinement 6.#

../../../../../_images/grid2.png

Fig. 213 Density field with detail of the mesh.#

The maximum vertical velocity is plotted against \(\phi_1\) in Fig. 214 and is found to match analytical results.

../../../../../_images/plot.svg

Fig. 214 Maximum velocity for three values of the \phi_1 parameter.#