An Analysis of Two-Dimensional Stratified Gravity Current Flow using Open FOAM

Direct numerical simulations (DNSs) of two-dimensional stratified gravity-current are simulated using OpenFOAM. Three different aspect ratio, h 0 /l 0 (where h 0 is the height of the dense fluid and l 0 is the length of the dense fluid) are simulated with stratification ranging from 0 (homogenous ambient) to 0.2 with a constant Reynolds number ( R e) of 4000. The stratification of the ambient air is determined by the density difference between the bottom and the top walls of the channel ( ρ b - ρ 0 , where ρ b is the density at the bottom of the domain and ρ 0 is the density at the top). The magnitude of the stratification ( S=ԑ b /ԑ ) can be determined by calculating the reduced density differences of the bottom fluid with the ambient fluid ( ԑ b = ( ρ b - ρ 0 )/ ρ 0 ) and the dense fluid with the ambient fluid ( ԑ = ( ρ c -ρ 0 )/ ρ 0 , where ρ c represents the density of the dense fluid). The configuration of the simulation is validated with a test case from Birman, Meiburg & Ungraish and the contour and front velocity (propagation speed) were in good agreement. The gravity current flow in the stratified ambient is analyzed qualitatively and compared with the gravity current in the homogenous ambient. Gravity current in homogenous ambient (S=0) and weak stratification (S=0.2) are supercritical flow where the flow is turbulent and Kelvin-Helmholtz (K-H) billow formed behind the gravity current head. The front location of the gravity is reduced as the stratification increase and denotes that the front velocity of the gravity current is reduced by the stratification.


Introduction
Gravity current or density current is a type of fluid flow that occurs mainly in the horizontal domain due to density difference of two fluids. Gravity current is observed in a natural phenomenon such as sandstorms, powder-snow avalanches, pyroclastic flow and haboob. Simpson [1] has reviewed a wide range of experimental studies on gravity current flow. Gravity current can be categorized into two types which is compositional gravity current (can be Boussinesq, where the density differences between two fluids are less than 5 % or non-Boussinesq) and particledriven gravity current (presence of sediments in the fluid). Gravity current flow has been widely studied with a variety of lock configurations such as lock shape [2] (elliptical release [3], [4] and circular release [5]), lock depth and initial aspect ratio [6]. Gravity current flow is determined by the Froude number (Fr) which is a dimensionless parameter defined as the ratio of inertial forces relative to the gravitational forces (Fr = u f / √h max , where u f is the front velocity and the h max is the maximum height of the current). When the Froude number is less than 1/π (Fr < 1/π), the gravity current resides in the subcritical flow regime. The flow is slow, and the formation of the gravitational waves which propagates faster than the front of the gravity current occurs. As the Froude number increases and becomes greater than 1/π (Fr > 1/π), the gravity current flow reaches the supercritical regime, where the flow becomes turbulent. The critical flow regime occurs in between the subcritical and supercritical regime, when the Froude number is 1/π (Fr = 1/π).
The propagation of the gravity current into a stratified ambient is the concern of this study. An experimental and numerical study has been conducted by Maxworthy et al. [7] to study the relationship between the Froude number (Fr), density ratio (R) and the effects of the internal waves on the dynamics of the gravity current flow. He found that the propagation of the wave is faster than the gravity current in the subcritical regime. Ungarish & Huppert [8] conducted a similar study using shallow-water approximation (SW) and the result shows excellent agreement with the experimental results presented by Maxworthy et al. [7]. He stated that the stratification will reduced the propagation speed of the gravity current compared to a homogenous ambient where the density is constant in wall normal direction. Ungarish (2005) [9] extended the stratified gravity current study with the 'Dam-break' release characteristic. In that study, several extensions of the stratified gravity current has been studied, such as shallow-water (SW) equation extended to non-linear stratification, gravity current released from an elliptical reservoir and axisymmetric configuration. Ungarish (2006) [10] developed the theory for the propagation of a steady gravity current speed into a linearly stratified ambient by generalizing the Benjamin's (1968) [11] theory on the steady state propagation of the gravity current into a homogeneous ambient. Birman et al. [12] then examined the theory on the dependency of the front velocity on the stratification (S) with varies depth ratio (ɸ=h 0 /L y ) by manipulating the height of the channel (L y ). It was found to work well for the weak stratification (S ≤ 0.5). Ungarish (2008) [13] also analyzed the energy-balances of an axisymmetric gravity current in the homogenous ambient and linearly stratified ambient with a two-dimensional geometry.
The purpose of the study is to investigate the effects of the initial aspect ratio and the stratification (where S=0, homogenous ambient and S=0.2, weak stratified ambient) to the dynamics of the gravity current flow using OpenFOAM. Two-dimensional stratified gravity current has been simulated with varies aspect ratio (h 0 /l 0 ) and stratification (S). The analysis of contour will be used to compare the flow structure of weak stratified (S = 0.2) gravity current with non-stratified (S = 0) gravity current.

