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

Authors: Jin Woo Lee; He-Rin Cho; Yong-Sik Cho
Journal: Abstract and Applied Analysis
 (2015)

---

## 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 inx- 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.


---

## Body

## 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 inx-axis and y-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 Figure1.Figure 1
A definitive sketch of the still water level and the free surface elevation.


## 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:(1)∂u∂x+∂v∂y+∂w∂z=0,(2)∂u∂t+∂uu∂x+∂uv∂y+∂uw∂z=-g+Δgz=h∂h∂x-∫zh∂Δg∂xdγ-1ρ0∂q∂x+∂∂xνT∂u∂x+∂∂yνT∂u∂y+∂∂zνT∂u∂z,(3)∂v∂t+∂uv∂x+∂vv∂y+∂vw∂z=-g+Δgz=h∂h∂y-∫zh∂Δg∂ydγ-1ρ0∂q∂y+∂∂xνT∂v∂x+∂∂yνT∂v∂y+∂∂zνT∂v∂z,(4)∂w∂t+∂uw∂x+∂vw∂y+∂ww∂z=-1ρ0∂q∂z+∂∂xνT∂w∂x+∂∂yνT∂w∂y+∂∂zνT∂w∂z,where x, y, and z are the Cartesian x-, y-, and z-coordinates and ρ0 is the reference density; Δρ is the deviation component from the reference density; u, v, and w are velocities in the x-, y-, and z-directions; t is the time; g is gravitational acceleration; νT is the kinematic viscosity.The kinematic boundary conditions at the bottomzb and free surface h, leading to(5)wz=zb=u∂zb∂x+v∂zb∂yz=zb,wz=h=∂h∂t+u∂h∂x+v∂h∂yz=h.For a free surface problem, the continuity equation (1) is integrated from the bottom to the surface leading to the free surface evolution equation:(6)∂h∂t+∂u¯∂x+∂v¯∂y=0.Transport of passive scalars such as temperature, salinity, and passive constituent is governed by the advection-diffusion equation [6]:(7)∂T∂t+∂uT∂x+∂vT∂y+∂wT∂z=Ts+∂∂xΓ∂T∂x+∂∂yΓ∂T∂y+∂∂zΓ∂T∂z,where T is the Reynolds averaged passive scalar and Ts 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]:(8)ρTw,S=+999.842597+6.793952×10-2×Tw-9.095290×10-3×Tw2+1.001685×10-2×Tw3-1.120083×10-3×Tw4+6.536332×10-2×Tw5+8.24493×10-3×S-4.0899×10-2×Tw×S-7.6438×10-3×Tw4×S-8.2467×10-2×Tw3×S-5.3875×10-3×Tw4×S+5.72466×10-2×S3/2-1.0227×10-3×Tw4×S,where Tw is the temperature of water and S 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 (Figure2). In the staggered grid system, velocities (ui+1/2,j,k,  vi,j+1/2,k, and wi,j,k+1/2) are assigned at half intervals and scalar quantities (q, h, and T) are defined at grid points where 1≤i≤N1, 1≤j≤N2, and 1≤k≤N3 are the grid points and N1, N2, and N3 are the number of the maximum grid points with the spatial increments Δx, Δy, and Δz in the x-, y-, and z-axes directions.Figure 2
A definitive sketch of a staggered grid system.Using this staggered grid scale oscillations for the hydrodynamic pressure or the free surface elevation [8] can be removed effectively.For the vertical systemz-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 x-axis directional momentum equation is discretized as(9)ui+1/2,j,k∗-ui+1/2,j,knΔt+ui+1/2,j,knui+1,j,kTVD3-ui,j,kTVD3nΔx+vi+1/2,j,knui+1/2,j+1/2,kTVD3-ui+1/2,j-1/2,kTVD3nΔy+wi+1/2,j,knui+1/2,j,k+1/2TVD3-ui+1/2,j,k-1/2TVD3nΔz=-Lh,u+LΔg,u+Lq,ui+1/2,j,kn+Du,xn+Du,yn+Du,zn+1/3i+1/2,j,k,(10a)Lh,ui+1/2,j,k=g+Δgi+1/2,j,N3hi+1,j-hi,jΔx,(10b)LΔg,ui+1/2,j,k=∫zh∂Δg∂x+zx∂Δg∂zi+1/2,j,kdγ,(10c)Lq,ui+1/2,j,k=1ρ0qi+1,j,k-qi,j,kΔx+zxi+1/2,j,k·qi,j,k+1-qi,j,k-1+qi+1,j,k+1-qi+1,j,k-14Δz.In the numerical scheme, the advection terms are explicitly discretized by the upwind scheme. The diffusion operatorD(a,b) indicates the second derivative term of a with respect to the b-direction. For instance, applying central difference method, x-axis directional diffusion term can be written as(11)Du,xi+1/2,j,k=∂∂xνT∂u∂xi+1/2,j,k=1Δx2νTi+1,j,kui+3/2,j,k-ui+1/2,j,k·-νTi,j,kui+1/2,j,k-ui-1/2,j,k.The y-axis and z-axis directional diffusion terms are similarly written as given in (11).Substituting all diffusion terms into (9), the x-axis directional momentum equation can be represented as the following tridiagonal system:(12)AAui+1/2,j,knui+1/2,j,k-1∗+BBui+1/2,j,knui+1/2,j,k∗+CCui+1/2,j,knui+1/2,j,k+1∗=DDui+1/2,j,kn,where (13a)AAi+1/2,j,ku=-ΔtΔz2cvui+1/2,j,k-1/2n,(13b)BBi+1/2,j,ku=1+ΔtΔz2cvui+1/2,j,k-1/2n+cvui+1/2,j,k+1/2n,(13c)CCi+1/2,j,ku=-ΔtΔz2cvui+1/2,j,k+1/2n,(13d)DDi+1/2,j,ku=ui+1/2,j,kn-ui+1/2,j,knΔtui+1,j,kTVD3-ui,j,kTVD3nΔx-vi+1/2,j,knΔtui+1/2,j+1/2,kTVD3-ui+1/2,j-1/2,kTVD3nΔy-Wi+1/2,j,knΔtui+1/2,j,k+1/2TVD3-ui+1/2,j,k-1/2TVD3nΔz-ΔtLh,u+LΔg,u+Lq,ui+1/2,j,kn+ΔtDu,x+Du,y+HOTui+1/2,j,kn.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. They-axis and z-axis directional momentum equation is similarly written as the x-axis directional equation (12).


