2DV FEM MODELING OF ESTUARINE TIDAL FLOW

S.S. Li and Z. Shi

Department of Harbour & Coastal Engineering,

Shanghai Jiao Tong University,

1954 Hua Shan Road, Shanghai 200030, China

Tel: +86+21+62933027, Fax: +86+21+62933160, E-mail: shsli@263.net

Abstract: -coordinate transformation with the finite difference method (FDM) has been widely applied to two-dimensional vertical (2DV) numerical model of estuarine tidal flow. However, after its transformation, several nonlinear terms are added to the width-integrated Navier-Stokes equations on the regular grids. As a result, the corresponding model becomes complicated and time consuming. This study has established an efficient 2DV numerical model of estuarine tidal flow, in which the governing equations have been solved by the moving grid finite element method (FEM) instead of -coordinate transformation. Verifications of the model suggest that there is a good agreement between theoretical and calculated results.

Keywords: 2DV FEM estuarine tidal flow

1  INTRODUCTION

Tidal flow is the main hydrodynamic force in a tidal estuary. Both scientists and engineers have been interested in tidal flow because of its significance in estuarine hydraulics, harbour construction, channel dredging and pollutant monitoring. To this end, numerical models have been developed and provided an effective method to simulate estuarine tidal flow field. Two principal types of 2D numerical models have been used, i.e., two-dimensional horizontal (2DH) (Chu & Yeh, 1980; Maa, 1990; Barer et al., 1995; and Mead, 1999) and two-dimensional vertical (2DV) models (Blumberg, 1977; Boericke & Hogan, 1977; Li et al., 1994; Kuo et al., 1996; and Mead, 1999). 2DH models, based on the depth-integrated equations, could simulate the horizontal distributions of tidal flows, but vertical variations are parameterised or treated in terms of analytic functions. 2DV models, based on the width-integrated equations, could predict the vertical distributions of estuarine flow in narrow estuaries. In those models above, FDM and FEM have been commonly used [e.g., FDM (Swift et al, 1979; Alfrink & Van Rijn, 1983; and Li et al., 1994); FEM (Malcherek, 2000)].

-coordinate transformation with FDM has also been used because of significant variations in estuarine bottom topography and tidal elevation when setting-up the computing point grids (Boericke & Hogan, 1977; Celik & Rodi, 1985; and Lin & Mehta, 1989). Flows throughout the computing domains are resolved by a constant number of sub-layers. The computation in the -direction can be performed on a fixed grid. After -transformation, several non- linear terms are added to the governing equations and cause more difficulties to be solved.

FEM, based on variational principles or weighted residual method (Chung, 1978), is an approximate approach to solving partial differential equations in engineering applications. It can better deal with estuarine complex bathymetry and land boundaries. According to whether the computing area and grid are varying or not, it can be grouped into three types: 1) traditional FEM with any fixed computing area and grid; 2) auto-adaptive FEM with a fixed computing area and varying grids; and 3) moving grid FEM with varying computing areas and grids. In 2DV numerical models of estuarine tidal flow, computing areas and grids must vary with the tidal elevation. By applying the moving grid FEM to the derivations of the governing equations, -coordinate transformation is no longer necessary.

Despite of those studies above, it is still necessary to develop an efficient numerical model to predict estuarine tidal flows because of their complexes. The primary objective of the present paper is to develop and verify a two-dimensional vertical (2DV) estuarine tidal model being derived by the moving grid FEM, coupled with a turbulence closure sub-model.

2  NUMERICAL MODEL

Governing equations

Considering a right-hand Cartesian coordinate system (Fig. 1), the governing equations can be obtained by integrating Navier-Stokes equations over estuarine width ( direction, Fig. 2). The following assumptions are made: 1) hydrostatic approximation of the vertical momentum equation; 2) ignoring lateral variations of the momentum distributions; 3) the Coriolis acceleration is omitted; 4) Boussinesq approximation; and 5) the non-slip conditions on the two sides of the estuary (where ). Then we can obtain the width-averaged equations for conservation of estuarine tidal water mass and momentum:

Continuity equation：

                                                         (1)

Tidal elevation equation：

                                                      (2)

Momentum equation：

                                           (3)

Hydraulic pressure：

                                                     (4)

Fig. 1  Cartesian coordinate system

Fig. 2  Schematic estuarine cross-section

where  denotes time;  the tidal elevation above the mean water level;  and  components of eddy viscosities in  and  directions; pressure; intensity of pressure;  the fluid density;  the gravity acceleration;  and  components of tidal flow velocities in  and  directions； the estuarine width;  the estuarine depth (Fig. 2);  the total water depth; and  and  horizontal and vertical Cartesian coordinates.

