Origin of lobe and cleft at the gravity current front

TL;DR

This paper investigates the formation of lobes and clefts at gravity current fronts, demonstrating that Rayleigh-Taylor instability explains their origin and confirming theoretical predictions with numerical simulations.

Contribution

It provides a detailed analysis linking RTI to the formation of lobes and clefts, validating semi-infinite RTI model predictions through 3D simulations.

Findings

RTI determines the formation of lobes and clefts.

Spanwise wavenumber depends on Prandtl number.

Wavenumber scales with Grashof number as 1/3.

Abstract

In this paper the Rayleigh-Taylor instability (RTI) along the density interfaces of gravity-current fronts is analyzed. Both the location and the spanwise wavenumber of the most unstable mode determined by the local dispersion relation agree with those of the strongest perturbation obtained from numerical simulations, suggesting that the original formation mechanism of lobes and clefts at the current front is RTI. Furthermore, the predictions of the semi-infinite RTI model, i.e. the original dominating spanwise wavenumber of the Boussinesq current substantially depends on the Prandtl number and has a 1/3 scaling law with the Grashof number, are confirmed by the three-dimensional numerical simulations.

Tables

Table 1. Table 1: The domain size ( L x , L z ) subscript 𝐿 𝑥 subscript 𝐿 𝑧 (L_{x},L_{z}) used in the simulations

	${10}^4$	${10}^5$	${10}^6$	${10}^7$	${10}^8$
$0.2$	(48, 83.1)	(48, 42.4)	(48, 19.6)	(48, 9.1)	(48, 4.2)
$1$	(48, 43.4)	(48, 19.9)	(48, 10.0)	(48, 4.4)	(48, 2.0)
$5$	(48, 28.5)	(48, 12.7)	(48, 5.9)	(48, 2.8)	(24, 1.3)
$650$	(48, 21.0)	(48, 9.4)	(48, 4.3)	(24, 2.0)	(12, 0.9)

Equations

