Two-dimensional shallow water equations with porosity and their numerical scheme on unstructured grids

Zhi-li WANG, Yan-fen GENG*

1. State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing Hydraulic Research Institute, Nanjing 210024, P. R. China

2. School of Transportation, Southeast University, Nanjing 210096, P. R. China

Two-dimensional shallow water equations with porosity and their numerical scheme on unstructured grids

Zhi-li WANG1, Yan-fen GENG*2

1. State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing Hydraulic Research Institute, Nanjing 210024, P. R. China

2. School of Transportation, Southeast University, Nanjing 210096, P. R. China

In this study, porosity was introduced into two-dimensional shallow water equations to reflect the effects of obstructions, leading to the modification of the expressions for the flux and source terms. An extra porosity source term appears in the momentum equation. The numerical model of the shallow water equations with porosity is presented with the finite volume method on unstructured grids and the modified Roe-type approximate Riemann solver. The source terms of the bed slope and porosity are both decomposed in the characteristic direction so that the numerical scheme can exactly satisfy the conservative property. The present model was tested with a dam break with discontinuous porosity and a flash flood in the Toce River Valley. The results show that the model can simulate the influence of obstructions, and the numerical scheme can maintain the flux balance at the interface with high efficiency and resolution.

shallow water equations with porosity; source term; Roe-type Riemann solver; finite volume method; unstructured grid

1 Introduction

Free surface flows exist in estuaries, along coasts, and in river and lake regions. In most free surface flows, the hydrostatic pressure assumption is generally valid. Situations where the hydrostatic pressure assumption may be questionable have been discussed in several previous studies (Heggelund et al. 2004; Yuan and Wu 2004; Lee et al. 2006). With the hydrostatic pressure assumption and Boussinesq approximation, the shallow water equations can be obtained from the Navier-Stokes (NS) equations. The shallow water equations have been widely applied in ocean and hydraulic engineering, including simulations of tidal flows in estuary and coastal water regions, wave propagation, stationary hydraulic jumps, and open channel flows (Wang et al. 2005b; Lu et al. 2005; Ding et al. 2004 ).

The shallow water equations are more simple than the NS equations, and they reduce the numerical calculations for free surface flows enormously, so the shallow water equations aresuitable for the numerical simulation of water flow in large-scale oceans, estuaries, rivers, and lakes (Lu et al. 2005; Wang 2005; Zhou et al. 2001). Numerical models of shallow water equations are effective for simulation of free surface flows, but when there are islands, buildings, or other structures in the numerical region, denser grids, which increase the numerical computation, are needed. Furthermore, the time step is often determined by the smallest grid size, which further reduces the efficiency of the model. In order to overcome these problems, we introduce the porosity to reflect the effects of obstructions and derive two-dimensional porous shallow water equations based on the laws of conservation of mass and momentum.

Recently, the Godunov-type numerical models have been widely used in computational aerodynamics (Harten et al. 1983; Li 2008) and computational hydraulics (Zoppou and Roberts 2000; Wang et al. 2005a; Geng et al. 2005), since they have several desirable properties, including being monotone and conservative, and having good shock-capturing capabilities with a correct shock speed value and an inherent upwind property. In the well-known work of Godunov, in order to reduce computational time, the exact solution of the Riemann problem was replaced by an approximate solution such as the Roe-type solver (Roe 1981), Osher-type solver (Osher and Solomon 1982), and HLL-type solver (Harten et al. 1983). In this study, the finite volume method and Roe-type approximate Riemann solver were used for the discretization of the two-dimensional porous shallow water equations, in which the source terms of the bed slope and porosity were treated by the local characteristic decomposition and upwind fashion to satisfy equilibrium and steady-state conditions.

2 Two-dimensional shallow water equations with porosity

The continuity and momentum equations with porosity are derived based on the laws of conservation of mass and momentum. There is an infinitesimally small control volume with a length of Δx and a width of Δy in the horizontal plane and a water depth h in the vertical direction. If there is a solid structure with a length of Δx′ and a width of Δy′ within the control volume, we can define the porosity as

The value of φ lies between 0 and 1, φ=1 means no solid structures in the control volume, and φ=0 means no water in the control volume. In this study, we assume that the porosity on the sides of the control volume is equal to φ.

2.1 Continuity equation

The total mass m of water in the control volume is

where ρ is the water density.

