UNCORRECTED PROOF
2. MATHEMATICAL MODEL17
surface in their modeling.
In this paper, a semi-3D layer-integrated numerical model based on interface continuity of 3
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 5
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.
7
[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 9
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.
11
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 13
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 15
capabilities of the proposed model.
2. MATHEMATICAL MODEL 17
2.1. Governing equations
The governing equations of unsteady flows in open channels generally based on the 3D Navier–
19
Stokes equations under hydrostatic pressure assumption can be expressed in the conservative form as follows:
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 25
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 27
zBis the elevation of channel bed, respectively.
2.2. Layer-integrated equations 29
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 31
interfaces are inclined in the longitudinal and horizontal in the transverse. A sketch of the layers
Copyrightq 2007 John Wiley & Sons, Ltd. Commun. Numer. Meth. Engng (2007)
UNCORRECTED PROOF
4
CNM 1061
M.-C. HUNG ET AL.
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, 1
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 3
topography, respectively. In contrast, the thickness of internal layers can be either uniform or non-uniform. Letting
where represents the physical variables such as u, v, and p, k+1/2represents its layer-average of 7
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;
9
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 11
and neglecting high-order terms, the kth layer-averaged momentum equations become
*
UNCORRECTED PROOF
CNM 1061
APPROACH FOR SHALLOW WATER FREE-SURFACE FLOW COMPUTATION 5 where the subscript(k +1/2) denotes the layer-averaged quantities; subscripts (k +1) and k denote 1
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
3
where wk is the vertical velocity at interface between (k −1)th and kth layers. Integration of 5
continuity equations and application of impermeable bottom boundary condition gives vertical velocity at the interface between kth and(k +1)th layer:
7
where m is the dummy index for summation. From Equations (6), (7), and (9), one can solve the 9
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 11
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
13
wheret is the time interval and K is the total number of layers.
15
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 17
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 19
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 21
with a quadratic shape function is described in the following paragraph.
2.3. Formulation of layer interface constraint relation 23
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 25
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 27
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 29
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) 31
in the interval[zk;zk+1] can be represented as
(z)=a +bz +cz2 (11)
33
Copyrightq 2007 John Wiley & Sons, Ltd. Commun. Numer. Meth. Engng (2007)
UNCORRECTED PROOF
6
CNM 1061
M.-C. HUNG ET AL.
Using the two nodal values of at the boundary points in a layer, namely, 1
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 5
interface of the kth layer can be expressed as top interface of the kth layer
*
bottom interface of the kth layer
*
Substituting Equations (18) and (19) into Equations (6) and (7) produces the layer-averaged momentum equations
UNCORRECTED PROOF
CNM 1061
APPROACH FOR SHALLOW WATER FREE-SURFACE FLOW COMPUTATION 7 The quadratic velocity distribution in one layer introduces two more unknowns which are 1
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 3
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 5
equal to the velocity gradient at the top of the (k −1)th layer multiplying the viscosity of the (k −1)th layer:
7
(v,k−1)
*
*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