$${\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\nabla p+\frac{1}{\sqrt{Gr}}{\nabla }^2\mathbf{u}+T{\mathbf{e}}_{\mathbf{y}},$$

$${\partial }_{t}T+(\mathbf{u}\cdot \nabla )T=\frac{1}{\sqrt{Gr}Pr}{\nabla }^{2}T,$$

$$\nabla \cdot \mathbf{u}=0,$$

$$(\gamma +\frac{1}{\sqrt{Gr}Pr}D)T_1=-v_1\frac{d{T}_0}{d\bar{y}},$$

$$\displaystyle(\gamma+\frac{1}{\sqrt{Gr}Pr}D)(\gamma+\frac{1}{\sqrt{Gr}}D)Dv_{1}=-\frac{Gr_{e}}{Gr}\cos(\theta)k^{2}v_{1}\frac{dT_{0}}{d\bar{y}},\ \ \$$

$$D=k^2-\frac{d^2}{d{\bar{y}}^2},$$

$$T_1=\left\{\begin{array}{c}{A}_{11}{e}^{{\alpha }_1\bar{y}}+A_{12}{e}^{-{\alpha }_1\bar{y}},-\bar{h}\le \bar{y}\le 0\\ A_{21}{e}^{-{\alpha }_1\bar{y}}+A_{22}{e}^{{\alpha }_1\bar{y}},0<\bar{y}\end{array}\right\},$$

$${\bar{v}}_1=\frac{1}{\sqrt{Gr}}\left\{\begin{array}{c}{B}_{11}{e}^{k\bar{y}}+B_{12}{e}^{-k\bar{y}}+C_{11}{e}^{{\alpha }_2\bar{y}}+C_{12}{e}^{-{\alpha }_2\bar{y}}\\ +A_{11}^{{}^{\prime }}{e}^{{\alpha }_1\bar{y}}+A_{12}^{{}^{\prime }}{e}^{-{\alpha }_1\bar{y}},-\bar{h}\le \bar{y}\le 0\\ B_{21}{e}^{-k\bar{y}}+B_{22}{e}^{k\bar{y}}+C_{21}{e}^{-{\alpha }_2\bar{y}}+C_{22}{e}^{{\alpha }_2\bar{y}}\\ +A_{21}^{{}^{\prime }}{e}^{-{\alpha }_1\bar{y}}+A_{22}^{{}^{\prime }}{e}^{{\alpha }_1\bar{y}},0<\bar{y}\end{array}\right\},$$

$$\gamma =\sqrt{\frac{k}{2}+\frac{k^4}{Gr}}-\frac{1}{\sqrt{Gr}}(1+\frac{1}{Pr})k^2.$$

$$k_{max}=(\frac{Gr}{2}{)}^{\frac{1}{3}}(\frac{1}{[2(1+\frac{1}{Pr})+\sqrt{4(1+\frac{1}{Pr}{)}^2-3}{]}^2-1}{)}^{\frac{1}{3}},$$

Full text

I Introduction

Lobes and clefts are the primary characteristics formed at the fronts of gravity currents Simpson et al. (1982, 1997); Huppert et al. (2006); Meiburg et al. (2010); Constantinescu et al. (2014); Hughes et al. (2016), dominating the local mass and momentum transportation Cantero et al. (2007); Espath et al. (2007); Jackson et al. (2007). Experimental study shown that these three-dimensional structures originated from the instability produced by the overrunning of light fluid by the dense current Simpson et al. (1972). However, it was argued recently that the overturning of ambient fluid or the density stratification might not be necessary and there might be some other Reynolds number-dependent formation mechanisms for lobes and clefts McElwaine et al. (2004). Linear stability analyses were carried out based on the flow fields of the current fronts obtained numerically, illustrating that the fronts were unstable to three-dimensional modes, and the mode growth rate was a weak function of the Schmidt number $Sc$ unless $Sc$ was very small Härtel et al. (2000). Three-dimensional numerical simulations were also carried out, and it was concluded that the lobe and cleft patterns were independent of $Sc$ by comparing the mean spanwise wavelength of lobes obtained numerically and the wavelength of the most unstable mode Bonometti et al. (2008). However, it should be noted that the mean spanwise wavelengths were estimated by counting the number of lobes observed in experiments and simulations Härtel et al. (2000); Bonometti et al. (2008), where the nonlinear cleft merging and lobe growing have occurred Simpson et al. (1972); Cantero et al. (2007); Espath et al. (2007), and hence the estimated wavelengths don’t correspond strictly to the linearly unstable mode.

Since the stratification is influenced by the diffusion property (strong diffusion corresponds to weak stratification), which is represented by the Schmidt number for mass diffusion or the Prandtl number (Pr) for thermal diffusion, a stratification related phenomenon should depend on $Sc$ or $Pr$. In fact, the discrepancy of the spanwise wavenumber between the linear result for $Sc=1$ and the experimental data for $Sc$ of $O(10^{3})$ Härtel et al. (2000) has implied $Sc$ effect. In addition, there are plenty of vortical structures in the viscous flow field of the current front Simpson et al. (1972); Härtel et al. (2000); Peng et al. (2010), and the previous linear stability analyses required detailed information of both the velocity field and the fluid’s properties Härtel et al. (2000). However, does the velocity field make much contribution to determining the spanwise wavenumber of the most unstable mode? It is shown in this paper by stability analyses and numerical simulations that the original spanwise wavenumber of lobes and clefts is intrinsically determined by the Rayleigh-Taylor instability, which depends on the fluid properties but not the velocity field, and its relations with the Grashof number and the Prantl number are obtained explicitly.

II Simulations and stability analyses

We consider gravity currents propagating along horizontal adiabatic bottom walls, which are governed by the incompressible Navier-Stokes equations. The density difference is caused by the temperature difference between the current fluid and the ambient fluid $\Delta\Theta$, which is assumed to be small enough to use the Boussinesq approximation. The Grashof number is defined as Gr$=\frac{g\beta\Delta\Theta h^{3}}{\nu^{2}}=\frac{2gAh^{3}}{\nu^{2}}$, where $g$, $\beta$, $\nu$, $h$ and $A$ are the gravitational acceleration, the thermal expansion coefficient, the kinematic viscosity, the characteristic length scale, and the Atwood number, respectively. The diffusion ratio is described by the Prandtl number, which is equivalent to the Schmidt number and is defined as Pr$=\nu/\kappa$, where $\kappa$ is the thermal diffusivity. By using the characteristic temperature $\Delta\Theta$, velocity $\nu\sqrt{Gr}/h$ and length $h$, the dimensionless governing equations are obtained as

[TABLE]

where $\mathbf{e_{y}}$ is the unit vector contrary to the gravity direction. These equations are solved with a pseudo-spectral method Chevalier et al. (2007) to study the lock-exchange flows in a cuboid box of $(L_{x},L_{y},L_{z})$ with periodic boundary conditions in the streamwise direction $x$ and the spanwise direction $z$, and no-slip conditions in the vertical or transverse direction $y$. The flow field is expanded in Fourier modes in the $x$ and $z$ directions and Chebyshev polynomials in $y$ direction. The fluids are initially at rest and the dense fluid fills two ends of the domain, where $0\leq x\leq L_{x}/4$ and $3L_{x}/4\leq x\leq L_{x}$.

After a short acceleration phase, the gravity current develops and reaches its slumping phase Benjamin et al. (1968); Rottman et al. (1983); Klemp et al. (1994); Shin et al. (2004); Scotti et al. (2008); Borden et al. (2013), where the front speed maintains constant. In this paper, we will focus on the three-dimensional instability at the beginning of the slumping phase. In order to examine the original dominating wavenumber $k_{max}$, the spanwise wavenumber of the linear unstable mode with the largest growth rate, random temperature disturbances are added to the whole field near the beginning of the slumping phase with the form of $\theta^{\prime}(x,y,z)=\delta\Sigma_{n=1}^{n_{dist}}\cos[2\pi nz/L_{z}+\phi(n)]$, where the amplitude $\delta$ should be small in order to trigger the linearly unstable mode and is set as 0.00033, $n_{dist}$ is the disturbing mode number, and $\phi(n)$ is a random number in the range of $(0,2\pi]$. For each couple of $Pr$ and $Gr$, preliminary simulations were carried out first to obtain the approximate dominating spanwise wavelength of lobes and clefts, and then $L_{z}$ is chosen to be at least twenty times larger than the approximation. The method determining the dominating wavelength or wavenumber is introduced next.

In the present simulations, $L_{y}=2$, spectral modes $(n_{x},n_{y},n_{z})=(2048,129,1024)$, and ($L_{x},L_{z}$) used for different parameters are listed in Table 1. Both the domain width and the mode resolution are larger and higher than those used in the previous three-dimensional simulations (e.g. Bonometti et al. (2008)) respectively, and it was checked that the dominating wavenumbers obtained with half modes agreed with the present results to better than $3\%$. Since vortex rolls will be induced by Kelvin-Helmholtz instability at high Grs and Prs Härtel et al. (2000), shorter $L_{x}$ and then finer resolutions are used as shown in Table 1. The front velocities and the dominating wavenumbers obtained from the test simulations within domains of $3L_{x}/4$ were checked to be consistent with the present results, indicating that $L_{x}$ is long enough to study the spanwise instability at the beginning stage of the slumping phase. In addition, the slumping velocity shown in Fig.1(b) and the dominating spanwise wavenumbers shown in Fig. 4(a) for $Pr=1$ agree well with the previous data (Fig. 3 and Fig. 6 of Härtel et al. (2000)) respectively, confirming that the present computational configurations are applicable.

It is shown in Fig. 2(a) that the spanwise modulation of the current front appears soon after the introduction of random disturbances, and develops quickly to the lobe and cleft structures [Fig. 2(b)]. The temperature perturbation amplitude $T^{*}_{max}$, representing the density amplitude, is shown in Fig. 2(c). During the slumping phase, the foremost part of the front in a translational frame moving with the front velocity is almost stationary. The strongest perturbation indicated by the empty square is very close to the interface of $\bar{T}=0.5$ and between the stagnation point of the mean flow and the maximum curvature of the mean interface averaged in the spanwise direction, agreeing with the previous observation that the maximum amplitude of the most unstable mode prevails between the stagnation point and the nose Härtel et al. (2000). According to the previous linear stability analyses, the growth rate or amplification rate is a continuous function of the wavenumber and has a single maximum at the dominating wavenumber $k_{max}$ for the unstable modes (Fig. 5 of Härtel et al. (2000)). Therefore, at the linear evolution stage where the perturbation amplitudes are small, the energy spectra of temperature perturbations caused by random disturbances with the same initial amplitude $E_{T}$, which satisfies $\frac{1}{L_{z}}\int^{L_{z}}_{0}T^{*2}dz=\int E_{T}(k)dk$, should be a smooth curve and has only one peak at $k_{max}$. Consequently, the original dominating spanwise wavenumber can be determined strictly by finding the wavenumber with the maximum spectral density at the initial stage of the instability. As illustrated in Fig. 2(d), the dominating wavenumber $k_{max}=18.2$. $n_{dist}=60$ corresponds to a wavenumber of $37.7>k_{max}$, and hence the wavenumber range of the introduced disturbances is wide enough to cover the dominating wavenumber. It is also checked that the dominating wavenumbers obtained with $n_{dist}=60$ and 100 are the same. During the nonlinear evolution stage, some lobes grow larger than their neighbors and some small lobes may disappear due to the cleft merging, and hence the spectral curve becomes irregular and has some peaks at low wavenumber range as the curves with $t>7$ shown in Fig. 2(d).

In order to clarify the contribution of the velocity field to the lobe-and-cleft formation, a new reference frame is built as shown in Fig. 1(a), where $\bar{x}$, $\bar{y}$ and $\bar{z}$ are in the tangent, normal and spanwise directions of the interface of $\bar{T}=0.5$, respectively. $\bar{y}$ deviates from $y$-direction with an angle of $\theta$, and the origin of the $\bar{x}\bar{y}\bar{z}$ frame has a distance of $\bar{h}$ from the bottom wall in the $\bar{y}$ direction. Now we study the Rayleigh-Taylor instability (RTI) only based on the temperature field in the $\bar{y}-\bar{z}$ plane, where the acceleration is $\bar{g}=g\cos\theta$. Because of the diffusivity the density interface has a finite thickness, whose stabilising effect on the RTI can be estimated by an effective Atwood number $\frac{A}{1+AkL}$ Lelevier et al. (1955); Cherfils et al. (1996); Mikaelian et al. (1997) or an effective Grashof number Gr${}_{e}=\frac{Gr}{1+AkL}$ and the temperature gradient length $L=(\bar{T}/\frac{d\bar{T}}{d\bar{y}})_{\bar{T}=0.5}$.

Assuming that the normal mode is in the form of $f^{\prime}=f_{1}(\bar{y})exp[ik\bar{z}+\gamma t]$, where $f^{\prime}$ represents the disturbing $\bar{y}$-direction velocity $v^{\prime}$, temperature $T^{\prime}$, and pressure $p^{\prime}$ respectively, and substituting them into the linearized governing equations in the $\bar{y}-\bar{z}$ plane, we have

[TABLE]

where $T_{0}$ is the temperature field obtained in simulations when the current flow is still two dimensional. These equations can be solved with a method similar as used in Mikaelian et al. (1996). The eigenfunctions $T_{1}$ and $v_{1}$ can be expressed as

[TABLE]

[TABLE]

where $A_{ij}^{{}^{\prime}}=\frac{Gr_{e}\cos(\theta)k^{2}}{GrPr(Pr-1)\gamma^{2}}A_{ij}$ with $i$ and $j=1,2$, $\alpha_{1}=\sqrt{k^{2}+Pr\sqrt{Gr}\gamma}$, and $\alpha_{2}=\sqrt{k^{2}+\sqrt{Gr}\gamma}$.

It is known that the influence of the unstable RT mode is limited in the range about one wavelength from the interface Drazin (2002). Since the dominating spanwise wavelength is smaller than the height of the current front as shown in Figs. 2(a) and 2(b), the unstable modes decay exponentially in the dense fluid $(\bar{y}>0)$ and vanish at the bottom of the light fluid where $\bar{y}=-\bar{h}$. By applying the boundary conditions and the continue conditions for velocity, stress, and temperature at the interface, Eqs. (4)-(7) can be written in a form of $MV=0$, where $M$ is an $12\times 12$ matrix and $V$ is a vector composed by the twelve constants, $A_{ij}$, $B_{ij}$, and $C_{ij}$. A nontrivial solution requires that the determinant of $M$ equals zero, and hence the dispersion relation $D(Gr,Pr,\theta,L,k,\gamma)=0$ is obtained. For details of the solving method, we refer to Mikaelian et al. (1996).

The maximum growth rate $\gamma_{m}$ and the corresponding spanwise wavenumber $k_{m}$ are calculated according to the local dispersion relation for given positions ($\theta$) along the interface. It is shown in Fig. 3(a) that the maximum of $\gamma_{m}$, which is referred as $\gamma_{max}$, lies at $\theta=0.54$ for Pr=1 and Gr=$10^{6}$. Stability analyses were carried out based on the temperature fields obtained at different instants (t = 5 and t=6) respectively, and $\gamma_{m}$ and $k_{m}$ curves almost coincide with each other. As shown in Fig. 3(b), $k_{max}$ corresponding to the maximum of $\gamma_{m}$ is 21.3, which is close to the dominating spanwise wavenumber 18.2 obtained from the numerical simulations [Fig. 2(c) and 2(d)], indicating that the origin of the lobe and cleft is determined mainly by the Rayleigh-Taylor instability at the front lower interface. This result is surprising at its first glance because the present stability analysis is solely based on the density field, while the previous gravitational stability analyses Härtel et al. (2000) required detailed information of both the velocity field and the density field of the current front. The underlying rationality is based on two features of the mean flow. First, the mean flow around the position of the strongest perturbation in the translational frame is weak, because it is close to the stagnation point, whose mean velocity is zero. It is noted that the front velocity field without perturbations is two-dimensional and has no spanwise velocity. Second, as shown by the streamlines in Fig. 2(c), the mean flow around the position of the strongest perturbation is nearly parallel to the interface or in the $\bar{x}$ direction, and then has little contribution to the RTI in the $\bar{y}\bar{z}$ plane.

The acceleration and the interface thickness can change evidently the growth rate of RT unstable mode but not the dominating wavenumber, which is mainly determined by the momentum and the density diffusion effects Duff et al. (1962). The agreement between the simulation and the RTI model shown in Fig. 3(b) encourages us to consider a further simplification, where the exact density and velocity fields of the current front are not required. Following the previous studies Duff et al. (1962); Bellman et al. (1954); Renaud et al. (1997), the dispersion relation governed by Eqs.(4)-(7) can be simplified for the semi-infinite fluids with a planar and infinitely thin interface, and the approximate but analytical solution is solved as

[TABLE]

The most unstable mode corresponds to $d\gamma/dk=0$, and it is easy to obtain that

[TABLE]

indicating that the original dominating spanwise wavenumber should scale as Gr1/3. For Gr$=10^{6}$ and Pr=1, Eq. (9) suggests $k_{max}=20.1$, which is close to the simulation value 18.2. By analyzing the simulation data with the same methods as discussed in Fig. 2 and Fig. 3, the dominating wavenumbers are obtained and agree well with the theoretical predictions for different Grs and Prs as shown in Fig. 4(a), indicating clearly the 1/3 scaling law and confirming that RTI is the dominating onset mechanism for lobes and clefts. This mechanism explains why $k_{max}$ was found numerically not sensitive to the length of the computational domain $L_{x}$: Short domain or small $L_{x}$ may affect the flow topology and the front speed due to the periodic boundary condition, but $Pr$ and $Gr$ remain the same and hence the wavenumber determined by RTI does not change.

When Pr$\rightarrow\infty$, it is easy to obtain the asymptotic value of $k_{max}=Gr^{1/3}/2^{4/3}$ from Eq.(9). For Boussinesq fluids, such a Pr dependence is reasonable because small Prandtl number corresponds to strong density diffusion, which suppresses high wavenumber components and lower the dominant wavenumber. As shown in Fig. 4(b), $k_{max}$ decreases with the decrease of Pr and has a strong dependence on Pr for Pr$<50$ as predicted by the semi-infinite RTI model. In fact, for the cases with small $Pr$ and $Gr$, the strong density and momentum diffusion effects thicken the density interface and shorten the light fluid layer beneath the dense one at the current front, and hence undermine the prerequisite of RTI. According to the simulations for $Pr$=0.2, the current front becomes stable to spanwise perturbations when $Gr$ is as low as $10^{4}$.

The spanwise wavenumber of the most unstable mode obtained previously by analyzing a rectangular domain Härtel et al. (2000) are shown as a dotted curve in Fig. 4(a), agreeing well with the present simulations. When lobes are clearly recognized, the nonlinear growth of lobes and cleft merging have occurred, leading to notable growth of low wavenumber components as shown in Fig. 2(d). Consequently, the mean wavenumber is smaller than the original dominating value, especially when it is estimated by counting the lobe’s number. For example, the estimated spanwise wavenumber of Simpson’s experiments Simpson et al. (1972) is shown in Fig. 4(a) as the filled square Härtel et al. (2000), which is lower than the simulation data.

III Conclusions

Though the velocity field and the density field are coupled for the gravity currents, the present simulations and stability analyses show that the local Rayleigh-Taylor instability at the density interface can determine the position and the original spanwise wavenumber of the strongest perturbations, illustrating that the density stratification is the crucial factor dominating the onset of lobes and clefts. Without any detailed information of the flow field but only the fluids’ properties ($Gr$ and $Pr$), the semi-infinite RTI model predicts successfully the original dominating spanwise wavenumber of the current front and the 1/3 scaling law, indicating that RTI is the original formation mechanism for lobes and clefts. In addition, it is shown theoretically and numerically that the Prandtl number has substantial effect on the spanwise wavenumber selection.

The authors thank Paul Linden for enlightening discussions on the current flows. Simulations code SIMSON from KTH and the help from Philipp Schlatter, Luca Brandt and Dan Henningson are gratefully acknowledged. The simulations were performed on TianHe-1(A) and this work has been supported by the National Natural Science Foundation of China (Grants No. 91752203, No. 11490553, and No.11521091).