Three-Dimensional Numerical Modeling of Large-Scale Free Surface Flows with Heat and Contaminant Transfer in Reservoirs

Diao, Wei; Yuan, Hao; Chen, Liang; Zhang, Xujin; Zhang, Cunze

doi:https://doi.org/10.1155/2021/9336868

Advances in Civil Engineering

On this page

Abstract Introduction Model Verification Conclusion Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2021 | Article ID 9336868 | https://doi.org/10.1155/2021/9336868

Three-Dimensional Numerical Modeling of Large-Scale Free Surface Flows with Heat and Contaminant Transfer in Reservoirs

Wei Diao,¹Hao Yuan,²Liang Chen,³Xujin Zhang,²and Cunze Zhang²

Academic Editor: Minghui Li

Received29 May 2021

Revised10 Aug 2021

Accepted27 Aug 2021

Published10 Sept 2021

Abstract

The temperature distribution and pollutant distribution in large reservoirs have always been a hotspot in the field of hydraulics and environmentology, and the three-dimensional numerical modeling that can effectively simulate the interactions between the temperature fields, concentration fields, and flow fields needs to be proposed. The double-diffusive convection lattice Boltzmann method is coupled with a single-phase volume of fluid model for simulating heat and contaminant transfer in large-scale free surface flows. The coupling model is used to simulate the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. The mechanism of convection-diffusion, gravity sinking flow, and the complexity of the temperature and the pollutant redistribution process are analyzed. Good agreements between the simulated results and the reference data validate the accuracy and effectiveness of the proposed coupling model in studying free surface flows with heat and contaminant transfer. At last, the temporal and spatial variations of flow state, water temperature stratification, and pollutant transport in the up-reservoir of a pumped-storage power station are simulated and analyzed by the proposed model. The obtained variations of the flow field agree well with the observations in the physical model test and in practical engineering. In addition, the simulated temperature field and concentration field are also consistent with the general rules, which demonstrates the feasibility of the coupling model in simulating temperature and pollutant distribution problems in realistic reservoirs and shows its good prospects in engineering application.

1. Introduction

A great number of large reservoirs have been built in China. By 2017, there are 732 reservoirs with a water capacity of over 100 million cubic meters [1]. The reservoirs play a crucial role in hydroelectric power, flood control, navigation, fisheries, and other comprehensive benefits. However, the water temperature stratification and pollutant accumulation exist widely in these huge reservoirs, which cause a detrimental effect on the ecological environment [2]. Therefore, to study the temperature and the pollutant redistributions in reservoirs, an effective and accurate numerical model for simulating heat and contaminant transfer in large-scale free surface (LS-FS) flows needs to be proposed, which is beneficial to the harmonious hydraulic development and ecological protection [3].

At present, the general methods for studying heat and contaminant transfer in LS-FS flows are mainly two-dimensional (2D) numerical models, especially shallow water equations coupling convection-diffusion equations for simulating planar flow problems [4, 5]. The 2D numerical simulations have been one of the most commonly used methods to analyze large-scale flows with heat and contaminant transfer, such as flows in reservoirs [6], lakes [7], and oceans [8]. However, the 2D method could not analyze the vertical temperature and pollutant distribution, and the plane simulation results are also not accurate enough because of the lack of vertical velocity effect, which evidently limits the applications of 2D models in engineering practice. The flows with an obvious three-dimensional (3D) convection-diffusion mechanism shall be simulated by more promising 3D methods. Typically, 3D numerical methods for simulating LS-FS flows are mainly based on gas-liquid two-phase flow models [9], interface tracking models [10–12], and interface capturing models [13–15]. Nevertheless, the flow problems in the field of hydraulic projects are usually complicated because of the large-scale and complex boundary conditions, which not only results in a dramatic amount of computing resource consumption but also makes the stable and accurate simulation challengeable [16].

Single-phase free surface lattice Boltzmann (SPFS-LB) model was originally proposed by Thürey in 2003 [17]. Compared with the traditional two-phase flow models in which both gas and water are simulated, the SPFS-LB model is more powerful to simulate 3D LS-FS flows, because it only simulates the water flows and ignores the movement of air, so the consumption of computational resources is lower [17, 18]. The lattice Boltzmann models for solving convection-diffusion problems can be mainly classified in the multidistribution function approach [19], the multispeed approach [20], and the combined approach [21]. Among the above approaches, the multidistribution function approach is the most mature and has been used successfully in many fields [19] but has seldom been applied to heat and contaminant transfer in LS-FS flows.

To study the 3D temperature and pollutant distribution in large reservoirs, the double-diffusive convection LB method is coupled with the SPFS-LB model. The rest of the paper is organized as follows. Firstly, the basic ideas and algorithms of methods are briefly described, and the treatments of the coupling scheme for improving stability and precision are addressed. Then, the accuracy and effectiveness of the proposed method are verified by two benchmark cases: the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. At last, as an attempt of simulating engineering practice, the process of the flow state, water temperature, and pollutant variations in a practical reservoir is simulated and analyzed.

2. Method

2.1. The Single-Phase Free Surface LB Model

