Numerical Simulation of Wind-Driven Circulation in a Thermally Stratified Flow

and Applied Analysis 3


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 threedimensional 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 twodimensional models have been developed to simulate winddriven 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.

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 0 is the reference density; Δ is the deviation component from the reference density; , V, 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]: ( , ) = + 999.842597 + 6.793952 × 10 −2 × − 9.095290 × 10 −3 × 2 + 1.001685 where is the temperature of water and is the salinity.

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.
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 * ] . (10c) In the numerical scheme, the advection terms are explicitly discretized by the upwind scheme. The diffusion operator D( , ) indicates the second derivative term of with respect to the -direction. For instance, applying central difference method, -axis directional diffusion term can be written as The -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 +1/2, , Abstract and Applied Analysis 5 +1/2, , 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 andaxis directional momentum equation is similarly written as the -axis directional equation (12).

Free Surface Correction
Step. In this step, the velocity fields ( * , V * , 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: * * whereh +1 =h +1 − ℎ andh 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] as Equation (15) is unconditionally stable for gravitational wave propagation when ≥ 0.5 [10].
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, * , , +1/2 is transferred into * * , , +1/2 without loss of generality.

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 as +1 +1/2, , Assuming the hydrostatic pressure distribution at the top free surface layer ( , , 3 +1/2 = +1 , , 3 = 0), (19) may be discretized with a new free surface elevation ℎ +1 as (23) A new hydrodynamic pressure is then obtained by combining the right hand sides of (19), (20), and (23) leading to

Scalar Transport Equation.
The scalar transport equation (7) is discretized as The discretization of advection terms is considered to be explicit. As | | > 1 case, the several steps, which were divided calculating time (Δ ), can persist in the condition of | | ≤ 1.
In this study, the explicit method is absolutely satisfied with the stability for | | ≤ 1. 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 as Abstract and Applied Analysis 7 The -axis and -axis directional diffusive flux is similarly written as shown in (26). The transport equation can be represented as the following tridiagonal system:

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 = 10 m and with a still water depth of 0 = 10 m. A zero initial velocity is assumed and the initial free surface elevation is given by where is the amplitude and the wave number is = 2 / . 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 Δ = 0.01 sec.
The computational domain uses a constant grid spacing of Δ = 0.25 m in the horizontal directions and Δ = 0.5 m in the vertical direction. A time interval of Δ = 0.01 sec is used to obtain solutions. A small time step of Δ = 0.01 sec is used for higher accuracy. The expected solution consists of a standing wave of length = 2 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 = /8, /2, 5 /8, and 2 sec. As seen from the results, there is a very good agreement between the numerical results and the analytical solution for short period waves. Figures 4 and 5 show the numerical predictions of the free surface elevations with hydrostatic and hydrodynamic pressure approximation at the boundaries = 0 m and = 10 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 = 2.02 sec, while for the inclusion of the hydrodynamic pressure, the period is = 3.59 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.

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/m 3 , 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 = 0.3506 hr, 0.5 hr, 1.0 hr, 1.5 hr, and 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.

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.