Computational Setup
The two-dimensional stratified gravity currents were simulated using OpenFOAM with the solver twoLiquidMixingFoam. The non-dimensional system of equation employed in the study takes the form where ρ, u, and p are the density of the fluid, velocity vector, and pressure respectively. e g is the unit vector for the direction of the gravity. The variables with an asterisk denote dimensional variables. The dimensionless term is normalised by the velocity scale (U * ), characteristic length (L * ) and time scale (T * ). The dimensionless density (ρ) is defined as Sc is the Schmidt number (where ν * is the kinematic viscosity and κ * is the molecular diffusivity of the current). The value of Sc=1, is used in the study. Generally, the saline (salt water) has Sc=700 but it is found that when Sc is in the order of 1, there is a weak scaling with the dynamics of the gravity current [6]. The velocity scale (U * ), characteristic length (L * ) and time scale (T * ) is defined as where g * is the gravitational acceleration, * * * is the initial volume per spanwise unit of the dense fluids (ρ c ). The magnitude of the stratification (S) can be determined by calculating the reduced density differences of the bottom fluid with the ambient fluid (ԑ b ) and the dense fluid with the ambient fluid where g ' is the reduced gravity and g is the gravitational acceleration. The reduced density differences of the bottom fluid (ρ b ) with the ambient fluid and the dense fluid with the ambient fluid (ρ 0 ) is defined as where ρ b is the density at the bottom, ρ 0 is the density at the top and ρ c is the density of the dense fluid. The equation then takes the form of where ρ a is the density respect to the height, y and L y is the height of domain in wall-normal direction. The second-order backward implicit time scheme was used for integration of the time and the divergence vanLeer scheme used to discretise the convection terms in the transport equation (3). The twoLiquidMixingFoam solver has been validated with a test case from Dai [6] who used pseudo-spectral code to conduct the simulation. The front velocity and energy budgets result from OpenFOAM have good agreement with the case slope angle, θ of 9° and depth ratio, ϕ = h 0 /L y of 0.16 where L y is the wall-normal length of the domain. The stratified ambient is added to the domain by using the funkySetFields utility. Figure 1 shows the sketch of the computational domain with the streamwise length L x =20, the wall-normal length L y =4. No-slip boundary condition is employed for top and bottom (streamwise) of the walls. However, the vertical side walls (wall-normal) allow for the slip condition. The dense fluid with a density of ρ c inside a rectangular channel containing a linearly stratified ambient with an aspect ratio of h 0 /l 0 (where h 0 is the height of the dense fluid and l 0 is the length of the dense fluid). Three aspect ratios of the dense fluid, h 0 /l 0 = 0.25, 1, and 4 with the stratification (S) of 0 and 0.2 is simulated in this study. A large range of aspect ratio is chosen to simulate the gravity current flow in different flow regines. The largest aspect ratio of 4 is 8 times larger than the smallest aspect ratio 0.25. This allow the gravity current flow can be investigated ranging from the subcritical regime to the supercritical regime. Weak stratification strength (S=0.2) is used for comparison with homogeneous ambient (S=0) in the study to determine the accuracy OpenFOAM before proceed to strong stratification (S>0.5). The Reynolds number (Re) of 4000 is used for all simulations. The size in the computational domain is equally spaced and has a streamwise length of Δx ≈ 0.0078 and wall-normal length of Δy ≈ 0.0156. This grid resolution is sufficient to resolve all the turbulent length scales and is smaller than requirement of Δx ≈ (ReSc) -1/2 .

