A two-dimensional coupled flow-mass transport model based on an improved unstructured finite volume algorithm.

Environmental Research ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Environmental Research journal homepage: www.elsevier.com/locate/envres

A two-dimensional coupled flow-mass transport model based on an improved unstructured finite volume algorithm Jianzhong Zhou a,n, Lixiang Song a,b, Suncana Kursan c, Yi Liu a a

School of Hydropower and Information Engineering, Huazhong University of Science and Technology, Wuhan 430074, PR China Pearl River Hydraulic Research Institute, Guangzhou 510611, PR China c Hrvatske Vode-Legal Entity for Water Management, Ulica grada Vukovara 220, Zagreb 10000, Croatia b

art ic l e i nf o

a b s t r a c t

Article history: Received 28 September 2014 Received in revised form 31 December 2014 Accepted 27 January 2015

A two-dimensional coupled water quality model is developed for modeling the flow-mass transport in shallow water. To simulate shallow flows on complex topography with wetting and drying, an unstructured grid, well-balanced, finite volume algorithm is proposed for numerical resolution of a modified formulation of two-dimensional shallow water equations. The slope-limited linear reconstruction method is used to achieve second-order accuracy in space. The algorithm adopts a HLLC-based integrated solver to compute the flow and mass transport fluxes simultaneously, and uses Hancock’s predictor– corrector scheme for efficient time stepping as well as second-order temporal accuracy. The continuity and momentum equations are updated in both wet and dry cells. A new hybrid method, which can preserve the well-balanced property of the algorithm for simulations involving flooding and recession, is proposed for bed slope terms approximation. The effectiveness and robustness of the proposed algorithm are validated by the reasonable good agreement between numerical and reference results of several benchmark test cases. Results show that the proposed coupled flow-mass transport model can simulate complex flows and mass transport in shallow water. & 2015 Elsevier Inc. All rights reserved.

Keywords: Mass transport Coupled model Advection-diffusion Well-balanced Finite volume algorithm

1. Introduction Over the past two decades, numerical models based on the non-linear shallow water equations have been widely used for applications in ocean and hydraulic engineering, such as tidal flows in estuaries, bore wave propagation, river, lake, and open channel flows (Audusse et al., 2004; Benkhaldoun and Elmahi (2007); Xing et al., 2011; Canestrelli et al., 2010; Cai et al., 2007; Lai et al., 2010). Moreover, with the rapid development of hydrodynamic models, there have been many developments in the field of flow-mass transport modeling, as the prediction of pollutant transport in flows is an important topic in many industrial and environmental projects. Majority of researchers have focused on the development of various flow-transport models. The numerical simulation of flow-mass transport in shallow water is not an easy task, since it involves steep gradients, wetting and drying processes, complex geometry, and so on. Therefore, research on robust, consistent and stable numerical schemes for pollutant transport simulation has attracted tremendous attention n

Corresponding author. E-mail addresses: [email protected] (J. Zhou), [email protected] (L. Song), [email protected] (S. Kursan), [email protected] (Y. Liu).

in the last years (Cai et al., 2007; Cea and Vázquez-Cendón, 2012; Issa et al., 2010; Li et al., 2012; Li and Duffy, 2012; Liang et al., 2010). Obviously, simulating the pollutant transport is not only determined by the characteristics of the pollutant but also the properties of fluid flow. So numerical resolution of the shallow water equations is as important as solving the pollutant transport equation. In recent years, many numerical models have been developed to resolve the shallow water equations in conservative form by using finite volume methods. In particular, the upwind Godunovtype schemes are very popular because of their robustness and accuracy in capturing discontinuities. Although the upwind Godunov-type schemes have been researched for several years, two persistent problems dog the schemes in practice. First, the shallow water equations with non-flat topography have stationary solutions in which the flux gradients are non-zero but should be exactly balanced by the source terms. This wellbalanced concept is also known as the exact conservation property (C-property). Since Bermudez and Vazquez (1994) proposed the ideal of the “exact C-property”, the development of well-balanced Godunov-type methods for numerical resolution of the shallow water equations with complex source terms have become one of the most important research topics, and several well-balanced methods have been developed (Audusse et al., 2004; Benkhaldoun

http://dx.doi.org/10.1016/j.envres.2015.01.017 0013-9351/& 2015 Elsevier Inc. All rights reserved.

Please cite this article as: Zhou, J., et al., A two-dimensional coupled flow-mass transport model based on an improved unstructured finite volume algorithm. Environ. Res. (2015), http://dx.doi.org/10.1016/j.envres.2015.01.017i

2

J. Zhou et al. / Environmental Research ∎ (∎∎∎∎) ∎∎∎–∎∎∎

and Elmahi, 2007; Canestrelli et al., 2010; Song et al., 2011a; Song et al., 2011b; Hou et al., 2013; Lai et al., 2010). Among these researches, a set of pre-balanced shallow water equations was derived by Liang and Borthwick (2009), and a well-balanced scheme was developed based on the pre-balanced equations. However, numerical instability would be produced near the wet/dry fronts due to the negative water depth. In order to resolve this problem, Liang and Borthwick (2009) proposed the additional source terms, which were introduced to balance spurious fluxes near the wet/ dry fronts. Song et al., (2011a) proposed a new formulation of the shallow water equations, which can essentially eliminate the need for momentum flux corrections, thus providing a robust alternative to the classical shallow water equations. Based on the new formulation, a robust and well-balanced scheme was proposed for shallow flows simulation. As for the second problem, how to design a robust approach to handle wet/dry fronts is one of the central challenges that researchers are tackling. Zia and Banihashemi, 2008 employed an explicit Godunov-type scheme, which can deal with wet/dry fronts and uneven beds. Song et al., (2011a) adopted a combination of the sloping bottom model and the modified shallow water equations to provide a robust technique for wet/dry fronts tracking, and used a semi-implicit scheme to resolve the stiff problem from friction source terms, which would introduce unphysical large velocity near wet/dry fronts where the water depth is vanishing. Hou et al., (2013) introduced a 2D well-balanced shallow water flow model, which was based on an unstructured cell-centered finite volume scheme, to simulate complex flows involving wet/dry fronts over arbitrary beds, and adopted a hydrostatic reconstruction method to reconstruct non-negative water depths at wet/dry interfaces. When considering the numerical solution of pollutant transport equation, advection and diffusion are the two fundamental transport elements for solutes in water flows, and advection is the predominant factor in most practical problems. The advection is governed by hyperbolic partial differential equation, thus numerical difficulties such as contact or shocks discontinuities often occurred in the process of discretization. These difficulties must be dealt with both mathematically and computationally. Traditional method would produce large numerical damping and unphysical oscillations near the discontinuities or regions with steep gradients due to the shocks. Actually, if a numerical method uses interpolation approach to relate the values on the characteristic lines, the numerical damping and oscillations would rise badly. To deal with these problems, Benkhaldoun and Elmahi (2007) proposed an adaptive finite volume method for the numerical simulation of pollutant transport. The method used adaptive unstructured meshes, incorporates upwind numerical fluxes and slope limiters to provide sharp resolution of steep bathymetric gradients. This adaptive scheme is non-oscillatory and mass conservative. However, a large amount of computer memory is needed to store the mesh topology information, which may reduce the computational efficiency. Cai et al., (2007) proposed an upwind scheme based on the central weighted essentially non-oscillatory (CWENO) reconstructions for computing the steady and unsteady pollutant transport, and this scheme have the non-oscillatory property. Liang et al., (2010) used an alternating operatorsplitting technique to separate the spatial discretization of the advection and diffusion terms in each step of the time marching procedure. In this model, a five-point Total Variation Diminishing (TVD) modification of the standard MacCormack scheme is used in the advection stage. This measure demonstrates an attractive balance between numerical accuracy and stability. No matter what numerical method is used to deal with the advection term, a robust scheme should reduce the numerical damping and prevent numerical oscillations near the shocks. This paper aims to present a two-dimensional coupled flow-

