Numerical Simulation of Density Current Evolution in a Diverging Channel

Abstract

When a buoyant inflow of higher density enters a reservoir, it sinks below the ambient water and forms an underflow. Downstream of the plunge point, the flow becomes progressively diluted due to the fluid entrainment. This study seeks to explore the ability of 2D width-averaged unsteady Reynolds-averaged Navier-Stokes (RANS) simulation approach for resolving density currents in an inclined diverging channel. 2D width-averaged unsteady RANS equations closed by a buoyancy-modified 𝑘−𝜀 turbulence model are integrated in time with a second-order fractional step approach coupled with a direct implicit method and discretized in space on a staggered mesh using a second-order accurate finite volume approach incorporating a high-resolution semi-Lagrangian technique for the convective terms. A series of 2D width-averaged unsteady simulations is carried out for density currents. Comparisons with the experimental measurements and the other numerical simulations show that the predictions of velocity and density field are with reasonable accuracy.

1. Introduction

When a river enters the relatively quiescent water of a lake or reservoir, it meets water of slightly different temperature, salinity, or turbidity. Three configurations may occur. First, if the river is lighter than the surface water, it will form an overflow. Such a current has been observed in Lake Kootenay in Canada [1]. Second, if the river reaches a depth where its density equals that of the ambient lake water, the plume separates from the bottom and spreads horizontally along surfaces of constant density [2]. Third, when the river is heavier than the water of the lake, it sinks and flows as an inclined plume along its morphological channel. For instance, such a current has been reported by Huppert and Simpson [3] in the Wellington Reservoir in Australia. The analysis of density currents is of great relevance for water quality management in reservoirs as they carry suspended matters and dissolved solids across the lake. Consequently, density currents often determine the distribution of pollutant substances in lakes. Density currents have been observed in situ [4–9], and several laboratory experiments have been performed to study density currents [10–13]. Some studies have essentially dealt with inclined channels of constant width [14–17], and the others have been considered channels of variable widths [18, 19] or unbounded 3D plumes [20]. With the rapid development of numerical methods and advancements in computer technology, CFD has been widely used to study density currents. Unsteady RANS approach with the statistical turbulence model provides a viable alternative for engineering computations of gravity currents at practical Reynolds numbers, and several such models have been proposed in the literature. The 𝑘−𝜀 two-equation model has been widely adopted by researchers to simulate 2D conservative gravity currents with constant width channels [21–28]. Bournet et al. [29] simulated plunging gravity currents in an inclined channel of constant width and then in a diverging channel. The commercially available code PHOENICS was used for calculations by Bournet et al. [29]. It is based on a coupling between pressure and momentum (SIMPLEST algorithm) and an iterative solution procedure. Kassem et al. [30] simulated negatively buoyant flow in diverging channels by the computational fluid dynamic solver FLUENT which solves the three-dimensional Reynolds-averaged Navier–Stokes equations, along with the heat transport equation using an implicit finite volume method. Simulations of Bournet et al. [29] and Kassem et al. [30] were presented after the head of the current has passed and the under current reaches steady state. Plunging flow is intrinsically unsteady because of continuous entrainment of ambient water [21]. On the other hand, existing commercial codes have some limitations. For instance, the memory requirements and execution time are often excessive for solving hydraulic engineering problems. The boundary conditions to be used are rigidly specified in the software. The CFD software is also very expensive. It is therefore beneficial to develop and validate a new numerical model that can speedily solve the specific problems in hydraulic engineering.

In this paper, the model employing the Boussinesq form of the width-averaged unsteady RANS equations in conjunction with a buoyancy-extended 𝑘−𝜀 closure for the turbulence [32] is developed for unsteady simulation of density currents in an inclined diverging channel. We carry out numerical simulations of various continuous density currents in an inclined diverging channel, and then the numerical results are compared with the experimental measurements and the numerical simulation published by Johnson and Stefan [31] and Bournet et al. [29], respectively. We also present the 2D width-averaged numerical simulation of the density current in an inclined diverging channel at selected time instants.

2. Governing Equations

