Non-hydrostatic versus hydrostatic modelings of free surface flows*

ZHANG Jing-xin (张景新)

MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China

Key Laboratory of Estuarine and Coastal Engineering, Ministry of Transport, Shanghai 201201, China

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China, E-mail: zhangjingxin@sjtu.edu.cn SUKHODOLOV Alexander N.

Department of Ecohydrology, Institute of Fresh-water Ecology and Inland Fisheries, Berlin 12587, Germany

LIU Hua (刘桦)

MOE Key Laboratory of Hydrodynamics, Shanghai Jiaotong University, Shanghai 200240, China

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

Non-hydrostatic versus hydrostatic modelings of free surface flows*

ZHANG Jing-xin (张景新)

MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China

Key Laboratory of Estuarine and Coastal Engineering, Ministry of Transport, Shanghai 201201, China

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China, E-mail: zhangjingxin@sjtu.edu.cn SUKHODOLOV Alexander N.

Department of Ecohydrology, Institute of Fresh-water Ecology and Inland Fisheries, Berlin 12587, Germany

LIU Hua (刘桦)

MOE Key Laboratory of Hydrodynamics, Shanghai Jiaotong University, Shanghai 200240, China

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

(Received December 7, 2012, Revised May 15, 2014)

The hydrodynamics of geophysical flows in oceanic shelves, estuaries, and rivers are often studied by solving shallow water equations under either hydrostatic or non-hydrostatic assumptions. Although the hydrostatic models are quite accurate and cost-efficient for many practical applications, there are situations when the fully hydrodynamic models are preferred despite a larger cost for computations. The present numerical model is implemented by the finite volume method (FVM) based on unstructured grids. The model can be efficiently switched between hydrostatic and non-hydrostatic modules. The case study shows that for waves propagating along the bar a criterion with respect to the shallowness alone, the ratio between the depth and the wave length, is insufficient to warrant the performance of shallow flow equations with a hydrostatic approach and the nonlinearity in wave dynamics can be better accounted with a hydrodynamic approach. Besides the prediction of the flows over complex bathymetries, for instance, over asymmetrical dunes, by a hydrodynamic approach is shown to be superior in accuracy to the hydrostatic simulation.

hydrostatic, non-hydrostatic, TVD scheme, non-orthogonal grid

Introduction

The flows in natural environments are often characterized by a vertical confinement between a solid boundary beneath and a free surface above the flow. These kinds of flows are common in oceanic shelves, estuaries, and rivers and they are generally referred to as shallow flows. The prediction of shallow flows is important for practical applications in coastal and civil engineering in the designs of nautical constructions, the management of waterways and for the prevention of natural disasters as floods and accidental spills.

Numerical modeling is one of the primary tools used in engineering practice for the prediction of shallow flows. The shallow flow models are mostly based on the hydrostatic assumption for the pressure distribution which allows considerable simplifications[1]. An advantage of these simplifications is the high efficiency of the models because of the neglecting of the vertical momentum equation and the relatively low computation costs. However, there are situations when the assumption about the hydrostatic pressure distribution might not be valid for physical processes and the computations would involve large errors.

To improve the predictions of shallow flows, models were developed to account for a non-hydrostatic pressure. Casulli[1,2]proposed a model with full hydro-dynamic pressure terms, in which the pressure is represented as a sum of hydrostatic and non-hydrostatic constituents. Jankowski[3]presented the details of a predictor-corrector method for calculation of the pressure field. In this method, the predictor step is used to calculate the hydrostatic pressure while the non-hydrostatic pressure is calculated in the following corrector step. Li[4]simulated the water wave by directly solving Navier-Stokes equations using a similar decomposition scheme for the pressure and adopting a fractional method to implement the numerical model. Kocyigit et al.[5]and Chen[6]solved the three-dimensional non-hydrostatic equations in the Cartesian coordinate system. Fringer et al.[7]presented a non-hydrostatic model for ocean flows based on unstructured grids with the finite volume method (FVM) and implemented it by highly efficient parallel computational codes. Because of these improvements the non-hydrostatic modelings become feasible for large scale characteristics of flows in natural environments.