The mass fluxesQwxandQsyacross the western and southern sides flowing into the control volume are, respectively,

whereuandvare the depth-averaged velocities in thexandydirections, respectively.

The mass fluxesQexandQnyacross the eastern and northern sides flowing out of the control volume are, respectively,

Based on the law of conservation of mass, we have

Substituting Eqs. (2) through (6) into Eq. (7), and assuming that water is not compressible, which means thatρis constant, the conservative form of the continuity equation can be obtained:

2.2 Momentum equation

To save space, the momentum equation in thexdirection is derived in detail only. The total momentum of water in thexdirection in the control volume is

The mass forces across the western, eastern, southern, and northern sides of the control volume are, respectively,

The hydrostatic pressures exerted on the western and eastern sides of the control volume are, respectively,

Bottom pressure in the x direction due to the bottom slope can be written as

The friction resistance at the bottom is accounted for by a classical Stickler law, and can be written as

The resistance caused by the obstructions, which includes the configuration resistance and friction resistance, is assumed to be identical over the entire flow region and to be proportional to the square of the velocities (Wang et al. 2008):

Based on the law of conservation of momentum (Newton’s second law) in the control volume, we obtain the following expression:

Substituting Eqs. (9) through (16) into Eq. (17) and dividing by ρ, the momentum equation in the x direction can be obtained:

The momentum equation in the y direction can be obtained in a similar way:

The continuity equation (Eq. (8)), and the momentum equations (Eqs. (18) and (19)), constitute a closed system of a two-dimensional depth-averaged shallow water model with porosity. Generally, a vector form of the model can be expressed as

3 Numerical solution

Eq. (20) is discretized using the finite volume method on unstructured grids. The average variables are stored at the center of each grid cell, and the edges of each grid cell are defined as the faces of control volume.

The domain is paved with a set of non-overlapping polygonal cells: {Ωi,i=1,2,…,I}, where i is the serial number of cells and I is the total number of cells. A polygonal cell Ωιis built with vertices labeled as mi,k, where k is the serial number of vertices, and k=1,2,…,Ei; and Ειdenotes the number of vertices of the cell Ωι. The two cells that share the jth side of the grid are identified by the indices j1and j2. If the jth side is the boundary of the computational region, we set j2=?1.

Integrating Eq. (20) over the cell Ωιyields

where F=(E, H), Αιis the area of the cell Ωι, and s is the integration variable. Using the Green formula, we obtain

The second term of Eq. (23) can be rewritten as

wherelijis the length of thejth side of the cellΩι;Fn=F·n=Enx+Hny; andFnijis the flux through thejth side of the cellΩι, which is discretized with a Roe-type approximate Riemann solver.

3.1 Flux computation

The fluxFndepends on the conservational variableUand porosityφ. With the variables stored at the cell center,Uandφare discontinuous at the cell boundary, which is known as the Riemann problem:

wheren0is the coordinate along the outward normal direction of the cell boundary oriented from the center of the cell boundary,t0is the initial time, and the variables with the subscripts L and R denote the values of the variables on the left and right sides of the interface between two cell, respectively. At thejth side of the cellΩι, ifUL=Uj1andUR=Uj2(whereUj1andUj2are the values of theconservational variableUat the cellsj1andj2that share thejth side, respectively), the numerical scheme is only first order in space. If we assume thatUis linearly varied in space, the second-order total variation diminishing (TVD) numerical scheme can be obtained with monotonic upstream schemes for conservation laws (van Leer 1979; Wang et al. 2005b).

Using the matrix theory, the matrixAcan be diagonalized as

For the Riemann problem of Eq. (26), the flux through the left side of the interface can be written as (Tan 1998)

where Δ(?) denotes the difference between the values of variables on the left and right sides of the interface. Substituting Eqs. (28) and (32) into (30), the flux through the left side of the interface can be obtained:

3.2 Computation of bed slope and porosity source term

The bed slope and porosity source termS0is a function of water depthh, the bed elevationzb, and porosityφ.h,zb, andφare discontinuous at the cell boundary and can be respectively expressed as

where

The upwind method is used for the discretization of the bed slope and porosity source terms (Wang et al. 2008; Bermudez et al. 1998). The bed slope and porosity source terms on the left and right sides of the cell boundary are, respectively,

whereIis the unit matrix, andγ=(γ1,γ2,γ3).