Unlike traditional two-phase flow methods that set up different distribution functions (DFs) for water and air, the SPFS-LB model only simulates water flows (the denser and more viscous phase) by the single-phase LB equations and tracks the free surface by a constructed boundary treatment. The following assumptions should be made before the simulation by the SPFS-LB model [17]: (1) the effects of air flows (lower density and viscosity) on water flows can be neglected; (2) the air can reach the equilibrium state immediately after the state of water is changed; (3) the simulated water is separated by a closed layer of the interface cells (the free surface) from air. Absolutely, the large-scale water flows in reservoirs are in line with the above assumptions. Three types of computational cells are introduced to capture the interface, including the empty cells, interface cells, and filled cells, as shown in Figure 1. The empty cells do not contain any water, so there is no need to define physical quantity. The filled cells are full of simulated water, in which the definition of physical quantity and the evolution processes are the same as those in the traditional single-phase LB models. The interface cells are tracked by computing the water volume fraction ε, which can be defined as ε = m/ρ, with m and ρ being the mass and density of cells, respectively. And the values of empty cells, interface cells, and filled cells are 0, 0-1, and 1, respectively.

2.1.1. Computation of Mass Flux

The movement of the interface in the SPFS-LB model is realized by the transformation of cell types. And the cell type is judged by the volume fraction ε which can be updated by calculating the inflow and outflow of cell mass. By conducting discretization in time and on uniform lattices and applying the divergence theorem to the convective term, the continuity equation in discrete form can be derived as m (x, t + Δt) − m (x, t) = Δm_i, where Δm_i denotes the mass variation on ith direction along with the particle velocity [22]. In the LB method, the mass flux can be simply computed through the two antiparallel particle distribution functions f_i and f_inv(i) (see equations (1) and (2)), where e_inv(i) = −e_i [23].

Therefore, the mass m of the interface cell at lattice node x for the time step t + Δt can be updated the following equation, which indicates that the updated m is obtained by summing the mass flux in all velocity directions.where n denotes the number of discrete velocities. It can be proved that the mass updated by equation (3) is conserved [22].

2.1.2. Reconstruction of Interface Cells

In the SPFS-LB model, the propagation and collision of LB equations can be normally performed in filled cells for that they are not adjacent to empty cells. However, the interface cells are always surrounded by empty cells, where physical quantity and DFs are not defined (shown in Figure 1). According to assumption 2, the missing DFs of the interface cells can be reconstructed based on the macroscopic boundary conditions and the force balance for opposite lattice directions [17]. The following equation shows the reconstructed DFs of interface cells at x with surrounding empty cells at x + e_iΔt, where u and denote the velocity and the pressure of interface cells at position x, respectively.

It should be pointed out that not only the missing DFs but also all the DFs whose discrete velocity e_i satisfies e_in < 0 need to be reconstructed through equation (4) [17]. The vector n refers to the surface normal direction, which is obtained by the second-order central difference approximation .

2.1.3. Mass Allocation

Interface cells need special treatments after the updating of their mass m and volume fraction ε. The interface cells whose updated ε is greater than 1 shall be converted to filled cells, and the excessive mass needs to be allocated to the adjacent interface cells and empty cells. After getting mass, the surrounding empty cells shall be turned into interface cells, and their macroscopic variables can be obtained by averaged velocity and density of the surrounding nonempty cells. The DFs of new emerging interface cells can be initialized by the equilibrium distribution functions of the LB method [18]. Similarly, the interface cells whose ε is less than 0 shall be converted to empty cells and the surrounding filled cells turn to interface cells. The negative m needs to be compensated to 0 by the neighboring nonempty cells. The proportion of the excess mass distribution and negative mass compensation can be obtained by the dot product between the interface normal direction and the distribution direction. The equations can be expressed as [17]where m^ex and υ refer to the excess mass and proportion, respectively. The interface moves as the completion of mass allocation and cell type conversion.

2.2. Double-Diffusive Convection Lattice Boltzmann Method

From a large-scale point of view, the water flows in reservoirs are normally considered incompressible, whose viscosity, density, thermal diffusivity, and contaminant diffusivity can also be assumed to be constant except for the buoyant. Therefore, the governing equations for simulating heat and contaminant transfer in reservoirs can be expressed with the Boussinesq assumption [21] aswhere F, ν, α, and D denote body force term, viscosity, thermal diffusivity, and contaminant diffusivity, respectively. In [24], the double-diffusive convection lattice Boltzmann (DDC-LB) method is used to solve equation (6). The DDC-LB method simulates velocity, temperature, and concentration fields by three sets of respective DFs, and the model couples the three fields by a body forcing in the velocity field. Here, the method is extended to 3D form on the basis of the 2D DDC-LB method in [24], which simulates velocity field by D3Q19 model and uses D3Q6 model for temperature fields and concentration fields simulation. The LB equation of the 3D model can be described as

f_i, , and h_i refer to the DFs for the velocity field, temperature field, and concentration field. Variable τ denotes the relaxation factor, which can be determined by the viscosity and diffusivity, and Prandtl number (Pr) can be obtained by τ (see the following equations):

The equilibrium distribution functions in equations (10)–(12) can be expressed as

To derivate governing equation (6) by Chapman–Enskog expansion, coefficients σ, λ, and γ need to satisfy λ + 2γ = σ and λ + 2γ = 1/3, which are chosen as 1/3, 1/6, and 1/12 in this work, respectively. s_i (u) is a function of the discrete velocity u:

ω_i is the weight coefficient, ω₀ = 1/3, ω_1∼6 = 1/18, and ω_7∼18 = 1/36 in the D3Q19 model. The macroscopic quantities u, P, T, and C are the macroscopic velocity, pressure, temperature, and concentration; they can be obtained by the first or the secondary moment of the corresponding DFs as