An important aspect of the hydrodynamic modeling is the accurate representation of the flow domain with computational grids. Although Cartesian grids used in early models are highly efficient, they are unsuitable for the accurate fitting of complex boundaries with a stair-stepped resolution. The fitting of boundaries is more accurate with curvilinear orthogonal structured grids. Most flexible are fully unstructured grids which were implemented in more recent models[2,6-8].

This short outline of the state of art in the modeling of shallow flows indicates that a significant progress has been achieved in modeling complex environmental flows over the last decade. However, the applications of advanced non-hydrostatic models are still few and far between due to a huge amount of computational resources[7]. Quite often in solving practical problems engineers have to find a compromise between the accuracy and the costs of modeling. There is evidently a knowledge gap, when a more advanced but costly non-hydrostatic models have to be applied.

The present study examines two typical shallow flows with both the hydrostatic and the non-hydrostatic models and compares the results in order to clarify the essential features in or not in favor of the costly non-hydrostatic methods. For this purpose, a previously coded in-house numerical model[9]is expanded to a fully hydrodynamic model with an unstructured horizontal grid and a vertical sigma-level grid. This expansion of the model and the application of the same basic codes suit best the aims of comparison between the hydrostatic and non-hydrostatic approaches. Both models have been used for examining the two flow cases: the non-linear wave propagation over a bar and a flow over an asymmetrical dune.

1. Numerical model

1.1 Governing equationsIn the above equations,g is the gravity acceleration, f=2ωsinφis the Coriolis term,φis the latitude, ωis the Earth’s angular velocity, andζis the free surface elevation. When the non-hydrostatic pressure pnis ignored, the vertical momentum equation should also be neglected, and the system (1)-(4) then degenerates to the conventional shallow water equations with two horizontal momentum equations and the continuity equation.

1.2 Turbulence model

To close the system the eddy viscosity concept is employed in the model and the vertical eddy viscosity coefficient νtVis determined by solving the one-equation Spalart-Allmaras (SA) model[11]. The transport equation forν～is

and dis the distance to the solid wall. The model constants include

The horizontal eddy viscosity coefficient is determined by a Smagorinsky-type formulation

where chis an arbitrary Smagorinsky constant, varying in the range of 0.01-0.5[12].

1.3 Numerical scheme

1.3.1 Outline of the numerical scheme

In the numerical method, a semi-implicit scheme is used, with a parameterθ. The importance of the value of θwas shown in previous papers[5,13]. The value of θis taken as 0.5 in the following computations. The computational domain is meshed by the unstructured grids and the flow variables are taken at the cell centers thus giving the method a name of the CC scheme. For the advection terms, the TVD (Total Variation Diminishing) scheme is used to calculate the cell face flux. An arbitrary variable ϕfat the cell face is calculated with a second order TVD scheme

where f denotes the cell face,ϕDand ϕCare, respectively, the cell-centered variables in the downwind and upwind nodes around the facef. Using a flux limiter functionψ(rf), which is simply a linear function ofrf, one can obtain a higher order TVD scheme. The ratiorfis calculated by the method extended by Darwish[14].

On the cell faces a non-orthogonal horizontal local coordinate(ξ,η)is used to replace the Cartisan coordinate(x,y )for the derivative calculations. From the chain rule of differentiations, one obtains

where J=xξyη-xηyξ, andξis directed from a control cell center to a neighbor cell center across the common face, andηis directed from one vertex to another on the same face anticlockwisely around the control cell center (Fig.1).

1.3.2 Predictor step with hydrostatic pressure

Firstly, the continuity equation is discretized as

and the discretized momentum equations are written as

whereF is an operator that includes the explicit discretization of the convective terms, the horizontal viscosity terms and the Coriolis term.

Setting q∗,q∗,q∗,q∗,η∗as the temporary

xyzσ flow variables controlled by the hydrostatic pressure ph, the discretized continuity and momentum equations are written as:

Furthermore, discretizing the above equations, one obtains

Hence, the governing equations are summarized in matrix notation as:

The details of the matrixes are

Aiis a series of tri-diagonal matrixes.