3.3 C-property relative to a stationary solution

The centered discretization of bed slope source terms gives rise to spurious waves (Wang et al. 2008). Bermudez and Vazquez (1994) proposed a conservation property (C-property) which prevents the appearance of spurious waves. The C-property characterizes the accuracy of a numerical scheme used for approximating a steady state solution representing water at rest. The steady state solution is characterized by

whereηis the water level.

To keep water at rest, the discharge flux must be zero:

From this we can get the average porosity at the cell face:

Under the rest water conditions (38), Eq. (20) can be simplified as

where the terms on the left side of Eq. (42) are the hydrostatic pressure, and the first and second terms on the right side are the bed slope source term and additional porosity source term, respectively. If the numerical scheme satisfies Eq. (42), we can say that it satisfies the C-property. Since the bed elevation is constant, from a mathematical point of view, if the numerical model satisfies the C-property, the momentum flux must equal the hydrostatic pressure flux:

3.4 Wet/dry fronts

3.5 Stability constraint

Since the numerical scheme is explicit in time, the time step is limited by the following Courant-Friedrichs-Lewy (CFL) condition (Cea et al. 2006):

4 Model test

Two applications of the model are presented for validation purposes. First, a dam break in a channel with discontinuous porosity was simulated, and the calculated results were compared with analytical results. Second, a flash flood in the Toce River Valley was simulated, and the calculated results were compared with the experimental results.

4.1 Dam break with discontinuous porosity

Initially, a gate was placed atx= 100 m and water was kept at rest. The water depths upstream and downstream of the gate wereh1= 10 m andh2= 1 m, respectively (Fig. 1). The gate was suddenly removed and a large volume of water was released. The numerical model used uniform rectangle grids with 200 computational cells and a time steptΔ of 0.04 s. Fig. 2 depicts the numerical and analytical water depths and velocities att= 0.3 s. The agreement between the computed and analytical solutions is seen to be quite good. The propagation speeds of the various waves are computed correctly.

Fig. 1Dam break with discontinuous porosity (Unit: m)

Fig. 2Numerical and analytical results of dam break with discontinuous porosity att= 0.3 s

4.2 Flash flood in Toce River Valley

Testa et al. (2007) presented a physical model study of a flash flood in the Toce River Valley, in Italy, which impacted an idealized urban district composed of an array of square blocks (Fig. 3). The aim of this test was to show that the shallow water model on fine grids could be advantageously replaced with the porous shallow water model on much coarser grids, where the effects of the blocks were modeled using porosity. There were two block configurations: aligned and staggered as in the physical experiments performed by Testa et al. (2007), but only the staggered layout (Fig. 3) was used for testing here.

Fig. 3Gauge stations with staggered arrangement

Two numerical models were built to simulate the physical model. The first model was a classical shallow water (CSW) model (Wang et al. 2005) on fine grids (Fig. 4(a)), in which the blocks were treated as the impermeable boundaries. The second model was the proposed porous shallow water (PSW) model on coarse grids, in which the effects of the blocks were modeled by the porosity (Fig. 4(b)). The CSW model and PSW model were made of 4 567 and 1 696 triangle cells, respectively. For wet and dry fronts, the values ofht1andht2were 0.000 1 m and 0.001 m, respectively. A uniform Manning’s roughness coefficientn=0.016 2 m-1/3?s was assigned to all cells to model bottom shear. In the PSW model, the approach for vegetative resistance modeling (Nepf 1999) was used to calculate the drag force of the blocks. The drag forceRsexerted by the blocks on the fluid is

whereRs=(Rsx,Rsy), withRsxandRsybeing thexandycomponents ofRs, respectively;u=(u,v);dis the projected length of the blocks, i.e., the length of the blocks as seen by an observer moving in the direction of flow; andCDis a bulk drag coefficient, which is a function of blocks’ density (Nepf 1999). The value ofCDfor square-shaped obstructions is tabulated in textbooks asCD=2.0 (Munson et al. 2006). Based on these considerations, the head loss coefficientCsin Eq. (21) can be calculated as follows:

Fig. 4Numerical meshes