Boundary and initial conditions

(a) Tidal elevation and velocity in each vertical layer can be specified at the open boundaries.

Upperstream open boundary:

, ,                                                      (5)

Lowerstream open boundary:

, ,                                                    (6)

Boundary conditions of tidal elevation and velocities should be given. However, the open boundary for tidal velocities in each vertical layer cannot be given in many cases. According to the St. Venant's principle, the influence will be limited to the nearby of the open boundary.

(b) At the bottom:

 ( is the bottom friction shear stress)                                      (7)

(c) At the free water surface:

  ( is the surface wind shear stress)                                       (8)

Initial values for both tidal elevation and velocity distribution must be specified. The former is obtained by linear interpolation between the upper and lower boundaries values, while the latter is taken as a logarithmic profile or a constant value: .

Eddy viscosities

The horizontal eddy viscosity is defined as (Sinha et al., 1999):

                                                 (9)

where ;  the time increment;  the average length of the main flow line;  the direct distance between the calculating point and the side of the estuary.

The vertical eddy viscosity is obtained from the Prandtl's mixing length assumption (Sinha et al., 1999):

                                                        (10)

where  is the Richardson number, accounting for the effect of vertical density gradients on the vertical momentum exchange; , ；  the coefficient； the mixing length ,  von Karman’s constant;  and  roughness length.

3  MOVING GRID FEM DERIVATIONS OF THE GOVERNING EQUATIONS

Derivations of the continuity and momentum equations

Derivations of the following equations are restricted within any local element (Fig. 3). Considering an arbitrarily shaped quadrilateral element (Fig. 3a), the isoparametric coordinate [ ] (Fig. 3b), ranging from 0 to +/-1, is established at the centroid of the element. The reference Cartesian coordinate [ ] (Fig. 3a) is related to the isoparametric coordinate as , . Where  is the isoparametric function (Ferziger & Peric, 1996):

.

Let ，and introduce a set of functions: , ， ， ，in which  are the numerical values of variables at the element node points. Because these functions are approximate ones, if being substituted into Eqs. (1) and (3), they will not exactly satisfy the governing differential equations and may result in an error or residual . Using the Galerkin weighted residual method, the inner product is set to zero [ ] (Chung, 1978).

 Fig. 3  Four-node isoparametric element

 Thus, the local finite element equations can be written as:

                                        (11)

                               (12)

where ;

; ; ;  representing matrices of the size  can be written as (the integration may be performed by means of the Gaussian quadrature): ， ， ， ， ， ， ， ； is Jacobian matrix； is the determinant of Jacobian ; and  is the element.

Derivation of the tidal elevation equation

Considering a one-dimensional element (two-node system) as shown in Fig. 4, the shape functions have the forms (Ferziger & Peric, 1996):   , in which, ，  is the length of the element, and  the reference Cartesian coordinate of the element-center (Fig. 4).

Let  and introduce a set of functions in the form , , .  are the numerical values of variables at the element node points. Because these functions are approximate ones too, if being substituted into Eq. (2), it may result in an error or residual . Using the Galerkin weighted residual method, the inner product is set to zero [ ] (Ferziger & Peric, 1996).

Thus, the local finite element equation has the form:

                                        (13)

where ； ； ； ； .

 Fig. 4  One-dimensional two-node element

For unsteady flow, the time derivative will be replaced by various forms of finite differences. In general, the explicit scheme is conditionally stable while the implicit scheme is unconditionally stable for any time increment and needs more computational resource due to non-linear matrix algebraic equations. Therefore, the explicit scheme is often used. However, any attempt to perform stability analysis is a formidable task and so the time increment is mostly acceptable throughout numerical modeling.

Calculating process and verification of model

Calculating process

Local element coefficients are assembled according to the Boolean matrices, which locate the appropriate local nodal contributions into the corresponding global system. Global finite element equations have the same forms with local ones. As shown in Fig. 5 (a), the computing areas and grids and the variables [ ] at time [ ] are all known. Tidal elevation  at time [ ] can be obtained from Eqs. (13), (7) and (8). Due to the varying free surface, the computing areas and grids should be updated according to  [Fig. 5 (b)]. Because of the varying areas and the grids, the velocity distributions [ or ] should be obtained by linear interpolation on the velocity distributions [ or ]. If the time increment  is small enough, this step can be omitted. Velocities [ ] and [ ] can be obtained from Eqs. (12) and (11). Then the same steps above will be repeated in the following cycle. The same topologies of two neighbouring grids keep the section locations fixed, but the  coordinate changes.

Fig. 5  Meshes of two neighboring times
        (a) Time   (b) Time