Results and Discussions
The configuration of the simulation is validated with a test case from Birman et al. [12] and is discussed in the first section. The qualitative analysis of the contour of the simulation study is discussed in Section 3.2.

Validation of the configuration of the Simulation
The configuration of the simulation is validated with a test case from Birman et al. [12] who investigated the prediction of the front velocity by theoretical analysis by Ungarish [10]. The paper is chosen for validation is to determine the accuracy of OpenFOAM with the results of the in house developed code by Birman et al. [12]. The countour plots and the front velocity results were agreeable with current simulation for the case, S = 0.2 and depth ratio, ɸ=h 0 /L y . Birman et al.   Figure 2 shows the contours of density at different time intervals for the validation case. Figure 2(a) is the initial condition for the validation study. Figure 2(b) is the contour of the gravity current flow at t = 5. The height of the dense fluid is reduced from 1 to approximately 0.5 as the dense fluid slumps due to gravity. The head of the gravity current approaches x ≈ 5, and three coherent vortices emerge behind the head of the current. At t = 10, the front location of the gravity current head is at x ≈ 7.8 and the height of the dense fluid is reduced to approximately 0.25. By comparing the figure 2(c) and figure 3(c), front location of the gravity current in figure 2(c), is slightly different with the front location in figure  3(c) where x ≈ 8.3, but the internal structure behind the head of the current looks about the same. It is worth to note that the waves of the stratification layer in figure 2(b) and 2(c) is the same with the wave of the stratification layer in figure 3(b) and 3(c). Figure 4 shows the comparison of the propagation velocity, U H , and it is clear that good agreement is obtained compared to Birman et al. [12]. The velocity of the current remains approximately constant after the initial acceleration phase. The propagation velocity in the current simulation is slightly slower by approximately 1.7% compared to Birman et al. [12]. This may be due to the different numerical method used by Birman et al. [12].  6 cases were simulated for the aspect ratio, h 0 /l 0 = 0.25, 1 and 4 with stratification (S) of 0 (homogeneous ambient) and 0.2. The gravity current flow regime in the current study is specified in Table 1. The concern of this study is the qualitative analysis of the contours to study the evolution of the flow regime of the gravity current. The numerical solution to determine the front location and the front velocity of the gravity current flow will not be described in this paper. Table 1 shows the flow regimes of the gravity current according to different aspect ratio (AR = h 0 /l 0 ) of the dense fluid and stratification (S). The flow regime can be determined by the Froude number (Fr) and the explanations were mentioned in the previous sections. The Froude number (Fr) can be calculated using the front velocity of the current, u f and the maximum height of the current, h max where it can be obtained directly from the simulation. According to Table 1, the flow regime remains supercritical flow for the aspect ratio, h 0 /l 0 = 0.25 1, and 4 with the stratification (S) increases from 0 (homogeneous ambient) to 0.2. Despite there being no changes in the flow regime for the gravity current flow, the internal structures of the gravity current flow are different from one another. The qualitative analysis of the contour is compared according to the aspect ratio in the sub-sections below.