For the discretizated Eqs.(23) and (24), the intermediate variablesare introduced, so

Rewriting the governing Eq.(22)

Integrating (31) in the horizontal control volume based on Gauss theorem

where the symbolf represents the cell face,Δsiis the horizontal area ofithcell,Δlis the length of is the sthface (NS is the total number of faces), and (cosαis,sinαis)is the x -and y-projections of the face normal unit. The derivative calculation on the cell face in the Cartisan coordinate(x,y)is transformed to the local coordinate (ξ,η). The final discretized equation of (32) is

Equation (33) is simplified in a compact form as

The matrix system (34) is solved using the Bi-CGSTAB method, and,can be calculated by

solving the tri-diagonal matrix systems (23) and (24). If the non-hydrostatic pressure is neglected, the new variables are justly the final flow variables at time step n+1. Before the further step for the nonhydrostatic flow, the tri-diagonal system (25) is firstly solved to obtain the variable.

1.3.3 Corrector step with non-hydrostatic pressure

The fully discretized equations at time stepn+1 following the resolved hydrostatic fluid flow can be written as:

By subtracting Eqs.(16), (17) and (18) from (35), (36) and (37), one obtains the following discretized equations in the corrector step for the flow induced by the non-hydrostatic pressure pn

Although the predicted flow variables satisfy the depth-averaged continuity equation, it does not necessarily satisfy the local continuity condition. Therefore, the non-hydrostatic pressure enforces the corrected flow to satisfy the continuity condition. In theσtransformed frame, the continuity Eq.(1) contains no vertical velocity qz, though another form of this equation containing qzcan be deduced as

A system of Poisson-type equations for the nonhydrostatic pressure pnn+1is derived by substituting Eqs.(38), (39) and (40) into (41),

where the cross-derivative terms of pnn+1are ignored. The final discretized matrix system for Eq.(42) is summarized as

At the solid wall and the bottom, a zero normal gradient is imposed while at the free surface a zero value of the non-hydrostatic pressure is specified. When the Poisson-type equation is solved, the updated flow velocities are calculated.

Finally the updated water surface is calculated at the corrector step. Casulli[2]presented a method of water surface updating after the second corrector procedure, and Zhang[9]implemented a third computational step similarly to the first for a more accurate update of the water surface and the velocities.

2. Case study and analysis

Two test cases are considered in this computational study with the focus on: (1) calibration of the numerical models, (2) validation of the pressure terms, (3) comparison between results of hydrostatic and non-hydrostatic modelings, and (4) analysis of the necessity of the non-hydrostatic modeling. The test cases include a nonlinear wave propagation over a bar, and a flow over an asymmetric dune.

2.1 Nonlinear wave propagation over a bar

In the numerical experiments the wave is generated by adding source terms in the governing equation. To absorb the incident waves in the simulations a dampening layer is set on the beach side of the computational domain. The quadrilateral horizontal mesh resolution is 0.025 m, and the vertical first grid pointnearest to the bottom is adjusted to be about 1.5(=yu/ν), which is required by the non-slip wall condition. Then the vertical grid is stretched with a ratio 1.15 upward until about 0.025 m equal to the horizontal mesh size. The integration time step is 0.0001s, and the total simulation takes about30T (wave peroid).

Measured and computed time series of the water surface elevation are compared for seven characteristic locations, as shown in Fig.3. Both models simulate the approaching waves quite accurately (Location 1) even when the waves become slightly asymmetrical due to the interaction with the bar (Location 2). As the shape of the waves becomes more sharp-crested (Locations 3-7), the computations of the hydrostatic model near the crests start to develop spurious oscillations while the non-hydrostatic model is still capable of the accurate simulation of the wave’s pattern.

2.2 Flow over dune

Subaqueous dunes are the micro-scale features of the natural environments which are common on marine beaches, estuarine bars, and on the riverbeds. Dunes are asymmetric morphological structures with steep lee sides. They affect the flow similarly to steps and the flow structure downstream their crests often exhibits a separation of the flow, a recirculation zone, and a reattachment of the flow on the stoss side. Therefore, our second case study concerns the modeling of the flow structure over a dune in the experiments of Balachandar et al.[17]In these experiments a train of dunes identical in geometry are physically modeled in a laboratory flume and measurements are performed over the 17th dune in the sequence. These experiments become popular for modeling and were already examined in several modeling studies[18-20].