Verification of model

Estuarine steady flow (bottom slope 1/1000 and water depth 9m) is used to test the model above. Fig. 6 shows theoretical and calculated results, while Fig. 7 shows calculating convergence. From these figures, we can conclude that the model is useful and steady with good convergence.

Fig. 6  Vertical distributions showing comparison between
theoretical and calculated results

Fig. 7  Calculating convergences

4  CONCLUSIONS

An efficient two-dimensional vertical (2DV) model, coupled with a turbulence closure sub-model, is established to simulate estuarine tidal flow. In the mean time, moving grid FEM (isoparametric element) is applied to the model without -coordinate transformation. The model is verified and should be applicable in engineering scheme.

The following conclusions can be drawn from this study: 1) -coordinate transformation is associated with the finite difference method (FDM) in the estuarine tidal model; 2) the FEM derivations of the governing equations with -coordinate transformation become more complicated than without it; and 3) if the governing equations are derived by the FEM, -coordinate transformation is no longer necessary.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (Marine Sciences Program Grant 498060005) and the State Education Ministry of China.

References

[1]  Alfrink, B.K. and Van Rijn, L.C. 1983. Two-equation turbulence model for flow in trenches. Journal of Hydraulic Engineering 109(7): 941-957.

[2]  Barer, R.W., Pearson, R.V., Roberts, A.P. and McKee, W.I. 1995. Numerical simulation of the Humber estuary using a non-orthogonal curvilinear coordinate system. In: C.A. Brebbia, L. Traveresoni and L.C. Wrobel (Eds.), Computer Modeling of Seas and Coastal Regions II, Computational Mechanics Publications, Southampton Boston: 29-37.

[3]  Blumberg, A.F. 1977. Numerical model of estuarine circulation. Journal of the Hydraulic Division 103 (HY3): 295-309.

[4]  Boericke, R.R. and Hogan, J.M. 1977. An X-Z hydraulic/thermal model for estuaries. Journal of the Hydraulics Division 103(HY1): 19-37.

[5]  Celik, I. and Rodi, W. 1985. Calculation of wave-induced turbulence flows in estuaries. Ocean Engineering 12(6): 531-542.

[6]  Chu, W.S. and Yeh, W.G. 1980. Two-dimensional tidally averaged estuarine model. Journal of the Hydraulics Division 106(HY4): 501-518.

[7]  Chung, T.J. 1978. Finite element analysis in fluid dynamics. McGraw-Hill International Book Company, New York: 1-378.

[8]  Ferziger, J.H. and Peric, M. 1996. Computational methods for fluid dynamics. Springer-Verlag Berlin Heidelberg, Germany: 321-340.

[9]  Johnson, B.H., Kim, K.W., Heath, R.E., Hsieh, B.B. and Butler, H.L. 1993. Validation of three-dimensional hydrodynamic model of Chesapeake Bay. Journal of Hydraulic Engineering 119 (1): 2-20.

[10]  Kuo, A.Y., Shen, J. and Hamrick, J.M. 1996. Effect of acceleration on bottom shear stress in tidal estuaries. Journal of Waterway, Port, Coastal, and Ocean Engineering 122 (2): 75-83.

[11]  Li, Z.H., Nguyen, K.D., Brun-Cottan, J.C. and Martin, J.M. 1994. Numerical simulation of the turbidity maximum transport in the Gironde estuary (France). Oceanologica ACTA 17 (5): 479-500.

[12]  Lin, C.P. and Mehta, A.J. 1989. Turbidity-induced sedimentation in closed-up channels. Journal of Coastal Research 5(3): 391-401.

[13]  Maa, J.P.-Y. 1990. An efficient horizontal two-dimensional hydrodynamic model. Coastal Engineering 14:1-18.

[14]  Malcherek, A. 2000. Application of TELEMAC-2D in a narrow estuarine tributary. Hydrological Processes 14: 2293-2300.

[15]  Mead, C.T. 1999. An investigation of the suitability of two-dimensional mathematical models for predicting sand deposition in dredged trenches across estuaries. Journal of Hydraulic Research 37 (4): 447-464.

[16]  Sinha, P.C., Rao, Y.R., Dube, S.K., Dube, S.K., Murthy, C.R. and Chatterjee, A.K. 1999. Application of two turbulence closure schemes in the modeling of tidal currents and salinity in the Hooghly Estuary. Estuarine, Coastal and Shelf Science 48: 649-663.

[17]  Swift, M.R., Reichard, R. and Celikkol, B. 1979. Stress and tidal current in a well-mixed estuary. Journal of the Hydraulic Division 105(HY7): 785-799.