Capturing the baroclinic effect in non-Boussinesq gravity currents

Shengqi Zhang a , Zhenhua Xia b , ?

a State Key Laboratory for Turbulence and Complex Systems, College of Engineering, Peking University, Beijing 100871, China b Department of Engineering Mechanics, Zhejiang University, Hangzhou 310027, China

Keywords:Baroclinic effect Gravity current

ABSTRACT Direct numerical simulations of two-dimensional gravity currents with small and medium density variations are performed using different non-Boussinesq buoyancy approximations. Taking the full low-Machnumber approximation as the reference, the accuracy of several buoyancy terms are examined. It is found that all considered buoyancy terms performed well in the cases with small density variation. In the cases with medium density variation, the classical gravitational Boussinesq’s buoyancy term showed the lack of accuracy, and a simple correction did not make any improvement. In contrast, the recently introduced second-order buoyancy term showed a significantly higher accuracy. The present results and our previous derivations indicate that simple algebraic buoyancy approximations extended from the Boussinesq’s gravitational buoyancy are unlikely to achieve an accuracy beyond first order. Instead, it seems necessary to solve at least one extra Poisson equation for buoyancy terms to capture the higher-order baroclinic effect. An approximate analysis is also provided to show the leading term of the non-Boussinesq effect corresponding to gravity.

The baroclinic effect is an important mechanism of vorticity generation, and could be regarded as the essence of the buoyancy effect which converts potential energy to kinetic energy [1] . Such effect is due to the nonalignment of pressure gradient and density gradient, and it appears in the vorticity equation as the curl of the pressure term ? ×(?p/ρ), namely the baroclinic torque.Due to the complexity caused by the coupling of pressure and density in the pressure term ?p/ρ, it is widely preferred to substituteρin the denominator by the constant reference densityρ0≈ρ. Apparently, such simplification completely eliminates the curl of the pressure term in the vorticity equation and requires a buoyancy termBto compensate the loss of the baroclinic torque.Therefore the curl of the buoyancy term ? ×Bshould be equal to or approaching ? ×(?p/ρ)in order to recover the true evolution of vorticity [1] . The earliest and simplest model of the buoyancy effect (baroclinic effect) is the classical Boussinesq’s gravitational buoyancy termρ′g/ρ0[2] , withρ′=ρ?ρ0being the density fluctuation. Following the pioneering work of Boussinesq, various buoyancy terms are proposed, including those related to the centrifugal force [3–5] , convection term [6] , and Coriolis force [7] .The gravitational and centrifugal buoyancy terms are proved to be accurate whenρ→ρ0, and are widely used in the theoretical analysis and numerical simulations of various flow problems,especially buoyancy-driven flows [8–15] . The latter two buoyancy terms are claimed to be corrections to the previous ones, and have been adopted in recent works [16–20] . However, as pointed out by Zhang et al. [1] , these two corrections are invalid in low-Machnumber (LMN) flows, because the errors they bring in have the same order as the corresponding parts of the baroclinic torque. Instead, based on the LMN equations [21] , Zhang et al. [1] introduced two types of buoyancy terms constructed from the perturbation solutions of the baroclinic torque, which can be made arbitrarily accurate by increasing the order of the perturbation expansion. Although the theoretical derivations are complete and rigorous [1] ,the supporting numerical simulations were not sufficient for people to fully understand the consequence of invalid buoyancy terms and the accuracy of the newly proposed buoyancy terms.

Non-Boussinesq gravity currents, which are important in various safety and environmental problems [22,23] , are potential flow problems for the examination of different buoyancy terms. There are various experimental and numerical works investigating non-Boussinesq gravity currents [22,24–28] , and numerical simulations with LMN equations have obtained consistent results with experiments [26] , indicating that the LMN equations are appropriate for the simulation of non-Boussinesq gravity currents.

In the present work, direct numerical simulations of gravity currents are performed in order to examine the accuracy of the classical buoyancy term, a simple correction corresponding to the convection term from Lopez et al. [6] , and the simplest non-trivial buoyancy term introduced in Ref. [1] . The numerical results are compared to the reference ones using the LMN equations, and we have shown more details about the mechanism of how the errors of buoyancy terms turn into the errors of flow fields in gravity currents. We hope that the present work, which is a follow-up one of Ref. [1] , can provide a better chance for people to understand the consequence of invalid buoyancy terms and the importance of the higher-order expansion in capturing the baroclinic effect.

