Abstract
The closed water bodies, such as reservoirs and lakes, can be polluted by an inflow of pollutants in the upstream as well as a stratification caused by seasonal natural phenomena. The vertical circulation particularly plays an important role in reducing environmental pollutants. The factors of the vertical circulation are the temperature, wind, thermal diffusivity, sunlight, and so on. The wind is the most significant factor among all possible factors causing the vertical circulation. Thus, it is necessary to describe the validation and application of a three-dimensional numerical model of wind-driven circulation in a thermally stratified flow. In this study, the numerical model is conducted in three steps to calculate the velocity components from the momentum equations in x- and y-directions, the elevations from the free surface equation, and the temperature from the scalar transport equation. The present model was applied to two tests for verification of the numerical accuracy. Numerical results are compared with analytical solutions of the sloshing free surface movement in a rectangular basin and the model is applied to the circulation for the wind-driven flow in a thermal stratification. Consequently, the developed model is validated by two verifications and phenomena of the internal flow.
1. Introduction
The closed water bodies, such as lakes and reservoirs, are major water sources to provide drinking water, power generation, irrigation for paddy fields, and recreation. But they may be polluted by an inflow of pollutants in the upstream as well as a stratification caused by seasonal natural phenomena. To maintain an environmental preservation, the flow pattern due to stratification and the condition of biological characteristics for a temperature must be analyzed. For instance, the stable stratification prevents the circulation from the bottom to the surface causing anaerobic conditions at the bottom, providing favorable conditions for excessive algal growth. Therefore, the vertical circulation particularly plays an important role in the flow to reduce serious environmental problems in closed water bodies. The factors of the vertical circulation are the temperature, wind, thermal diffusive, sunlight, and so on. The wind would be the most significant factor among all possible factors causing the vertical circulation in closed water bodies because the wind-driven flow crucially affects the water quality in the closed water area.
The purpose of this paper is to describe thermal circulation and stratification with a three-dimensional numerical model of wind-driven circulation and thermally stratified flow in closed water bodies. During the last two decades, two-dimensional and three-dimensional numerical models have been extensively refined and used in simulating free surface flows. Using the hydrostatic pressure approximation, the Navier-Stokes equations can be integrated over depth, resulting in the shallow water equations. Although the shallow water equations are easier to solve than the full Navier-Stokes equations, vertical variations of the flow characteristics are not available in the depth-averaged equations. Also, researchers have recognized that the vertical variation of the flow characteristics is very important for environmental problems in the late 1980s and the early 1990s [1], and three-dimensional numerical models have been developed with the hydrostatic pressure approximation but without the depth integration [2, 3]. In order to solve a three-dimensional model, the velocity projection method is used to consider free surface elevation in depth-averaged equation [3, 4]. The projection method is defined where the momentum equations were projected on to the free surface equation, obtaining the predictor-correction of free surface elevations. Also, to resolve baroclinic gradients, many numerical models are introduced by several authors. Huang and Spaulding [5] developed a hydrostatic model and applied it to environmental problems which are described for coastal circulation and pollutant transport.
To analyze the phenomenon of internal flow, many two-dimensional models have been developed to simulate wind-driven currents. Recently, many three-dimensional numerical models have been developed to simulate the wind-driven circulation. In this study, numerical models of closed water bodies for wind-driven circulation in the thermally stratified flow are developed to calculate in three steps the velocity components from the momentum equations in -axis and -axis directions, the free surface elevation, and the scalar transport for the temperature.
The numerical verifications are conducted to simulate the sloshing free surface movement and to compare with analytical solutions in a closed basin. In this study, the numerical model is applied to the circulation driven by the wind in the thermally stratified flow. Numerical simulation will describe the phenomenon of the internal flow and thermal circulation driven by the wind in a domain as shown in Figure 1.
2. Governing Equations
Assuming the hydrostatic pressure and the Boussinesq approximation, the governing equations for the closed water bodies become incompressible Reynolds averaged Navier-Stokes equations. The continuity and momentum equations are expressed in a conservative form as follows:where , , and are the Cartesian -, -, and -coordinates and is the reference density; is the deviation component from the reference density; , , and are velocities in the -, -, and -directions; is the time; is gravitational acceleration; is the kinematic viscosity.
The kinematic boundary conditions at the bottom and free surface , leading to
For a free surface problem, the continuity equation (1) is integrated from the bottom to the surface leading to the free surface evolution equation:Transport of passive scalars such as temperature, salinity, and passive constituent is governed by the advection-diffusion equation [6]:where is the Reynolds averaged passive scalar and is the source or sink term. is the eddy diffusivity. After solving the scalar transport equation (7) for temperature and salinity, the density variation due to these scalars can be considered using the international equation of state [7]:where is the temperature of water and is the salinity.
3. Numerical Model
When the wavelength is relatively small compared to the depth, vertical acceleration is important, so that the effect of the hydrodynamic pressure cannot be ignored. However, in most cases, a direct relation between pressure and density is not available, because of the incompressibility of the fluid. To solve the hydrodynamic pressure, we consider the velocity projection method.
The velocity projection method uses fractional or predictor-corrector steps to solve the hydrodynamic pressure. Initially, the momentum equations are solved using values of hydrodynamic pressure calculated at the previous time step. For this reason, the velocities from the first fractional step (which is called the hydrostatic pressure step, since the effect of the hydrodynamic pressure is not considered fully in this step) may not satisfy the continuity equation. At the subsequent step, or the hydrodynamic pressure correction step, the velocity field from the previous step is updated to ensure local mass conservation.
3.1. Hydrostatic Pressure Step
In this step, numerical gird allocation uses a staggered grid system (Figure 2). In the staggered grid system, velocities (, , and ) are assigned at half intervals and scalar quantities (, , and ) are defined at grid points where , , and are the grid points and , , and are the number of the maximum grid points with the spatial increments , , and in the -, -, and -axes directions.
Using this staggered grid scale oscillations for the hydrodynamic pressure or the free surface elevation [8] can be removed effectively.
For the vertical system -axis direction, an implicit method for the diffusion terms is employed to reduce numerical instability caused by relatively small grid increments in the vertical direction compared to the horizontal directions. For example, the -axis directional momentum equation is discretized as
In the numerical scheme, the advection terms are explicitly discretized by the upwind scheme. The diffusion operator indicates the second derivative term of with respect to the -direction. For instance, applying central difference method, -axis directional diffusion term can be written asThe -axis and -axis directional diffusion terms are similarly written as given in (11).
Substituting all diffusion terms into (9), the -axis directional momentum equation can be represented as the following tridiagonal system:where
It should be noted that this approximation may decrease numerical accuracy when the off diagonal elements are in the same magnitude of the diagonal elements. The -axis and -axis directional momentum equation is similarly written as the -axis directional equation (12).
3.2. Free Surface Correction Step
In this step, the velocity fields (, , and ) obtained from the hydrostatic pressure step are updated by considering free surface movements. According to Chen [4, 9], velocities are updated by considering the following equations:where and is the provisional free surface elevation which is corrected in the next step. Note that combining (14a) with (12) leads to a second-order time discretization for the free surface gradient terms in the horizontal momentum equation.
In order to solve (14a), (14b) for the horizontal velocities and the free surface elevation together, this equation is projected into the free surface equation which is discretized using a semi-implicit method similar to that introduced by Casulli and Cheng [3] asEquation (15) is unconditionally stable for gravitational wave propagation when [10].
For velocity projection in (15), (14a) and (14b) are written in terms of the transformed velocities and then integrated over the depth leading to
Substituting (16a), (16b) into (15), the unknown velocities and can be eliminated leading to the free surface correction equation. After some algebraic manipulation, this can be written aswhere and . In order to solve (17), boundary conditions for the free surface elevation are needed. If the surface elevation is explicitly given at the boundaries, it is implemented directly into (17); otherwise, normal velocity conditions are applied. For example, assuming is a solid wall, the normal velocity at the wall is set to be zero: that is, , which removes most of the free surface elevations outside the boundary from (16a).
The free surface correction equation leads to a symmetric system only when the water depth is constant which is not the case for free surface flows. For this reason, it is solved using the biconjugate stabilized (BICGSTAB) method. When the free surface elevation is obtained by solving (17), the horizontal velocities are updated using (14a), (14b). At this stage, velocities satisfy the global mass conservation represented by the depth averaged continuity (or free surface) equation. Note that the vertical velocity does not affect the global mass conservation because of depth integration so that the vertical momentum equation is independent on the free surface. Therefore, is transferred into without loss of generality.
3.3. Hydrodynamic Pressure Correction Step
As the vertical velocity is not taken into account at the previous step, the velocity field may not satisfy the local mass conservation at each computational grid cell. In this step, the velocities and the free surface elevation are updated to conserve the local mass by considering the effect of hydrodynamic pressure correction. To this end, velocity fractional equations are considered aswhere and is the provisional hydrodynamic pressure. With the provisional free surface elevation, this hydrodynamic pressure can be related to the total pressure asor with the new free surface elevation [11],In transformed form, the vertical fraction equation can be written asIf the continuity equation is discretized using a central difference method as
by substituting (18a), (18b), and (21) into the discrete continuity equation (22) except for the top surface, the discrete Poisson equation can be obtained. Similar to the free surface correction equation, it is solved using the biconjugate stabilized (BI-CGSTAB) method. Similarly to the free surface correction equation, the normal velocity condition is applied at solid horizontal boundary. At open boundaries, the hydrostatic pressure distribution is assumed: that is, . At the bottom, can be applied which is equivalent to because . However, even after applying this normal velocity condition, extra values outside the boundary remained that is caused from grid transformation. These values are extrapolated using values at two inner grid points. For instance, .
Assuming the hydrostatic pressure distribution at the top free surface layer (), (19) may be discretized with a new free surface elevation as A new hydrodynamic pressure is then obtained by combining the right hand sides of (19), (20), and (23) leading to
3.4. Scalar Transport Equation
The scalar transport equation (7) is discretized asThe discretization of advection terms is considered to be explicit. As case, the several steps, which were divided calculating time (, can persist in the condition of . In this study, the explicit method is absolutely satisfied with the stability for . For this reason, the velocity of an internal wave is the less than that of external wave and internal wave is influenced significantly in the flow.
The diffusive flux at a cell face is approximated asThe -axis and -axis directional diffusive flux is similarly written as shown in (26). The transport equation can be represented as the following tridiagonal system:where
4. Verification and Application
4.1. Standing Wave in a Closed Basin
The developed model is applied to free surface problems and verified against analytic solutions and experimental data. The first case deals with nonbreaking waves resulting from a relatively large ratio of total depth to the wave length . In this case, an inviscid fluid of constant density is confined in a closed basin with a square closed basin of length m and with a still water depth of m. A zero initial velocity is assumed and the initial free surface elevation is given bywhere is the amplitude and the wave number is . The amplitude of the wave is taken to be 0.1 m, 1% of the water depth so that small amplitude wave theory applies. Applying the small amplitude wave theory [12], the phase speed of the wave is obtained from the dispersion relation: where is the wave phase speed. By neglecting bottom friction and horizontal and vertical viscosity, the calculation is carried out with a time step sec.
The computational domain uses a constant grid spacing of m in the horizontal directions and m in the vertical direction. A time interval of sec is used to obtain solutions. A small time step of sec is used for higher accuracy. The expected solution consists of a standing wave of length and frequency of wave , where is given by the above dispersion relation. Comparisons of the free surface elevations obtained from the analytical solution and numerical model predictions are given in Figure 3 for the , , , and sec. As seen from the results, there is a very good agreement between the numerical results and the analytical solution for short period waves.
(a)
(b)
(c)
(d)
Figures 4 and 5 show the numerical predictions of the free surface elevations with hydrostatic and hydrodynamic pressure approximation at the boundaries m and m, respectively. A reasonable agreement is observed between the analytical solution and those of the hydrodynamic pressure model. For the case of the hydrostatic pressure assumption, the period of the waves sec, while for the inclusion of the hydrodynamic pressure, the period is sec. Thus, the wave obtained from the hydrostatic pressure model will propagate with a faster speed than that of the hydrodynamic pressure model. This is because the hydrodynamic pressure may play a role like the frequency dispersion.
4.2. Wind-Driven Circulation in a Closed Basin
In this study, the wind-driven circulation for a thermal stratified flow is applied by using the verified eddy viscosity, the scalar transport equation, the internal equation of state, and the baroclinic terms, and then we analyzed phenomenon of the flow circulation.
In the present model, a channel measuring 250 m by 50 m and with a depth of 6 m is considered, with a strong wind blowing over the surface in the -axis direction. The density is a function of the temperature. The density variation on the temperature is assumed to be linear. Using the international equation of states, the density of the fluid is obtained. For initial temperature, a linear distribution is used as a function of depth :
The computational domain employs a constant grid spacing system of 5 m in the horizontal directions, with 20 equally divided layers in the -axis direction. The constant horizontal and vertical eddy viscosities are 0.01 m/s and 0.001 m/sec. At the free surface, a constant wind speed of 10 m/sec in the -axis direction is assumed, and if the blowing time reaches 2,000 sec, the wind should stop. In relation to wind stress according to the empirical formulation, 1.2 kg/m3, and = 1.5 10−3 [6] are assumed. The calculated time is a total 3 hr with the time step of 2.5 sec.
Figure 6 shows that the wind-driven circulation is caused in thermal stratified flow at hr, 0.5 hr, 1.0 hr, 1.5 hr, and 2.0 hr.
(a) 0.3506 hr
(b) 0.5 hr
(c) 1.0 hr
(d) 1.5 hr
(e) 2.0 hr
While the wind is blowing, Figures 6(a) and 6(b) show that the thermocline interface in the thermal stratified flow tends to have the steep slope. The circulation is generated by the wind in the upper and down layers. Upwelling and mixing are observed by the induced-wind external force.
As shown in Figures 6(c) to 6(e), once wind shuts off suddenly, the oscillation of free surface occurred and the slope of thermocline interface will be fluent. After the time passed, the oscillation and an interface of the sloping thermocline will be achieved with original thickness of equilibrium state. Therefore, this indicates that small movements of the free surface produce huge internal waves. The scalar transport equation for the temperature and the numerical approximations are validated and phenomenon of the internal flow is described well.
5. Concluding Remarks
In this study, to maintain an environmental preservation in the closed water bodies, the flow pattern and the conditions of biological characteristics for a temperature must be analyzed. The vertical circulation plays a significant role in the flow to reduce an environmental pollution. The wind is considered to be the most significant factor among all possible factors causing the vertical circulation in the closed water bodies, because of the wind-induced circulation significantly affects water quality in the closed water area.
In the numerical application, when the wind suddenly stops blowing, the present model is applied to the circulation for inducing the wind in the stratified flow. The oscillation of free surface occurred and the thermocline interface has the steep slope. If time of the numerical simulation is large enough, the free surface and thermocline will reach an equilibrium state due to the total energy dissipation by eddy viscosity. The developed model is verified by two numerical tests and phenomena of the internal flow.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was a part of the project titled “Countermeasure System against Typhoons and Tsunamis in Harbor Zones” funded by the Ministry of Oceans and Fisheries, Korea.