Simulating 3D temperature fields and concentration fields in large-scale flows by the D3Q6 model can reduce calculation amount and save computer memory usage. However, when precisely analyzing the small-scale flow problems with complex boundary conditions, using the D3Q19 model instead of the D3Q6 model can retain better accuracy and stability.

The number of computational grids for simulating the model reservoir (in Section 3.2) is over 10 million (2440 × 46 × 92), and there are over 37 million (1024 × 1024 × 36) for simulating the practical reservoir (in Section 4). The computer memory usage and computation time are estimated to be increased by 45% and 23%, respectively, if the temperature fields and concentration fields are simulated by the D3Q19 model instead of by the D3Q6 model. Due to the constraint of computer performance, the D3Q6 model is used for the present work. Good agreements between the simulated results and the reference data in Section 3 show that the results simulated by the D3Q6 model are accurate enough.

2.3. Key Treatments of the Coupling Scheme

The SPFS-LB model is firstly coupled to the DDC-LB method in the present work, and several key treatments need to be adopted to make the coupling scheme compatible and to obtain reasonable simulated results.

2.3.1. Reconstruction and Initialization of Temperature Fields and Concentration Fields in Interface Cells

As described in Section 2.1.2, the undefined DFs of empty cells would enter into interface cells, so the DFs of temperature fields and concentration fields in interface cells need to be reconstructed after the LB propagation. According to assumption 2 (introduced in Section 2.1) and referring to the reconstruction of flow fields, the missing DFs in the temperature field and concentration field can be calculated by the following equations:

Some empty cells would convert to interface cells after the update of cell types, and the temperature T, concentration C, DFs , and DFs h_i should be initialized. This coupling model neglects the exchange (absorption or release) of heat and contaminant between water and air. Therefore, the macroscopic T and C of new emerging interface cells originate from the surrounding excessive cells, and they can be obtained by the weighted average of the variables of surrounding new filled cells. The weight can be computed according to the proportion of the assigned excess mass. Referring to the initialization of DFs f_i, DFs , and DFs h_i in new converted interface cells can be initialized by equations (10)–(12).

2.3.2. The Force Term for the Buoyancy Flows

The flow velocity dominates the transfer of heat and contaminant, while the buoyancy induced by the nonuniform distribution of temperature and concentration has an impact on the water flows. To simulate the buoyancy flows and realize the two-way coupling between the three fields, an additional body forcing term Fi needs to be attached to equation (7). The present model adopts the approach proposed by Cheng and Li [25] on account of its stability and accuracy for treating nonuniform, unsteady source terms. The forcing term Fi can be calculated by the following equation according to [25]:wherein which A and B denote the source term in the continuity equation and momentum equation, respectively. Considering equation (6), the present model sets A = 0 and B = − (β_T (T − T₀) + β_C (C − C₀)).

2.3.3. Large Eddy Simulation for Turbulence

Turbulence exists in practical large-scale water flows with heat and contaminant transfer. The subgrid-scale model based on large eddy simulation is introduced to the coupling scheme. The physical variables are separated into large-scale variables and small-scale variables by a certain filter function in the subgrid-scale model [26, 27]. Take φ = + φ′, for example. The large-scale variable can be expressed as follows. φ′ presents the small-scale pulsation, which is modeled in a simplified way.

The following equation should be obtained through filtering equation (7) and assuming .

The relaxation time τ_e is not a constant. It varies with space and time, which can be calculated by through equation (8). The variable ν_e can be understood as the total viscosity ν_e = ν₁ + ν_t, in which ν₁ and ν_t denote the molecular viscosity and eddy viscosity (or turbulence viscosity), respectively. ν₁ is determined by the physical property of water, and ν_t can be computed by , where C, Δ, and refer to Smagorinsky constant, filter width, and strain-rate tensor, respectively. According to [27], the strain-rate tensor S can be easily computed by using the second-order moments of the filtered particle distributions.

2.3.4. The Coupling Procedure of the Coupling Scheme

This section makes a summary of the calculation procedure of the coupling scheme proposed in the paper. Before the iteration, the geometry needs to be initialized, and the reasonable boundary condition should be adopted in accordance with the simulated flow problems. Then given all the variables at time step n, the procedure for marching to step n + 1 is as follows:(i)Perform LB propagation with the boundary conditions(ii)Compute the mass flux, and update the mass of cells as described in Section 2.1.1(iii)Calculate the volume fraction of nonempty cells(iv)Reconstruct the DFs of the flow field, temperature field, and concentration field by equations (4) and (15)∼(16)(v)Compute the macroscopic velocity, pressure, temperature, and concentration through equation (14)(vi)Obtain relaxation time τ_e (Section 2.3.3), add body force term (Section 2.3.2), and perform LB collision through equation (7)(vii)Allocate the excessive or negative mass through equation (5), and update cell types as detailed in Section 2.1.3(viii)Initialize the macroscopic quantity and distribution functions of new emerging interface cells as described in Sections 2.1.3 and 2.3.1(ix)Perform the analysis of convergence

3. Model Verification

To validate the feasibility and accuracy of the proposed coupling scheme, the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir are simulated, and the detailed comparison and analysis are made between the simulated results and the reference data. If not otherwise stated, the acceleration of gravity, water density, and kinematic viscosity are set as = −9.8 m/s², ρ = 10³ kg/m³, and ν = 10⁻⁶ m²/s, respectively.