transport model based on an incremental modification of the previously published model (Song et al., 2011a). The Hancock's predictor–corrector method is adopted for efficient time stepping because fluxes and gradients need only be evaluated once per time step. The continuity and momentum equations are updated in both wet and dry cells. The bed slope terms in a wet cell are exactly discretized based on the water surface elevation and bed slopes of the cell, whereas the DFB (divergence form for bed slope source term) technique is adopted for bed slope terms approximation in dry cell. This new hybrid method for bed slope terms approximation can preserve the well-balanced property of the algorithm for simulations involving flooding and recession. To simulate the flow-mass transport in shallow water with complex topography and geometry, a HLLC-based integrated solver is proposed to compute the flow and mass transport fluxes simultaneously, and then a two-dimensional coupled flow-mass transport model based on Godunov-type finite volume scheme is developed. The proposed coupled finite volume model has good performance on flow-mass transport simulations.

2. Governing equations The two-dimensional depth-averaged shallow water equations are applicable to many engineering hydrodynamic problems including river flows, shallow lake currents, estuarine and coastal flows, etc. Different formulations of shallow water equations are proposed and adopted by many researchers for numerical modeling (Begnudelli and Sanders, 2006; Nguyen et al., 2006; Liang and Marche, 2009; Liang, 2010; Song et al., 2011a; Lai and Khan 2012). In their paper Song et al., (2011a) proposed a new formulation of the two-dimensional shallow water equations, which can essentially eliminate the need for momentum flux corrections, and thus provides a robust alternative to the classical form of shallow water equations (Song et al., 2011a). In this paper, we adopt the new formulation combined with the convection–diffusion equation as the governing equations, which are given by

∂E ∂G ∂U + + =S ∂t ∂x ∂y

(1)

where t is the time; x and y are Cartesian coordinates; U , E , G and S are vectors of the conserved variables, fluxes in the x - and y -directions, and source terms, respectively; in which ⎡h ⎤ ⎢ ⎥ hu U=⎢ ⎥ ⋯ ⎢ hv ⎥ ⎢⎣ hc ⎥⎦ ⎤ ⎡ hu ⎥ ⎢ 2 2 2 ⎥ ⎢ E= ⎢ hu + g (h − b )/2⎥ ⋯ ⎥ ⎢ huv ⎦ ⎣ huc ⎤ ⎡ hv ⎥ ⎢ ⎥ ⎢ huv G= ⎢ 2 2 2 hv + g (h − b )/2⎥ ⎥ ⎢ ⎦ ⎣ hvc ⎤ ⎡ qin ⎥ ⎢ ⎥ ⎢ g (h + b) α − ghS fx ⎥ ⎢ S=⎢ g (h + b) β − ghS fy ⎥ ⎥ ⎢∂ ∂c ∂ ∂c ⎢ (D x h ) + (D y h ) − khc + qin c in ⎥ ∂x ∂y ∂y ⎦ ⎣ ∂x

(2)

where h is the water depth; u and v are the depth-averaged velocity components in the x - and y -directions, respectively; c is the depth-averaged contaminating material concentration; b is the bed elevation; g = 9.81 m/s2 is the acceleration due to gravity; α



and β are bed slopes in the x - and y -directions, respectively; S fx and S fy are friction slopes in the x - and y -directions, respectively; qin is the inflow discharge of the point source; cin is the contaminating material concentration of the point source; Dx and Dy are the mass diffusion coefficients in the x - and y -directions, respectively; k is the degradation coefficient. Assume that the bed is fixed, i.e. b = b (x, y) , then α = − ∂b /∂x and β = − ∂b /∂y . The friction slopes are estimated by Manning formulae