The discharge hydrograph lasting 60 seconds, obtained from the experiments performed by Testa et al. (2007), was used as the upstream boundary condition (Fig. 5). A zero-order extrapolation or so-called soft boundary condition was implemented at the downstream boundary (Wang et al. 2005). Fig. 6 shows time series of water depth at gauge stations 3, 5, 7, and 10. Predictions of the CSW model (Wang et al. 2005) and PSW model are shown, along with water depth measurements reported by Testa et al. (2007). Fig. 6 shows that the results of both the CSW model and the PSW model match the experimental results. At gauge stations 3 and 7, the two models’ predictions are both close to the experimental results. At gauge station 5, compared with the experimental results, predictions of the two models are both underestimates. At gauge station 10, the CSW model overestimates the water depth while the PSW model underestimates the water depth.

Fig. 5Discharge hydrograph at upstream boundary

Fig. 6 Predicted and measured water depths at different gauge stations

The PSW model predictions match the CSW model predictions, but it should be noted that the CPU time required by the 60-second simulation with the Pentium 4 processor is 18 seconds for the PSW model and 176 seconds for the CSW model. This example shows that the PSW model can simulate the influence of blocks with high efficiency.

5 Conclusions

Through introducing the porosity to reflect the effects of obstructions, the two-dimensional porous shallow water equations including the continuity and momentum equations were derived based on the laws of conservation of mass and momentum. The unstructured finite volume method and modified Roe-type approximate Riemann solver were used for the solution of the two-dimensional porous shallow water equations. The bed slope source term and the additional porosity source term were both decomposed in the characteristic direction. It has been shown that the numerical scheme exactly satisfies the conservative property. Numerical results of a dam break in a channel with discontinuous porosity and a flash flood in the Toce River Valley with an urban district show that the porous shallow water equations and numerical scheme can simulate the influence of blocks with high efficiency and resolution.

Alcrudo, F., and Benkhaldoun, F. 2001. Exact solutions to the Riemann problem of the shallow water equations with a bottom step. Computers and Fluids, 30(6), 643-671. [doi:10.1016/S0045-7930 (01)00013-5]

Bermudez, A., Dervieux, A., Desideri, J. A., and Vazquez, M. E. 1998. Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes. Computer Methodsin Applied Mechanics and Engineering, 155(1-2), 49-72. [doi:10.1016/S0045-7825(97)85625-3]

Bermudez, A., and Vazquez, M. E. 1994. Upwind methods for hyperbolic conservation laws with source terms. Computers and Fluids, 23(8), 1049-1071. [doi:10.1016/0045-7930(94)90004-3]

Cea, L., French, J. R., and Vazquez-Cendon, M. E. 2006. Numerical modelling of tidal flows in complex estuaries including turbulence: An unstructured finite volume solver and experimental validation. International Journal for Numerical Methods in Engineering, 67(13), 1909-1932. [doi:10.1002/ nme.1702]

Ding, L., Pang, Y., Zhao, D. H., Wu, J. Q., and Lu, J. 2004. Analysis of applicability of flux difference splitting scheme on 2D flow-pollutant calculation. Advanced in Water Science, 15(5), 561-565. (in Chinese)

Geng, Y. F., Wang, Z. L., and Jin, S. 2005. A high resolution Godunov-type scheme for one dimensional shallow water flow. Journal of Hydrodynamics, Ser. A, 20(4), 507-512 (in Chinese)

Harten, A., Lax, P. D., and van Leer, B. 1983. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25(1), 35-61. [doi:10.1137/1025002]

Heggelund, Y., Vikebo, F., Berntsen, J., and Furnes, G. 2004. Hydrostatic and non-hydrostatic studies of gravitational adjustment over a slope. Continental Shelf Research, 24(18), 2133-2148. [doi: 10.1016/j.csr.2004.07.005]

Lee, J. W., Teubner, M. D., Nixon, J. B., and Gill, P. M. 2006. A 3-D non-hydrostatic pressure model for small amplitude free surface flows. International Journal for Numerical Methods in Fluids, 50(6), 649-672. [doi:10.1002/fld.1054]

Li, Z. W. 2008. Study on the dissipative effect of approximate Riemann solver on hypersonic heatflux simulation. Chinese Journal of Theoretical and Applied Mechanics, 40(1), 19-25. (in Chinese)

Lu, Y. J., Zhuo, L. Q., Shao, X. J., Wang, H. C., and Li, H. L. 2005. A 2D mathematical model for sediment transport by waves and tidal currents. China Ocean Engineering, 19(4), 571-586.