3.1. Double-Diffusive Natural Convection in a Cubic Cavity

In the heat and contaminant transfer flows, the double-diffusive natural convection in a cubic cavity reflects the essential problem of two-way coupling between temperature field/concentration field and flow field. At present, there have been many classic physical experiments or numerical research results [28, 29]. In a cubic cavity of L × L × L, the front (wall 5) and back (wall 6) are smooth walls, and the rest (wall 1∼4) are nonslip walls (see Figure 2). The temperature and concentration at the left and right boundaries are kept constant at T₀/C₀ and T₁/C₁, respectively (T₀ < T₁; C₀ < C₁); the top and bottom are adiabatic for heat transfer and impermeable for contaminant transfer; that is, the derivatives of temperature fields and concentration fields are always 0 along the normal direction of the boundary. Initially, the cubic cavity is filled with homogeneous fluid with velocity u = 0, temperature T = (T₀ + T₁)/2, concentration C = (C₀ + C₁)/2, and Pr = 0.71. The flow characteristics of double-diffusive natural convection can be estimated by the following three dimensionless numbers: temperature Rayleigh number , concentration Rayleigh number , and Lewis number . The initial values of temperature T₀/T₁ and concentration C₀/C₁ can be assigned to arbitrary numbers that satisfy the formulas above, such as T₀/C₀ = −1 and T₁/C₁ = 1 in the present simulation.

In order to quantify the simulated results, several dimensionless numbers along the right hot wall are calculated below, which are local Nusselt number, average Nusselt number, local Sherwood number, and average Sherwood number, respectively, as shown in the following equations:

Firstly, the Rayleigh number and Lewis number are set as Ra_T = 10⁴, Ra_C = 10⁴, and Le = 1, respectively. The fluid in the cavity is driven by both temperature field and concentration field. The added driving flows are generated if the driving directions of temperature and concentration are the same. Otherwise, the opposed driving flows are generated. Both flows are simulated in this section. The computational grid is set as 128 × 128 × 128. It should be noted that the simulated distribution of temperature field and concentration field is exactly the same, for the same convection-diffusion equation is applied to simulate heat transfer and contaminant transfer, and the same calculation parameters are set in the present simulation in order to compare with the research results in [29].

Figures 3 and 4, respectively, show the streamline and isothermal surfaces (or isoconcentration surfaces) after the convergence of simulation. For the added driving flows, the isothermal surfaces move farther on the upper and lower wall compared with the single-diffusive convection under the same Rayleigh number, which shows that the heat and contaminant spread more fully in the cubic cavity. The reason is that the flow field is driven by the buoyancy from temperature difference and concentration difference simultaneously, and the convection that leads to the spread of heat and contaminant is more significant. As shown in Figures 3(a) and 4, driven by the downward buoyancy around the left wall and the upward buoyancy around the right wall, the thin boundary layers are formed on the two walls. The flow forms a stable cycle, while a stagnation area appears in the center of the cubic cavity.

(a)

(b)

For the opposed driving flows, the buoyancy caused by temperature difference and concentration difference is exactly equal in magnitude and opposite in direction, and their driving effect on the flow field also cancels each other out. The fluid in the cavity does not flow at all (see Figure 3(b)). In this case, the propagation of heat and contaminant is completely dependent on diffusion, which can be proved by the distribution of equidistant vertical isothermal surfaces/isoconcentration surfaces in Figure 3(b). The above results are consistent with the research in [28, 29].

The local Nusselt number of the added driving flows along the right wall (hot wall) is calculated and plotted in Figure 5. In the region of Z/L = 0∼0.13, the isothermal surfaces are vertically distributed, and the temperature gradient along the Y direction is basically unchanged (see Figure 3(a)). The Nu rises slightly from 4.18 to 4.43 with the increase of boundary layer velocity. After reaching Z/L = 0.13∼0.96 regions, the velocity increases sharply and the transverse convection enhances, leading to a significant decrease of the temperature negative gradient −∂T/∂y along the Y direction. Therefore, Nu plummets from 4.43 to 0.64. Subsequently, influenced by the upper wall, the velocity and temperature gradient remain relatively stable, and Nu also stabilizes at about 0.64. It can be seen from Figure 5 that the simulated results are very close to the data in [29], except for slight differences on the upper and lower corner of the right wall (the maximum relative error is 0.72%). The relative error is inevitable, and it is enlarged especially at the corner, for different numerical methods and boundary treatments are adopted in this paper and the reference.

Additionally, the added driving double-diffusive convection with the temperature Rayleigh Ra_T = 10⁷ and concentration Rayleigh number Ra_C = 10⁶, 5 × 10⁶, 1.5 × 10⁷, and 5 × 10⁷ is simulated. The average Nusselt number and average Sherwood number are calculated along the hot wall. According to Table 1, the dimensionless numbers from the present simulations agree well with the results in [29], and the errors are all within 1% (the literature results only kept a decimal place, so it is difficult to measure the errors more accurately). These analyses demonstrate the accuracy of the proposed coupling model in simulating the theoretical heat and contaminant transfer flows.

3.2. Temperature Distribution in a Model Reservoir