We consider the LMN variable-density fluid in a twodimensional rectangular domainV= [ 0,L?x] ×[ 0,L?y] , with no-slip horizontal walls and free-slip non-penetrative vertical walls. The fluid has uniform kinematic viscosityν?and density diffusion coefficientκ?, and is subject to a uniform body forcef?=f?(exsinβ?eycosβ). Initially, the heavy fluid with densityρ?Hoccupied the domain [0,L?0]×[0,H?0], and was surrounded by the light fluid with densityρ?L=γ ρ?H. TakingL?r=H?0 as the reference length scale,T?r= [H?0/((1 ?γ)f?)]1/ 2 as the reference time scale, andρ?r=ρ?Has the reference density scale, the LMN governing equations [21] ,boundary conditions and initial conditions are written as

here, the Reynolds numberRe=(L?r)2(T?r)?1 /ν?= 40 0 0 ; the Schmidt numberSc=ν?/κ?= 1 ; the tilt angle of the body forceβ= 9?; the non-dimensionalized sizes of the fluid domain areLx= 12.5 andLy= 6.25 , and the initial sizes of the heavy fluid domain areL0= 1.25 andH0= 1 . The present set-up is sketched in Fig. 1 . It should be noted that the difference between the LMN equations and fully compressible Navier-Stokes equations is that, in the LMN momentum equation the thermodynamic pressure is considered to be spatially homogeneous and the hydrodynamic pressureΠshould be specified by an second-order elliptic equation in order to make the momentum equation consistent with the continuity equation [1] . The physical parameters and the basic set-ups are the same as Dai and Huang [26] , except that their continuity equation was slightly inconsistent with their density diffusion term and their viscous term was simplified.

In order to examine the performance of buoyancy terms, the widely used simplified momentum equation is considered:

The simulations of LMN Eq. (1a) –(1c) are performed with a second-order finite-difference code on a homogeneousNx×Ny=10 0 0 ×10 0 0 grid, and the pressureΠis computed with the successive over-relaxation (SOR) method. The LMN equations with the simplified momentum equation and three buoyancy terms are solved with the same grid and code as described above, except that～Πin Eq. (4) and ?Π(1)in Eq. (7) are both computed with a highly efficient discrete cosine transform because they are governed by Poisson equations. Density ratiosγ∈ { 0.998,0.98,0.9,0.8,0.7,0.6 }are considered to study the performance of buoyancy terms at small and medium density variation. Eachγrequires four simulation cases, including the reference case corresponding to original LMN equations, and three cases corresponding to the buoyancy terms. All variables from the reference cases are marked with superscript ‘ ? ’, and others are marked with the same superscripts as the corresponding buoyancy terms.

A normalized density is defined asθ=(ρ?γ)/(1 ?γ), and the equivalent heights of gravity currents are defined as [26]

Figure 2 shows the equivalent heights of gravity currents atγ=0.998 , indicating that the present results are in good consistency with those from Dai and Huang [26] (their validation case). Therefore the solvers for LMN equations and simplified equations are reliable, and the grid resolution is sufficient. In addition, Fig. 2 indicates that at very small density variation, the classical Boussinesq’s buoyancy is a satisfactory approximation of the baroclinic effect.

Fig. 1. The sketch of the present set-up.

Fig. 2. Equivalent heights h of gravity currents at γ = 0 . 998 . Here, h? refers to the results from the LMN equations, h B denotes the results using the classical Boussinesq’s buoyancy term. The results from Dai and Huang [26] are also shown for comparison.

Fig. 3. Relative errors of velocity fields at ( a ) γ = 0 . 998 and ( b ) γ = 0 . 8 .

