- Email: [email protected]
A coupled morphodynamic model for applications involving wetting and drying
273
2011,23(3):273-281 DOI: 10.1016/S1001-6058(10)60113-8
A COUPLED MORPHODYNAMIC MODEL FOR APPLICATIONS INVOLVING WETTING AND DRYING* LIANG Qiuhua School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne NE1 7RU, England, UK, E-mail: [email protected] (Received April 22, 2011, Revised May 23, 2011) Abstract: This work presents a new finite volume Godunov-type model for predicting morphological changes under the rapidly varying flood conditions with wetting and drying. The model solves the coupled shallow water and Exner equations, with the interface fluxes evaluated by an Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver. Well-balanced solution is achieved using the surface gradient method and wetting and drying are handled by a non-negative reconstruction approach. The new model is validated against several theoretical benchmark tests and promising results are obtained. Key words: morphodynamic model, Exner equation, shallow water equations, Godunov-type method, wetting and drying, wellbalanced scheme
Introduction Morphodynamics concerns the evolution of the channel bathymetry in response to sediment activities driven by flows. It may have great impacts on local habitats, affect river conveyance and hence have implications for flooding during storm events. An understanding of morphodynamics and the relevant processes is therefore vital for river management systems. In the past few decades, with the rapid development of computer technology, computer modelling has become a common practice in assessing river morphodynamics although it may not entirely replace the laboratory studies for understanding the relevant processes[1]. A morphodynamic modelling system is typically consisted of a hydrodynamic component that describes the flow hydrodynamics and a sediment transport/morphological component that accounts for the bed evolution[2]. Conventionally, morphodynamic modelling is carried out in a decoupled way where the flow hydrodynamics is first obtained and then used to drive sediment transport and update the channel bed * Project supported by the UK Engineering and Physical Sciences Research Council (Grant No. EP/F030177/1). Biography: LIANG Qiu-hua (1974-), Male, Ph. D., Lecturer
change[2-4]. Ignoring the unsteady hydrodynamic effects, this approach has been recognized to be more appropriate for slow-varying flows or quasi-steady processes[5]. For a rapidly-varying flow (e.g., dam breaks, floods, etc.) where the flow unsteadiness can strongly impact the sediment transport, a coupled model, in which the flow hydrodynamics and morphodynamics are solved instantaneously, is desirable and may lead to more reliable predictions[5]. A number of coupled models that solve the integrated shallow water and Exner equations have been reported in recent years for simulating morphodynamics. For example, Liu et al.[6] presented a Roe’s approximate Riemann solver based finite volume Godunov-type morphodynamic model on unstructured grids with applications to coastal waters. They also compared the results with those produced by the quasi-steady approach and confirmed the necessity of using a coupled approach for a problem involving strong fluid-sediment interactions. Benkhaldoun et al.[7] reported an alternative finite volume Godunov-type model combined with the Roe’s approximate Riemann solver on adaptive triangular grids. The authors also compared the performance of the coupled and decoupled methods and indicated that the coupled approach out-performed the quasi-steady approach when the bed-load rapidly interacts with the hydrodynamics.
274
Murillo and García-Navarro[8] formulated the sediment fluxes in the Exner equation in a universal form for a number of different empirical formulae and solved the coupled system using a first-order finite volume Roe-type scheme. Soares-Frazão and Zech[9] presented a first-order finite volume Godunov-type morphodynamic model using an Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver, where new wave-speed estimates were derived based on the novel approximate eigenvalues to the coupled governing equations. Attempts have also been made to solve the coupled system using Essentially NonOscillatory (ENO) and Weighted Essentially NonOscillatory (WENO) schemes[10], higher-order Central Weighted Essentially Non-Oscillatory (CWENO) scheme[11] and finite element discontinuous Galerkin method[12,13]. However, most of these existing models concerns only wet-bed applications, i.e., the problem domains are covered entirely by water. A morphodynamic model that can handle applications involving wetting and drying may lead to wider applications. Therefore, this work aims to develop a new 2-D shallow flow model to include morphological change under the rapidly varying flood conditions that involves wetting and drying. 1. Governing equations The coupled governing equations, consisting of the 2-D SWEs and the Exner equation, describing shallow flows over erodible bed may be written in a matrix form as wq wf wg + + =s (1) wt wx wy where t , x and y are the time and the two Cartesian coordinates, q denotes the vector of flow variables, f and g are the flux vectors in the x and y -direction, and s is the vector containing the source terms. The vectors are given by ª qy º ª qx º « qq » « » 2 2 ªhº x y « » « qx + gh » «q » « » « » h h 2 x» « , g=« 2 q= , f =« », » 2 «qy » qx q y » « q y + gh » « « » «h « » h 2 » ¬z¼ « » « q » «¬ qsy »¼ sx ¬ ¼
ª « 2 « C f qx (qx « h3 s=« 2 « C f q y (qx « h3 « «¬
º » 2 + q y ) ghwz » wx » » 2 + q y ) ghwz » wy »» »¼ 0
where h is the water depth, qx = uh and qv = vh are the unit-width discharge in the x and y -direction with u and v being the corresponding depth-averaged velocities, z represents the bed elevation above an arbitrary datum, g = 9.81 m / s 2
denotes the acceleration due to gravity, qs x = ] Ag qx <
q
2 x
+ q 2y / h3 and qs y = ] Ag q y qx2 + q 2y / h3 are the
sediment fluxes, ] = 1 / (1 p ) is a fixed parameter with p being the bed porosity, Ag (0 Ag d 1) is an empirical constant that indicates a stronger flowsediment interaction as it reaches 1, C f = gn 2 / h1/3 is the bed roughness coefficient with n being the Manning coefficient, wz / wx and wz / wy are the two bed slopes. In order to construct a proper eigenstructure, the coupled System (1) and (2) are rewritten as follows by moving the bed slope source terms to the left hand side[9] wq wf (q ) wq wg (q ) wq + + H (q ) + + K (q ) = s f wt wx wx wy wy
(3)
where
ª0 «0 H (q ) = « «0 « ¬0 ª « « « C f qx sf = « « « C f q y « «¬
0 0º 0 c 2 »» , K (q ) 0 0» » 0 0¼
0 0 0 0
ª0 «0 « «0 « ¬0
0 0 0 0
0 0º 0 0 »» , 0 c2 » » 0 0¼
º » q + q » » h3 » qx2 + q y2 »» h3 » »¼ 0 0
2 x
2 y
(4)
where c = gh . Then the coupled system may be written in a quasi-linear form as wq wq wq + A(q ) + B (q ) = sf wt wx wy
(5)
where A (q ) and B (q ) are the pseudo-Jacobian matrices given by[9]
0
(2)
§ 0 ¨ 2 c u2 wf (q ) + H (q ) = ¨ A (q) = ¨ uv wq ¨ © M
1 2u v N
0 0· ¸ 0 c2 ¸ u 0 ¸, ¸ R 0¹
275
§ 0 ¨ uv wg (q ) + K (q ) = ¨ 2 B (q ) = ¨ c v2 wq ¨ © D
0 1 0· ¸ v u 0¸ 0 2v c 2 ¸ ¸ R E 0¹
(6)
ximate Riemann solver is chosen in this work to solve the Riemann problems and it is well-documented in literature[14] for calculating the hydrodynamic fluxes in Eq.(1). Herein, the HLLC approximate Riemann solver is briefly introduced herein for completion and modified for the current application. For the first three components of f E (i.e., the hydrodynamic fluxes),
where M=
R=
E=
3] Ag u (u 2 + v 2 ) h
2] Ag uv h
, D=
, N=
] Ag (3u 2 + v 2 ) h
3] Ag v(u 2 + v 2 ) h
f E = f L if ,
f E = f*L if
,
] Ag (u 2 + 3v 2 ) h
Considering A (q ) , a set of approximate eigenvalues may be derived as[9]
O1 = 0.5 ªu c (u c)2 + 4c 2 N º , ¬ ¼ O2 = 0.5 ªu c (u c)2 + 4c 2 N º , ¬
¼
O3 = u , O4 = u + c
(7)
With these eigenvalues, we may derive the wave speeds based on which an HLLC approximate Riemann solver can be implemented. 2. Numerical model In this work, the coupled matrix System (1) is solved by a 1st order finite volume Godunov-type scheme, which updates the flow variables to a new time step using the following formula k +1 i, j
q
§ f fW g N gS · = q 't ¨ E + si , j ¸ 'y © 'x ¹ k i, j
(8)
where the superscript k represents time level, 't is the time interval, 'x and 'y are the cell dimensions, f E , fW , g N and gS are the fluxes across the four cell interfaces. The interface fluxes are evaluated by solving local Riemann problems, e.g., f E = F (qEL , qER ) in which qEL and qER are the Riemann states on either side of the cell interface under consideration and they are reconstructed from the cell-centred flow variables. The HLLC appro-
0 d SL SL d 0 d SM
(9a) (9b)
f E = f*R if S M d 0 d S R
(9c)
f E = f R if 0 t S R
(9d)
where f L = f (qEL ) , f R = f (qER ) and f*L and f*R are the HLLC fluxes in the star region of the solution structure[14]. The wave speeds are estimated by S L = u ER 2 ghER if hEL = 0
(10a)
S L = min O1L , O1* , O2 L , O2* if hEL ! 0
(10b)
S R = u EL + 2 ghEL if hER = 0
(10c)
S R = max O4 L , O4* if hER ! 0
(10d)
SM =
S L hER (u ER S R ) S R hEL (u EL S L ) hER (u ER S R ) hEL (uEL S L )
(10e)
with h* =
v* =
hEL uEL + hER u ER hEL + hER , , u* = 2 hEL + hER hEL vEL + hER vER hEL + hER
(11)
The sediment component of the fluxes may be estimated as[9] f E (4) = qsxL if 0 d S sL
(12a)
f E (4) = q*sx if S sL d 0 d S sR
(12b)
f E (4) = qsxR if 0 t S sR
(12c)
276
where the sediment flux in the star region is given by[9] q*s x =
S sR qsLx S sL qsRx + S sL S sR hEL ( z ER z EL ) S sR S sL
(13)
with the wave speeds approximated by S sL = min O1L , O1 , S sR = max O2 R , O2
(14)
In order to calculate the interface fluxes using the above HLLC approximate Riemann solver, the Riemann states must be properly reconstructed. For a first-order accurate scheme, the face values of flow variables are assumed to be the same as those at the cell centres, i.e.,
q EL = qi , j , q ER = qi +1, j , KEL = Ki , j , KER = Ki +1, j
(15)
where h = (hEL + hER ) / 2 . For the friction source terms, a point-implicit scheme proposed by Liang[16] is used to ensure non-negative water depth for simulations involving wetting and drying over complex domain topographies. The aforementioned finite volume Godunov-type morphodynamic model is overall explicit and its numerical stability is controlled by the CourantFriedrichs-Lewy (CFL) criterion. Either transmissive or reflective (slip) boundary conditions are used in the test cases considered in this work. 3. Results In this section, the aforementioned coupled morphological change model is validated against idealized tests. In all of the following test cases, the constants p = 0.4 and Ag = 0.01 are used. The
CFL number is set to 0.75.
where K = h + z is the water level, which is defined to assist the implementation of the surface gradient method[15]. In order to define the effective water depth at the cell interface and hence calculate the hydrodynamic fluxes, a single value of bed elevation is first defined z E = max( zEL , zER )
(16)
The Riemann states are then reconstructed as Fig.1 Exner’s dune: initial bed elevation L L hEL = max(0,KEL z E ) , qxE = u EL hEL , q yE = vEL hEL ,
zEL = zEL
(17)
The Riemann states on the right hand side of the eastern interface are defined accordingly. For applications involving wetting and drying over irregular topography, the following step is essential 'z = max 0, z E KEL
(18)
§ 2Sx · z0 = z ( x, 0) = A0 + A1 cos ¨ ¸ © O ¹
(21)
where A0 = A1 = 1 , O = 20 are constants. Denoting c( z ) { wqsx / wz
3] Ag qx3 / (K z )4 and assuming a
constant water level and discharge, the analytical solution is given by
and z E m z E 'z , zEL m z EL 'z , zER m z ER 'z
(19)
In the context of the above surface gradient method, the slope source terms are directly discretized by a central-differencing scheme, e.g., in the x -direction gh
3.1 Analytical test of Exner’s dune We first consider the analytical test of Exner’s dune. As shown in Fig.1, the initial bed profile is described by[17]
wz § z zW · = gh ¨ E ¸ wx © 'x ¹
(20)
ª 2S º z ( x, t ) = A0 + A1 cos « ( x c( z )t ) » ¬O ¼
(22)
which can be solved using the Newton-Raphson method. Simulations occur in a 40 m long domain for 50 s. The water level and discharge are respectively K { 3 m and qx { 1m 2 / s . Figure 2 shows the
277
numerical predictions obtained on grids of different resolution, which agree closely with the analytical solution and appear to be convergent.
analytical solution in the original paper in Ref.[17].
Fig.3 Analytical 1-D test with shock: numerical predictions at t = 3.5 s Fig.2 Exner’s dune: numerical predictions at t = 50 s where (b) is the zoom-in view
3.2 Analytical test with shock The second analytical case involves a so-called Riemann problem for the Exner equation[17]. In order to simplify the procedure of finding the analytical solution, the 1-D Exner equation takes the form of ( z )t + (D u ) x = 0 with D = ] Ag and the derivative terms M , N , R , D and E in Eqs.(6) and (7) are changed accordingly. For the initial bed profile z0 = z L , x 0
(23a)
z0 = z R , x ! 0
(23b)
the analytical solution for the case of z L z R may be derived to be[17]
z ( x, t ) = z L , x c( z L )t z ( x, t ) = K D q x
t , c( z L )t d x d c( z R )t x
z ( x, t ) = z R , x ! c( z L )t
(24a) (24b) (24c)
where c( z ) = D qx / (K z ) 2 for the simplified Exner equation. Please note the typing mistake of the
Fig.4 Dam break over three humps: initial conditions
Simulation occurs at a [– 1 m, 5 m] domain with z L = 1 m , zR = 2 m , K { 3 m and qx { 1 m 2 / s .
278
Fig.5 Dam break over three humps: results at t = 6 s
Fig.6 Dam break over three humps: results at t = 12 s
279
Fig.7 Dam break over three humps: results at t = 30 s
The numerical results obtained after t = 3.5 s on uniform meshes of different resolution are plotted in Fig.3, comparing with the analytical solution. The original discontinuous bed profile is found to be washed out by the current, where the numerical predictions again agree closely with the analytical solution and seem to be convergent with increasing grid refinement. Therefore, the capability of the current coupled model in modelling discontinuous sediment flux is confirmed. 3.3 Dam break over three humps This test case is about a dam-break flow running over three humps in the originally dry floodplain. The initial flow conditions are indicated in Fig.4 and the detailed settings can be found in Ref.[14]. Figures 5-8 present the numerical predi- ctions at different output times in terms of 3-D surface view and depth contours for both fixed and movable bed tests. The contour lines of eroded bed topography for the movable bed test are also illustrated. It is evident that the dam-break flow induces substantial change to the floodplain topography meanwhile the morphological changes also feedback and modify the flow hydrodynamics. As shown for t = 6 s , the dam- break wave rushes downstream after the dam has been removed, carrying great momentum to reshape the landform. After t = 12 s , the two small humps have
been substantially eroded after the dam-break wave passing through and sediment has been transported to downstream with the wave front. The eroding and depositing processes continue as the dam-break wave propagates further downstream and interacts with the topography, as shown in Fig.3 for t = 30 s . Apart from the erosion to the small humps, erosion also occurs to the part of the big mount that is affected by the flow. Substantial deposit of sediment is found at the downstream end of the domain where the flow stops and reverses its direction. As the simulation time increases, the momentum of the flow is dissipated by bed friction as well as the interaction with the erodible bed, which causes both flow and sedimentation processes to slow down, as for t = 60 s . The numerical simulation reproduces the complex physical process associated with the dam break over erodible bed. The numerical solutions are found to capture smoothly the wet-dry fronts and the sediment moves according to the flow hydrodynamics, which is physically correct. Since this test case involves most of the features that are found in a realistic situation, e.g. complex domain topography, flow discontinuity, wetting and drying and bed roughness, successfully handling of this challenging test essentially indicates the potential of the current model for real-world applications.
280
Fig.8 Dam break over three humps: results at t = 60 s
4. Conclusions This paper presents a new morphodynamic model that is capable of simulating complex flows over movable topographies involving wetting and drying. The model solves a coupled morphodynamic system consisting of the shallow water and Exner equations using a finite volume Godunov-type scheme. An HLLC approximate Riemann solver is used to evaluate the interface fluxes for both hydrodynamics and sediment transport. The wave speeds that facilitate the implementation of the Riemann solver are estimated according to a set of approximate eigenvalues to the flux Jacobians of the coupled governing equations. Wetting and drying are treated by a nonnegative reconstruction technique that predicts nonnegative water depth over complex domain topographies. As a whole, the new model is able to predict bed load transport and morphological changes caused by complex flow over irregular topographies involving wetting and drying. After being validated against two analytical test cases, the model is applied to simulate the morphodynamics caused by a dam-break flow over three humps. Since this test involves most of the essential factors in practical applications such as nonuniform topography, complex flow hydrodynamics, wetting and drying and bed roughness, etc., the model
is therefore suited for more complicated real-world simulations. Acknowledgements The author thanks Dr. Soares-Frazão Sandra from Université Catholique de Louvain for providing useful information on sediment transport modelling and the eigenstructure of the coupled governing equations. References [1]
[2]
[3]
[4] [5]
YAN Jun, CAO Zhi-xian and LIU Huai-han et al. Experimental study of landslide dam-break flood over erodible bed in open channels[J]. Journal of Hydrodynamics, 2009, 21(1): 124-130. KUBATKO E. J., WESTERINK J. J. and DAWSON C. An unstructured grid morphodynamic model with a discontinuous Galerkin method for bed evolution[J]. Ocean Modelling, 2006, 15(1-2): 71-89. CUNGE J. A., HOLLY F. M. and VERWAY A. Practical aspects of computational river hydraulics[M]. London, UK: Pitman Publishing Limited, 1980. DE VRIEND H. J., ZYSERMAN J. and NICHOLSON J. et al. Medium-term 2DH coastal area modelling[J]. Coastal Engineering, 1993, 21(1-3): 193-224. HUDSON J., DAMGAARD J. and DODD N. et al. Numerical approaches for 1D morphodynamic
281
[6]
[7]
[8]
[9]
[10]
[11]
modelling[J]. Coastal Engineering, 2005, 52(8): 691707. LIU X., LANDRY B. J. and GARCÍA M. H. Twodimensional scour simulations based on coupled model of shallow water equations and sediment transport on unstructured meshes[J]. Coastal Engineering, 2008, 55(10): 800-810. BENKHALDOUN F., SAHMIM S. and SEAÏD M. A two-dimensional finite volume morphodynamic model on unstructured triangular grids[J]. International Journal for Numerical Methods in Fluids, 2010, 63(11): 1296-1327. MURILLO J., GARCÍA-NAVARRO J. An Exner-based coupled model for two-dimensional transient flow over erodible bed[J]. Journal of Computational Physics, 2010, 229(23): 8704-8732. SOARES-FRAZÃO S., ZECH Y. HLLC scheme with novel wave-speed estimators appropriate for twodimensional shallow-water flow on erodible bed[J]. International Journal for Numerical Methods in Fluids, 2010, 66(8): 1019-1036. ýRNJARIû-ŽIC N., VUKOVIû S. and SOPTA L. Extension of ENO and WENO schemes to one- dimensional sediment transport equations[J]. Computers and Fluids, 2004, 33(1): 31-56. CALEFFI V., VALIANI A. and BERNINI A. Highorder balanced CWENO scheme for movable bed shallow water equations[J]. Advances in Water Resources, 2007, 30(4): 730-741.
[12]
[13]
[14]
[15]
[16] [17]
TASSI P. A., RHEBERGEN S. and VIONNET C. A. et al. A discontinuous Galerkin finite element model for river bed evolution under shallow flows[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(33-40): 2930-2947. MIRABITO C., DAWSON C. and KUBATKO E. J. et al. Implementation of a discontinuous Galerkin morphological model on two-dimensional unstructured meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(1-4): 189-207. LIANG Q., BORTHWICK A. G. L. Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography[J]. Computers and Fluids, 2009, 38(2): 221-234. ZHOU J. G., CAUSON D. M. and MINGHAM C. G. et al. The surface gradient method for the treatment of source terms in the shallow-water equations[J]. Journal of Computational Physics, 2001, 168(1): 1-25. LIANG Q. Flood simulation using a well-balanced shallow flow model[J]. Journal of Hydraulic Engineering, ASCE, 2010, 136(9): 669-675. KUBATKO E. J., WESTERINK J. J. Exact discontinuous solutions of exner’s bed evolution model: Simple theory for sediment bores[J]. Journal of Hydraulic Engineering, ASCE, 2007, 133(3): 305-311.
2011,23(3):273-281 DOI: 10.1016/S1001-6058(10)60113-8
A COUPLED MORPHODYNAMIC MODEL FOR APPLICATIONS INVOLVING WETTING AND DRYING* LIANG Qiuhua School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne NE1 7RU, England, UK, E-mail: [email protected] (Received April 22, 2011, Revised May 23, 2011) Abstract: This work presents a new finite volume Godunov-type model for predicting morphological changes under the rapidly varying flood conditions with wetting and drying. The model solves the coupled shallow water and Exner equations, with the interface fluxes evaluated by an Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver. Well-balanced solution is achieved using the surface gradient method and wetting and drying are handled by a non-negative reconstruction approach. The new model is validated against several theoretical benchmark tests and promising results are obtained. Key words: morphodynamic model, Exner equation, shallow water equations, Godunov-type method, wetting and drying, wellbalanced scheme
Introduction Morphodynamics concerns the evolution of the channel bathymetry in response to sediment activities driven by flows. It may have great impacts on local habitats, affect river conveyance and hence have implications for flooding during storm events. An understanding of morphodynamics and the relevant processes is therefore vital for river management systems. In the past few decades, with the rapid development of computer technology, computer modelling has become a common practice in assessing river morphodynamics although it may not entirely replace the laboratory studies for understanding the relevant processes[1]. A morphodynamic modelling system is typically consisted of a hydrodynamic component that describes the flow hydrodynamics and a sediment transport/morphological component that accounts for the bed evolution[2]. Conventionally, morphodynamic modelling is carried out in a decoupled way where the flow hydrodynamics is first obtained and then used to drive sediment transport and update the channel bed * Project supported by the UK Engineering and Physical Sciences Research Council (Grant No. EP/F030177/1). Biography: LIANG Qiu-hua (1974-), Male, Ph. D., Lecturer
change[2-4]. Ignoring the unsteady hydrodynamic effects, this approach has been recognized to be more appropriate for slow-varying flows or quasi-steady processes[5]. For a rapidly-varying flow (e.g., dam breaks, floods, etc.) where the flow unsteadiness can strongly impact the sediment transport, a coupled model, in which the flow hydrodynamics and morphodynamics are solved instantaneously, is desirable and may lead to more reliable predictions[5]. A number of coupled models that solve the integrated shallow water and Exner equations have been reported in recent years for simulating morphodynamics. For example, Liu et al.[6] presented a Roe’s approximate Riemann solver based finite volume Godunov-type morphodynamic model on unstructured grids with applications to coastal waters. They also compared the results with those produced by the quasi-steady approach and confirmed the necessity of using a coupled approach for a problem involving strong fluid-sediment interactions. Benkhaldoun et al.[7] reported an alternative finite volume Godunov-type model combined with the Roe’s approximate Riemann solver on adaptive triangular grids. The authors also compared the performance of the coupled and decoupled methods and indicated that the coupled approach out-performed the quasi-steady approach when the bed-load rapidly interacts with the hydrodynamics.
274
Murillo and García-Navarro[8] formulated the sediment fluxes in the Exner equation in a universal form for a number of different empirical formulae and solved the coupled system using a first-order finite volume Roe-type scheme. Soares-Frazão and Zech[9] presented a first-order finite volume Godunov-type morphodynamic model using an Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver, where new wave-speed estimates were derived based on the novel approximate eigenvalues to the coupled governing equations. Attempts have also been made to solve the coupled system using Essentially NonOscillatory (ENO) and Weighted Essentially NonOscillatory (WENO) schemes[10], higher-order Central Weighted Essentially Non-Oscillatory (CWENO) scheme[11] and finite element discontinuous Galerkin method[12,13]. However, most of these existing models concerns only wet-bed applications, i.e., the problem domains are covered entirely by water. A morphodynamic model that can handle applications involving wetting and drying may lead to wider applications. Therefore, this work aims to develop a new 2-D shallow flow model to include morphological change under the rapidly varying flood conditions that involves wetting and drying. 1. Governing equations The coupled governing equations, consisting of the 2-D SWEs and the Exner equation, describing shallow flows over erodible bed may be written in a matrix form as wq wf wg + + =s (1) wt wx wy where t , x and y are the time and the two Cartesian coordinates, q denotes the vector of flow variables, f and g are the flux vectors in the x and y -direction, and s is the vector containing the source terms. The vectors are given by ª qy º ª qx º « qq » « » 2 2 ªhº x y « » « qx + gh » «q » « » « » h h 2 x» « , g=« 2 q= , f =« », » 2 «qy » qx q y » « q y + gh » « « » «h « » h 2 » ¬z¼ « » « q » «¬ qsy »¼ sx ¬ ¼
ª « 2 « C f qx (qx « h3 s=« 2 « C f q y (qx « h3 « «¬
º » 2 + q y ) ghwz » wx » » 2 + q y ) ghwz » wy »» »¼ 0
where h is the water depth, qx = uh and qv = vh are the unit-width discharge in the x and y -direction with u and v being the corresponding depth-averaged velocities, z represents the bed elevation above an arbitrary datum, g = 9.81 m / s 2
denotes the acceleration due to gravity, qs x = ] Ag qx <
q
2 x
+ q 2y / h3 and qs y = ] Ag q y qx2 + q 2y / h3 are the
sediment fluxes, ] = 1 / (1 p ) is a fixed parameter with p being the bed porosity, Ag (0 Ag d 1) is an empirical constant that indicates a stronger flowsediment interaction as it reaches 1, C f = gn 2 / h1/3 is the bed roughness coefficient with n being the Manning coefficient, wz / wx and wz / wy are the two bed slopes. In order to construct a proper eigenstructure, the coupled System (1) and (2) are rewritten as follows by moving the bed slope source terms to the left hand side[9] wq wf (q ) wq wg (q ) wq + + H (q ) + + K (q ) = s f wt wx wx wy wy
(3)
where
ª0 «0 H (q ) = « «0 « ¬0 ª « « « C f qx sf = « « « C f q y « «¬
0 0º 0 c 2 »» , K (q ) 0 0» » 0 0¼
0 0 0 0
ª0 «0 « «0 « ¬0
0 0 0 0
0 0º 0 0 »» , 0 c2 » » 0 0¼
º » q + q » » h3 » qx2 + q y2 »» h3 » »¼ 0 0
2 x
2 y
(4)
where c = gh . Then the coupled system may be written in a quasi-linear form as wq wq wq + A(q ) + B (q ) = sf wt wx wy
(5)
where A (q ) and B (q ) are the pseudo-Jacobian matrices given by[9]
0
(2)
§ 0 ¨ 2 c u2 wf (q ) + H (q ) = ¨ A (q) = ¨ uv wq ¨ © M
1 2u v N
0 0· ¸ 0 c2 ¸ u 0 ¸, ¸ R 0¹
275
§ 0 ¨ uv wg (q ) + K (q ) = ¨ 2 B (q ) = ¨ c v2 wq ¨ © D
0 1 0· ¸ v u 0¸ 0 2v c 2 ¸ ¸ R E 0¹
(6)
ximate Riemann solver is chosen in this work to solve the Riemann problems and it is well-documented in literature[14] for calculating the hydrodynamic fluxes in Eq.(1). Herein, the HLLC approximate Riemann solver is briefly introduced herein for completion and modified for the current application. For the first three components of f E (i.e., the hydrodynamic fluxes),
where M=
R=
E=
3] Ag u (u 2 + v 2 ) h
2] Ag uv h
, D=
, N=
] Ag (3u 2 + v 2 ) h
3] Ag v(u 2 + v 2 ) h
f E = f L if ,
f E = f*L if
,
] Ag (u 2 + 3v 2 ) h
Considering A (q ) , a set of approximate eigenvalues may be derived as[9]
O1 = 0.5 ªu c (u c)2 + 4c 2 N º , ¬ ¼ O2 = 0.5 ªu c (u c)2 + 4c 2 N º , ¬
¼
O3 = u , O4 = u + c
(7)
With these eigenvalues, we may derive the wave speeds based on which an HLLC approximate Riemann solver can be implemented. 2. Numerical model In this work, the coupled matrix System (1) is solved by a 1st order finite volume Godunov-type scheme, which updates the flow variables to a new time step using the following formula k +1 i, j
q
§ f fW g N gS · = q 't ¨ E + si , j ¸ 'y © 'x ¹ k i, j
(8)
where the superscript k represents time level, 't is the time interval, 'x and 'y are the cell dimensions, f E , fW , g N and gS are the fluxes across the four cell interfaces. The interface fluxes are evaluated by solving local Riemann problems, e.g., f E = F (qEL , qER ) in which qEL and qER are the Riemann states on either side of the cell interface under consideration and they are reconstructed from the cell-centred flow variables. The HLLC appro-
0 d SL SL d 0 d SM
(9a) (9b)
f E = f*R if S M d 0 d S R
(9c)
f E = f R if 0 t S R
(9d)
where f L = f (qEL ) , f R = f (qER ) and f*L and f*R are the HLLC fluxes in the star region of the solution structure[14]. The wave speeds are estimated by S L = u ER 2 ghER if hEL = 0
(10a)
S L = min O1L , O1* , O2 L , O2* if hEL ! 0
(10b)
S R = u EL + 2 ghEL if hER = 0
(10c)
S R = max O4 L , O4* if hER ! 0
(10d)
SM =
S L hER (u ER S R ) S R hEL (u EL S L ) hER (u ER S R ) hEL (uEL S L )
(10e)
with h* =
v* =
hEL uEL + hER u ER hEL + hER , , u* = 2 hEL + hER hEL vEL + hER vER hEL + hER
(11)
The sediment component of the fluxes may be estimated as[9] f E (4) = qsxL if 0 d S sL
(12a)
f E (4) = q*sx if S sL d 0 d S sR
(12b)
f E (4) = qsxR if 0 t S sR
(12c)
276
where the sediment flux in the star region is given by[9] q*s x =
S sR qsLx S sL qsRx + S sL S sR hEL ( z ER z EL ) S sR S sL
(13)
with the wave speeds approximated by S sL = min O1L , O1 , S sR = max O2 R , O2
(14)
In order to calculate the interface fluxes using the above HLLC approximate Riemann solver, the Riemann states must be properly reconstructed. For a first-order accurate scheme, the face values of flow variables are assumed to be the same as those at the cell centres, i.e.,
q EL = qi , j , q ER = qi +1, j , KEL = Ki , j , KER = Ki +1, j
(15)
where h = (hEL + hER ) / 2 . For the friction source terms, a point-implicit scheme proposed by Liang[16] is used to ensure non-negative water depth for simulations involving wetting and drying over complex domain topographies. The aforementioned finite volume Godunov-type morphodynamic model is overall explicit and its numerical stability is controlled by the CourantFriedrichs-Lewy (CFL) criterion. Either transmissive or reflective (slip) boundary conditions are used in the test cases considered in this work. 3. Results In this section, the aforementioned coupled morphological change model is validated against idealized tests. In all of the following test cases, the constants p = 0.4 and Ag = 0.01 are used. The
CFL number is set to 0.75.
where K = h + z is the water level, which is defined to assist the implementation of the surface gradient method[15]. In order to define the effective water depth at the cell interface and hence calculate the hydrodynamic fluxes, a single value of bed elevation is first defined z E = max( zEL , zER )
(16)
The Riemann states are then reconstructed as Fig.1 Exner’s dune: initial bed elevation L L hEL = max(0,KEL z E ) , qxE = u EL hEL , q yE = vEL hEL ,
zEL = zEL
(17)
The Riemann states on the right hand side of the eastern interface are defined accordingly. For applications involving wetting and drying over irregular topography, the following step is essential 'z = max 0, z E KEL
(18)
§ 2Sx · z0 = z ( x, 0) = A0 + A1 cos ¨ ¸ © O ¹
(21)
where A0 = A1 = 1 , O = 20 are constants. Denoting c( z ) { wqsx / wz
3] Ag qx3 / (K z )4 and assuming a
constant water level and discharge, the analytical solution is given by
and z E m z E 'z , zEL m z EL 'z , zER m z ER 'z
(19)
In the context of the above surface gradient method, the slope source terms are directly discretized by a central-differencing scheme, e.g., in the x -direction gh
3.1 Analytical test of Exner’s dune We first consider the analytical test of Exner’s dune. As shown in Fig.1, the initial bed profile is described by[17]
wz § z zW · = gh ¨ E ¸ wx © 'x ¹
(20)
ª 2S º z ( x, t ) = A0 + A1 cos « ( x c( z )t ) » ¬O ¼
(22)
which can be solved using the Newton-Raphson method. Simulations occur in a 40 m long domain for 50 s. The water level and discharge are respectively K { 3 m and qx { 1m 2 / s . Figure 2 shows the
277
numerical predictions obtained on grids of different resolution, which agree closely with the analytical solution and appear to be convergent.
analytical solution in the original paper in Ref.[17].
Fig.3 Analytical 1-D test with shock: numerical predictions at t = 3.5 s Fig.2 Exner’s dune: numerical predictions at t = 50 s where (b) is the zoom-in view
3.2 Analytical test with shock The second analytical case involves a so-called Riemann problem for the Exner equation[17]. In order to simplify the procedure of finding the analytical solution, the 1-D Exner equation takes the form of ( z )t + (D u ) x = 0 with D = ] Ag and the derivative terms M , N , R , D and E in Eqs.(6) and (7) are changed accordingly. For the initial bed profile z0 = z L , x 0
(23a)
z0 = z R , x ! 0
(23b)
the analytical solution for the case of z L z R may be derived to be[17]
z ( x, t ) = z L , x c( z L )t z ( x, t ) = K D q x
t , c( z L )t d x d c( z R )t x
z ( x, t ) = z R , x ! c( z L )t
(24a) (24b) (24c)
where c( z ) = D qx / (K z ) 2 for the simplified Exner equation. Please note the typing mistake of the
Fig.4 Dam break over three humps: initial conditions
Simulation occurs at a [– 1 m, 5 m] domain with z L = 1 m , zR = 2 m , K { 3 m and qx { 1 m 2 / s .
278
Fig.5 Dam break over three humps: results at t = 6 s
Fig.6 Dam break over three humps: results at t = 12 s
279
Fig.7 Dam break over three humps: results at t = 30 s
The numerical results obtained after t = 3.5 s on uniform meshes of different resolution are plotted in Fig.3, comparing with the analytical solution. The original discontinuous bed profile is found to be washed out by the current, where the numerical predictions again agree closely with the analytical solution and seem to be convergent with increasing grid refinement. Therefore, the capability of the current coupled model in modelling discontinuous sediment flux is confirmed. 3.3 Dam break over three humps This test case is about a dam-break flow running over three humps in the originally dry floodplain. The initial flow conditions are indicated in Fig.4 and the detailed settings can be found in Ref.[14]. Figures 5-8 present the numerical predi- ctions at different output times in terms of 3-D surface view and depth contours for both fixed and movable bed tests. The contour lines of eroded bed topography for the movable bed test are also illustrated. It is evident that the dam-break flow induces substantial change to the floodplain topography meanwhile the morphological changes also feedback and modify the flow hydrodynamics. As shown for t = 6 s , the dam- break wave rushes downstream after the dam has been removed, carrying great momentum to reshape the landform. After t = 12 s , the two small humps have
been substantially eroded after the dam-break wave passing through and sediment has been transported to downstream with the wave front. The eroding and depositing processes continue as the dam-break wave propagates further downstream and interacts with the topography, as shown in Fig.3 for t = 30 s . Apart from the erosion to the small humps, erosion also occurs to the part of the big mount that is affected by the flow. Substantial deposit of sediment is found at the downstream end of the domain where the flow stops and reverses its direction. As the simulation time increases, the momentum of the flow is dissipated by bed friction as well as the interaction with the erodible bed, which causes both flow and sedimentation processes to slow down, as for t = 60 s . The numerical simulation reproduces the complex physical process associated with the dam break over erodible bed. The numerical solutions are found to capture smoothly the wet-dry fronts and the sediment moves according to the flow hydrodynamics, which is physically correct. Since this test case involves most of the features that are found in a realistic situation, e.g. complex domain topography, flow discontinuity, wetting and drying and bed roughness, successfully handling of this challenging test essentially indicates the potential of the current model for real-world applications.
280
Fig.8 Dam break over three humps: results at t = 60 s
4. Conclusions This paper presents a new morphodynamic model that is capable of simulating complex flows over movable topographies involving wetting and drying. The model solves a coupled morphodynamic system consisting of the shallow water and Exner equations using a finite volume Godunov-type scheme. An HLLC approximate Riemann solver is used to evaluate the interface fluxes for both hydrodynamics and sediment transport. The wave speeds that facilitate the implementation of the Riemann solver are estimated according to a set of approximate eigenvalues to the flux Jacobians of the coupled governing equations. Wetting and drying are treated by a nonnegative reconstruction technique that predicts nonnegative water depth over complex domain topographies. As a whole, the new model is able to predict bed load transport and morphological changes caused by complex flow over irregular topographies involving wetting and drying. After being validated against two analytical test cases, the model is applied to simulate the morphodynamics caused by a dam-break flow over three humps. Since this test involves most of the essential factors in practical applications such as nonuniform topography, complex flow hydrodynamics, wetting and drying and bed roughness, etc., the model
is therefore suited for more complicated real-world simulations. Acknowledgements The author thanks Dr. Soares-Frazão Sandra from Université Catholique de Louvain for providing useful information on sediment transport modelling and the eigenstructure of the coupled governing equations. References [1]
[2]
[3]
[4] [5]
YAN Jun, CAO Zhi-xian and LIU Huai-han et al. Experimental study of landslide dam-break flood over erodible bed in open channels[J]. Journal of Hydrodynamics, 2009, 21(1): 124-130. KUBATKO E. J., WESTERINK J. J. and DAWSON C. An unstructured grid morphodynamic model with a discontinuous Galerkin method for bed evolution[J]. Ocean Modelling, 2006, 15(1-2): 71-89. CUNGE J. A., HOLLY F. M. and VERWAY A. Practical aspects of computational river hydraulics[M]. London, UK: Pitman Publishing Limited, 1980. DE VRIEND H. J., ZYSERMAN J. and NICHOLSON J. et al. Medium-term 2DH coastal area modelling[J]. Coastal Engineering, 1993, 21(1-3): 193-224. HUDSON J., DAMGAARD J. and DODD N. et al. Numerical approaches for 1D morphodynamic
281
[6]
[7]
[8]
[9]
[10]
[11]
modelling[J]. Coastal Engineering, 2005, 52(8): 691707. LIU X., LANDRY B. J. and GARCÍA M. H. Twodimensional scour simulations based on coupled model of shallow water equations and sediment transport on unstructured meshes[J]. Coastal Engineering, 2008, 55(10): 800-810. BENKHALDOUN F., SAHMIM S. and SEAÏD M. A two-dimensional finite volume morphodynamic model on unstructured triangular grids[J]. International Journal for Numerical Methods in Fluids, 2010, 63(11): 1296-1327. MURILLO J., GARCÍA-NAVARRO J. An Exner-based coupled model for two-dimensional transient flow over erodible bed[J]. Journal of Computational Physics, 2010, 229(23): 8704-8732. SOARES-FRAZÃO S., ZECH Y. HLLC scheme with novel wave-speed estimators appropriate for twodimensional shallow-water flow on erodible bed[J]. International Journal for Numerical Methods in Fluids, 2010, 66(8): 1019-1036. ýRNJARIû-ŽIC N., VUKOVIû S. and SOPTA L. Extension of ENO and WENO schemes to one- dimensional sediment transport equations[J]. Computers and Fluids, 2004, 33(1): 31-56. CALEFFI V., VALIANI A. and BERNINI A. Highorder balanced CWENO scheme for movable bed shallow water equations[J]. Advances in Water Resources, 2007, 30(4): 730-741.
[12]
[13]
[14]
[15]
[16] [17]
TASSI P. A., RHEBERGEN S. and VIONNET C. A. et al. A discontinuous Galerkin finite element model for river bed evolution under shallow flows[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(33-40): 2930-2947. MIRABITO C., DAWSON C. and KUBATKO E. J. et al. Implementation of a discontinuous Galerkin morphological model on two-dimensional unstructured meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(1-4): 189-207. LIANG Q., BORTHWICK A. G. L. Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography[J]. Computers and Fluids, 2009, 38(2): 221-234. ZHOU J. G., CAUSON D. M. and MINGHAM C. G. et al. The surface gradient method for the treatment of source terms in the shallow-water equations[J]. Journal of Computational Physics, 2001, 168(1): 1-25. LIANG Q. Flood simulation using a well-balanced shallow flow model[J]. Journal of Hydraulic Engineering, ASCE, 2010, 136(9): 669-675. KUBATKO E. J., WESTERINK J. J. Exact discontinuous solutions of exner’s bed evolution model: Simple theory for sediment bores[J]. Journal of Hydraulic Engineering, ASCE, 2007, 133(3): 305-311.