The Aspect Ratio (h 0 /l 0 ) of 0.25
In this section, the contour of the gravity current flow with initial aspect ratio, h 0 /l 0 = 0.25 is compared with the stratification S = 0 (homogeneous ambient) and 0.2. The contour line of the stratification (solid black lines) have an interval of approximately 0.063 to clearly visualise the flow disturbance in the stratified region. This is to generate equivalent amount of the stratification layer showed in figure 5 in [12]. The intensity of the contour increasing from white at the top of the domain to light grey at the bottom of the domain due to the stratification. The time between each frames is t = 5.  The flow regime of the gravity current in figures 5 and 6 is supercritical flow. By comparing the contour of the gravity current flow, the size of the gravity current head and the internal structure (vortices) in the gravity current head is reducing as the stratification increases. The propagation of the waves is formed in figure 6. The motion of the wave is weak, mainly caused by the height of the gravity current head. The front location of the gravity current is decreasing as the stratification increases. This statement can be validated by observing the front location of the gravity current in figure 5(e) to figure 6(e). The front location of the gravity current flow in homogeneous ambient is approaching x ≈ 9.5 and reduced to x ≈ 8.8 as the stratification increased from 0 to 0.2 in figure 6(e).

The Aspect Ratio (h 0 /l 0 ) of 1
The contour of the gravity current flow with initial aspect ratio, h 0 /l 0 = 1 is compared with the stratification S = 0 (homogeneous ambient) and 0.2. The intensity of the colour at the bottom of the domain is increasing from light-grey to grey as the stratification is increased. The time between each frames is t = 5.

The Aspect Ratio (h 0 /l 0 ) of 4
The contour of the gravity current flow with initial aspect ratio, h 0 /l 0 = 4 is compared with the stratification S = 0 (homogeneous ambient) and 0.2. The intensity of the colour at the bottom of the domain is increasing from light-grey to grey as the stratification is increased. The time between each frames is t = 5. The flow regime for both cases is supercritical flow. Kelvin-Helmholtz (K-H) able to see clearly at the back of the head in both cases. The motion of the waves is strong and is presented in figure  10. The sinusoidal wave is formed between x ≈ 0 and x ≈ 8 at y ≈ 3 with maximum amplitude of 2.9 and minimum amplitude of 2.5 in figures 10(d) and 10(e). A complete sinusoidal wave cycle is difficult to be observed in other simulation cases. The study about the internal waves in the stratified ambient can be found in [7], [14]. The front location of the gravity current flow in the homogeneous ambient is approximately 10.5 and reduced to x ≈ 10 as the stratification (S) increases from 0 to 0.2 in figure 10(e).

Conclusion
Direct numerical simulations (DNSs) of a two-dimensional stratified gravity-current is simulated with varying aspect ratios and the stratification, S = 0 (homogeneous ambient) and 0.2 using Open-FOAM. A test case from Birman et al. [12] is simulated to validate the configuration of the simulation. The contour of the gravity current flow and the propagation velocity plot shows a good agreement with current simulations. The propagation velocity in the current simulation is slightly slower by approximately 1.7% compared to Birman et al. [12]. All the gravity current flow in the study is supercritical flow where the flow is turbulent with Kelvin-Helmholtz (K-H) behind the gravity current head. According to Ungarish (2002) [8], the propagation speed of the gravity current is reduced by the stratification compared with the motion in homogenous ambient where the density is remained constant at y = H. This phenomenon can be obtained in the current study where the front location of the gravity is reduced as the stratification is increased. Despite the analysis of the front velocity is not included in current study, this can be supported by observing the front location of the gravity current in current study. The front location of the gravity current, x is reduced from approximately 9.5 to 8.8, 10.5 to 9.8 and 10.5 to 10 at t = 20 in figure 5(e) to figure 10(e) respectively as the stratification (S) is increased from 0 to 0.2. A complete cycle of sinusoidal waves can be seen clearly between x ≈ 0 and x ≈ 8 at y ≈ 3 with maximum amplitude of 2.9 and minimum amplitude of 2.5 in figures 10(d) and 10(e). It is difficult to observe in other simulation cases. The structure of the gravity current head is smaller and vortices behind the gravity current head is reduced as the stratification increases. It is worth to point out that the aspect ratio, h 0 /l 0 of the dense fluid will affect the motion of the gravity current flow and affect the motion of the wave. There has been limited work on two-dimensional simula-tions of stratified gravity current where the aspect ratio of the lock is altered.