Figure 3 a shows the time evolutions of the relative errors ofuatγ= 0.998 , defined with the vectorL?2 norm. It can be seen that the relative error ofuBcorresponding to the Boussinesq’s buoyancyρ′f/ρ0is small but increasing with time. The correction term ?ρ′u·?u/ρ0introduced by Lopez et al. [6] does not make any difference at the beginning because the convection term is very small initially. However, withtincreasing, the correction makes the error even larger than that of the original Boussinesq’s buoyancy term. This is because such correction term has an error ofO(?)[1] , where?is almost equal to the relative density differenceχ=(1 ?γ)/γduringt∈ [0,5] . Since the error ofBBis alreadyO(?), the simple correction apparently will not provide any significant improvement. Figure 3 a also shows that the error ofu(2)is much smaller, which can be explained by the fact that the error of ?B(2)isO(?2)convergent. At a largerγ= 0.8 , the relative errors of all buoyancy terms are larger, as shown in Fig. 3 b. Although the the error ofu(2)is about 104times larger than that atγ= 0.998 , it is still much smaller than the errors ofuBanduL, and remains below 4% withint∈ [0,3] . In contrast, the relative errors of bothuBanduLexceed 25% att= 3 , and the extra relative error caused by the correction term ?ρ′u·?u/ρ0is around 5% .

Fig. 4. Contours of the normalized density at γ = 0 . 8 and t = 3 . The solid isolines show an isovalue of 0.5: ( a ) θ? ; ( b ) θB ; ( c ) θL ; ( d ) θ(2) . The purple dashed lines in ( b,c,d )represent the isoline θ? = 0 . 5 in panel ( a ).

Fig. 5. The maximum front velocities U f, max of gravity currents as a function of the relative density difference χ = (1 ?γ) /γ. The experimental results of Dai [25] with Reynolds number ranging from Re ≈90 0 0 for γ = 0 . 998 to Re ≈260 0 0 for γ = 0 . 85 and numerical results from Dai and Huang [26] with Re = 40 0 0 are also shown for comparison.

Fig. 6. Relative errors of ( a ) velocity fields and ( b ) density fields with different ?= (1 ?γ) /γ at t = 3 .

Therefore, the order of convergence ofδutwith respect to?should be the same as ? ×(B?B?). Since the velocity appears in the convection term of density of Eq. (1c) , the errors of normalized density fields should have the same convergence rate as that of the velocity fields, which is verified by Fig. 6 b.

For numerical applications, the applicable range of buoyancy terms should be considered, and can be indicated from Fig. 6 . For example, if it is required that the relative errors of both velocity and density fields are below 1% , the Boussinesq’s buoyancy termBBand the modified buoyancy termBLwith the simple correction?ρ′u·?u/ρ0are applicable atγ≥ 0.998 , while the new buoyancy term ?B(2)is applicable atγ≥ 0.9 . A wider applicable range requires a type-II buoyancy term ?B(n)with higher order, which involves more Poisson equations to be solved. It should be mentioned that sinceρ0is chosen to be 1, aγwhich is smaller than 0.5 will cause?> 1 , at which the convergence of ?B(n)withn→ ∞has not been proved. However, takingρ0= 2γ/(1 +γ)will make?≤(1 ?γ)/(1 +γ)< 1 , ensuring that

The computational efficiencies of different treatments of the baroclinic torque is also worth mentioning. The present simulations are performed with MATLAB R2016a and CUDA 7.5 on an NVIDIA Quadro K40 0 0 graphics card. Computing the second-order type-II buoyancy term (one iteration) only requires 0.075 s with the highly efficient discrete cosine transform, while solving the pressure equation with SOR and optimized relaxation parameter requires more than 200 iterations, which cost over 2.2 s. Therefore, taking the final projection operation into account, computing the simplified LMN equations with the second-order type-II buoyancy term using the highly efficient discrete cosine transform only takes around 7% of the computational time in solving the original LMN equations with SOR pressure solver.

In conclusion, the present work used the simulations of twodimensional gravity currents to examine the accuracy of three different buoyancy terms. The results show that the classical Boussinesq’s buoyancy term is only reliable at very small density variation and a simple correction by Lopez et al. [6] will possibly worsen the performance, while a new buoyancy term proposed by Zhang et al. [1] , which requires solving an extra Poisson equation,is applicable at relatively larger density variations.

Declaration of Competing Interest

No potential conflict of interest was reported by the authors.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grants 11822208 , 92152101 , 11772297 , and 91852205 ).