In this section, the variation process of temperature field and flow field in a model reservoir is simulated by the proposed coupling model, so as to demonstrate its effectiveness in simulating practical physical problems. The simulation is based on the data of gravity undercurrent experiments, tested by the US Army Corps of Engineers, Waterways Experiment Station [30, 31]. The geometry is shown in Figure 6. The inlet is a 0.15 m high hole at the bottom of the left-hand wall, while the outlet is a banded hole that is 0.04 m high and 0.15 m away from the bottom of the right-hand wall. Initially, the reservoir water is stationary with a uniform temperature of 21.44°C. After the simulation starts, the 16.67°C cold water is discharged at a flow rate of 0.00063 m³/s, while the water flows out from the outlet at the same flow rate.

Due to the lateral symmetry of the simulation, only a half is selected for calculation and results show (see Figure 6). Correspondingly, the computing grid is 2440 × 46 × 92, and the symmetric plane is treated as the symmetric boundary condition. The inlet and outlet are set as a velocity boundary condition of 0.014 m/s and 0.0173 m/s, respectively, and the rest are adiabatic and nonslipping solid wall boundaries.

Figures 7 and 8 show the distribution of water temperature and flow velocity at different times, respectively. It can be seen that the entire variation process of temperature in the model reservoir can be divided into three stages:(1)Jet penetrating stage: once the cold water flows into the reservoir, it sinks and forms an undercurrent that moves along the bottom of the reservoir due to its higher density. And a large temperature gradient zone appears near the inlet, as shown in Figures 7(a) and 7(b). Under the influence of the shear force formed by the undercurrent, reverse flows occur in the surface waters and at the forward position of cold current, resulting in the tongue-shaped uplift at the front of the isothermal surface. See Figures 8(a) and 8(b) for details.(2)Tongue-shape isothermal surface stage: when the cold current forward reaches the vertical expansion section of the reservoir along the bottom slop, it accelerates and keeps a high speed due to the gravity, and the tongue-shaped isothermal surfaces form (see Figures 7(c) and 7(d)). By this time, the temperature difference between the cold current and the upper warm water inhibits the transfer of vertical momentum and heat, which maintains a large vertical velocity gradient and temperature gradient. Gradually, the tongue-shaped region becomes thicker and local thermal stratification appears with the further advance of the cold current into the deeper reservoir. The velocity of cold current forward decreases and temperature rises slowly during the process (see in Figures 8(c) and 8(d)).(3)Horizontal thermal stratification stage: when the tongue-shaped cold current comes to the right wall of the reservoir, the thermal stratification formally formed. The backflow area at the top of the reservoir stretches across the entire reservoir and reaches the longest, as shown in Figure 8(e). As the cold inflow water forward flows out of the reservoir, its disturbance to the upper warm water decreases and the surface backflow intensity declines gradually. The major heat-transfer mechanism changes from convection to diffusion, which is also the reason why the water temperature in the reservoir ultimately presents horizontal stratification. At last (T = 50 min), the 20°C isothermal surface expands to the whole reservoir, and the 18°C and 19°C isothermal surfaces are also in parallel with the bottom. The flow condition of the whole reservoir is composed of warm backflow at the upper part and cold current at the lower part, as shown in Figures 7(f) and 8(f).

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(c)

(d)

(e)

(f)

Based on the above analyses, the flow field and the temperature field are strongly coupled. The accurate simulation of buoyant flows is the key to study the thermal stratification of the reservoir water. Caused by the temperature difference, the inflow water sinks and the transfer of vertical heat is restricted during the whole moving process of the cold current, which strongly affects the formation, development, and stabilization of flow state and temperature distribution. It also explains the mechanism of horizontal thermal stratification in reservoirs.

Figure 9 compares the velocity profile along a vertical line at 11.43 m from the reservoir inlet at t = 11 min (see Figure 6 for the position of the line), and Figure 10 shows the histories of the discharge water temperature. The maximum velocity of the cold current and upper warm water is about 0.024 m/s and −0.0065 m/s, respectively. By comparing the results simulated by the proposed coupling model with the experimental results [30, 31] and the results obtained by Fluent software [32], the current layer simulated in this paper is thicker, and the velocity of backflow is smaller. Some errors exist because it cannot sensitively capture the sharp velocity gradient of the boundary layer. Besides, no boundary layer model is adopted to simulate the turbulence near the reservoir bottom. In terms of the temperature of the outflow water, the results from the present coupling model are superior to the results calculated by commercial software [32]. In summary, the above simulation shows that coupling the DDC-LB method and SPFS-LB model can simulate the large-scale thermal flows in model reservoirs effectively.

4. Temporal and Spatial Variations of Temperature and Concentration in a Practical Reservoir

An upper reservoir of a pumped-storage power station is about 700 m in length and width, and the maximum depth is about 25 m (shown in Figure 11). The inlet, located at the bottom of the right-hand wall, is a rectangular region within the coordinate range of X = 233∼265 m, Y = 700 m, and Z = 0∼5 m. Under the pumping operation modes, there is a continuous point pollution source with the unit concentration (C₀ = 1) in the center of the reservoir bottom (X = 350 m, Y = 350 m, and Z = 3.5 m). The simulation is carried out with a lattice of 1024 × 1024 × 36 (Δx = 0.68 m). The bottom and walls of the reservoir are set as nonslipping solid wall boundaries that are adiabatic and permeable.