In this study, we focus on flows with relatively small density differences for which the usual Boussinesq approximation can be assumed to be valid—that is, all variations in density can be neglected except for the buoyancy term. Incorporating the Boussinesq hypothesis to relate the Reynolds stresses with the mean rate of strain tensor via an eddy viscosity, 2D width-averaged unsteady RANS equations for incompressible, stratified flow read as𝜕𝐵𝑢+𝜕𝑥𝜕𝐵𝑤𝜕𝑧=0,𝜕𝐵𝑢+𝜕𝑡𝜕𝐵𝑢²+𝜕𝑥𝜕𝐵𝑢𝑤𝐵𝜕𝑧=−𝜌_∘𝜕𝑃^∗+𝜕𝜕𝑥2𝜕𝑥𝜈+𝜈_𝑡𝐵𝜕𝑢+𝜕𝜕𝑥𝜕𝑧𝜈+𝜈_𝑡𝐵𝜕𝑢+𝜕𝑧𝜕𝑤,𝜕𝑥𝜕𝐵𝑤+𝜕𝑡𝜕𝐵𝑤²+𝜕𝑧𝜕𝐵𝑢𝑤𝜕𝑥=−𝜌−𝜌_∘𝜌_∘𝐵𝑔𝐵−𝜌_∘𝜕𝑃^∗+𝜕𝜕𝑧𝜕𝑥𝜈+𝜈_𝑡𝐵𝜕𝑢+𝜕𝑧𝜕𝑤+𝜕𝜕𝑥2𝜕𝑧𝜈+𝜈_𝑡𝐵𝜕𝑤,𝜕𝑧 (1)(2)(3) where 𝐵 is the width of channel at 𝑥; 𝑢 and 𝑤 are the width-averaged velocities along 𝑥 and 𝑧 directions, respectively; 𝜌_∘ is the density of the ambient fluid; 𝑔 is the acceleration due to gravity; 𝜈 and 𝜈_𝑡 are the molecular and eddy viscosity, respectively; 𝜌 is the local mixture density; 𝑡 is time. 𝑃^∗=𝑃−𝜌_∘𝑔ℎ is the modified width-averaged pressure including the gravity terms, where 𝑃 is total width-averaged pressure and ℎ is distance from the reference in 𝑧 direction. The following assumptions are made in (2) and (3).(1)Side shear stress terms that appear on the boundaries, due to the width averaging, are equal to zero, and slipping conditions are assumed.(2)The dispersion terms due to lateral nonuniformities of the flow quantities, which appear in the width-averaging of the 3D equations, are neglected [33].

The eddy viscosity is modeled by the buoyancy-modified 𝑘−𝜀 turbulence model [32] as𝜈_𝑡=𝐶_𝜇𝑘²𝜀̃, (4) where 𝑘and 𝜀̃ are the width-averaged values of the turbulent kinetic energy and dissipation and 𝐶_𝜇 is a model constant. Quantities of 𝑘 and 𝜀̃ are obtained from the solution of the transport equations:𝐷𝑘=𝜕𝐷𝑡𝜈𝜕𝑥𝜈+^𝑡𝜎_𝑘𝐵𝜕𝑘+𝜕𝜕𝑥𝜈𝜕𝑧𝜈+^𝑡𝜎_𝑘𝐵𝜕𝑘𝑃𝜕𝑧+𝐵_𝑟+𝑃_𝑏−𝜀,𝐷𝜀̃=𝜕𝐷𝑡𝜈𝜕𝑥𝜈+^𝑡𝜎_𝜀𝐵𝜕𝜀̃+𝜕𝜕𝑥𝜈𝜕𝑧𝜈+^𝑡𝜎_𝜀𝐵𝜕𝜀̃𝜕𝑧+𝐵𝐶_𝜀1𝜀̃𝑘𝑃_𝑟+𝐶_𝜀3𝑃_𝑏−𝐵𝐶_𝜀2𝜀²̃𝑘. (5)(6)

In the above equations, 𝐷/𝐷𝑡=𝜕𝐵()/𝜕𝑡+𝜕𝐵𝑢()/𝜕𝑥+𝜕𝐵𝑤()/𝜕𝑧 denotes the material derivative, 𝑃_𝑟 and 𝑃_𝑏 are production terms of turbulent kinetic energy due to the mean velocity gradient and the buoyancy and are estimated, respectively, as𝑃_𝑟=𝜈_𝑡2𝜕𝑢𝜕𝑥²+𝜕𝑢+𝜕𝑧𝜕𝑤𝜕𝑥²+2𝜕𝑤𝜕𝑧²,𝑃_𝑏=𝜈_𝑡𝑔𝜌_∘𝜎_𝑡𝜕𝜌.𝜕𝑧 (7)(8)