### 3.2. Free Surface Correction Step

In this step, the velocity fields (u∗, v∗, and w∗) 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:(14a)ui+1/2,j,k∗∗=ui+1/2,j,k∗-θΔtg∂∂xδh~n+1i+1/2,j(14b)ui,j+1/2,k∗∗=ui,j+1/2,k∗-θΔtg∂∂yδh~n+1i,j+1/2,where δh~n+1=h~n+1-hn and h~ 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(15)δh~i,jn+1+θΔt1Δxu˘i+1/2,j∗∗-u˘i-1/2,j∗∗+1Δyv˘i,j+1/2∗∗-v˘i,j-1/2∗∗=-1-θ·Δt1Δxu˘i+1/2,jn-u˘i-1/2,jn+1Δyv˘i,j+1/2n-v˘i,j-1/2n.Equation (15) is unconditionally stable for gravitational wave propagation when θ≥0.5 [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(16a)u˘i+1/2,j∗∗=u˘i+1/2,j∗-Hi+1/2,jngθΔt∂∂xδh~n+1i+1/2,j,(16b)v˘i,j+1/2∗∗=v˘i,j+1/2∗-Hi,j+1/2ngθΔt∂∂yδh~n+1i,j+1/2.Substituting (16a), (16b) into (15), the unknown velocities u˘∗∗ and v˘∗∗ can be eliminated leading to the free surface correction equation. After some algebraic manipulation, this can be written as(17)δh~i,jn+1-gθΔt2α1+Δx∂∂xδh~n+1i+1/2,j-α1-Δx·∂∂xδh~n+1i-1/2,j+α2+Δy∂∂yδh~n+1i,j+1/2-α2-Δy∂∂yδh~n+1i,j-1/2=-θΔt1Δxu˘i+1/2,j∗-u˘i-1/2,j∗+1Δyv˘i,j+1/2∗-v˘i,j-1/2∗-1-θ·Δt1Δxu˘i+1/2,jn-u˘i-1/2,jn+1Δyv˘i,j+1/2n-v˘i,j-1/2n,where α1±=Hi±1/2,jn and α2±=Hi,j±1/2n. 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 i=1/2 is a solid wall, the normal velocity at the wall is set to be zero: that is, u˘1/2,j, 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, wi,j,k+1/2∗ is transferred into wi,j,k+1/2∗∗ 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 as(18a)ui+1/2,j,kn+1=ui+1/2,j,k∗∗-Δtρ0∂∂xδq~n+1+zxn∂∂zδq~n+1i+1/2,j,k,(18b)vi,j+1/2,kn+1=vi,j+1/2,k∗∗-Δtρ0∂∂yδq~n+1+zyn∂∂yδq~n+1i,j+1/2,k,(18c)wi,j,k+1/2n+1=wi,j,k+1/2∗∗-Δtρ0zzn∂∂zδq~n+1i,j,k+1/2,where δq~n+1=q~n+1-qn and q~n+1 is the provisional hydrodynamic pressure. With the provisional free surface elevation, this hydrodynamic pressure can be related to the total pressure as(19)pi,j,kn+1ρ0=gh~i,jn+1-zi,j,k+qi,j,knρ0+q~i,j,kn+1ρ0or with the new free surface elevation [11],(20)pi,j,kn+1ρ0=ghi,jn+1-zi,j,k+qi,j,knρ0.In transformed form, the vertical fraction equation can be written as(21)W^i,j,k+1/2n+1=W^i,j,k+1/2∗∗-Δtρ0zxn∂∂xδq~n+1+zyn∂∂yδq~n+1+zxzx+zyzy+zzzzn∂∂zδq~n+1i,j,k+1/2.If the continuity equation is discretized using a central difference method as(22)1Δxun+1Jni+1/2,j,k-un+1Jni-1/2,j,k+1Δyvn+1Jni,j+1/2,k-vn+1Jni,j-1/2,k+1ΔzW^n+1Jni,j,k+1/2-W^n+1Jni,j,k-1/2=0by 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, q=δq=0. At the bottom, Wi,j,1/2=0 can be applied which is equivalent to W^i,j,1/2=0 because (zt)i,j,1/2=0. 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, δq~i,j,0=2δq~i,j,1-δq~i,j,2.Assuming the hydrostatic pressure distribution at the top free surface layer (qi,j,N3+1/2n=qi,j,N3n+1=0), (19) may be discretized with a new free surface elevation hn+1 as (23)hi,jn+1=h~i,jn+1+δq~i,j,N3n+1ρ0g.A new hydrodynamic pressure is then obtained by combining the right hand sides of (19), (20), and (23) leading to(24)qi,j,kn+1=δq~i,j,kn+1-δq~i,j,N3n+1.


### 3.4. Scalar Transport Equation

The scalar transport equation (7) is discretized as(25)1ΔtTi,j,kn+1-Ti,j,kn+1Δxui,j,kn+1Ti+1/2,j,kTVD3-ui,j,kn+1Ti-1/2,j,kTVD3+1Δyvi,j,kn+1Ti,j+1/2,kTVD3-vi,j,kn+1Ti,j-1/2,kTVD3+1Δzwi,j,kn+1Ti,j,k+1/2TVD3-wi,j,kn+1Ti,j,k-1/2TVD3=DT,xn+DT,yn+DT,zn.The discretization of advection terms is considered to be explicit. As cr>1 case, the several steps, which were divided calculating time (Δt), can persist in the condition of cr≤1. In this study, the explicit method is absolutely satisfied with the stability for cr≤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(26)DT,xi+1/2,j,k=∂∂xνh∂T∂xi+1/2,j,k=1Δx2νhi+1,j,kTi+3/2,j,k-Ti+1/2,j,k-νhi,j,kTi+1/2,j,k-Ti-1/2,j,k.The y-axis and z-axis directional diffusive flux is similarly written as shown in (26). The transport equation can be represented as the following tridiagonal system:(27)AATi,j,knTi,j,k-1n+1+BBTi,j,knTi,j,kn+1+CCTi,j,knTi,j,k+1n+1=DDTi,j,kn,where(28)AAi,j,kT=-ΔtΔz2cvTi,j,k-1/2n,BBi,j,kT=1+ΔtΔz2cvTi,j,k-1/2n+cvTi,j,k+1/2n,CCi,j,kT=-ΔtΔz2cvTi,j,k+1/2n,DDi,j,kT=Ti,j,kn-ui,j,kn+1ΔtTi+1/2,j,kTVD3-Ti-1/2,j,kTVD3nΔx-vi,j,kn+1ΔtTi,j+1/2,kTVD3-Ti,j-1/2,kTVD3nΔy-Wi,j,kn+1Δtwi,j,k+1/2TVD3-wi,j,k-1/2TVD3nΔz+ΔtTsi,j,kn+ΔtDT,z+DT,yi,j,kn,cvTi,j,k±1/2n=Γi,j,k±1/2Δzi,j,k±1/2.


## 3.1. Hydrostatic Pressure Step

In this step, numerical gird allocation uses a staggered grid system (Figure2). In the staggered grid system, velocities (ui+1/2,j,k,  vi,j+1/2,k, and wi,j,k+1/2) are assigned at half intervals and scalar quantities (q, h, and T) are defined at grid points where 1≤i≤N1, 1≤j≤N2, and 1≤k≤N3 are the grid points and N1, N2, and N3 are the number of the maximum grid points with the spatial increments Δx, Δy, and Δz in the x-, y-, and z-axes directions.Figure 2
A definitive sketch of a staggered grid system.Using this staggered grid scale oscillations for the hydrodynamic pressure or the free surface elevation [8] can be removed effectively.For the vertical systemz-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 x-axis directional momentum equation is discretized as(9)ui+1/2,j,k∗-ui+1/2,j,knΔt+ui+1/2,j,knui+1,j,kTVD3-ui,j,kTVD3nΔx+vi+1/2,j,knui+1/2,j+1/2,kTVD3-ui+1/2,j-1/2,kTVD3nΔy+wi+1/2,j,knui+1/2,j,k+1/2TVD3-ui+1/2,j,k-1/2TVD3nΔz=-Lh,u+LΔg,u+Lq,ui+1/2,j,kn+Du,xn+Du,yn+Du,zn+1/3i+1/2,j,k,(10a)Lh,ui+1/2,j,k=g+Δgi+1/2,j,N3hi+1,j-hi,jΔx,(10b)LΔg,ui+1/2,j,k=∫zh∂Δg∂x+zx∂Δg∂zi+1/2,j,kdγ,(10c)Lq,ui+1/2,j,k=1ρ0qi+1,j,k-qi,j,kΔx+zxi+1/2,j,k·qi,j,k+1-qi,j,k-1+qi+1,j,k+1-qi+1,j,k-14Δz.In the numerical scheme, the advection terms are explicitly discretized by the upwind scheme. The diffusion operatorD(a,b) indicates the second derivative term of a with respect to the b-direction. For instance, applying central difference method, x-axis directional diffusion term can be written as(11)Du,xi+1/2,j,k=∂∂xνT∂u∂xi+1/2,j,k=1Δx2νTi+1,j,kui+3/2,j,k-ui+1/2,j,k·-νTi,j,kui+1/2,j,k-ui-1/2,j,k.The y-axis and z-axis directional diffusion terms are similarly written as given in (11).Substituting all diffusion terms into (9), the x-axis directional momentum equation can be represented as the following tridiagonal system:(12)AAui+1/2,j,knui+1/2,j,k-1∗+BBui+1/2,j,knui+1/2,j,k∗+CCui+1/2,j,knui+1/2,j,k+1∗=DDui+1/2,j,kn,where (13a)AAi+1/2,j,ku=-ΔtΔz2cvui+1/2,j,k-1/2n,(13b)BBi+1/2,j,ku=1+ΔtΔz2cvui+1/2,j,k-1/2n+cvui+1/2,j,k+1/2n,(13c)CCi+1/2,j,ku=-ΔtΔz2cvui+1/2,j,k+1/2n,(13d)DDi+1/2,j,ku=ui+1/2,j,kn-ui+1/2,j,knΔtui+1,j,kTVD3-ui,j,kTVD3nΔx-vi+1/2,j,knΔtui+1/2,j+1/2,kTVD3-ui+1/2,j-1/2,kTVD3nΔy-Wi+1/2,j,knΔtui+1/2,j,k+1/2TVD3-ui+1/2,j,k-1/2TVD3nΔz-ΔtLh,u+LΔg,u+Lq,ui+1/2,j,kn+ΔtDu,x+Du,y+HOTui+1/2,j,kn.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. They-axis and z-axis directional momentum equation is similarly written as the x-axis directional equation (12).


## 3.2. Free Surface Correction Step

In this step, the velocity fields (u∗, v∗, and w∗) 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:(14a)ui+1/2,j,k∗∗=ui+1/2,j,k∗-θΔtg∂∂xδh~n+1i+1/2,j(14b)ui,j+1/2,k∗∗=ui,j+1/2,k∗-θΔtg∂∂yδh~n+1i,j+1/2,where δh~n+1=h~n+1-hn and h~ 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(15)δh~i,jn+1+θΔt1Δxu˘i+1/2,j∗∗-u˘i-1/2,j∗∗+1Δyv˘i,j+1/2∗∗-v˘i,j-1/2∗∗=-1-θ·Δt1Δxu˘i+1/2,jn-u˘i-1/2,jn+1Δyv˘i,j+1/2n-v˘i,j-1/2n.Equation (15) is unconditionally stable for gravitational wave propagation when θ≥0.5 [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(16a)u˘i+1/2,j∗∗=u˘i+1/2,j∗-Hi+1/2,jngθΔt∂∂xδh~n+1i+1/2,j,(16b)v˘i,j+1/2∗∗=v˘i,j+1/2∗-Hi,j+1/2ngθΔt∂∂yδh~n+1i,j+1/2.Substituting (16a), (16b) into (15), the unknown velocities u˘∗∗ and v˘∗∗ can be eliminated leading to the free surface correction equation. After some algebraic manipulation, this can be written as(17)δh~i,jn+1-gθΔt2α1+Δx∂∂xδh~n+1i+1/2,j-α1-Δx·∂∂xδh~n+1i-1/2,j+α2+Δy∂∂yδh~n+1i,j+1/2-α2-Δy∂∂yδh~n+1i,j-1/2=-θΔt1Δxu˘i+1/2,j∗-u˘i-1/2,j∗+1Δyv˘i,j+1/2∗-v˘i,j-1/2∗-1-θ·Δt1Δxu˘i+1/2,jn-u˘i-1/2,jn+1Δyv˘i,j+1/2n-v˘i,j-1/2n,where α1±=Hi±1/2,jn and α2±=Hi,j±1/2n. 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 i=1/2 is a solid wall, the normal velocity at the wall is set to be zero: that is, u˘1/2,j, 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, wi,j,k+1/2∗ is transferred into wi,j,k+1/2∗∗ 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 as(18a)ui+1/2,j,kn+1=ui+1/2,j,k∗∗-Δtρ0∂∂xδq~n+1+zxn∂∂zδq~n+1i+1/2,j,k,(18b)vi,j+1/2,kn+1=vi,j+1/2,k∗∗-Δtρ0∂∂yδq~n+1+zyn∂∂yδq~n+1i,j+1/2,k,(18c)wi,j,k+1/2n+1=wi,j,k+1/2∗∗-Δtρ0zzn∂∂zδq~n+1i,j,k+1/2,where δq~n+1=q~n+1-qn and q~n+1 is the provisional hydrodynamic pressure. With the provisional free surface elevation, this hydrodynamic pressure can be related to the total pressure as(19)pi,j,kn+1ρ0=gh~i,jn+1-zi,j,k+qi,j,knρ0+q~i,j,kn+1ρ0or with the new free surface elevation [11],(20)pi,j,kn+1ρ0=ghi,jn+1-zi,j,k+qi,j,knρ0.In transformed form, the vertical fraction equation can be written as(21)W^i,j,k+1/2n+1=W^i,j,k+1/2∗∗-Δtρ0zxn∂∂xδq~n+1+zyn∂∂yδq~n+1+zxzx+zyzy+zzzzn∂∂zδq~n+1i,j,k+1/2.If the continuity equation is discretized using a central difference method as(22)1Δxun+1Jni+1/2,j,k-un+1Jni-1/2,j,k+1Δyvn+1Jni,j+1/2,k-vn+1Jni,j-1/2,k+1ΔzW^n+1Jni,j,k+1/2-W^n+1Jni,j,k-1/2=0by 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, q=δq=0. At the bottom, Wi,j,1/2=0 can be applied which is equivalent to W^i,j,1/2=0 because (zt)i,j,1/2=0. 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, δq~i,j,0=2δq~i,j,1-δq~i,j,2.Assuming the hydrostatic pressure distribution at the top free surface layer (qi,j,N3+1/2n=qi,j,N3n+1=0), (19) may be discretized with a new free surface elevation hn+1 as (23)hi,jn+1=h~i,jn+1+δq~i,j,N3n+1ρ0g.A new hydrodynamic pressure is then obtained by combining the right hand sides of (19), (20), and (23) leading to(24)qi,j,kn+1=δq~i,j,kn+1-δq~i,j,N3n+1.


## 3.4. Scalar Transport Equation

The scalar transport equation (7) is discretized as(25)1ΔtTi,j,kn+1-Ti,j,kn+1Δxui,j,kn+1Ti+1/2,j,kTVD3-ui,j,kn+1Ti-1/2,j,kTVD3+1Δyvi,j,kn+1Ti,j+1/2,kTVD3-vi,j,kn+1Ti,j-1/2,kTVD3+1Δzwi,j,kn+1Ti,j,k+1/2TVD3-wi,j,kn+1Ti,j,k-1/2TVD3=DT,xn+DT,yn+DT,zn.The discretization of advection terms is considered to be explicit. As cr>1 case, the several steps, which were divided calculating time (Δt), can persist in the condition of cr≤1. In this study, the explicit method is absolutely satisfied with the stability for cr≤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(26)DT,xi+1/2,j,k=∂∂xνh∂T∂xi+1/2,j,k=1Δx2νhi+1,j,kTi+3/2,j,k-Ti+1/2,j,k-νhi,j,kTi+1/2,j,k-Ti-1/2,j,k.The y-axis and z-axis directional diffusive flux is similarly written as shown in (26). The transport equation can be represented as the following tridiagonal system:(27)AATi,j,knTi,j,k-1n+1+BBTi,j,knTi,j,kn+1+CCTi,j,knTi,j,k+1n+1=DDTi,j,kn,where(28)AAi,j,kT=-ΔtΔz2cvTi,j,k-1/2n,BBi,j,kT=1+ΔtΔz2cvTi,j,k-1/2n+cvTi,j,k+1/2n,CCi,j,kT=-ΔtΔz2cvTi,j,k+1/2n,DDi,j,kT=Ti,j,kn-ui,j,kn+1ΔtTi+1/2,j,kTVD3-Ti-1/2,j,kTVD3nΔx-vi,j,kn+1ΔtTi,j+1/2,kTVD3-Ti,j-1/2,kTVD3nΔy-Wi,j,kn+1Δtwi,j,k+1/2TVD3-wi,j,k-1/2TVD3nΔz+ΔtTsi,j,kn+ΔtDT,z+DT,yi,j,kn,cvTi,j,k±1/2n=Γi,j,k±1/2Δzi,j,k±1/2.


## 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 depthH=h+ζ to the wave length L. In this case, an inviscid fluid of constant density is confined in a closed basin with a square closed basin of length l=10 m and with a still water depth of H0=10 m. A zero initial velocity is assumed and the initial free surface elevation is given by(29)ζx=acos⁡kx+H0with0≤xl,where a is the amplitude and the wave number is k=2π/L. 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: (30)c=gL2πtanh⁡2πLH01/2,where c is the wave phase speed. By neglecting bottom friction and horizontal and vertical viscosity, the calculation is carried out with a time step Δt=0.01 sec.The computational domain uses a constant grid spacing ofΔx=0.25 m in the horizontal directions and Δz=0.5 m in the vertical direction. A time interval of Δt=0.01 sec is used to obtain solutions. A small time step of Δt=0.01 sec is used for higher accuracy. The expected solution consists of a standing wave of length L=2l and frequency of wave f=c/L, where c 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 t=T/8, T/2, 5T/8, and 2T sec. As seen from the results, there is a very good agreement between the numerical results and the analytical solution for short period waves.Figure 3
Comparison of free surface elevations obtained from the numerical model and analytical solutions fort=T/8, T/2, 5T/8, and 2T with T being the period.
(a)
t
=
T
/
8
(b)
t
=
T
/
2
(c)
t
=
5
T
/
8
(d)
t
=
2
TFigures4 and 5 show the numerical predictions of the free surface elevations with hydrostatic and hydrodynamic pressure approximation at the boundaries x=0 m and x=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 T=2.02 sec, while for the inclusion of the hydrodynamic pressure, the period is T=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.Figure 4
Comparison of analytical solutions (circle) and numerical predictions (dashed line: hydrostatic pressure model; solid line: hydrodynamic pressure model) for the free surface elevation at boundary ofx=0 m.Figure 5
Comparison of analytical solutions (circle) and numerical predictions (dashed line: hydrostatic pressure model; solid line: hydrodynamic pressure model) for the free surface elevation at boundary ofx=10 m.


### 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 thex-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 d:(31)T=27.5,d<2.027.5-4.02.0d-2.0,2≤d≤422.5,d>4.0.The computational domain employs a constant grid spacing system ofΔx=Δy= 5 m in the horizontal directions, with 20 equally divided layers in the z-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 z-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, ρa= 1.2 kg/m3, and cf = 1.5 × 10−3 [6] are assumed. The calculated time is a total 3 hr with the time step of 2.5 sec.Figure6 shows that the wind-driven circulation is caused in thermal stratified flow at t=0.3506 hr, 0.5 hr, 1.0 hr, 1.5 hr, and 2.0 hr.Figure 6
Wind-driven circulation for thermal stratified flow in a closed basin aty=25m.
(a)
t
= 0.3506 hr
(b)
t
= 0.5 hr
(c)
t
= 1.0 hr
(d)
t
= 1.5 hr
(e)
t
= 2.0 hrWhile the wind is blowing, Figures6(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 Figures6(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.


## 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 depthH=h+ζ to the wave length L. In this case, an inviscid fluid of constant density is confined in a closed basin with a square closed basin of length l=10 m and with a still water depth of H0=10 m. A zero initial velocity is assumed and the initial free surface elevation is given by(29)ζx=acos⁡kx+H0with0≤xl,where a is the amplitude and the wave number is k=2π/L. 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: (30)c=gL2πtanh⁡2πLH01/2,where c is the wave phase speed. By neglecting bottom friction and horizontal and vertical viscosity, the calculation is carried out with a time step Δt=0.01 sec.The computational domain uses a constant grid spacing ofΔx=0.25 m in the horizontal directions and Δz=0.5 m in the vertical direction. A time interval of Δt=0.01 sec is used to obtain solutions. A small time step of Δt=0.01 sec is used for higher accuracy. The expected solution consists of a standing wave of length L=2l and frequency of wave f=c/L, where c 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 t=T/8, T/2, 5T/8, and 2T sec. As seen from the results, there is a very good agreement between the numerical results and the analytical solution for short period waves.Figure 3
Comparison of free surface elevations obtained from the numerical model and analytical solutions fort=T/8, T/2, 5T/8, and 2T with T being the period.
(a)
t
=
T
/
8
(b)
t
=
T
/
2
(c)
t
=
5
T
/
8
(d)
t
=
2
TFigures4 and 5 show the numerical predictions of the free surface elevations with hydrostatic and hydrodynamic pressure approximation at the boundaries x=0 m and x=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 T=2.02 sec, while for the inclusion of the hydrodynamic pressure, the period is T=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.Figure 4
Comparison of analytical solutions (circle) and numerical predictions (dashed line: hydrostatic pressure model; solid line: hydrodynamic pressure model) for the free surface elevation at boundary ofx=0 m.Figure 5
Comparison of analytical solutions (circle) and numerical predictions (dashed line: hydrostatic pressure model; solid line: hydrodynamic pressure model) for the free surface elevation at boundary ofx=10 m.


## 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 thex-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 d:(31)T=27.5,d<2.027.5-4.02.0d-2.0,2≤d≤422.5,d>4.0.The computational domain employs a constant grid spacing system ofΔx=Δy= 5 m in the horizontal directions, with 20 equally divided layers in the z-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 z-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, ρa= 1.2 kg/m3, and cf = 1.5 × 10−3 [6] are assumed. The calculated time is a total 3 hr with the time step of 2.5 sec.Figure6 shows that the wind-driven circulation is caused in thermal stratified flow at t=0.3506 hr, 0.5 hr, 1.0 hr, 1.5 hr, and 2.0 hr.Figure 6
Wind-driven circulation for thermal stratified flow in a closed basin aty=25m.
(a)
t
= 0.3506 hr
(b)
t
= 0.5 hr
(c)
t
= 1.0 hr
(d)
t
= 1.5 hr
(e)
t
= 2.0 hrWhile the wind is blowing, Figures6(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 Figures6(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.



---

*Source: 727315-2015-07-06.xml*