Munson, B. R., Young, D. F., and Okiishi, T. H. 2006. Fundamentals of Fluid Mechanics. 2nd ed. New York: John Wiley & Sons.

Nepf, H. M. 1999. Drag, turbulence and diffusion in flow through emergent vegetation. Water Resources Research, 35(2), 479-489. [doi:10.1029/1998WR900069]

Osher, S., and Solomon, F. 1982. Upwind difference schemes for hyperbolic system of conservation laws. Mathematics of Computation, 38(158), 339-374.

Roe, P. L. 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43(2), 357-372. [doi:10.1016/0021-9991(81)90128-5]

Shi, H. D., and Liu, Z. 2005. A finite volume method with unstructured triangular grids for numerical modeling of tidal current. China Ocean Engineering, 19(4), 693-700.

Sleigh, P. A., and Gaskell, P. H., Berzins, M., and Wright, N. G. 1998. An unstructured finite-volume algorithm for predicting flow in rivers and estuaries. Computers and Fluids, 27(4), 479-508.

Stoker, J. J. 1957. Water Waves: The Mathematical Theory with Applications. New York: Wiley-Interscience.

Tan, W. Y. 1998. Computational Hydraulics-Finite Volume Method. Beijing: Tsinghua University Press. (in Chinese)

Testa, G., Zuccala, D., Alcrudo, F., Mulet, J., and Soares-Frazao, S. 2007. Flash flood flow experiment in a simplified urban district. Journal of Hydraulic Research, 45(s1), 37-44. [doi:10.1080/00221686. 2007.9521831]

van Leer, B. 1979. Towards the ultimate conservative difference scheme, V: A second-order sequel to Godunov’s method. Journal of Computational Physics, 32(1), 101-136. [doi:10.1016/0021-9991 (79)90145-1]

Wang, J. S., Ni, H. G., and He, Y. S. 2000. Finite-difference TVD scheme for computation of dam-break problems. Journal of Hydraulic Engineering, 126(4), 253-262. [doi:10.1061/(ASCE)0733-9429 (2000)126:4(253)]

Wang, Z. L. 2005. The Unstructured 2D and 3D Shallow Water Model Study Based on Godunov and Semi-Lagrangian Method. Ph. D. dissertation. Dalian: Dalian University of Technology. (in Chinese)

Wang, Z. L., Geng, Y. F., and Jin, S. 2005a. An unstructured finite volume algorithm for nonlinear two-dimensional shallow water equation. Journal of Hydrodynamics, Ser. B, 17(3), 306-312.

Wang, Z. L., Geng, Y. F., and Jin, S. 2005b. Flux balance method for shallow water equation with source terms. Advanced in Water Science, 16(3), 373-379. (in Chinese).

Wang, Z. L., Lu, Y. J., and Geng, Y. F. 2008. One dimensional shallow water equations with porosity and their numerical discretization schemes. Chinese Journal of Theoretical and Applied Mechanics, 40(5), 585-592. (in Chinese).

Yuan, H. L., and Wu, C. H. 2004. A two-dimensional vertical non-hydrostatic σ model with an implicit method for free-surface flows. International Journal for Numerical Methods in Fluids, 44(8), 811-835. [doi: 10.1002/fld.670]

Zhao, D. H., Shen, H. W., Tabios, III, G. Q., Lai, J. S., and Tan, W. Y. 1994. Finite-volume two-dimensional unsteady-flow model for river basins. Journal of Hydraulic Engineering, 120(7), 863-883. [doi: 10.1061/(ASCE)0733-9429(1994)120:7(863)]

Zhou, J. G., Causon, D. M., Mingham, C. G., and Ingram, D. M. 2001. The surface gradient method for the treatment of source terms in the shallow water equation. Journal of Computational Physics, 168(1), 1-25. [doi:10.1006/jcph.2000.6670]

Zoppou, C., and Roberts, S. 2000. Numerical solution of the two-dimensional unsteady dam break. Applied Mathematical Modelling, 24(7), 457-475. [doi:10.1016/S0307-904X(99)00056-6]

(Edited by Yan LEI)

This work was supported by the National Natural Science Foundation of China (Grants No. 50909065 and 51109039) and the National Basic Research Program of China (973 Program, Grant No. 2012CB417002).

*Corresponding author (e-mail: yfgeng@seu.edu.cn)

Received Dec. 5, 2011; accepted May 9, 2012