In the above equations, 𝜎_𝑘,𝜎_𝜀, and 𝜎_𝑡, turbulent Prandtl number for 𝑘, 𝜀̃ and the source of density current (in this study, width-averaged temperature 𝑆), respectively. 𝐶_𝜀1, 𝐶_𝜀2, and 𝐶_𝜀3are turbulence model constants. The following standard values are specified for the turbulence model constants: 𝐶_𝜇=0.09,𝐶_𝜀1=1.44,𝐶_𝜀2𝜎=1.92,_𝑘=1.0,𝜎_𝜀=1.3. (9)

The turbulent Prandtl number 𝜎_𝑡, the ratio of the turbulent viscosity and the turbulent diffusivity set equal to 0.9. The empirical coefficient 𝐶_𝜀3 is directly related to the buoyancy effect, and its optimal value does not have a firm basis, unlike other model coefficients. Rodi [32] suggested that the coefficient should vary between 0 for buoyant horizontal shear flows and 1 for buoyant vertical shear flow. That is, the optimal value of 𝐶_𝜀3 is dependent upon the flow conditions and its choice is controversial. 𝐶_𝜀3 are suggested values ranging from 0 to 0.4 to give a good agreement with experimental results on density currents [34–37]. In this study, we use a constant value of 𝐶_𝜀3=0.2. The 2D width-averaged unsteady RANS and turbulence closure equations are solved simultaneously with the following transport equations for the width-averaged temperature𝑆, which is used to determine scalar transport and the variation of the fluid density:𝐷𝑆=𝜕𝐷𝑡𝜈𝜕𝑥𝜈+^𝑡𝜎_𝑡𝐵𝜕𝑆+𝜕𝜕𝑥𝜈𝜕𝑧𝜈+^𝑡𝜎_𝑡𝐵𝜕𝑆𝜕𝑧. (10)

Density is assumed to be linearly related to the mean volumetric concentration through the equation of state as 𝜌=𝜌_∘(1.0+𝑚𝑆), where 𝑚 is a constant.

3. Numerical Method