In this study the computational domain is restricted to a single dune because of the periodicity of the flow structure in a train. The geometry of the dune is shown in Fig.4. The coordinate system withx in the main stream direction,y in the spanwise direction, andz in the upward direct4ion is adopted here. The Reynolds number is 5.7×10 calculated based on the water depth and the free-surface velocity at the inlet of the domain. The Froude number is 0.44 based on the same flow parameters. The three-dimensional domain covers the dune within an area of 0.4 m long and 0.2 m wide. The horizontal grid is uniform with a scale of 0.0025 m while the vertical grid is stretched from the bottom with a ratio 1.15. The point nearest to the bottom is set within the viscosity boundary layer of about 1.6z+. At the inlet, the discharge and the eddy viscosity are provided by another simulation in a long channel with a flat bottom to get a fully developed turbulent flow. A zero water level boundary condition is set at the outlet. Simulations are carried out with both hydrostatic and non-hydrostatic models. The results of modeling are compared for six characteristic locations in vertical profiles of the stream velocities, as shown in Fig.4.

The comparison between simulations and measurements shows that the performance of the hydrostatic modeling is especially poor in the area of flow separation and recirculation＜12,-0.2＜z/L＜0.1, as shown in Fig.5. On the other hand, the results of the non-hydrostatic modeling see a good agreement with the measurements everywhere in the computational domain. It is worth noticing that the non-hydrostatic modeling is capable of accurately resolving the separation flow with the recirculation. The patterns of streamlines simulated with hydrostatic and non-hydrostatic models are compared in Fig.5. The streamlines of the flow in the hydrostatic modeling remain parallel to the riverbed and indicate neither the flow separation nor recirculation. The pattern of recirculation is well captured by the non-hydrostatic model and the free surface sees an increase in pressure, which is commonly reported for flows over dunes. It may be concluded that the accuracy of the model with the nonhydrostatic pressure is superior in comparison with the hydrostatic modeling. However, in the downstream part on the stoss side of the wave (＞12)the results of modeling with both approaches agree well with the measurements, as shown in Fig.6.

Results of this case study suggest that a switch from the hydrostatic modeling to the non-hydrostatic one is necessary for micro-scale bedforms with steep angles of lee sides. The precondition for the non-hydrostatic modeling is the necessity to resolve the lower part of the water column with complex patterns of flow separation and recirculation. The knowledge of these complex features in the flow structure is vital for the modeling of morphological processes.

3. Conclusions

A three-dimensional numerical model based on the hydrostatic assumption is extended to a fully hydrodynamic model with additional non-hydrostatic pressure terms. Two cases are investigated in order to validate both hydrostatic and non-hydrostatic models with the focus on the impact of the non-hydrostatic pressure for the numerical simulation of shallow flows.

The results of the study show that for waves propagating along the bar a criterion of shallowness alone, the ratio between the depth and the wave length, is insufficient to warrant the performance of the shallow flow equations with a hydrostatic approach and the dispersion in the wave dynamics can be better accounted with the non-hydrostatic modeling approach.

Beyond the laboratory scale, the hydrostatic model is widely used in the simulation of natural flows in oceans, estuaries and rivers because of the lower computation cost and the higher efficiency. However, the additional calculation of the Poisson equation for the non-hydrostatic pressure costs much more CPU time, especially, for natural flows in a large domain (commonly on the kilometer scale). A further development should probably focus on the development of a hybrid hydrostatic and non-hydrostatic model based on a series of practical and credible switch conditions.

Acknowledgement

This work was supported by the Deutsche Forschungsgemeinschaft (Grant No. DFG SU 405/3, SU 405/4).

[1] CASULLI V. A semi-implicit finite fifference method for non-hydrostatic, free-surface flows[J]. International Journal for Numerical Methods in Fluids, 1999, 30(4): 425-440.