The changing processes of the flow field, temperature field, and concentration field in the upper reservoir under the pumping condition are simulated. The simulation conditions are as follows: ① at the time of t = 0 min, there is no water in the reservoir, and the warm water with the temperature T₀ = 0°C flows into the reservoir at a flow rate of 640 m³/s; ② at t = 70 min, the reservoir has about 1/2 capacity left. At this time, adjust the inflow into the cold water of T₁ = −1°C, and remain the flow rate unchanged; and ③ at t = 100 min, only about 1/4 of the reservoir capacity is left, and the pumping operation ends.

4.1. Variations of Flow Fields

Figure 12 shows the streamlines in the reservoir and the velocity distribution on the slice of Z = 7 m after pumping for 5 min, 15 min, 85 min, and 100 min, among which the blue area is the water area.

(a)

(b)

Under the pumping operation modes, the water bypasses Mountain A and Mountain B (Figure 11) and rushes to the deepest part of the reservoir basin once it flows into the reservoir. Under the influence of the bottom slope, the water flow keeps a high velocity, up to 4 m/s (Figure 12(a)), which will scour the reservoir basin and disturb the deposited sediment. As shown in Figure 12(b), the water flow at the inlet is deflected by the right bank slop of the reservoir, and it impacts Mountain A from the left side, which will cause bank erosion to Mountain A and thus affect its stability. The water flows into the open area of the reservoir basin after passing through Mountain A, forming a clockwise circulation with a diameter of over 100 m. After 15 minutes, the inflow water covers the whole basin, and the water surface calms down and gradually rises. The incoming water accounts for 65% of the reservoir capacity at the time of t = 85 min. As the flow rate at the right side of Mountain A is much larger than that at the left side, a counterclockwise circulation across the entire reservoir forms with a diameter of over 300 m, as shown in Figure 12(a). Under the disturbance of Mountain A and Mountain B, a small counterclockwise circulation (about 50 m in diameter) appears between the two mountains, and this place will be a gathering area of floating objects (such as garbage, vegetation, and others). The flow velocity decreases gradually with the rise of the water level. At t = 100 min, the flow state of recirculation zones becomes stable, and the streamlines become smooth; see Figure 12(d). The above-obtained variation laws of the flow field agree with the observations in the physical model test and in practical engineering [33].

(a)

(b)

(c)

(d)

4.2. Variations of Concentration Fields

Figure 13 shows the isoconcentration surfaces (C = 0.4, 0.2, 0.1, and 0.05) in the reservoir after pumping for 5, 15, 85, and 100 minutes.

(a)

(b)

(c)

(d)

The point pollution source is located between the two mountains. According to the analysis of flow fields in Section 4.1, the water flows from right to left at a fast speed under the action of Mountain A, and the pollutants also rapidly spread through the water before the time of t = 15 min. The spread speed of isoconcentration surface (C = 0.05) is equal to the flow velocity, and the pollutants almost do not propagate in the reverse flow direction (Figures 11(a), 11(b), 13(a), and 13(b)). It is because the velocity here is pretty high. Compared with diffusion, convection plays a dominant role in pollutants transfer. Pumping water should be stopped and measures should be taken to prevent the dispersal of pollutants immediately. At the time of t = 85 min, affected by the counterclockwise circulation, the pollutants bypass Mountain A and spread to the deepest part of the reservoir basin, while their propagation distance in the clockwise direction is relatively close due to the significant inhibition effect of circulation, as shown in Figures 11(c) and 13(c). In order to control the transfer pollution, effective engineering measures need to be taken to stop the counterclockwise circulation of reservoir water. As the flow velocity in the reservoir slows down, the convection that causes the spread of pollutants gradually weakens and the diffusion effect on pollutants propagation appears. Finally, at 100 min after pumping, the pollutants almost diffuse to the entire bottom area, as shown in Figure 13(d). The above analysis of the pollutant propagation process is of certain guiding significance to the formulation of remedial measures and emergency treatment plans for pollution leakage accidents in reservoirs.

4.3. Variations of Temperature Fields

At the time of t = 70 min, the cold water is discharged into the reservoir. The isothermal surfaces of T = −0.2, −0.1, −0.05, and −0.02 at t = 71, 75, 85, and 100 minutes are shown in Figure 14, in which the light blue area is the water surface.

(a)

(b)

(c)

(d)

When the cold water just enters the reservoir (t = 71 min), it sinks to the bottom due to its relatively high density. The cold water is not fully mixed with the warm water in the reservoir, and a large vertical temperature gradient is maintained at the inlet (Figure 14(a)). Five minutes later, the cold water forward reaches Mountain B, forming an undercurrent that slides into the depth of the reservoir basin at a speed of 1 m/s along the bottom slope of the reservoir. As the vertical temperature difference inhibits the transfer of vertical momentum and heat, the isothermal surfaces present a flat tongue shape, as shown in Figure 14(b). Subsequently, the advance speed of the undercurrent slows down as the bottom slope of the reservoir decreases. As shown in Figure 14(c), cold water reaches the opposite reservoir bank at t = 85 min. At the time of t = 100 min, the cold water bypasses Mountain A and reaches the deepest part of the reservoir basin. In this process, due to the slowing down of the undercurrent velocity, the dominated heat-transfer mechanism changes from convection to diffusion, and the isothermal surfaces also develop from tongue shape to pie shape, as shown in Figure 14(d). The development process and formation mechanism of vertical water temperature distribution in this reservoir are basically consistent with the analysis results in Section 3.2.

