Published online 16 November 2007 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/cnm.1061
A layer-integrated approach for shallow water free-surface
flow computation
Meng-Chi Hung
1,∗,†,‡, Te-Yung Hsieh
2,§, Tung-Lin Tsai
3,¶and Jinn-Chuang Yang
4,1Natural Hazards Mitigation Research Center, National Chiao Tung University, 1001 Ta Hsueh Rd., Hsinchu,
30010 Taiwan, Republic of China
2Energy and Environment Research Laboratories, Industrial Technology Research Institute, Bldg. 64,
195, Sec. 4, Chung-Hsing Rd., Chutung, Hsinchu, 310 Taiwan, Republic of China
3Department of Civil and Water Resources Engineering, National Chiayi University, No. 300 Syuefu Rd.,
Chiayi City, 60004 Taiwan, Republic of China
4Department of Civil Engineering, Natural Hazards Mitigation Research Center, National Chiao Tung University,
1001 Ta Hsueh Rd., Hsinchu, 30010 Taiwan, Republic of China
SUMMARY
A novel semi-three-dimensional (semi-3D) layer-integrated approach was proposed in this study for the shallow water free-surface flow computation. A quadratic shape function proposed to approximate the velocity distribution in the layer ensures the continuity of velocities and shear stresses at interfaces. To discretize the governing equations, the finite volume formulation with staggered grid is used. As the two-dimensional (2D) sub-model for locating the free surface is not needed in this approach, the computational time consumption has been dramatically reduced. The model was verified through the wind-induced circulation in a closed basin by comparing the velocity profiles to the analytical solution. The formation of the velocity profiles due to change of viscosity distribution and the evolution process of water surface slope were also investigated. Further, the eddy viscosity is estimated by a function related to discharge, water depth, and Manning’s n under the assumption of parabolic distribution of longitudinal velocity along the vertical direction for open channel flows. Two designed hypothetical cases were conducted to examine the proposed eddy viscosity relation and to demonstrate the capabilities of the proposed model. Copyrightq 2007 John Wiley & Sons, Ltd.
Received 14 September 2006; Revised 2 June 2007; Accepted 5 September 2007 KEY WORDS: layer integrated; quadratic; shape function
∗Correspondence to: Meng-Chi Hung, Natural Hazards Mitigation Research Center, National Chiao Tung University,
1001 Ta Hsueh Rd., Hsinchu, 30010 Taiwan, Republic of China. †E-mail: [email protected]
‡Postdoctoral Fellow.
§Researcher. ¶Assistant Professor. Professor.
Contract/grant sponsor: National Science Council of Taiwan, Republic of China; contract/grant number: NSC-89-2211-E-009-091
1. INTRODUCTION
Over the past decades, two-dimensional (2D) depth-averaged numerical models have been widely used to simulate unsteady free-surface flows in open channel, reservoir, and coastal water. The depth-averaged approach is computationally efficient with trade-off of discarding depth-related information but incapable of understanding the mechanism of the natural phenomenon in which vertical variations of velocities are the crucial information, such as flows in channel bends, density currents in reservoirs, or seawater intrusions in river estuaries.
With the rapid increase in computer power in recent years, several three-dimensional (3D) numerical models solving the Reynolds-averaged equations have been developed to simulate flows in natural rivers. Demuren and Rodi [1] developed a 3D model using the ‘rigid-lid’ approach, which takes the free surface into account as a symmetry plane. Hervouet and Van Haren [2] proposed a 3D Navier–Stokes solver in a domain with a moving boundary method, incorporating a 2D shallow water equations solver assuming hydrostatic pressure. Because of the pressure field calculated by 3D Navier–Stokes solver that may locally differ from the hydrostatic pressure, a procedure is employed to correct the free surface given by 2D shallow water equations according to the 3D pressure field. Meselhe and Sotiropoulos[3] presented a new code solving Reynold-averaged Navier–Stokes equations, in which the free-surface elevation is determined by allowing the computational mesh to deform during the iterative solution procedure so that the proper kinematic and dynamic conditions are satisfied at convergence. The model developed by Meselhe and Sotiropoulos [3] will be severely computationally time demanding for unsteady simulations because the mesh has to be regenerated at each iteration. As volume of fluid (VOF) [4] method has been widely used to simulate free-surface flow in computational fluid dynamics (CFD) and chemical engineering, the CFX-4 [5], reserving some modules can be extended by users for their own modeling purposes, is one of the most popular commercial CFD packages. Faure et al.[6] embedded the treatment of inlet/outlet boundary conditions in using the CFX-4 to simulate flow with free surface. This treatment is indispensable to the application of common 3D CFD codes in situ due to the lack of field velocity investigations to set up the boundary conditions.
As the fully 3D numerical models mentioned above are severely computationally time demanding and hard to be applied to the field problem due to its large demand of boundary conditions not being investigated in situ, the layer-integrated concept with hydrostatic pressure assumption was proposed to simplify the 3D shallow water free-surface flow computation and has been widely used in coastal water simulations. Li and Zhan [7] presented a multi-level finite element model with vertical integration for each water layer and applied to the tidal current flow in Tokyo Bay. Kim and Lee[8] applied the finite difference model BACHOM-3 based on the multi-level concept to simulate the tidal and wind-induced current in Suyoung Bay, Korea. Shankar et al.[9] developed a finite difference model using the same concept to simulate flow patterns for the sudden contraction and expansion flows and show reasonable flow characteristics downstream of an opening. Li and Gu[10] applied the layer-integrated concept to the mass exchange and tidal flushing problem in the harbor. The multi-level or multi-layer approach allows free surface on top layer to vary with time or iteration. The mesh does not need to be regenerated in each iteration. It significantly reduces the computational time consumption. In large-scale coastal simulation, this approach discards some local phenomena but gains practicability. Lin and Falconer[11] extended the use of layer-integrated approach to the modeling of tidal circulation with flooding and drying in Hummer Estuary along the northeast coast of England. Owing to the lack of mechanism on balancing the water depth
due to change of flow discharges, a 2D depth-averaged model is still needed to locate the water surface in their modeling.
In this paper, a semi-3D layer-integrated numerical model based on interface continuity of velocity and shear stress was proposed for the shallow water free-surface flow computation. A quadratic shape function, being applied as an auxiliary function for deriving the velocities distribution within a layer, is to assure the continuity of velocity and shear stress at interfaces. This quadratic shape function was examined through advection–diffusion equation by Tsai et al. [12]. The proposed model was verified through the wind-induced circulation in a closed basin by comparing the simulated velocity profiles to the analytical solution. The developing process of water surface and velocity profile was also investigated. Further, the flexibility of the proposed approach was tested through the formation of the velocity profiles due to change of viscosity distribution. As the turbulence model is not embedded, an estimation of eddy viscosity, related to discharge, water depth, and Manning’s roughness n under the assumption of uniform distribution of viscosity along the vertical direction, was also derived for open channel flows. Two designed hypothetical cases were conducted to examine the proposed eddy viscosity relation and to demonstrate the capabilities of the proposed model.
2. MATHEMATICAL MODEL
2.1. Governing equations
The governing equations of unsteady flows in open channels generally based on the 3D Navier– Stokes equations under hydrostatic pressure assumption can be expressed in the conservative form as follows: *u *x+ *v *y+ *w *z = 0 (1) *u *t + *(uu) *x + *(vu) *y + *(wu) *z = − 1 *p *x+ 1 *x x *x + *yx *y + *zx *z (2) *v *t+ *(uv) *x + *(vv) *y + *(wv) *z = − 1 *p *y+ 1 *x y *x + *yy *y + *zy *z (3) p= g(H +zB−z) (4)
where t is the time; x, y, z are the Cartesian co-ordinates in the horizontal plane (x, y) and vertical axis(z); u, v, w are components of velocity in the x, y, z directions, respectively; p is the pressure; is the density of water; g is the gravitational acceleration; and x x,yx,zx,x y,
yy,zy are components of the stress tensor in the x–z and y–z planes; H is the flow depth, and
zBis the elevation of channel bed, respectively.
2.2. Layer-integrated equations
Under the assumption of incompressible fluid and hydrostatic pressure, the layer-averaged equations can be obtained from integrating Navier–Stokes equations along the flow depth by layers. The layer interfaces are inclined in the longitudinal and horizontal in the transverse. A sketch of the layers
Figure 1. Schematic illustration of layer-integrated concept.
and the relative position of the variables in the x–z plane is illustrated in Figure 1. As illustrated, three layer types exist, including top layer, bottom layer, and internal layer. The thickness of top and bottom layers varies with the x, y co-ordinates to prescribe both the free surface and bottom topography, respectively. In contrast, the thickness of internal layers can be either uniform or non-uniform. Letting k+1/2= 1 h zk+1(x,y) zk(x,y) (x, y,z,t)dz (5)
where represents the physical variables such as u, v, and p, k+1/2represents its layer-average of
kth layer; zk+1refers to the elevations of the layer interfaces between the(k +1)th and kth layers, and zkrefers to that between kth and(k −1)th layers, respectively; K is the total number of layers; and h is the thickness of the layer. The governing equations for mass and momentum are first integrated over the kth layer, where k=1,2,3,..., K . Applying Leibniz’s rule to the layer-integral and neglecting high-order terms, the kth layer-averaged momentum equations become
* *t(huk+1/2)+ * *x(huk+1/2uk+1/2)+ * *y(hvk+1/2uk+1/2)+(Fu)k+1−(Fu)k =− gh*(H+zB−z) *x k+1/2 + * *x hh*uk+1/2 *x + * *y hh*uk+1/2 *y + zx k+1 − zx k (6) * *t(hvk+1/2)+ * *x(huk+1/vk+1/2)+ * *y(hvk+1/2vk+1/2)+(Fv)k+1−(Fv)k =− gh*(H+zB−z) *y k+1/2 + * *x hh*vk+1/2 *x + * *y hh*vk+1/2 *y + zy k+1 − zy k (7)
where the subscript(k +1/2) denotes the layer-averaged quantities; subscripts (k +1) and k denote the quantities at top and bottom layer interfaces, respectively; h is the layer thickness; h is the horizontal eddy viscosity; and the layer across flux F can be expressed as
Fk=wk− u*z *x k − v*z *y k (8) where wk is the vertical velocity at interface between (k −1)th and kth layers. Integration of continuity equations and application of impermeable bottom boundary condition gives vertical velocity at the interface between kth and(k +1)th layer:
wk+1= u*z *x k+1 + v*z *y k+1 − k m=1 * *x(hum+1/2)+ * *y(hvm+1/2) (9)
where m is the dummy index for summation. From Equations (6), (7), and (9), one can solve the velocity components u,v, and w for each layer. The depth variation will be solved only for the top layer from continuity equation. Re-applying Equation (9) with the kinematic boundary condition on free surface and integrating along the time step, the change of water depth,H = Ht+t− Ht, in a time step can be determined by
Ht+t− Ht t + K m=1 * *x(hum+1/2)+ * *y(hvm+1/2) =0 (10)
wheret is the time interval and K is the total number of layers.
However, there exist two more unknowns for momentum equation in the x and y directions, respectively, which are the shear stress,(/)k,(/)k+1appearing at the top and bottom interfaces of each layer. To evaluate the shear stress in terms of the velocity gradient, one may need a function to be able to present the velocity variation in the layer. Obviously, three unknowns including the velocity at center, top, and bottom of the layer exist in one layer. A quadratic function might be a reasonable and suitable choice. The formulation of layer interface constraint relation integrated with a quadratic shape function is described in the following paragraph.
2.3. Formulation of layer interface constraint relation
The shear stress at interfaces which remained in the layer-integrated momentum equations (6) and (7) can be evaluated by transforming the vertical velocity gradient into a combination of the nearby components of velocity. The transformation can be conducted by applying a quadratic shape function. A thorough analysis for the use of quadratic shape function for the dispersion equation has been performed by Tsai et al.[12]. This quadratic shape function is applied as an auxiliary function either to improve the velocity resolution or to assure the continuity of velocity and shear stress at interfaces. The following shows the derivation of first derivatives at interfaces in each layer. By adopting the quadratic polynomial interpolation function in each layer, (can be u or v) in the interval[zk;zk+1] can be represented as
Using the two nodal values of at the boundary points in a layer, namely,
(z =0) = a =k (12)
(z =z) = a +bz +cz2=k
+1 (13)
and the definition ofk+1/2, i.e.
k+1/2= 1 z z 0 (a +bz +cz 2)dz (14)
The coefficients a, b, and c in Equation (11) can be specified as follows:
a= k (15)
b= 1
z(−2k+1−4k+6k+1/2) (16)
c= 1
z(3k+1+3k−6k+1/2) (17)
where the node zk is taken as the origin and z =zk+1−zk. Thus, the velocity gradient at the interface of the kth layer can be expressed as top interface of the kth layer
* *z zt,k = 1 z(4k+1+2k−6k+1/2) (18)
bottom interface of the kth layer * *z zb,k = 1 z(−2k+1−4k+6k+1/2) (19)
Substituting Equations (18) and (19) into Equations (6) and (7) produces the layer-averaged momentum equations * *t(huk+1/2)+ * *x(huk+1/2uk+1/2)+ * *y(hvk+1/2uk+1/2)+(Fu)k+1−(Fu)k =− gh*(H +zB−z) *x k+1/2 + * *x hh*uk+1/2 *x + * *y hh*uk+1/2 *y
+v,k+1/2(6uk+1−12uk+1/2+6uk) (20)
* *t(hvk+1/2)+ * *x(huk+1/2vk+1/2)+ * *y(hvk+1/2vk+1/2)+(Fv)k+1−(Fv)k =− gh*(H +zB−z) *y k+1/2 + * *x hh*vk+1/2 *x + * *y hh*vk+1/2 *y +v,k+1/2(6vk+1−12vk+1/2+6vk) (21)
The quadratic velocity distribution in one layer introduces two more unknowns which are the velocities at the top and bottom layer interfaces. Therefore, one needs the internal interface boundary condition to close the system of equations for the flow field. By following that shear stress at the interface between two adjacent, (k −1)th and kth layers should be continuous, the velocity gradient at the bottom of the kth layer multiplying the viscosity of kth layer should be equal to the velocity gradient at the top of the (k −1)th layer multiplying the viscosity of the
(k −1)th layer: (v,k−1) * *z zt,k−1 =(v,k) * *z zb,k (22)
Then the layer interface restriction equations are obtained
2uk−1−6uk−1/2+4(1+rk)uk−6rkuk+1/2+2rkuk+1= 0 (23) 2vk−1−6vk−1/2+4(1+rk)vk−6rkvk+1/2+2rkvk+1= 0 (24) where rk=v,k+1/2/v,k−1/2is the ratio of the vertical viscosities of two adjacent layers.
2.4. Boundary conditions
The treatment of boundary conditions for the model is described as follows: • The boundary condition for the lateral walls can be slip or non-slip.
• At the free surface, the wind shear and the kinematic free-surface boundary conditions were used. At the bottom boundary, the no slip condition and impermeable condition were used. • At the layer interfaces, the continuity of velocity and shear stress is automatically satisfied
by using Equations (23) and (24) under the assumption of quadratic distribution in the layer. • At the upstream end, assuming the vertical and lateral velocities are zero, the longitudinal velocity assumed uniform along the depth can be set as the discharges per unit width divided by wetted area obtained from water depth in every iteration.
• At the downstream end, the normal depth is used and velocities were assumed homogeneous in the longitudinal direction, i.e.*u/*x =0.
3. NUMERICAL ALGORITHM
3.1. Finite volume formulation
The control volume-based finite volume formulation implies that the integral conservation of mass and momentum is exactly satisfied over any group of control volume and, of course, over the whole computational domain. By following the control volume concept in a staggered grid (plane view) shown in Figure 2, as the momentum equation has been integrated along the layer thickness, the finite volume formulation was applied only in the plane coordinates. Integrating the layer-averaged momentum equations (20) and (21) along the control volume and time interval, the conservative
Figure 2. The location of variables and their control volumes in the x–y plane.
form can be derived and given as
[(huk+1/2xy)t+t−(huk+1/2xy)t]/t
+(huk+1/2uk+1/2y)e−w+(hvk+1/2uk+1/2x)n−s+((Fu)k+1−(Fu)k)xy = hh *uk +1/2 *x y e−w + hh *uk +1/2 *y x n−s
+v(6uk+1−12uk+1/2+6uk)xy −[gh(H +zB−z)k+1/2]e−wy (25)
[(hvk+1/2xy)t+t−(hvk+1/2xy)t]/t
+(huk+1/2vk+1/2y)e−w+(hvk+1/2vk+1/2x)n−s+((Fv)k+1−(Fv)k)xy = hh *vk+1/2 *x y e−w + hh *vk+1/2 *y x n−s +v(6vk+1−12vk+1/2+6vk)xy −[gh(H +zB−z)k+1/2]n−sx (26)
where the superscripts t and t+t indicate time, and the subscripts (e,w), (n,s) indicate the faces of the CV on the x- and y-axis, respectively.
3.2. Difference schemes
The approximation to the surface integral in Equations (25) and (26) requires the values of variables at locations other than computational nodes. The upwind difference scheme is applied to the convective flux due to its dependence on the flow direction:
e=
P if (v·n)e>0
E if (v·n)e<0
(27) where variable is either u or v, subscripts e, P, E indicate locations of variables as shown in Figure 2. The upwind difference scheme has been widely used in the practice. The scheme is only 1st order accurate but can be easily extended to the 3D problem. The central difference scheme is applied to the diffusive flux due to its physical nature.
hh* *xy e ∼ =(hhy)e 1 x(E−P) (28)
where variable is either u or v, subscript P indicates the center of the CV, and subscript E indicates the center of the neighboring CVs as shown in Figure 2.
3.3. Algebraic equation system for the water column
By substituting Equations (27) and (28) into Equations (25) and (26) and then summing all the flux approximations and source terms, one can produce an algebraic equation that relates the variable value at the center of that CV to the values at the surrounding CVs.
Aututi, j,k+1+t + AuPuti, j,k+1/2+t + Aubuti, j,k+t = Au Wuit−1, j,k+1/2+ AuEuti+1,k+1/2+ AuNuti, j+1,k+1/2+ AuSuti, j−1,k+1/2+ Su (29) Avtvi, j,k+1t+t + AvPvi, j,k+1/2t+t + Avbvi, j,kt+t = Av Wvit−1, j,k+1/2+ AvEvit+1,k+1/2+ AvNvit, j+1,k+1/2+ AvSvit, j−1,k+1/2+ Sv (30) where subscripts i , j , and k are the indices of the variable along the x, y, and z axes, and superscripts
u andv indicate the coefficients for u and v momentum, respectively, the subscript P indicates the
center of the CV, and subscripts N , E, W , S indicate the center of the neighboring CVs as shown in Figure 2, S is the summation of the remaining terms including the pressure gradient and previous time step contribution. The coefficients AuP, AvP,... can be computed from Equations (25) and (26). Applying Equation (29) to the water column from bottom up to water surface layer by layer,
K layer-averaged equations will be established for the K layers. Furthermore, taking boundary
conditions on the bottom bed, water free surface, and(K −1) layer interface constraint (23) into account, a series of x-momentum equations consisting of(2K +1) equations can be established for a water column. Consequently the(2K +1) unknowns can be solved by the (2K +1) equations
system. Similarly, a series of y-momentum equations consisting of (2K +1) equations can be established for a water column from Equations (30) and (24) in the same way.
The layer-integrated equations and interface shear continuous conditions constitute the algebraic equation system of a water column. Consequently, the flow field can be solved implicitly in the vertical to satisfy the boundary conditions at the bed, interface, and the free surface simultaneously.
3.4. Computational procedure
A fully implicit time integration method is applied in the vertical direction and leads to a penta-diagonal system of equations with respect to the vertical direction. The vertical coupling of the discretized momentum equations is solved implicitly in the vertical direction by the Thomas algorithm and column by column iteratively implicitly in the horizontal direction. The main procedure can be summarized as follows and presented as a flowchart in Figure 3:
(a) Set the initial depth and velocities. (b) Compute the eddy viscosity.
(c) Solve the (2K +1) equations of x-momentum equation system ((23) and (29)) to obtain velocity u for each water column by using the Thomas algorithm.
(d) Solve the (2K +1) equations of y-momentum equation system ((24) and (30)) to obtain velocityv for each water column by using the Thomas algorithm.
(e) Repeat steps (c) and (d) if velocity field does not converge.
(f) The vertical velocity w and water depth H are determined by Equations (9) and (10), respectively.
(g) Repeat steps (b)–(f) until both velocity and depth converge, then march to the next time step.
3.5. Numerical stability limitations
The momentum equations were solved by using a combined implicit (vertical) and iteratively implicit (horizontal) schemes, but the continuity equation was solved explicitly. For the 3D shallow water equations, the horizontal velocity components are coupled in the vertical direction by the vertical advection and viscosity term. An explicit time integration of the vertical exchange in shallow areas would lead to the time step limitations[13]:
tz2 2v (31) and tz w (32)
4. ESTIMATION OF EDDY VISCOSITY FOR CHANNEL FLOWS
The law of flow resistance in open channel flow can be derived by balancing the retarding shear force at the boundary against the propulsive force acting in the flow direction. The water depth is determined by the interaction of the mean flow and bottom roughness in 1D or plane 2D
simulations. But from the point of view of 3D modeling, the bottom roughness only acts on the bottom layer. Hence, the velocity profile is determined not only by the mean flow discharge and bottom friction but also by the eddy viscosity, an improper use of eddy viscosity may result in the incorrect water depth. However, the eddy viscosity is indeed a function of the flow field, not a fluid property. As the turbulent model has not been embedded in the proposed model at the present stage, an estimation of eddy viscosities is still required. The Elder’s relation[14] between water depth, shear velocity, and eddy viscosity,=Hu∗, is frequently used in practice. Following the Elder’s relation, the determination of was derived as below.
In order to find a clear and simple relation for in Elder’s relation, uniform turbulent intensity (i.e. uniform eddy viscosity) along the depth is assumed and will lead to parabolic distribution of horizontal velocities in the vertical direction. Consider 3D Navier–Stokes equation in the mean flow direction *u *t +u *u *x+v *u *y+w *u *z=− 1 *P *x + 1 *x x *x + *yx *y + *zx *z +gx (33)
Assuming that flow is a steady 2D uniform flow and slope is small, the equation can be reduced to
v* 2
u
*z2=−gS0 (34)
Integrating Equation (34) along the water depth becomes *u
*z=−
gS0
v
z+C1
and applying bottom shear condition *u/*z|z=0=0/v yields C1=0/v, then the vertical gradient of velocity can be expressed as
*u *z=− gS0 v z+ 0 v (35) By applying the surface condition
*u *z z=H = s v Equation (35) becomes s v=− gS0 v H+ 0 v or −gS0H+0 = s
Therefore, bottom shear 0 can be determined by water depth, bottom slope, and surface (wind) shear
Integrating Equation (35) along water depth again becomes u=−gS0 2v z2+0 v z+C2
Applying no slip condition at the bottom bed, C2 will be diminished. Then the velocity profile can be expressed as u=−gS0 2v z2+0 v z (37)
For mass conservation, integration of the velocity along water depth should be consistent to the flow discharge, i.e.
H 0
u dz=q (38)
Substituting Equation (37) into Equation (38), then
gS0 3v
H3+ s
2vH
2=q (39)
The eddy viscosity cannot be determined by Equation (39) alone due to another remaining unknown H . If there were no wind shear acting on water surface, the energy slope is equal to the bottom slope because of uniformity in flow direction, then the kinematic eddy viscosity can be expressed as
v=
gS0 3q H
3 (40)
or in the form of Elder’s relationv=Hu∗, in which,
= u∗
3U
Consider Manning’s equation in a wide rectangular channel, the depth-averaged velocity,
U= 1 nH 2/3S1/2 0 = H1/6 n√gu∗ (41)
Then coefficient in Elder’s relation can also be expressed as
= u∗
3U=
n√g
3H1/6 (42)
Hence, the eddy viscosity can be determined by Manning’ roughness n, water depth h, and shear velocity u∗. Also, the bed slope in Equation (41) can be replaced by the friction slope so that the non-uniform flow situation can also be satisfied. Consequently, Equation (42) could be useful for the practical problems. Thus, can be estimated easily by Manning’ roughness n,
water depth H in situ if no investigated data were available. The value of is about 0.014–0.058 calculated by Equation (42) for natural river flows (about n=0.02–0.05, H =0.5–10 m). This is slightly smaller but of the same order compared with data investigated by Nezu and Rodi[15] due to the assumption of uniform turbulent intensity.
5. MODEL VERIFICATION
5.1. Wind-induced flow
Wind-induced flow in a closed basin in which wind shear acting on the water surface results in a drift current in the direction it blows and a bottom return flow in the opposite direction is chosen to verify the model’s accuracy and applicability. The analytical solution for steady laminar flow under no-slip condition, given by Heaps[16], for constant vertical eddy viscosity is
u(z)= s 4v z(3z −2H) H (43)
where u(z) is the velocity at any depth z; H the water depth; v the vertical eddy viscosity;sthe wind stress acting on the surface; and the water density.
This numerical experiment is carried out to simulate a steady wind-induced flow in a uniform rectangular channel of length 400 m, width 160 m, and depth 2 m. The water initially at rest was subject to a uniform wind stress, increasing linearly from 0 N/m2at t=0s to 0.5N/m2at t=10s, and maintaining 0.5N/m2at t>10s. The following parameters were adopted: x =20m, y =8m,
=1000kg/m3,
v=0.015m2/s.
A constant viscosity leads to a parabolic velocity profile in the case of laminar flow. Effect of number of layers was also carried out with four cases of 5, 10, 20, and 30 layers. Comparison of model predictions at the middle length of the basin to the analytical solution is shown in Figure 4.
Figure 5. Effect of number of layers on water depth prediction.
Figure 6. Variation of velocity in the wind-induced flow at middle length(x =200m).
The major improvement of predictions of velocity profile takes place between 5 and 10 layers, and no more significant difference takes place between predictions over 20 layers. Comparison of the u profile at middle length using 20 layers to the analytical solution is as expected. The upper layer currents are in the wind direction while the lower layer currents are in the opposite direction. It can be seen that the prediction has a slightly larger velocity. The predicted depth with zero velocity is at 13 of water depth under the free surface. It is consistent with the analytical solution. The wind setup of the water surface has the same tendency as shown in Figure 5. With increase
Figure 7. Variation of water depth in the wind-induced flow.
in the surface wind shear or decrease in the eddy viscosity, the maximum Reynolds number tested for limited stable condition can be raised up to around 300.
The evolution of the wind-induced circulation could be classified into two steps: (a) surface flow and water setup induced by wind shear and (b) bottom reverse flow induced by pressure difference. The transient process was also carried out to examine the physical mechanism. Figures 6 and 7 show the simulated evolution of velocity profile at middle length and water surface setup varying with time with 20 layers, respectively. From Figure 6 it is obvious that the shear-induced surface flow in the beginning and the pressure-induced reverse flow at middle length were successfully predicted by the proposed model. These simulation results include surface setup and pressure-induced circulation, two significant components of the physical mechanism.
5.2. Effect of viscosity distribution on wind-induced flow
Owing to the distribution of vertical eddy viscosities existing in natural flows not being uniform, three more types of viscosity distribution were considered to further investigate the capability of the proposed quadratic shape function for solving free-surface flows. The parabolic law of eddy viscosity distribution is the most frequently assumed type in natural flows [15]. It has the maximum viscosity at middle depth and the minimum near the surface and the bottom. Figure 8 shows the simulated velocity profile with parabolic viscosity distribution compared to the one with uniform viscosity. The larger velocity gradient is found near both surface and bottom due to its smaller viscosity. Contrarily, the smaller velocity gradient is found at middle depth due to its larger viscosity. The simulated profile behaves in a well-consistent tendency with its physical mechanism. To further examine the flexibility of the auxiliary shape function, the trapezoidal distribution of viscosity was also conducted. The comparison of the velocity profiles due to change of the viscosity distribution is shown in Figure 9. For the case of trapezoidal type I with increased viscosity from top to bottom, the simulated velocity profile is well adapted as its gradient decreases from top to
Figure 8. Velocity profiles corresponding to their own viscosities (uniform vs parabolic).
bottom by the quadratic shape function. By preserving the condition of zero net flow, the point of zero is risen a little compared with the uniform law. Figure 9 also shows that the quadratic shape function can also adapt the velocity profile to the trapezoidal type II of viscosity distribution as well.
5.3. Examination of proposed eddy viscosity relation
Because there is no net flow in the wind-induced flow in a closed basin, the non-uniformity of viscosity distribution does not affect the mean water depth but the velocity profile. As far as 1D or plane 2D is concerned, the normal depth is well known as a function of flow discharge and roughness in natural rivers. However, it is not sufficient in the 3D free-surface flow computations because the water depth and velocity profile are also affected by the eddy viscosity. The following case is designed hypothetically to investigate the effect of viscosity variation on channel flow. A steady flow is carried in a uniform rectangular channel with length 5000 m, width 100 m, slope 0.0005, Manning’s roughness n of 0.035, and discharge of 3.987 cm per unit width. The normal depth determined by Manning’s equation is 3.0 m. The slip condition was used at the side walls but no-slip condition at the channel bed. The water body was divided into four layers of equal thickness. The following parameters were specified: x =500m, y =10m, =1000kg/m3. Upstream, the lateral and vertical velocities v and w were set as zero, and the longitudinal velocity u was equal to the inflow discharge divided by the instantaneous flow area. Downstream, zero longitudinal velocities gradient*u/*x =0, *v/*x =0 were applied. Three types of viscosity variation distribution employed as shown in Figure 10 and the coefficient in Elder’s relation
Figure 11. Velocity profiles corresponding to the three types of viscosity variation distribution.
can be expressed as follows:
Type I: = v H u∗=0.0304 Type II: =0.0304 0.5+4 z H 1− z H Type III: = z H 1− z H , =0.41
where type I is calculated by Equation (42); type III is a parabolic distribution function proposed by Elder[14]; and type II is an arbitrary function.
Figure 11 shows the velocity profiles corresponding to the three types of viscosity variation distributions, where the shape of velocity profile apparently varies with the type of viscosity variation distribution employed. The shapes of velocity profiles for type I and type III are parabolic and logarithmic, respectively, while the shape for type II is in between.
The water depth in 1D and 2D shallow water flow computation is determined by the flow discharge, energy slope, and channel roughness. From Figure 11 one can see that the water depth is affected not only by these three factors but also by the viscosity variation distribution in the 3D free-surface flow computations. However, the flow regime is indeed turbulent because of their high roughness and irregular bed forms and the eddy viscosity is a function of flow field,
Table I. Comparison of computed water depth at mid-length from simulation and Manning’s equation.
n=0.03 n=0.035 n=0.04 Manning’s equation 2.74 3.00 3.25
Simulated 2.79 3.00 3.28
Figure 12. Simulation of water surface profiles with different n values.
not a fluid property. The difficulty of the 3D channel flow simulation is that the eddy viscosity changes anytime and anywhere so that manually calibrating eddy viscosity in the model is almost impossible. Before the turbulent model is considered in the proposed model, an estimation of eddy viscosity, related to Manning’s roughness n, water depth, and flow discharge is needed. The case designed hypothetically is the same as above, used to examine the eddy viscosity relation derived previously in this study. is estimated to be about 0.0304 by Equation (42), and the corresponding eddy viscosity is about 0.011m2/s.
This numerical experiment with same conditions as the previous case is carried out to verify the relation between flow depth and Manning’s roughness n. Three different values 0.03, 0.035, and 0.04 of Manning’s roughness n were considered in the numerical experiments. Table I shows the comparison of computed water depths with those calculated by Manning’s equation. Both the water depth and velocity variation obviously demonstrate the response of the change of Manning’s roughness n. The normal depths are 2.74, 3.0, and 3.25 m for n=0.03, 0.035, and 0.04, respectively. Computed water depths are a slightly larger than those from Manning’s equation. Comparison of the water depths under different n values is shown in Figure 12. The velocity profiles of inner cross section at middle length for various Manning’s n are also compared in Figure 13. The larger the value of roughness considered, the larger the depth and smaller the velocity obtained in the simulation.
Figure 13. Comparison of velocity profiles with different n values at x=2500m.
Figure 14. Simulation of water surface profile under different wind shears.
5.4. Channel flow with wind shear
The final case demonstrated in the following is to solve the wind shear effect on the flow patterns. All conditions were the same as the previous experiment with n=0.035 except that the wind shears acting on the surface are 0.0, 2.0, and 5.0N/m2directing to the upstream. Figure 14 shows
Figure 15. Comparison of velocity profiles under different wind shears at x=2500m.
the comparison of wind-induced setups under different shear stresses. The water surface profile behaves like the M2 curve under the action of wind shears because the rating curve was used as downstream boundary conditions. Increase of the water depth drives the flow against the additional resistance on the surface. The maximum velocity takes place at one-fourth of depth under free surface due to wind shear resistance acting on the water surface. The converged water depth is 3.0 m downstream and increases gradually to 3.17 m upstream under wind shear of 2.0N/m2due to wind-induced setup. The velocity profiles were also compared in Figure 15 under different values of wind shears.
6. CONCLUSION
There were a few 3D free-surface flow models developed in the past decades but least applied to field problem due to its heavy computational time consumptions. A novel interface continuity-based approach proposed in this study makes the layer approach applicable in the shallow water free-surface flow computation. Use of a quadratic shape function guarantees that there is no discon-tinuity situation taking place between two adjacent layers further. The modeling procedure adopted provides a simple and robust method for simulating 3D free-surface shallow water flows. The verification performed in the case of wind-induced flow in a closed rectangular basin by comparing with the analytical solution shows a reasonable agreement in the simulated steady-state profiles. The evolution of water surface slope and the velocity profiles at the middle length of the basin
were examined to investigate whether the model responding to the wind shear is consistent with the physical nature. Alternate shapes of velocity profile formed by their correspondent distribution of vertical viscosity also show the flexibility of the proposed approach either in wind-induced flows or in channel flows. The variations of velocity profile due to change of viscosity distribution and near bed velocity are also investigated. Further, a simplified relation between discharge, water depth, and eddy viscosity is derived as the closure of the governing equations. The proposed model applicability for the 3D shallow free-surface water flow simulations was examined through two hypothetical cases and shows convincing results. Further development of considering the turbulent model is needed for real-world application.
NOTATION h layer thickness n Manning’s roughness p hydrostatic pressure u x-component of velocity u∗ shear velocity v y-component of velocity
zB elevation of bottom bed
A coefficients of the algebraic equation for water columns
C arbitrary constant
F layer across flux
H water depth
K number of layers
S slope
Su, Sv source terms
U depth-averaged velocity
arbitrary constant for Elder’s relation
specific weight of water
dynamic eddy viscosity
kinematic eddy viscosity
water density
shear stress
s wind shear
Subscripts
e,s,w,n denote faces of the considered cell
h horizontal direction
k denote the layer interface quantity between kth and(k −1)th layers
k+1/2 denote the kth layer-averaged quantity
m dummy index used for integration along water depth
v vertical direction
E, S, W, N denote centroid of neighbors of the considered cell
ACKNOWLEDGEMENTS
Partial financial support of this study from National Science Council of Taiwan, Republic of China, through Contract NSC-89-2211-E-009-091 is greatly appreciated. The computer resources used in this study were provided by National Center for High-Performance Computing in Taiwan, Republic of China.
REFERENCES
1. Demuren AO, Rodi W. Side discharges into open channels: mathematical model. Journal of Hydraulic Engineering 1983; 109(12):1707–1722.
2. Hervouet JM, Van Haren L. Recent advances in numerical methods for fluid flows. In Floodplain Processes, Chapter 6, Anderson MG, Walling DE, Bates PD (eds). Wiley: New York, 1996.
3. Meselhe EA, Sotiropoulos F. Three-dimensional numerical model for open channel with free surface variations.
Journal of Hydraulic Research 2000; 38(2):115–121.
4. Hirt CW, Nicholls BD. Volume of fluid (VOF) method for dynamics of free boundaries. Journal of Computational
Physics 1981; 39:201–221.
5. CFX. CFX-4.2: Solver. AEA Technology: U.K., 1997.
6. Faure JB, Buil N, Gay B. 3-D modeling for unsteady free-surface flow in open channel. Journal of Hydraulic
Research 2004; 42(3):263–272.
7. Li YS, Zhan JM. An efficient three-dimensional semi-implicit finite element scheme for simulation of free surface flows. International Journal for Numerical Methods in Fluids 1993; 16:187–198.
8. Kim CK, Lee JS. A three-dimensional PC-based hydrodynamic model using an ADI scheme. Coastal Engineering 1994; 23:271–287.
9. Shankar NJ, Chan ES, Zhang QY. Three-dimensional numerical simulation for an open channel flow with a constriction. Journal of Hydraulic Research 2001; 39(2):187–201.
10. Li CW, Gu J. 3D layered-integrated modeling of mass exchange in semi-enclosed water bodies. Journal of
Hydraulic Research 2001; 39(4):403–412.
11. Lin B, Falconer RA. Three-dimensional layer-integrated modeling of estuarine flows with flooding and drying.
Estuarine, Coastal and Shelf Science 1997; 44:737–751.
12. Tsai TL, Yang JC, Huang LH. An accurate integral-based scheme for advection–diffusion equation.
Communications in Numerical Methods in Engineering 2001; 17:701–713.
13. Delft Hydraulics. User Manual of Delft3D-FLOW—Simulation of Multi-dimensional Hydrodynamic Flows and
Transport Phenomena, Including Sediments. Delft Hydraulics, 2003.
14. Elder JW. The dispersion of marked fluid in turbulent shear flow. Journal of Fluid Mechanics 1959; 5:544–560. 15. Nezu I, Rodi W. Open-channel flow measurements with a laser Doppler anemometer. Journal of Hydraulic
Engineering 1986; 112:332–355.
16. Heaps NS. Vertical structure of current in homogeneous stratified waters. Hydrodynamics of Lakes, Hutter K (ed.). Springer: New York, 1984; 154–202.