Nonuniform bathymetry of the computational domains is handled by the sigma-type body-fitted grids which fit the vertical direction of the physical domain. The finite-volume integration method is used for different terms in the governing equations. A fractional step [38], or operator splitting, method is employed to integrate the governing equations in time coupled with a projection method for satisfying the continuity equation [39, 40]. This method in which the time advancement is decomposed into a sequence of steps: the advection, the diffusion, and the projection. The fractional step method allows us deploy the most efficient numerical technique for each of the individual process. Dimensional splitting is also used to reduce the multidimensional problem to a series of one-dimensional problems. These fractional step and dimensional splitting approaches save considerable amounts of computational time and preserve overall accuracy of the temporal and spatial discretization schemes [39]. The fractional step approach for solving the momentum equations is also employed to integrate the transport equations for the turbulence quantities and the temperature. The later equation only the advection and diffusion steps are required. For the former equations, however, the diffusion step is followed by a step in which the source terms are taken into account explicitly. The discretization of the advection terms in the above iterative scheme is critical for spatial accuracy of the numerical method. We employ the semi-Lagrangian method ([39] based on the Fromm scheme [41] for discretization of the advection terms of the one-dimensional form isolated by the dimensional splitting approach. The ULTIMATE conservative difference scheme of Leonard [42] is employed to maintain monotonicity while preserving spatial accuracy. Using a monotone scheme is essential for numerical stability and the smoothness of the computed solutions in the discontinuities where sharp scalar gradients are present and dominate the dynamics of the flow. The diffusion terms are solved by the semi-implicit second-order Crank-Nicolson scheme to obtain the second intermediate velocity components. During the advection step of the fractional step procedure, the discrete equations are integrated in time explicitly. In the diffusion step, the Thomas tridiagonal matrix algorithm is used to solve the discrete equations implicitly. The pressure-Poisson equation is solved using a block tridiagonal algorithm, which can be formulated by writing the discrete Poisson equation for all the cells. After calculating the velocity components and pressure field in new time step, the transport equations for the turbulence kinetic energy, energy dissipation rate, and the temperature are solved using the same fractional step algorithm used to integrate the momentum equations. Finally, the density is calculated using the state equation.

4. Boundary Conditions

The boundaries of the computational domain are inlet, outlet, free-surface, and solid walls. In the absence of wind shear, the net fluxes of horizontal momentum and turbulent kinetic energy are set equal to zero at the free surface. Flat, slip boundary condition of zero velocity gradient normal to the surface is applied, and the pressure is set equal to the atmospheric value at the surface. The dissipation rate 𝜀 is calculated from the relation given by [32]:𝑘𝜀̃=^3/2(0.43𝐷)at𝑧=𝑧_𝑠, (11) where 𝐷 is the local water depth and𝑧_𝑠 is the coordinate of free surface.

The wall function approach is used to specify boundary conditions at the bottom of the channel in order to avoid the resolution of viscous sublayer. The first grid point off the wall (center of the control volume adjacent to the wall) is placed inside the logarithmic layer based on the buoyancy velocity𝑈_𝑏. In the wall-function approach, the resultant wall shear stress vector 𝜏_𝑤 is related to the velocity vector at the first grid point 𝑉₂ by the standard, smooth-wall, log law [43]: 𝜏_𝑤=−𝜆_𝑤𝑉₂with𝜆_𝑤=𝜌𝑐_𝜇^1/4𝑘₂^1/2𝜅ln𝐸𝑧⁺₂,𝑧⁺₂=𝑢_∗𝑧^′₂𝜈, (12) where the subscript 2 refers to the point at the first control volume center at the wall, 𝑧^′₂ is the normal distance to the wall 𝑢_∗=𝜏_𝑤/𝜌 is the bed shear velocity; 𝜅 is von Karman constant (= 0.41); 𝐸 is an integration constant that for smooth beds is equal to 8.43. Here, we employ an algorithm for computing 𝑢_∗ in (12). The near-wall values of turbulent kinetic energy 𝑘 and dissipation rate 𝜀 are given as𝑢𝑘=^2∗𝐶_𝜇^1/2𝑢,𝜀=^3∗𝜅𝑧^′₂. (13)

On the other hand, if the first grid node from the wall locates within the viscous layer (𝑧⁺₂<11.6) due to the evolution of complex vertical structures, these near-wall values are approximated as𝑧𝑘=^′2𝑧^𝜈𝑘_𝑝,𝜀=2𝜈𝑘𝑧^′₂². (14)

At the inlet, known quantities are specified for the inflow velocity, the concentration of the density current source, and the current thickness, while turbulent kinetic energy and dissipation rate at the inlet are estimated as follows [37]:𝑘_in=0.1𝑢_in²,𝜀̃_in=10𝑘^3/2_in𝐶_𝜇^3/4𝜅ℎ_in, (15) where 𝑢_in, 𝑘_in and ℎ_in are the averaged velocity, the turbulent kinetic energy, and the thickness of the gravity current at the inlet. At the outlet the normal gradients of all dependent variables are set equal to zero.

5. Results and Discussion

In this section, we report numerical results for density currents from continuous sources along with comparisons with the experimental and numerical results of Johnson and Stefan [31] and Bournet et al. [29], respectively. In Figure 1, plan and longitudinal section view of diverging channel with sloping bottom at the experiment of Johnson and Stefan [31] is shown. The parameter values used in the calculation for the numerical simulation of these cases are summarized in Table 1 which includes inflow rate(𝑞₀), inlet temperature(𝑇_inlet), ambient water temperature(𝑇₀), divergence half-angle of channel(𝛿), the bottom slope(𝑆), width of inflow channel (𝐵₀), and depth of inflow channel (𝐻₀).

Table 1 

The parameter values for the numerical simulations of the continuous density currents in an inclined diverging channel.

(a) Plan

(b) Longitudinal section

Figure 1 

Plan and longitudinal section view of diverging channel with sloping bottom at the experiment of Johnson and Stefan [31].

We first present the comparison of the numerical solutions and experimental measurements published by Johnson and Stefan [31]. A 21.8 m long computational domain is employed in order to avoid reflections from the outlet. The calculations consider 200×30uniform grids in the 𝑥- and 𝑧-directions, respectively, and a time step of 0.1 second was adopted. We have continued our unsteady simulations until the solutions in the region of interest (0≤𝑥≤6.0m) reach a quasiequilibrium state. In Figure 2, the simulated dimensionless temperature distribution of the Johnson and Stefan experiment [31] is shown at a quasiequilibrium state. The numerical simulation of dimensionless temperature is compared with experimental results of Johnson and Stefan [31] at two sections (𝑥=4.2m and 𝑥=6.0 m) in Figure 3 which reveals very good agreement between two sets (numerical and experimental) of results. The dimensionless temperature is expressed in the following form [31]:𝑇𝑇=⁰−𝑇𝑇₀−𝑇_inlet, (16) where 𝑇, 𝑇₀, 𝑇, and 𝑇_inlet are dimensionless temperature, ambient water temperature, temperature, and inlet temperature, respectively.

Figure 2 

Simulated dimensionless temperature distribution of the Johnson and Stefan experiment [31] at a quasiequilibrium state.

(a)

(b)

Figure 3 

Comparison of simulated dimensionless temperature and experimental results of Johnson and Stefan [31] at two sections.

In the second stage, we present the comparison of the solutions and numerical results published by Bournet et al. [29] after the head of the current has passed and the under current reaches steady state. In Figure 4, plan and longitudinal section view of diverging Channel with sloping bottom at the simulation of Bournet et al. [29] is shown. A 5000 m long computational domain is employed in order to avoid reflections from the outlet. The parameter values used in the calculation for the numerical simulation of this case are summarized in Table 1.

(a) Plan

(b) Longitudinal section

Figure 4 

Plan and longitudinal section view of diverging channel with sloping bottom at the simulation of Bournet et al. [29].

In the upstream region of the channel (0≤𝑥≤1500m), the calculations consider 70×20 nonuniform grids in the 𝑥- and 𝑧-direction, respectively, similar to that in Bournet et al. [29]. Cells are unevenly spaced, and the grid is squeezed upstream and near the bottom where shear stress occurs and velocity gradients develop. A time step of 0.5 second is adopted.

To compare numerical solutions with the results of Bournet et al. [29], we have continued our unsteady simulations until the solutions in the region of interest (0≤𝑥≤1500m) reach a quasiequilibrium state. In Figures 5 and 6, the simulated velocity and temperature field are compared with the results of Bournet et al. [29] which reveal very good agreement between two sets of results.

(a)

(b)

Figure 5 

Simulated Velocity Field (a) the results of Bournet et al. [29] and (b) the present results.

Figure 6 

Comparison of simulated iso-contours (dashed lines) of temperature distribution and the results of Bournet et al. [29] (solid lines) at a quasiequilibrium state.

Figure 7 shows the flow evolution of the density current computed on at several instants in time. Note that the 𝑥-axis scales are different in subfigures of Figure 7. The numerical solutions reproduce the variation of the longitudinal profile of the density current height and the propagation of the current head. At the first because of dominant inertia force, the total depth of density current is moving forward (Figure 7, 𝑡=1 min). At time about 2 minutes, due to dominant buoyancy force, density current plunges to the reservoir bottom and then moves near it. At the plunge point, intense turbulence causes too much mixing of dense fluid and less dense ambient fluid. A dense fluid propagates rightward along the bottom boundary, while the ambient fluid moves leftward along the surface boundary naturally due to the density difference. About 5 minutes, density current has been established and then the current head is seen to grow as it travels downstream as the dense fluid is pushed by the vortices from the rear part into the head region.

Figure 7 

Sequential snapshots of density distributions computed at selected time instants (note that the 𝑥-axis scales are different in subfigures).

6. Conclusions

A series of 2D width-averaged URANS simulations are carried out to resolve continuous density currents in an inclined diverging channel. Comparison of the numerical solutions with available experimental measurements and numerical results leads to the conclusion that 2D width-averaged unsteady RANS simulations using a buoyancy-extended 𝑘−𝜀 turbulence model can resolve reasonably well steady state flow in the rear part. The numerical solutions reproduce the variation of the longitudinal profile of the density current height and the propagation of the current head.

Detailed comparisons have shown in general that this two-dimensional width-averaged model can be applied to the simulation of the density current in reservoirs which the Froude number of the inlet dense flow and divergence half-angle of reservoir (𝛿)are less than 2 and 20°, respectively. Although the model presented in this paper is for density current, it can be extended to model two-dimensional turbidity current in an inclined diverging channel. Such a modification is currently being undertaken by the authors.