This paper summarizes the general conclusion that explains the spatial-temporal variation process of water temperature in the case of the injection of cold water into a warm reservoir. On the whole, the coupling model proposed in this paper is effective and feasible to simulate the large-scale free surface flows with the heat and contaminant transfer in reservoirs.

5. Conclusion

The SPFS-LB model is firstly coupled to the DDC-LB method in this paper to simulate large-scale flows with heat and contaminant transfer. To make the coupling scheme compatible and to obtain reasonable results, several key treatments described in Section 2.3 need to be adopted. The feasibility and accuracy of the proposed model are validated by two typical examples, namely, the double-diffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. Finally, as an attempt at the field of engineering applications, the coupling scheme is applied to study the temporal and spatial variations of temperature and concentration in a realistic up-reservoir of a pumped-storage power station in detail, and some general conclusions of heat and contaminant transfer in reservoirs are drawn.

Compared with the traditional 2D numerical models, 3D temperature, pollutant, and velocity distribution can be obtained in this work. In addition, the multifield coupling mechanism of this model makes the simulation more accurate. These are the main advantages of the proposed coupling model. It should be noted that the proposed model ignores the exchange (absorption or release) of heat and contaminant between water and air. The nonpoint source term will be introduced in equation (7) to solve the problems in the following work. The wind and drainage also have an effect on the flow field of reservoirs, and they should be taken into consideration in the future. Besides, all the numerical simulations in the present work are performed on uniform mesh, and the practical applications require a large number of grid nodes to reach acceptable accuracy. Therefore, in order to enhance the computing efficiency, other works in the future will focus on local mesh refinement to this coupling scheme and realize massive parallel computing.

Data Availability

Some or all data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

Wei Diao contributed to methodology, investigation, writing of the original draft; Hao Yuan contributed to resources and validation; Cunze Zhang contributed to funding acquisition; Liang Chen contributed to software; Xujin Zhang contributed to supervision and writing, review, and editing.

Acknowledgments

The authors are grateful for the financial support from the Science and Technology Research Program of Chongqing Municipal Education Commission, Grant no. KJQN201800712.

References

