UNCORRECTED PROOF
3. NUMERICAL ALGORITHM23
*z
zt,k−1
=(v,k)
*
*z
zb,k
(22) Then the layer interface restriction equations are obtained
9
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 11
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.
13
• 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.
15
• 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.
17
• 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 19
by wetted area obtained from water depth in every iteration.
• At the downstream end, the normal depth is used and velocities were assumed homogeneous 21
in the longitudinal direction, i.e.*u/*x =0.
3. NUMERICAL ALGORITHM 23
3.1. Finite volume formulation
The control volume-based finite volume formulation implies that the integral conservation of mass 25
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) 27
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 29
momentum equations (20) and (21) along the control volume and time interval, the conservative form can be derived and given as
31
[(huk+1/2xy)t+t−(huk+1/2xy)t]/t
Copyrightq 2007 John Wiley & Sons, Ltd. Commun. Numer. Meth. Engng (2007)
UNCORRECTED PROOF
8
CNM 1061
M.-C. HUNG ET AL.
Figure 2. The location of variables and their control volumes in the x–y plane.
+(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 1
of the CV on the x- and y-axis, respectively.
UNCORRECTED PROOF
CNM 1061
APPROACH FOR SHALLOW WATER FREE-SURFACE FLOW COMPUTATION 9 3.2. Difference schemes
1
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 3
convective flux due to its dependence on the flow direction:
e=
P if (v·n)e>0
E if (v·n)e<0 (27)
5
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 7
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.
9
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 11
indicates the center of the neighboring CVs as shown in Figure 2.
3.3. Algebraic equation system for the water column 13
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 15
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 + Aubut+ti, j,k
= AuWuit−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
= AvWvit−1, j,k+1/2+ AvEvit+1,k+1/2+ AvNvi, j+1,k+1/2t + AvSvi, j−1,k+1/2t + Sv (30) where subscripts i , j , and k are the indices of the variable along the x, y, and z axes, and superscripts 17
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 19
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 21
(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 23
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 25
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 27
established for a water column from Equations (30) and (24) in the same way.
Copyrightq 2007 John Wiley & Sons, Ltd. Commun. Numer. Meth. Engng (2007)
UNCORRECTED PROOF
10
CNM 1061
M.-C. HUNG ET AL.
The layer-integrated equations and interface shear continuous conditions constitute the algebraic 1
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
3.4. Computational procedure
A fully implicit time integration method is applied in the vertical direction and leads to a penta-5
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 7
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:
9
(a) Set the initial depth and velocities.
(b) Compute the eddy viscosity.
11
(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.
13
(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.
15
(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), 17
respectively.
(g) Repeat steps (b)–(f) until both velocity and depth converge, then march to the next time 19
step.
3.5. Numerical stability limitations 21
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 23
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 25
shallow areas would lead to the time step limitations[13]:
tz2 2v
27 (31) and
tz
w (32)
29