S fx =

n2u u2 + v2 h4/3

, S fy =

n2v u2 + v2 h4/3

(3)

where n is the empirical Manning roughness coefficient.

γ1 = bi,3 − 3bi,1 γ2 = 3hi bi,1 − 3hi bi,3 − bi,2 bi,3 + bi,1bi,2 + bi2, 1

3

(7)

and bi,1 ≤ bi,2 ≤ bi,3. It should be noted that the water depth is directly updated by numerical fluxes at each time step, and then the water surface elevation is updated by (6). Once the water surface elevation is updated, a cell can be classified as wet or dry. In this study, the horizontal free surface classification method proposed by Begnudelli and Sanders (2007) is adopted for wet/dry cell classification, i.e. a cell is classified as wet only when the water level is higher than the highest nodal bed elevation, and otherwise the cell is classified as dry. 3.2. Limited slope evaluation

3. Numerical model The unstructured finite volume scheme used in this study is an incremental modification of the previously published model (Song et al., 2011a). The present model applies the MUSCL-Hancock scheme to achieve second-order accuracy in space and time. A HLLC-based integrated solver is proposed to evaluate the flow and mass transport fluxes simultaneously. The sloping bottom grids which the bed elevations are defined at nodes are used to represent the topography with second-order accuracy. The major subroutines of the present model are briefly described in the sections below. 3.1. Sloping bottom grids The unstructured triangular grids are used in this study for computational domain discretization. The flow variables (i.e. water depth, water surface elevation, and velocities) and the contaminating material concentration are specified at cell centroid, whereas the bed elevations are defined at nodes to represent terrain with second-order spatial accuracy. The bed slopes of the cell Ci are given by

∂bi (x, y) ∂x 1 = [(yi,2 − yi,3 ) bi,1 + (yi,3 − yi,1 ) bi,2 + (yi,1 − yi,2 ) bi,3 ] 2Ω i

αi =

∂bi (x, y) ∂y 1 = [(xi,3 − xi,2 ) bi,1 + (xi,1 − xi,3 ) bi,2 + (xi,2 − xi,1) bi,3 ] 2Ω i

(4)

(5)

where Ωi denotes the area of cell Ci ; x i, k , yi, k , and bi, k are x -, y -Cartesian coordinates and bed elevation at vertex vk (k ¼ 1, 2, 3) in cell Ci , supposing that the vertices v1−v2−v3 is in the anticlockwise direction. Provided that the water depth is given and the free surface elevation is uniform in a cell, then the water surface elevation of the cell Ci can be computed by VFRs (Begnudelli and Sanders, 2006; Song et al., 2011a)

where

∇pi = ( − μ2 /μ3 ,

− μ1/μ3)

(8)

where p is a generic dependent variable (i.e. η , u, v , c ); ∇p is the unlimited gradient; μ = (μ1, μ2 , μ3) = (Ρ1 − Ρ3) × (Ρ2 − Ρ3); Ρk = (x ik , yik , pik )T (k = 1, 2, 3); pi1, pi2, and pi3 are cell central values of three neighboring cells; (x i1, yi1), (x i2, yi2 ), and (x i3, yi3 ) are the Cartesian coordinates of centroids of three neighboring cells, which are numbered in the counter-clockwise direction. To achieve the TVD property, the unlimited slopes is limited by using the approach proposed by Barth and Jespersen (1989), and limits at vertices:

∇pi = ϕi ∇pi

(9)

ϕi = min (ϕk ) k = 1,2,3