B. Zhang, “Review of new China's hydropower development in 70 years—to celebrate the 70th anniversary of the founding of the People's Republic of China,” Hydropower and Pumped Storage, vol. 005, no. 5, pp. 1–6, 2019, (in Chinese).
View at: Google Scholar
S. Y. Chi, J. X. Hu, M. Li, and J. X. Zheng, “Developing a baseline and tools for the future assessment of the ecological status of rivers within the Three Gorges Reservoir catchment, China,” Hydrobiologia, vol. 775, no. 1, pp. 185–196, 2016.
View at: Publisher Site | Google Scholar
W. X. Duan, G. W. Zhu, J. J. Liu, W. B. Yang, and P. C. Shi, “An evaluation method for ecological security of water resource reservoirs,” China Environmental Science, vol. 40, no. 9, pp. 4135–4145, 2020, (in Chinese).
View at: Google Scholar
Y. Xiong, “Coupling sediment transport and water quality models,” Mississippi State University, MS, USA, 2010, Ph.D. thesis.
View at: Google Scholar
S. Zhang, Z. Xu, and K. Kong, “Simulation and analysis on flow and pollutant of jinan queshan reservoir,” Advances in Computer Science, Intelligent System and Environment, vol. 105, no. 2, pp. 693–699, 2011.
View at: Publisher Site | Google Scholar
W. Li, J. H. Liu, Y. F. Zhang, and W. H. Xiao, “Research on reservoir water temperature simulation and fish response,” Sustainable Development of Water Resources and Hydraulic Engineering in China, vol. 42, pp. 327–339, 2019.
View at: Publisher Site | Google Scholar
M. Elshemy, M. Khadr, Y. Atta, and A. Ahmed, “Hydrodynamic and water quality modeling of lake manzala (Egypt) under data scarcity,” Environmental Earth Sciences, vol. 75, no. 19, 2016.
View at: Publisher Site | Google Scholar
T. Sun and J. H. Tao, “Numerical simulation of pollutant transport acted by wave for a shallow water sea bay,” International Journal for Numerical Methods in Fluids, vol. 51, no. 5, pp. 469–487, 2006.
View at: Google Scholar
S. Prayitno, R. A. Santoso, Deendarlianto, H. Thomas, and D. Lucas, “Counter current flow limitation of gas-liquid two-phase flow in nearly horizontal pipe,” Science and Technology and Nuclear Installations, vol. 2012, Article ID 513809, 9 pages, 2012.
View at: Publisher Site | Google Scholar
Y. Lang, L. Wang, X. Cai, and Z. Hu, “Research on hydrodynamics with water temperature characteristics and spring algal blooms in a typical tributary bay of three gorges reservoir,” Mathematical Problems in Engineering, vol. 2019, Article ID 7654543, 10 pages, 2019.
View at: Publisher Site | Google Scholar
Y. Dai, Z. Zhou, J. Lin, and J. Han, “Modeling of two-phase flow in rough-walled fracture using level set method,” Geofluids, vol. 2017, pp. 1–11, 2017.
View at: Publisher Site | Google Scholar
J. Shen and X. Yang, “A phase-field model and its numerical approximation for two-phase incompressible flows with different densities and viscosities,” SIAM Journal on Scientific Computing, vol. 32, no. 3, pp. 1159–1179, 2010.
View at: Publisher Site | Google Scholar
M. Dai and D. P. Schmidt, “Adaptive tetrahedral meshing in free-surface flow,” Journal of Computational Physics, vol. 208, no. 1, pp. 228–252, 2005.
View at: Publisher Site | Google Scholar
F. Harlow and J. Welch, “Numerical calculations of time-dependent viscous incompressible flow of fluid with free surfaces,” Physics of Fluids, vol. 10, pp. 2182–2189, 1967.
View at: Google Scholar
D. M. Kelly, Q. Chen, and J. Zang, “Picin: a particle-in-cell solver for incompressible free surface flows with two-way fluid-solid coupling,” SIAM Journal on Scientific Computing, vol. 37, no. 3, pp. B403–B424, 2015.
View at: Publisher Site | Google Scholar
X. H. Zhang, B. D. Li, W. J. Wei, J. Wang, and X. X. Zhang, “Numerical Study on the 3-D water temperature in the reservoir region adjacent to a dam,” Journal of Sichuan University (Engineering Science Edition), vol. 48, no. 2, pp. 14–20, 2016, (in Chinese).
View at: Google Scholar
N. Thürey, “A single-phase free-surface lattice boltzmann method,” Friedrich-Alexander-Universität of Erlangen-Nürnberg, Erlangen, Germany, 2003, Master's Thesis.
View at: Google Scholar
C. Körner, M. Thies, T. N. Hofmann, N. Thürey, and U. Rüde, “Lattice Boltzmann model for free surface flow for modeling foaming,” Journal of Statistical Physics, vol. 121, no. 1-2, pp. 179–196, 2005.
View at: Publisher Site | Google Scholar
Z. Guo, B. Shi, and C. Zheng, “A coupled lattice BGK model for the Boussinesq equations,” International Journal for Numerical Methods in Fluids, vol. 39, no. 4, pp. 325–342, 2002.
View at: Publisher Site | Google Scholar
X. Shan, X. F. Yuan, and H. Chen, “Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation,” Journal of Fluid Mechanics, vol. 550, no. 1, pp. 413–441, 2006.
View at: Publisher Site | Google Scholar
F. Verhaeghe, B. Blanpain, and P. Wollants, “Lattice Boltzmann method for double-diffusive natural convection,” Physical review. E, Statistical, nonlinear, and soft matter physics, vol. 75, no. 4, Article ID 046705, 2007.
View at: Publisher Site | Google Scholar
C. Janßen and M. Krafczyk, “Free surface flow simulations on GPGPUs using the LBM,” Computers & Mathematics with Applications, vol. 61, no. 12, pp. 3549–3563, 2011.
View at: Publisher Site | Google Scholar
N. Thürey and U. Rüde, “Stable free surface flows with the lattice Boltzmann method on adaptively coarsened grids,” Computing and Visualization in Science, vol. 12, no. 5, pp. 247–263, 2009.
View at: Publisher Site | Google Scholar
G. Z. Li, L. I. Qing, and Z. C. Guang, “Lattice Boltzmann simulation of double diffusive natural convection,” Chinese Journal of Computational Physics, vol. 19, no. 6, pp. 483–487, 2002, (in Chinese).
View at: Google Scholar
Y. G. Cheng and J. P. Li, “Introducing unsteady non-uniform source terms into the lattice Boltzmann model,” International Journal for Numerical Methods in Fluids, vol. 56, no. 6, pp. 629–641, 2007.
View at: Google Scholar
S. Hou, J. Sterling, S. Chen, and G. D. Doolen, “A lattice Boltzmann subgrid model for high Reynolds number flows,” Fields Institute Communications, vol. 6, no. 13, 1994.
View at: Google Scholar
G. Mayer, J. Pales, and G. Hazi, “Large eddy simulation of subchannels using the lattice Boltzmann method,” Annals of Nuclear Energy, vol. 34, no. 1, pp. 140–149, 2007.
View at: Publisher Site | Google Scholar
S. Ostrach, “Natural convection with combined driving forces,” Physicochemical Hydrodynamics, vol. 1, pp. 233–247, 1980.
View at: Google Scholar
C. Béghein, F. Haghighat, and F. Allard, “Numerical study of double-diffusive natural convection in a square cavity,” International Journal of Heat and Mass Transfer, vol. 35, no. 4, pp. 833–846, 1992.
View at: Publisher Site | Google Scholar
E. M. Buchak and J. E. Edinger, User guide for LARM2:A longitudinal-vertical time-varying hydrodynamic reservoir model, U.S.Amry Corps of Engineerings, Waterways Experiment Station, Vicksburg, MS, USA, 1982.
B. H. Johnson, “A review of multidimensional reservoir hydrodynamic modeling,” in Proceedings of the Symposiyum on Surface Water Impoundments, A4SCE, pp. 497–507, Minneapolis, Minnesota, USA, June 1980.
View at: Google Scholar
X. Tang, “Reservoir water temperature and quality prediction based on cfd software fluent,” Wuhan University, Wuhan, China, 2009, Master's Thesis.
View at: Google Scholar
J. D. Yang and X. T. Yang, Physical model test for upper and lower storage reservoirs of Shenzhen pumped-storage power station, State Key Laboratory of Water Resource and Hydropower Engineering Science, Wuhan University, Wuhan, China, 2005, (in Chinese).

Copyright

Copyright © 2021 Wei Diao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

206

Downloads

765

Citations