Research Article  Open Access
Wei Diao, Hao Yuan, Liang Chen, Xujin Zhang, Cunze Zhang, "ThreeDimensional Numerical Modeling of LargeScale Free Surface Flows with Heat and Contaminant Transfer in Reservoirs", Advances in Civil Engineering, vol. 2021, Article ID 9336868, 17 pages, 2021. https://doi.org/10.1155/2021/9336868
ThreeDimensional Numerical Modeling of LargeScale Free Surface Flows with Heat and Contaminant Transfer in Reservoirs
Abstract
The temperature distribution and pollutant distribution in large reservoirs have always been a hotspot in the field of hydraulics and environmentology, and the threedimensional numerical modeling that can effectively simulate the interactions between the temperature fields, concentration fields, and flow fields needs to be proposed. The doublediffusive convection lattice Boltzmann method is coupled with a singlephase volume of fluid model for simulating heat and contaminant transfer in largescale free surface flows. The coupling model is used to simulate the doublediffusive natural convection in a cubic cavity and the temperature distribution of a model reservoir. The mechanism of convectiondiffusion, 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 upreservoir of a pumpedstorage 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 largescale free surface (LSFS) 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 LSFS flows are mainly twodimensional (2D) numerical models, especially shallow water equations coupling convectiondiffusion equations for simulating planar flow problems [4, 5]. The 2D numerical simulations have been one of the most commonly used methods to analyze largescale 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 threedimensional (3D) convectiondiffusion mechanism shall be simulated by more promising 3D methods. Typically, 3D numerical methods for simulating LSFS flows are mainly based on gasliquid twophase 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 largescale 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].
Singlephase free surface lattice Boltzmann (SPFSLB) model was originally proposed by Thürey in 2003 [17]. Compared with the traditional twophase flow models in which both gas and water are simulated, the SPFSLB model is more powerful to simulate 3D LSFS 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 convectiondiffusion 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 LSFS flows.
To study the 3D temperature and pollutant distribution in large reservoirs, the doublediffusive convection LB method is coupled with the SPFSLB 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 doublediffusive 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 SinglePhase Free Surface LB Model
Unlike traditional twophase flow methods that set up different distribution functions (DFs) for water and air, the SPFSLB model only simulates water flows (the denser and more viscous phase) by the singlephase LB equations and tracks the free surface by a constructed boundary treatment. The following assumptions should be made before the simulation by the SPFSLB 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 largescale 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 singlephase 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, 01, and 1, respectively.
2.1.1. Computation of Mass Flux
The movement of the interface in the SPFSLB 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 SPFSLB 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_{i}n < 0 need to be reconstructed through equation (4) [17]. The vector n refers to the surface normal direction, which is obtained by the secondorder 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. DoubleDiffusive Convection Lattice Boltzmann Method
From a largescale 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 doublediffusive convection lattice Boltzmann (DDCLB) method is used to solve equation (6). The DDCLB 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 DDCLB 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, ω_{0} = 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 largescale flows by the D3Q6 model can reduce calculation amount and save computer memory usage. However, when precisely analyzing the smallscale 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 SPFSLB model is firstly coupled to the DDCLB 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 twoway 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_{0}) + β_{C} (C − C_{0})).
2.3.3. Large Eddy Simulation for Turbulence
Turbulence exists in practical largescale water flows with heat and contaminant transfer. The subgridscale model based on large eddy simulation is introduced to the coupling scheme. The physical variables are separated into largescale variables and smallscale variables by a certain filter function in the subgridscale model [26, 27]. Take φ = + φ′, for example. The largescale variable can be expressed as follows. φ′ presents the smallscale 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} = ν_{1} + ν_{t}, in which ν_{1} and ν_{t} denote the molecular viscosity and eddy viscosity (or turbulence viscosity), respectively. ν_{1} is determined by the physical property of water, and ν_{t} can be computed by , where C, Δ, and refer to Smagorinsky constant, filter width, and strainrate tensor, respectively. According to [27], the strainrate tensor S can be easily computed by using the secondorder 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 doublediffusive 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^{2}, ρ = 10^{3} kg/m^{3}, and ν = 10^{−6} m^{2}/s, respectively.
3.1. DoubleDiffusive Natural Convection in a Cubic Cavity
In the heat and contaminant transfer flows, the doublediffusive natural convection in a cubic cavity reflects the essential problem of twoway 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_{0}/C_{0} and T_{1}/C_{1}, respectively (T_{0} < T_{1}; C_{0} < C_{1}); 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_{0} + T_{1})/2, concentration C = (C_{0} + C_{1})/2, and Pr = 0.71. The flow characteristics of doublediffusive 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_{0}/T_{1} and concentration C_{0}/C_{1} can be assigned to arbitrary numbers that satisfy the formulas above, such as T_{0}/C_{0} = −1 and T_{1}/C_{1} = 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^{4}, Ra_{C} = 10^{4}, 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 convectiondiffusion 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 singlediffusive 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 doublediffusive convection with the temperature Rayleigh Ra_{T} = 10^{7} and concentration Rayleigh number Ra_{C} = 10^{6}, 5 × 10^{6}, 1.5 × 10^{7}, and 5 × 10^{7} 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 lefthand wall, while the outlet is a banded hole that is 0.04 m high and 0.15 m away from the bottom of the righthand 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^{3}/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 tongueshaped uplift at the front of the isothermal surface. See Figures 8(a) and 8(b) for details.(2)Tongueshape 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 tongueshaped 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 tongueshaped 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 tongueshaped 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 heattransfer 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 DDCLB method and SPFSLB model can simulate the largescale thermal flows in model reservoirs effectively.
4. Temporal and Spatial Variations of Temperature and Concentration in a Practical Reservoir
An upper reservoir of a pumpedstorage 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 righthand 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_{0} = 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} = 0°C flows into the reservoir at a flow rate of 640 m^{3}/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} = −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 aboveobtained 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 heattransfer 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 spatialtemporal 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 largescale free surface flows with the heat and contaminant transfer in reservoirs.
5. Conclusion
The SPFSLB model is firstly coupled to the DDCLB method in this paper to simulate largescale 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 doublediffusive 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 upreservoir of a pumpedstorage 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 gasliquid twophase 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 twophase flow in roughwalled fracture using level set method,” Geofluids, vol. 2017, pp. 1–11, 2017. View at: Publisher Site  Google Scholar
 J. Shen and X. Yang, “A phasefield model and its numerical approximation for twophase 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 freesurface 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 timedependent 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 particleincell solver for incompressible free surface flows with twoway fluidsolid 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 3D 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 singlephase freesurface lattice boltzmann method,” FriedrichAlexanderUniversität of ErlangenNü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. 12, 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 NavierStokes 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 doublediffusive 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 nonuniform 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 doublediffusive 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 longitudinalvertical timevarying 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 pumpedstorage 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.