⎧ p max − pi ⎪ min (1, i ) if pik − pi > 0 ⎪ pik − pi ⎪ ⎪ pimin − pi ϕk = ⎨ ) if pik − pi < 0 ⎪ min (1, pik − pi ⎪ ⎪ ⎪ 1 if pik − pi = 0 ⎩ pik = pi + ∇pi ⋅ri, k

(10)

where ϕi is the limit function; ∇pi is the limited gradient; pi is the

βi =

⎧ η 1 = b + 3 3h (b − b )(b − b ) if b < η 1 ≤ b i,1 i i,2 i,1 i,3 i,1 i,1 i,2 i ⎪ i ⎪ 1 2 2 2 ηi = ⎨ ηi = ( − γ1 + γ1 − 4γ2 ) if bi,2 < ηi ≤ bi,3 2 ⎪ ⎪ 3 if ηi 3 > bi,3 ⎩ ηi = hi + (bi,1 + bi,2 + bi,3 )/3

In this work the primitive variables η , u, v , and c are used for slope evaluation. The unlimited slope is computed by the differencing of three neighboring centroids data

variable value at the centroid of Ci; pik is the unlimited variable value at the vertex vk of Ci; pimax = max (pi , pi1, pi2 , pi3 ) and

pimin = min (pi , pi1, pi2 , pi3 ); ri,k is the position vector of vertex vk relative to the centroid of the Ci. 3.3. Predictor The primitive variables η , u, v , and c are advanced from t to t + Δt/2 in the predictor step. The updated solutions will be used by the reconstruction step. If the cell Ci is classified as dry, then the predictor results are given by

ηit + Δt /2 = ηit , uit + Δt /2 = uit ,

(6)

vit + Δt /2 = vit , cit + Δt /2 = cit

(11)



4

If the cell Ci is classified as wet, then the predictor results are given by

ηit + Δt /2 = ηit −

t Δt ¯ (h∂ xu + h∂¯ yv + u∂¯ xh + v∂¯ yh) 2 i

4m

(12)

1 + gn2h−4/3

⎡ ⎛ ⎞⎤ ⎢u − Δt ⎜u∂¯ xu + g ∂¯ xη + v∂¯ yu⎟ ⎥ ⎜ ⎟⎥ 2 ⎝ u2 + v2 Δt /2 ⎢⎣ ⎠⎦

0

t

y i

x ⎡ ⎛ ⎞⎤ ⎢v − Δt ⎜u∂¯ xv + g ∂¯ yη + v∂¯ yv⎟ ⎥ ⎜ ⎟⎥ 2 ⎝ u2 + v2 Δt /2 ⎢⎣ ⎠⎦

t Δt ¯ (u∂ xc + v∂¯ yc) 2 i

t

Fig. 1. Partial dam-break experiment: schematic view of the experimental set. i

(14)

(15)

where the superscript t + Δt/2 indicates the predictor results of half-time level solution; the superscript t indicates the base time level solution; ∂ x and ∂ y indicate the limited slopes in the x - and y -directions, respectively; Δt is the time step constrained by the Courant–Friedrichs–Lewy (CFL) condition for numerical stability. 3.4. Variables reconstruction In the reconstruction step, the primitive variables η , u, v , and c at the midpoints of edges are linearly extrapolated from the cellcentral values of the half-time level solution (Hu et al., 2000; Liang and Borthwick, 2009). The reconstructed variables will be used to compute numerical fluxes. The limited linear reconstruction is given by

piRec ,k

= pi + ∇pi ⋅mi, k

(16)

where piRec , k is the reconstructed variable (i.e. η , u, v , and c ) at the midpoints of edges; m is the position vector of the midpoints of edges to the centroid of the cell; ∇p is the limited gradient. The reconstructed depth is given by (Begnudelli and Sanders, 2007)

hRec

y (m) 1 .00 0 .40 1 .16 1 .00 1 .00

1m 1

cit + Δt /2 = cit −

Point x (m) -5A 0.18 C 0.48 4 1 .00 0 1 .00 8 A 1.72

C

(13)

vit + Δt /2 =

1+

8A 2m

1

gn2h−4/3

4

-5A

uit + Δt /2 =

⎧ 0 if η Rec ≤bmin ⎪ ⎪ (η Rec − bmin )2 ⎪ if bmin bmax ⎪ η Rec − ⎩ 2

invariants and boundary conditions (Song et al., 2011a). Assume that the water flow fluxes are denoted as F WF = (F 1, F 2, F 3), in which F 1 denotes the flux of water volume, F 2 and F 3 denote the momentum fluxes. Then the convection flux of mass transport is given by

⎧ 1 ⎪ F cL F adv = ⎨ ⎪ 1 ⎩ F cR

if

F1 ≥ 0

if

F1 < 0

(18)

where denotes the convection flux; cL and cR are the reconstructed values from the left-cell and the right-cell, respectively; F 1 > 0 means that water flows from the left-cell to the rightcell, and F 1 < 0 means that water flows from the right-cell to the left-cell. It should be noted that (18) is also used for boundary edges. If the boundary is inflow (i.e. F 1 < 0), then cR should be given by the time series of boundary condition. By using the divergence theorem, the diffusion flux of mass transport is given by

F adv

F dif = 0.5(hLRec + hRRec )(Dx ⋅∂ xc⋅n x + D y ⋅∂ yc⋅n y ) hLRec

(19) hRRec

where denotes the diffusion flux; and are the reconstructed water depth from the left-cell and the right-cell, respectively; n ¼(nx' ny)T is the unit outward normal vector; ∂ x and ∂ y are the unlimited slopes in the x - and y -directions, respectively.

F dif

3.6. Corrector In the corrector step, the solution is advanced from t to t + Δt . The conserved vector is first updated by using numerical fluxes and bed slope approximation, which is given by

(17) bmin

where the superscript Rec indicates the reconstructed value, and bmax are the minimum and maximum bed elevations of the two endpoints of edge, respectively. If hRec o 10 6 m, then the corresponding variables hRec , uRec , v Rec , and c Rec are all set as zero. 3.5. Fluxes computation Numerical fluxes of all edges are computed by the reconstructed variables at the midpoints of edges. If the edge is a boundary edge, the boundary conditions should also be used for fluxes computation. For inner edges, the HLLC approximate Riemann solver (Song et al., 2011a) is adopted to compute the water flow fluxes. Boundary edges are classified as different types, such as solid wall, free outflow, water level or discharge boundary. The fluxes of boundary edges could be computed by using Riemann

⎡ ˜ ti + Δt = Uti + Δt ⎢ − U Ωi ⎢⎣

⎤t + Δt /2

3

∑ Fi, k⋅ni, k li, k + Si,0 ⎥ ⎥⎦

k=1

(20)

˜ ˜ , hv ˜ ˜ , hc ˜ ˜] denotes the intermediate results; ˜ = [h˜ , hu where U the subscripts i and k indicate the index of cell and cell’s edge, respectively; l is the edge length; F⋅ n= (F 1, F 2, F 3, F adv + F dif ) is the outward normal flux vector for both water flows and mass transport; Si,0 = (0, Si,0x, Si,0y, 0)T is the bed slope estimation. The intermediate velocities are further updated by using a semi-implicit scheme for friction (Begnudelli and Sanders, 2006; Song et al., 2011a):

ut + Δt =

1 u˜ t + Δt 1 + Δtτ˜t + Δt

vt + Δt =

1 v˜t + Δt 1 + Δtτ˜t + Δt

(21)

t + Δt −4/3 where τ˜t + Δt = gn2 (u˜ t + Δt )2 + (v˜ t + Δt )2 (h˜ . So the conserved ) vector at the new time level t + Δt is given by t + Δt t + Δt t + Δt t + Δt t + Δt t + Δt t + Δt Ut + Δt = [h˜ , h˜ u , h˜ v , h˜ c˜ ].



0.8

0.8 Numerical Measured

0.4

Numerical Measured

0.6 h (m)

0.6 h (m)

5

0.4 0.2

0.2 0.0 0

2

4

6

8

0.0 0

10

2

4

t (s)

8

10

t (s)

0.8


0.4 0.2

Numerical Measured

0.6 h (m)

0.6 h (m)

6

0.4 0.2

0.0 0

2

4

6

8

10

0.0 0

2

4

t (s)

6

8

10

t (s)


h (m)

0.6 0.4 0.2 0.0 0

2

4

6

8

10

t (s) Fig. 2. Partial dam-break experiment: comparison of water depths at surveyed points.

If the cell Ci is classified as wet, then the bed slopes are estimated by

Si,0y = − =

Si,0x = −

∫C

i

g (h + b) αdΩ = − gηi αi Ωi

(22)

∫C

∫C

i

∑ k=1

Si,0y = −

∫C

i

g (h + b) βdΩ = − gηi βi Ωi

(23)

If the cell Ci is classified as dry, then the bed slopes are estimated by using the DFB technique (Valiani and Begnudelli, 2006)

Si,0x = −

=

∫C

∫C

i

g (h + b) αdΩ

gh

3

=

i

∑ k=1

∂h dΩ − ∂x

∫C

i

gb

∂b dΩ ∂x

1 x 2 g (hiRec , k ) l i , k ni , k − 2

3

∑ k=1

1 g (bi, k )2li, k nix, k 2

(24)

g (h + b) βdΩ

gh

3

=

i

∂h dΩ − ∂y

∫C

i

gb

∂b dΩ ∂y

1 y 2 g (hiRec , k ) l i , k ni , k − 2

3

∑ k=1

1 g (bi, k )2li, k niy, k 2

(25)

where bi, k is the bed elevation at the middle point of the kth edge of the ith cell; nx and n y are the components of the unit outward normal vector.

3.7. The well-balanced property The well-balanced property means that a numerical scheme could maintain the stationary state in the presence of a stationary flow field. The well-balanced property is approved for the proposed model as below. Considering the cell Ci which could be either wet or dry, the following relations stand in the stationary flow field situation for the kth edge (k ¼1, 2, 3):



6

t = 10 s

t =0s 30

c (mg/L) 1.0001 1.0001

20

1.0000 1.0000

10

0.9999 0.9999

b (m)

10 m/s

30

3.0 2.0 1.5

10

y (m)

y (m)

2.5

20

1.0 0.5

0

0

10

20

30

40

50

60

70

0

0.0

0

10

20

30

t = 20 s y (m)

2.5

20

2.0 1.5

10

y (m)

3.0

1.0 0.5

10

20

30

40

60

70

t = 10 s

b (m)

10 m/s

30

0

50

50

60

70

30

c (mg/L) 1.0001 1.0001

20

1.0000 1.0000

10

0.9999 0.9999

0

0.0

0

10

20

30

40

50

60

70

t = 30 s t = 30 s

b (m)

10 m/s

30

30

c (mg/L) 1.0001 1.0001

20

1.0000 1.0000

10

0.9999 0.9999

2.5 2.0 1.5

10

y (m)

3.0

20

y (m)

0.0000 0.0000

x (y)

x (m)

1.0 0.5

0

0.0000 0.0000

x (y)

x (m)

0

40

0

10

20

30

40

50

60

70

0.0

x (m) Fig. 3. Flow-mass transport with a uniform concentration: the simulated velocity fields.

0

0

10

20

30

40

50

60

70

0.0000 0.0000

x (y) Fig. 4. Flow-mass transport with a uniform concentration: the simulated concentration.

4. Numerical results and discussions

ηi Rec , k = η = h + b = Constant (if h > 0)

Rec uiRec , k = vi , k = 0

(26) 4.1. Test 1: Partial dam-break experiment with a symmetrical breach

Then we have

3

∑ Fi, k⋅ni, k li, k k=1

⎡ ⎢ Rec 3 ⎢ [(hi, k )2 − 1 = g∑⎢ Rec 2 k = 1 ⎢ [(hi, k )2 − ⎢⎣

⎤ ⎥ (bi, k )2] li, k nix, k ⎥ ⎥ (bi, k )2] li, k niy, k ⎥ ⎥⎦ 0 0

(27)

If the cell Ci is wet, then it is easy to achieve that 3

1 x 2 2 − g ∑ [(hiRec , k ) − (bi, k ) ] li, k ni , k = gηi αi Ωi 2 k=1

(28)

3

1 y 2 2 − g ∑ [(hiRec , k ) − (bi, k ) ] li, k ni , k = gηi βi Ωi 2 k=1

(29)

So we have 3

− ∑ Fi, k ⋅ni, k li, k + Si,0 = 0 k=1

(30)

Besides, from (24), (25), and (27) it is easy to ascertain that the equality of (30) also holds when the cell Ci is dry. So (30) approves the well-balanced property for both wet and dry cells, validating that the proposed model can preserve the well-balanced property for simulations involving flooding and recession.

Fraccarollo and Toro (1995) performed a laboratory experiment case of dam-break flood wave propagation on a horizontal bottom, which has been employed for the 2D validation of hydrodynamic models (Gottardi and Venutelli, 2004; Singh et al., 2011). Fig. 1 shows the schematic view of the experimental set and the location of five gauge points. The reservoir is 2 m long and 1 m wide with a downstream floodplain of 3 m length. The dam is 2 m wide with a 0.4 m breach symmetrically located in the middle. The particular case that the initially stationary water depth in the reservoir is 0.6 m and the downstream floodplain is initially dry is selected in this paper. The dam breaches instantaneously at t¼0. The reservoir boundaries are solid wall, whereas all sides of the downstream floodplain are free outflow. Some nodes are specified at the dam location before triangulating, and then the dam is defined by raising the bed elevation of these nodes. Fig. 2(a)–(e) compare the observed water depth at five gauge points with the numerical predictions. The simulated depth hydrographs at the upstream gauge points “ 5 A” and “C” are in good agreement with the experimental data, see Fig. 2(a) and (b). The model slightly underestimated the water depth at the dam-site gauge points “4” and “0” in the initial phase of dam-break flood propagation, but subsequently reproduced the water depth reasonably well, see Fig. 2 (c) and (d). The depth hydrograph at the downstream gauge point “8A” is predicted up to a reasonable accuracy; see Fig. 2(e). Overall, the numerical predictions and the measured data agree well and it



t =6s

30

c (mg/L) 1.0001

30

c (mg/L) 1.0001

20

1.0000

20

1.0000

10

0.9999

10

0.9999

0

0

10

20

30

40

50

60

70

y (m)

y (m)

t =0s

7

0

0.0000

0

10

20

30

x (m)

60

70

0.0000

t = 18 s

30

c (mg/L) 1.0001

30

c (mg/L) 1.0001

20

1.0000

20

1.0000

10

0.9999

10

0.9999

0

0

10

20

30

40

50

60

70

y (m)

y (m)

50

x (m)

t = 12 s

0.0000

0

0

10

20

30

x (m)

40

50

60

70

0.0000

x (m)

t = 24 s

t = 30 s

30

c (mg/L) 1.0001

30

c (mg/L) 1.0001

20

1.0000

20

1.0000

10

0.9999

10

0.9999

0

0.0000 0

10

20

30

40

50

60

70

y (m)

y (m)

40

0

0.0000 0

10

x (m)

20

30

40

50

60

70

x (m)

Fig. 5. Flow-mass transport with a concentration step: the simulated concentration.

c (mg/L)

t =0s

30

y (m)

1.00 20

0.75 0.50

10

0.25 0

0

10

20

30

40

50

60

70

0.00

x (m) c (mg/L)

t =8s

30

y (m)

1.00 20

0.75 0.50

10

0.25 0

0

10

20

30

40

50

60

70

0.00

x (m) c (mg/L)

t = 20 s

30

y (m)

1.00 20

0.75 0.50

10

0.25 0

can be concluded that the performance of the proposed model in simulation the two-dimensional, partial dam-break experimental case is good.

0

10

20

30

40

50

60

70

0.00

x (m) Fig. 6. Mass diffusion in a stationary flow: the simulated concentration.

4.2. Test 2: Flow-mass transport over three Humps in a closed channel To evaluate the performance of the model in problems involving wetting–drying process, in this test we consider the case of flow-mass transport over three humps in a closed channel. The channel is 75 m long and 30 m wide, with the bed topography defined by

⎡ 1 b (x, y) = max ⎢0, 1 − (x − 30)2 + (y − 6)2 , ⎣ 8 1 1− (x − 30)2 + (y − 24)2 , 8

3−

⎤ 3 (x − 47.5)2 + (y − 15)2 ⎥ ⎦ 10

(31)

Three different cases are considered in this test as follows. 4.2.1. Flow-mass transport with a uniform concentration A dam is located at x¼ 16 m that initially retains still water with constant surface elevation 1.875 m, while the rest of the channel is dry. The Manning coefficient n is set to 0.018. Solid slip conditions are applied at the lateral walls of the channel. The contaminating material concentration is uniformly 1 mg/L. The degradation coefficient k is set to 0. The mass diffusion coefficients Dx and Dy are set to 0.5 m2/s. At t¼0, the dam collapses instantaneously. Since the initial concentration is uniform, the computed concentration should be 1 mg/L for any time.



8

1.2

1 Numerical Analytical

1

c (mg/L)

c (mg/L)

0.8 0.6 0.4

Numerical Analytical

0.8 0.6 0.4

0.2 0

0.2

0

5000 x (m)

0

10000


10000

c (mg/L)

0.6


1

0.8 c (mg/L)

5000 x (m)

1.2

1

0.4

0.8 0.6 0.4

0.2 0

0

0.2 0

0

5000 x (m)

10000

0.3

0

5000 x (m)

10000

1.2 Numerical Analytical

1

c (mg/L)

c (mg/L)

0.2

0.1


0.8 0.6 0.4 0.2

0

0

5000 x (m)

10000

Fig. 7. Mass transport of Gaussian distribution in uniform flow: results comparison. (a) D ¼ 0 m2 /s, (b) D ¼ 2 m2 /s and (c) D ¼ 50 m2 /s.

0

0

5000 x (m)

10000

Fig. 8. Mass transport of step distribution in uniform flow: results comparison. (a) D ¼ 0 m2 /s, (b) D ¼ 2 m2 /s and (c) D ¼ 50 m2 /s.



The hydrodynamic simulation of this case of dam-break floods was reported by Song et al., (2011b), and in this paper the details of flooding process are not discussed for conciseness. Fig. 3 presents the simulated velocity fields at t¼10 s, 20 s, and 30 s, from which it can be seen to have the expected physical behavior of symmetrical flows. Fig. 4 presents the initial and simulated concentration at t ¼0 s, 10 s, and 30 s. Results show that the computed concentration is very close to 1 mg/L, validating the accuracy of the proposed model. 4.2.2. Flow-mass transport with a concentration step All computational conditions are the same as that of the former case, except that the initial concentration is 1 mg/L for xo 8 m, otherwise the concentration is 0. Since there is a concentration step at x¼ 8 m, the convection and diffusion process should be observed in this case. Fig. 5 presents the initial and simulated concentration at t¼ 0 s, 6 s, 12 s, 18 s, 24 s, and 30 s, from which it can be seen to have the expected flow-mass transport behavior. The water body with high concentration is diluted by the clean water. Besides, the simulated results are highly symmetrical along the central line of y ¼15 m, due to the symmetry of bed and initial conditions. 4.2.3. Mass diffusion in a stationary flow This case is used for testing the well-balanced property of the proposed model. All computational conditions are the same as that of the former case, except that the initial condition is given by

η = 1.875 m u=v=0 ⎧1 mg/L c=⎨ ⎩0

for 10m < x< 16m,

12m < y< 18m

9

c (x, y , t) =

⎛ x − x 0 − ut ⎞ ⎛ u (x − x 0 ) ⎞ ⎛ x − x 0 + ut ⎞ ⎞ 1⎛ ⎜ erfc ⎜ ⎟⎟ ⎟ + exp ⎜ ⎟ erfc ⎜ 2⎝ ⎝ 2 Dt ⎠ ⎝ D ⎠ ⎝ 2 Dt ⎠⎠

where erfc is the complementary error function. We run again the simulations until t¼9600 s for Dx = Dy ¼0 m2/s, 2 m2/s, and 50 m2/s, respectively. Fig. 8 shows the computed concentration distributions at t¼9600 s, compared with the analytical solutions: how it can be seen, a very good agreement is observed. The results validate the accuracy of the proposed coupled model for simulations with discontinuousness or large gradient in concentration distribution.

5. Conclusions A two-dimensional coupled flow-mass transport model, which is based on the modified shallow water equations and the convection–diffusion equation, has been developed and tested to simulate water quality in shallow water. Advantages of the present model include the computational effectiveness by using the MUSCL-Hancock’s predictor–corrector scheme, the numerical robustness in wet–dry treatment by using the new hybrid method for bed slope approximation, and the capacity of mass transport simulation with discontinuity by using an integrated HLLC-type flux solver. The well-balanced property for both wet and dry cells is theoretically approved for the proposed scheme. Numerical results show that the proposed model can simulate complex flows and mass transport with low numerical damping and oscillations.

Competing interests

otherwise

Fig. 6 presents the initial and simulated concentration at t¼ 0 s, 8 s, and 20 s, from which it can be seen to have the expected mass diffusion behavior in a stationary flow. Besides, the magnitude of resultant velocity was computed to be 0 during the simulation, validating the well-balanced property of the proposed model for cases involving flooding and recession.

The authors have no conflict of interest to disclose.

Acknowledgements

This test case is used to show the calculation accuracy of the proposed model. The domain is [0, 12800 m] [0, 1000 m]. The constant velocity of the uniform flow is u0 = 0.5 m /s . The initial

This work was supported by a grant from the State Key Program of National Natural Science Foundation of China (Project no. 51239004), a grant from the Open Research Foundation of PRHRI (Project no. 2013KJ01), a grant from the National Basic Research Program of China (Project no. 2007CB714107), a grant from the New Star of Pearl River on Science and Technology of Guangzhou (Project no. 2012J2200071), and a grant from the CRSRI Open Research Program (Program SN:CKWV2014220/KY).

concentration is c(x,y,0)¼exp[-(x 2000)2/(2σ02)], where σ0 = 264 m . The analytical solution of the concentration distribution is given by

References

4.3. Test 3: mass transport with uniform flow

c (x, y , t) = σ0/σ⋅exp ⎡⎣− (x − x´)2 /(2σ 2) ⎤⎦ 2

where s

¼ σ02 þ2Dt,

(32)

t

x´¼ 2000 þ ∫ u (τ)dτ . 0 We run the simulations until t¼9600 s for Dx = Dy ¼ 0 m2/s, 2 m2/s, and 50 m2/s, respectively. Fig. 7 presents the comparisons between the computed and analytical solutions at t¼9600 s. The results show that the computed concentration distributions match the exact solutions well. Furthermore, the initial concentration distribution is changed to be a step, i.e.

⎧1, c (x, y , 0) = ⎨ ⎩ 0, ⎪

x ≤ x0 x > x0

and the boundary condition is c (0, y, t) = 1, exact solution is given by,

t ≥ 0, then the

Audusse, E., Bouchut, F., Bristeau, M.O., Klein, R., Perthame, B., 2004. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput. 25 (6), 2050–2065. Barth, T.J., Jespersen, D.C., 1989. “The design and application of upwind schemes on unstructured meshes.”. AIAA, 890366. Begnudelli, L., Sanders, B.F., 2006. “Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying.”. J. Hydraul. Eng. 132 (4), 371–384. Begnudelli, L., Sanders, B.F., 2007. “Conservative wetting and drying methodology for quadrilateral grid finite-volume models.”. J. Hydraul. Eng. 133 (3), 312–322. Benkhaldoun, F., Elmahi, I., 2007. “Well-balanced finite volume schemes for pollutant transport by shallow water equations on unstructured meshes.”. J. Comput. Phys. 226 (1), 180–203. Bermudez, A., Vazquez, M.E., 1994. “Upwind methods for hyperbolic conservation laws with source terms.”. Comput. Fluids 23 (8), 1049–1071. Cai, L., Xie, W.X., Feng, J.H., Zhou, J., 2007. “Computations of transport of pollutant in shallow water.”. Appl. Math. Model. 31 (3), 490–498. Canestrelli, A., Dumbser, M., Siviglia, A., Toro, E.F., 2010. “Well-balanced high-order centered schemes on unstructured meshes for shallow water equations with fixed and mobile bed.”. Adv. Water Resour. 33 (3), 291–303.


10


Cea, L., Vázquez-Cendón, M.E., 2012. “Unstructured finite volume discretisation of bed friction and convective flux in solute transport models linked to the shallow water equations.”. J. Comput. Phys. 231 (8), 3317–3339. Fraccarollo, L., Toro, E.F., 1995. “Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems.”. J. Hydraul. Res. 33 (6), 843–864. Gottardi, G., Venutelli, M., 2004. “Central scheme for two-dimensional dam-break flow simulation.”. Adv. Water Resour. 27 (3), 259–268. Hou, J., Simons, F., Mahgoub, M., Hinkelmann, R., 2013. “A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography”. Comput. Methods Appl. Mech. Eng. 257, 126–149. Hu, K., Mingham, C.G., Causon, D.M., 2000. “Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations.”. Coast. Eng. 41 (4), 433–465. Issa, R., Rouge, D., Benoit, M., Violeau, D., Joly, A., 2010. “Modelling algae transport in coastal areas with a shallow water equation model including wave effects.”. J. Hydro-environ. Res. 3 (4), 215–223. Lai, J.S., Guo, W.D., Lin, G.F., Tan, Y.C., 2010. “A well-balanced upstream flux-splitting finite-volume scheme for shallow-water flow simulations with irregular bed topography.”. Int. J. Numer. Methods Fluids 62 (8), 927–944. Lai, W., Khan, A.A., 2012. “A discontinuous Galerkin method for two-dimensional shallow water flows.”. Int. J. Numer. Methods Fluids 70 (8), 939–960. Li, G., Gao, J., Liang, Q., 2012. “A well-balanced weighted essentially non-oscillatory scheme for pollutant transport in shallow water.”. Int. J. Numer. Methods Fluids 71, 1566–1587. Li, S., Duffy, C.J., 2012. “Fully-coupled modeling of shallow water flow and pollutant transport on unstructured grids.”. Proc. Environ. Sci. 13, 2098–2121. Liang, D., Wang, X., Falconer, R.A., Bockelmann, B.N., 2010. “Solving the depth-integrated solute transport equation with a TVD-MacCormack scheme.”. Environ.

Model. Softw. 25 (12), 1619–1629. Liang, Q., 2010. “Flood simulation using a well-balanced shallow flow model.”. J. Hydraul. Eng. 136 (9), 669–675. Liang, Q., Borthwick, A.G.L., 2009. “Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography.”. Comput. Fluids 38 (2), 221–234. Liang, Q., Marche, F., 2009. “Numerical resolution of well-balanced shallow water equations with complex source terms.”. Adv. Water Resour. 32 (6), 873–884. Nguyen, D.K., Shi, Y.E., Wang, S.S.Y., Nguyen, T.H., 2006. “2D shallow-water model using unstructured finite-volumes methods.”. J. Hydraul. Eng. 132 (3), 258–269. Singh, J., Altinakar, M.S., Ding, Y., 2011. “Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme.”. Adv. Water Resour. 34 (10), 1366–1375. Song, L., Zhou, J., Guo, J., Zou, Q., Liu, Y., 2011a. “A robust well-balanced finite volume model for shallow water flows with wetting and drying over irregular terrain.”. Adv. Water Resour. 34 (7), 915–932. Song, L., Zhou, J., Li, Q., Yang, X., Zhang, Y., 2011b. “An unstructured finite volume model for dam‐break floods with wet/dry fronts over complex topography.”. Int. J. Numer. Methods Fluids 67 (8), 960–980. Valiani, A., Begnudelli, L., 2006. “Divergence form for bed slope source term in shallow water equations.”. J. Hydraul. Eng. 132 (7), 652–665. Xing, Y., Shu, C.W., Noelle, S., 2011. “On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations.”. J. Sci. Comput. 48 (13), 339–349. Zia, A., Banihashemi, M.A., 2008. “Simple efficient algorithm (SEA) for shallow flows with shock wave on dry and irregular beds.”. Int. J. Numer. Methods Fluids 56 (11), 2021–2043.


Research on an Improved Medical Image Enhancement Algorithm Based on P-M Model.

Microcellular propagation prediction model based on an improved ray tracing algorithm.

An improved FCM medical image segmentation algorithm based on MMTD.

Use of an Improved Matching Algorithm to Select Scaffolds for Enzyme Design Based on a Complex Active Site Model.

An improved piecewise linear chaotic map based image encryption algorithm.

CT based on an improved watershed algorithm.

Development of a hip joint model for finite volume simulations.

An Improved Interacting Multiple Model Filtering Algorithm Based on the Cubature Kalman Filter for Maneuvering Target Tracking.

An Improved Multiobjective Optimization Evolutionary Algorithm Based on Decomposition for Complex Pareto Fronts.

An Improved Inertial Frame Alignment Algorithm Based on Horizontal Alignment Information for Marine SINS.

An improved fuzzy c-means clustering algorithm based on shadowed sets and PSO.

Parameters identification for photovoltaic module based on an improved artificial fish swarm algorithm.

An efficient technique for nuclei segmentation based on ellipse descriptor analysis and improved seed detection algorithm.

Visual Tracking Based on an Improved Online Multiple Instance Learning Algorithm.

Automatic ICD-10 coding algorithm using an improved longest common subsequence based on semantic similarity.

Construction of monitoring model and algorithm design on passenger security during shipping based on improved Bayesian network.

Improved algorithm for gradient vector flow based active contour model using global and local information.

Extending the MODPATH Algorithm to Rectangular Unstructured Grids.

Network intrusion detection based on a general regression neural network optimized by an improved artificial immune algorithm.

Improved random-starting method for the EM algorithm for finite mixtures of regressions.

A combined approach to cartographic displacement for buildings based on skeleton and improved elastic beam algorithm.

A coupled model of fast axonal transport of organelles and slow axonal transport of tau protein.

Registration of the Cone Beam CT and Blue-Ray Scanned Dental Model Based on the Improved ICP Algorithm.

An Improved PID Algorithm Based on Insulin-on-Board Estimate for Blood Glucose Control with Type 1 Diabetes.