[2] CASULLI V., ZANOLLI P. Semi-implicit numerical modeling of nonhydrostatic free-surface flows for environmental problems[J]. Mathematical and computer modeling, 2002; 36(10): 1131-1149.

[3] JANKOWSKI J. A. A non-hydrostatic model for free surface flows[D]. Doctoral Thesis, Hannover, Germany: University of Hannover, 1999.

[4] LI B., FLEMING C. A. Three-dimensional model of Navier-Stokes equations for water waves[J]. Journal of Waterway, Port, Coastal and Ocean Engineering, 2001, 127(1): 16-25.

[5] KOCYIGIT M. B., FALCONER R. A. and LIN B. Three-dimensional numerical modeling of free surface flows with non-hydrostatic pressure[J]. International Journal for Numerical Methods in Fluids, 2002, 40(9): 1145-1162.

[6] CHEN X. J. A fully hydrodynamic model for three-dimensional, free-surface flows[J]. International Journal for Numerical Methods in Fluids, 2003, 42(9): 929-952.

[7] FRINGER O. B., GERRITSEN M. and STREET R. L. An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator[J]. Ocean modeling, 2006, 14(3): 139-173.

[8] CASULLI V., ZANOLLI P. High resolution methods for multidimensional advection-diffusion problems infree-surface hydrodynamics[J]. Ocean modeling, 2005, 10(2): 137-151.

[9] ZHANG Jing-xin, LIU Hua and XUE Lei-ping. A vertical 2-D mathematical model for hydrodynamic flows with free surface in σcoordinatet[J]. Journal of Hy- drodynamics Ser. B, 2006, 18(1): 82-90.

[10] PHILLIPS N. A. A coordinate system having some special advantages for numerical forecasting[J]. Journal of Meteorology, 1957, 14: 184-185.

[11] SPALART P. R. Strategies for turbulence modelling and simulations[J]. International Journal of Heat and Fluid Flow, 2000, 21(3): 252-263.

[12] ZHANG Q. Y., CHAN E. S. Sensitivity studies with the three-dimensional multi-level model for tidal motion[J]. Ocean Engineering, 2003, 30(12): 1489-1505.

[13] CHEN C., LIU H. and BEARDSLEY R. C. An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[J]. Journal of Atmospheric and Oceanic Technology, 2003, 20(1): 159-186.

[14] DARWISH M. S., MOUKALLED F. TVD schemes for unstructured grids[J]. International Journal of heat and Mass Transfer, 2003, 46(4): 599-611.

[15] BEJI S., BATTJES J. A. Experimental investigation of wave propagation over a bar[J]. Coastal Engineering, 1993, 19(2): 151-162.

[16] OHYAMA T., BEJI S. and BATTJES J. A. Experimental verification of numerical model for nonlinear wave evolutions[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1994, 20(6): 637-644.

[17] BALACHANDAR R., POLATEL C. and HYUN B.-S. et al. LDV, PIV, and LES investigation of flow over a fixed dune[C]. Sedientation and Sediment Transport: Proceedings of the Symposium Held in Monte Verta. Monte Verita, Switzerland, 2002, 171-178.

[18] YUE W., LIN C. L. and PATEL V. C. Large eddy simulation of turbulent open-channel flow with free surface simulated by level set method[J]. Physics of fluids, 2005, 17(2): 1-12.

[19] YUE W., LIN C. L. and PATEL V. C. Large-eddy simulation of turbulent flow over a fixed two-dimensional dune[J]. Journal of Hydraulic Engineering, ASCE, 2006, 132(7): 643-651.

[20] BALEN W. V., UIJTTEWAAL W. S. J. and BLANCKAERT K. Large-eddy simulation of a curved openchannel flow over topography[J]. Physics of Fluids, 2010, 22(7): 1-17.

10.1016/S1001-6058(14)60058-5

* Project supported by the National Natural Science Foundation of China (Grant No.10702042), the Non-profit Industry Financial Program of MWR (Grant No. 201401027) and the National Key Basic Research Development Program of China (973 Program, Grant No. 2014CB046200).

Biography: ZHANG Jing-xin (1975-), Male, Ph. D.,

Associate Professor