- Email: [email protected]
Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help
Advances in Water Resources 34 (2011) 980–989
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help S. Cordier a,⇑, M.H. Le a,b, T. Morales de Luna c a
Laboratoire de Mathématiques - Analyse, Probabilités, Modélisation - Orléans (MAPMO), UMR CNRS 6628, Université d’Orléans , France Bureau de Recherches Géologiques et Minières (BRGM), France c Departamento de Matemáticas, Universidad de Córdoba, Spain b
a r t i c l e
i n f o
Article history: Received 23 December 2010 Received in revised form 8 May 2011 Accepted 8 May 2011 Available online 23 May 2011 Keywords: Shallow water system Exner equation Sediment transport Hyperbolicity Stability Splitting methods
a b s t r a c t In this paper, we are concerned with sediment transport models consisting of a shallow water system coupled with the so called Exner equation to describe the evolution of the topography. We show that, for some bedload transport models like the well-known Meyer-Peter and Müller model, the system is hyperbolic and, thus, linearly stable, only under some constraint on the velocity. In practical situations, this condition is hopefully fulfilled. Numerical approximations of such system are often based on a splitting method, solving first shallow water equation on a time step and, updating afterwards the topography. It is shown that this strategy can create spurious/unphysical oscillations which are related to the study of hyperbolicity. Using an upper bound of the largest eigenvector may improve the results although the instabilities cannot be always avoided, e.g. in supercritical regions. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Soil erosion is caused by the movement of sediment due to mechanical actions of flow. In the context of bedload transport, a mass conservation law also called Exner equation [1] is used to update the bed elevation. This equation is often coupled with the shallow water equations describing the overland flows (see [2] and references therein). The complete system of PDE may be written in the form
8 > < @ t h þ @ x ðhuÞ ¼ 0; 2 2 @ t ðhuÞ þ @ x ðhu þ gh =2Þ ¼ gh@ x zb þ gh@ x H; > : @ t zb þ n@ x qb ¼ 0;
ð1Þ
where h is the water depth, u is the flow velocity and zb is the thickness of sediment layer which can be modified by the fluid. The sediment layer is assumed to be located on a fixed bedrock layer at depth H which is not modified by the fluid (see Fig. 1). n = (1 c)1 with c the porosity of the sediment layer. qb is the volumetric bedload sediment transport rate per unit time and width and g is the acceleration due to gravity. The conservative variable hu is also called water discharge and noted by q. In what follows and for the sake of simplicity, we shall assume that H is constant ⇑ Corresponding author. E-mail addresses: [email protected] (S. Cordier), [email protected] (M.H. Le), [email protected] (T. Morales de Luna). 0309-1708/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2011.05.002
so that the term @ xH can be neglected. We shall also assume that n = 1. Let us first quote some of the expressions proposed for the bedload transport rate in the literature. Many researches have developed different formulae for predicting and estimating bedload transport rate. Their expressions are usually proposed for granular non-cohesive sediments and quantified empirically. One of the simplest expression is the one proposed by Grass [3] which considers qb as a function of the flow velocity and a constant coefficient Ag which depends on soil properties
Grass ð1981Þ :
qb ¼ Ag ujujmg 1 ;
1 6 mg 6 4:
ð2Þ
In practice, estimations of bedload transport rate are mainly based on the bottom shear stress sb, i.e. the force of water acting on the bed during its routing. In the context of laminar flows, the bottom shear stress is given as
sb ¼ qghSf ;
ð3Þ
where q is the water density and Sf is the friction term that can be quantified by different empirical laws such as the Darcy–Weisbach or Manning formulae (see [4]) Darcy–Weisbach:
Sf ¼
fujuj ; 8gh
where f is the Darcy–Weisbach’s coefficient.
ð4Þ
981
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
0.07
Bedload rate (m^2/s)
0.06
Meyer-Peter & Muller Fernandez Luque & Van Beek Nielsen Ribberink Camenen & Larson
0.05 0.04 0.03 0.02 0.01 0
0
0.05
0.1
0.15
n2 ujuj h
4=3
ð5Þ
;
where n is the Manning’s coefficient. The bottom shear stress is usually used in dimensionless form, noted sH b , which is also called Shields parameter. It is defined in terms of the ratio between drag forces and the submerged weight by
sHb ¼
0.25
0.3
0.35
0.4
0.45
0.5
Fig. 2. Bedload rate as function of Shields parameter with ds = 32 mm, q = 1000 kg/ m3, qs = 2600 kg/m3 and sH cr ¼ 0:05.
Manning:
Sf ¼
0.2
Shields parameter
Fig. 1. Sketch of shallow water over an erodible bed.
sb ; ðqs qÞgds
ð6Þ
where qs is the sediment density and ds is the diameter of sediment. H The main hypothesis is that sH b must exceed a threshold value scr in order to initiate motion. The threshold value sH depends on the cr physical properties of sediment and is usually computed experimentally. One of the first works on this topic was done by Shields [5] in which sH cr is determined in relation with the boundary Reynolds number. Bedload transport rate may be represented as a function of sH b via a non-dimensional parameter U by
s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qs 3 qb ¼ U 1 gds :
ð7Þ
q
The followings expressions, illustrated by Fig. 2, have been often applied [6–10]:
€ ller ð1948Þ : U ¼ 8 Meyer-Peter & Mu
sHb sHcr
ndez Luque&Van Beek ð1976Þ : U ¼ 5:7 Ferna qffiffiffiffiffiffi H H Nielsen ð1992Þ : U ¼ 12 sH b sb scr þ ; H 1:65 Ribberink ð1998Þ : U ¼ 11 sH b scr þ ;
3=2
þ
ð8Þ
;
sHb sHcr
3=2 þ
or one could consider a modification of the classical formulae by including gravity effects as it was done in [11]. As we have said and for the sake of simplicity we shall assume here that H(x) is constant. Remark that the general case where H(x) is non-constant can be considered by introducing the trivially satisfied equation
@t H ¼ 0
ð13Þ
and following the ideas presented in [12]. PDE system (1) can be written in vectorial form as
@ t W þ @ x FðWÞ ¼ BðWÞ@ x W;
ð14Þ
where W = (h, q, zb)T is the state vector in conservative form, t 2 2 FðWÞ ¼ hu; hu þ 2g h ; qb ðh; huÞ and B(W) = ghE23 where Eij denotes the canonical base of the 3 3 matrices. System (1) can also be written in non conservative form as
@ t W þ AðWÞ@ x W ¼ 0;
ð15Þ
where A(W) = DWF B(W) is the matrix of transport coefficients. An important property of such systems is hyperbolicity [13] which requires that the matrix A(W) is R diagonalizable (or strictly hyperbolic when eigenvalues are distinct). Let us quickly give an interpretation of this property as a stability condition for the linearized system. Assume that W is a small wave perturbation of a constant state, i.e. of the form
Wðx; tÞ ¼ W 0 þ W 1 eiðkxxtÞ ;
1:
ð16Þ
This is a solution of the linearized problem
;
ð9Þ ð10Þ ð11Þ
1:5 H Camenen and Larson ð2005Þ : U ¼ 12 sH exp 4:5sH b cr =sb : ð12Þ Remark 1. Classical formulae for bedload transport do not take into account gravity effects. As a consequence, particles located at the advancing front of a dune do not fall due to gravity as the motion of particles only depends on the hydrodynamical variables. This means that vertical profiles may be observed in numerical simulations which are not found in physical situations, for example at the front of an advancing dune. In order to obtain more realistic profiles, a diffusion term may be added to the third equation of (1)
ixW 1 þ AðW 0 ÞðikÞW 1 ¼ 0
ð17Þ
if and only if W1 is an eigenvector associated to the eigenvalue xk . In other words, if one can write any arbitrary perturbation W1 as a linear combination of such eigenvector associated with real velocity, the solution will propagate without amplification of the perturbation. On the contrary, if any of the eigenvalues is complex, this will lead to instability. Let us remark that these eigenvalues are also important when using explicit upwind schemes to ensure stability. Indeed, the time step, for a constant mesh size Dx, has to satisfy the so called Courant–Friedichs–Lewy condition (CFL condition in what follows) (see [13,14] and references therein)
jkmax jDt < Dx;
ð18Þ
where jkmaxj is the maximum of the modulus of the eigenvalues or in other words, kmax is the maximum velocity for propagation of
982
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
information. We recall that when using finite volume schemes technique for (14), this condition can be seen as the definition of a time step sufficiently small so that the different Riemann problems at each intercell do not interact between each other, that is, information of each Riemann problem does not cross more than one cell. This paper is organized as follow: in next section, we study the hyperbolicity of system (1) and establish the main result. Next, we compare two different strategies for numerically solving (1): a time-splitting strategy and a coupled scheme. While some works have already been done comparing coupled and uncoupled schemes, we show how information on the eigenvalues may give better results. Finally, we present some numerical tests to justify the arguments mentioned.
Let us find relations between a and b in the different models proposed before: (2) and (8)–(12). Grass: qb = Agujujm1
u ¼ hq ) @qb @q
2
0
6 AðWÞ ¼ 4 gh u2 a
1
0
3
7 2u gh 5; b
ð19Þ
0
b b with a ¼ @q and b ¼ @q . The characteristic polynomial of A(W) can @h @q be written as
pA ðkÞ ¼ k
k 1 k 1 gh a b u þ gh 2u k
@u ; @q
@u @h @qb @h
)
¼ hq2
¼
ðqb Þ0u
)
@u @h
@qb q @qb ¼ : h @q @h
@qb @q
sb ¼ ðqb Þ0sb @@q ¼ ðqb Þ0sb ðsb Þ0u @u @q
@qb @h
sb ¼ ðqb Þ0sb @@h ¼ ðqb Þ0sb ðsb Þ0u @u @h
)
)
@qb q @qb ¼ h @q @h
Manning : 2 ujuj qjqj sb ¼ qgh nh4=3 ¼ a hujuj 1=3 ¼ a 7=3 h @ sb @q
jqj ¼ 2ha7=3
@ sb @h
¼ 76
9 =
q 2ajqj h h7=3
¼ 76
q @ sb h @q
;
)
@qb 7 q @qb ¼ : 6 h @q @h
Thus, we can summarize that we have two types of relation: a = ub or a ¼ 76 ub. We shall assume that b > 0 since sediment rate increases with that of the flow. The line d(k) can be rewritten as d(k) = ghb(k ku) with k = 1 or k = 7/6 and we shall assume the slope ghb > 0. Proposition 2.1. Consider system (1) with qb = qb(h, q) such that
@qb q @qb ¼ k : h @q @h
ð21Þ
For a given state (h, q), the system is strictly hyperbolic if and only if
a < ku < aþ ;
ð22Þ
where a± is defined by expression (26). More explicitly, In the case k = 1 (or a = ub), the system is always strictly hyperbolic. In the case k = 7/6 (or a ¼ 76 ub), a sufficient condition for system (1) to be strictly hyperbolic is
pffiffiffiffiffiffi juj < 6 gh:
2
¼ k½ðu kÞ2 gh þ ghðbk þ aÞ:
ðqb Þ0u
¼ 1h ;
Meyer-Peter and Müller, Fernández Luque and Van Beek, Nielsen, Ribberink, Camenen and Larson: qb = qb(sb) Darcy–Weisbach : sb ¼ qgh fujuj ¼ aujuj ða ¼ cstÞ 8gh
2. Domain of hyperbolicity We are concerned here by the hyperbolicity of system (1) for the different models proposed (2)–(11) for the bedload transport rate qb. It is known that, for Grass model, the system is always hyperbolic [2]. Nevertheless, to our knowledge, the study the hyperbolicity of system (1) for more general bedload transport fluxes has not been done. The usual approach in such cases relies b n on the typically small magnitude of @q which permits a perturba@u h tion solution of the cubic characteristic equation. This approach allows to give approximations of the eigenvalues of the system in that framework, see for instance [15–17]. Nevertheless this does not answer the fact of whether the model is always hyperbolic or not. In fact, Castro et al. [2] have stated numerically that hyperbolicity may be lost in some cases for the Meyer-Peter and Müller bedload transport flux. Moreover, it would be interesting to have an easier way to know if one would find complex eigenvalues or not rather than computing directly the roots of characteristic polynomial. Proposition 2.1 will give an answer to this fact and states a necessary and sufficient condition for system (1) to be hyperbolic for the different bedload transport fluxes proposed before. Note that qb = qb(h, q) and thus the matrix of transport coefficients is given by
¼
@u @q
ð23Þ
ð20Þ
System (1) is thus strictly hyperbolic if and only if pA(k) has three different solutions noted by k1 < k2 < k3. In other words if the curve f(k) = k[(u k)2 gh] and the line d(k) = gh(bk + a) have three points of intersection. This is illustrated in Fig. 3 for the case of a subcritical flow.
Proof. For simplicity, let us consider the case u > 0 so a = kub < 0. The case u < 0 can be treated with the same arguments. We define the two tangents of the curve f(k) which are parallel to d(k). Their intersections with f(k) are characterized by f0 (k) = ghb which yields to two values of k of the form def
k ¼
2u
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ 3ðgh þ ghbÞ : 3
ð24Þ
The two tangents are such that d±(k±) = f(k±). This implies the equations for the tangents are given by
d ðkÞ ¼ ghbðk a Þ;
ð25Þ
with def
a ¼ k Fig. 3. Eigenvalues of the transport coefficients matrix.
f ðk Þ : ghb
ð26Þ
The roots of pA(k) which correspond to the eigenvalues of A(W) are given as the intersection of f(k) and d(k) (see Fig. 3). Recall that
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
n pffiffiffiffiffiffio f(k) is a third order polynomial with roots 0; u gh . The equation pA(k) = d(k) f(k) will have 3 distinct solutions if and only if the line d(k) lies in between d-(k) and d+(k). This can be equivalently written as a- < ku < a+. pffiffiffiffiffiffi Itpcan ffiffiffiffiffiffi be checked that we always have a < u gh < u < u þ gh < aþ so the systempis hyperbolic =ffiffiffiffiffi 1.ffi ffiffiffiffiffialways ffi pffiffiffiffiffiffiin the case kp In the case k = 7/6, if juj < 6 gh, we have u gh < 76 u < u þ gh and thus the hyperbolicity condition is verified which concludes the proof. Note that similar arguments have been used by one of the author in [18–20] to characterize the hyperbolicity domains for systems of moments equations arising in plasma physics. h
Remark 2. Condition (22) is interesting because it gives a necessary and sufficient condition to check if the model remains hyperbolic without computing the eigenvalues. Moreover, condition (22) shows that the model is indeed hyperbolic in most physical situations. For instance, Fig. 4 shows the region where condition (22) is satisfied for the Meyer-Peter and Müller formula with values sHcr ¼ 0:47, qs = 2612 and n = 0.0196. 3. Time-splitting method The strategy of time-splitting consists in solving separately the shallow water system during a first step for a fixed topography zb, and updating afterward the topography on a second step using the third equation. See, for example, the references [21,22] or Chapter 3 of [23]. It seems natural such approach as the variation of topography is slow compared to the characteristic time associated to hydrodynamical variables. Using a time-splitting strategy allows to minimize computational costs and makes easier to implement the codes for different bedload transport fluxes in a modular way. Moreover, time-splitting techniques allows to use robust and efficient numerical schemes already developed for shallow water system like, for example, those presented in the publication [24]. In particular this would be an interesting approach for the code FullSWOF_2D (licence DL 03434-01) developed in the framework of the multidisciplinary project METHODE ANR-07-BLAN-0232 (see [25] and http://www.univ-orleans.fr/mapmo/methode/). However, it is a known fact that such an approach is not a good one as this may re-
sult in numerical instabilities, see for instance [26]. The main reason is related to the bad approximation of the true eigenvalues of the system. We shall propose in next section different numerical tests in order to study how reducing CFL condion or the new upper bound (29) may help. More explicitly: 1. As explained in the introduction, the stability condition for shallow water equations only requires the CFL limitation associated with
juj þ
pffiffiffiffiffiffi gh Dt < Dx:
ð27Þ
Remark that the CFL condition (27) is not enough to ensure the stability condition (18) for the coupled system (1) because the maximum eigenvalue k3 of the matrix pffiffiffiffiffiffi of transport coefficients A(W) is always larger than juj þ gh (see Fig. 3). This may lead to numerical instabilities in some numerical simulations when a splitting technique is used. These instabilities are less frequent or do not exist if the numerical scheme used is sufficiently diffusive. For instance, when Lax–Friedrichs scheme is used it is very difficult to find instabilities while for Roe type schemes they are easily found. We recall that something similar has been studied in [27,28] where the splitting technique for the two-layer shallow water system was studied. In [26] a comparison between coupled and decoupled schemes have also been done. Instead of computing the exact eigenvalues of A(W), an upper bound k03 of k3 may be used, defined asffi the intersection of d(k) pffiffiffiffiffi and the tangent dT(k) of f(k) at u þ gh (see Fig. 5)
pffiffiffiffiffiffi pffiffiffiffiffiffi dT ðkÞ ¼ f 0 u þ gh k u gh :
ð28Þ
k03 is thus of the form
def k03 ¼
uþ
pffiffiffiffiffiffi pffiffiffiffiffiffi gh f 0 u þ gh þ gha : pffiffiffiffiffiffi f 0 u þ gh ghb
ð29Þ
This value can be used to impose a CFL condition associated to the full coupled system (1). With this upper-bound the true CFL condition is now granted. Nevertheless we remark that now the true maximum propagation speed is overestimated. It will be interesting to study the splitting technique when this estimation is used. 2. Another problem is related to the sign of eigenvalues of the matrix A(W). It is known that there exist always a negative eigenvalue even in the supercritical case (see [15,16,26]). Indeed, rewriting the characteristic polynomial pA(k) in form
pA ðkÞ ¼ ðk k1 Þðk k2 Þðk k3 Þ; the product of three eigenvalues, in the case u > 0, satisfies
k1 k2 k3 ¼ pA ð0Þ ¼ gha < 0:
Fig. 4. Hyperbolicity region (in red) for Meyer-Peter and Müller. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
983
Fig. 5. Upper bound first order of the maximum eigenvalue.
984
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
pffiffiffiffiffiffi This means that in a supercritical flow, i.e. juj > gh, k2 and k3 are positive and consequently k1 < 0. This can be interpreted as a wave propagating upstream. The using of splitting strategy in such situation cannot account this information as we would have pffiffiffiffiffiffi take into pffiffiffiffiffi ffi 0 < u gh < u þ gh. This will lead to numerical instabilities when using splitting method with upwind schemes. 4. Numerical results
where
is assumed to be a Roe matrix for (1) in the particular case that
@qb q @qb ¼ k : h @q @h
4.1. Description of the numerical schemes Let us describe briefly the numerical schemes used hereafter. As usual, we consider a set of computing cells Ii ¼ ½xi1=2 ; xiþ1=2 ; i 2 Z. We shall assume that these cells have a constant size Dx and that xi+1/2 = iDx. xi = (i 1/2)Dx is the center of the cell Ii. Let Dt be the time step and tn = nDt. We denote by W ni the approximation of the cell averages of the exact solution
W ni ffi
1 Dx
Z
xiþ1=2
Wðx; tn Þdx:
ð30Þ
xi1=2
n
biþ1=2
pffiffiffiffi pffiffiffiffiffiffiffiffiffi g 1 hi þ hiþ1 mX Ag ¼ n pffiffiffiffi pffiffiffiffiffiffiffiffiffi ðuiþ1 Þk ðui Þmg 1k : hi hiþ1 þ hiþ1 hi k¼0
i Dt h þ n Ai1=2 W i W ni1 þ Aiþ1=2 W niþ1 W ni ; Dx ð31Þ
where
ð38Þ
For other models considered in this work, it is not always possible to obtain an explicit formula and its implementation is very costly. So, in practice the following approximation is used
biþ1=2 ¼
Splitting scheme: We use a Roe scheme for the shallow water system plus an upwinding technique based on the sign of velocity for the bedload transport equation, that is,
W inþ1=2 ¼ W ni
Remark 3. As it was stated in [2] it is possible to write a Roe matrix for to Grass model with (35) and
n
We consider two different schemes for (1):
ð37Þ
@qb ðhiþ1=2 ; hiþ1=2 uiþ1=2 Þ: @ðhuÞ
ð39Þ
As a consequence, we are using an approximated Roe matrix and we cannot use directly schemes in the form (36) but rather equivalent expressions based on the physical flux function F(W) as it is done in [2]. Remark 4. In order to compute the eigenstructure of matrix Ai+1/2, one could use the Cardano–Vieta’s formula (see [31]). Once the eigenvalues are known, one could express the eigenvectors as functions of these eigenvalues (see [32]). Nevertheless, numerical simulations in next section have been carried out by computing numerically the eigenstructure of Ai+1/2. 4.2. Numerical tests
is a Roe matrix for the shallow water system with non-flat topography, and then nþ1
hi
nþ1=2
¼ hi
;
nþ1=2
unþ1 ¼ ui i
ðzb Þnþ1 ¼ ðzb Þni i
;
ð32Þ
i Dt h nþ1 n ðqb Þnþ1 iþ1=2 ðqb Þi1=2 ; Dx
ð33Þ
with
ðqb Þnþ1 iþ1=2
8 nþ1 nþ1 > ; if uiþ1=2 P 0; < qb hi ; qi ¼ > : q hnþ1 ; qnþ1 ; if u iþ1=2 < 0 b iþ1 iþ1
In what follows and unless said otherwise, the computational domain is discretized using 400 cells. Boundary conditions have been set by defining a ghost cell duplicated from the first and last cell, respectively for each time step. We shall take the usual value mg = 3 for Grass model. First, we shall show a case where the splitting technique gives the same result as the fully coupled scheme. Let us consider as
ð34Þ
and ui+1/2 the velocity at the interface of the cells i and i + 1, for instance the value used in the Roe’s matrix. In particular, the usual definitions
hiþ1=2 ¼
hi þ hiþ1=2 ; 2
uiþ1=2 ¼
pffiffiffiffi pffiffiffiffiffiffiffiffiffi hi ui þ hiþ1 uiþ1 pffiffiffiffi pffiffiffiffiffiffiffiffiffi hi þ hiþ1
ð35Þ
is taken, which corresponds to the choice of segments as the path connecting two different states. See [27,29,30] for further details. Coupled Scheme: We use a Roe scheme for the full system which can be written in the form
W nþ1 ¼ W ni i
Dt þ n Ai1=2 W i W ni1 þ Aiþ1=2 ðW niþ1 W ni Þ ; Dx ð36Þ
Fig. 6. Initial condition for (40).
985
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
Fig. 7. Solution for (40) with Ag = 0.005.
Table 1 Comparison at time t = 10 between numerical solutions computed with a Roe scheme for the full system and a splitting scheme and a reference solution computed on a grid with 3200 points. # Points
50
100
200
400
Coupled L1 error L1 error L1 error
0.1056 0.1026 0.0385
0.0534 0.0503 0.0204
0.0264 0.0244 0.0107
0.0126 0.0114 0.0051
0.1021 0.1067 0.0330
0.0503 0.0534 0.0167
0.0240 0.0269 0.0084
0.0111 0.0135 0.0038
scheme in h in hu in z
Splitting L1 error in h L1 error in hu L1 error in z
initial condition a subcritical steady state for the classical shallow water system given by
8 huðx; t ¼ 0Þ ¼ 0:5; > < 2 zb ðx; t ¼ 0Þ ¼ 0:1 þ 0:1eðx5Þ ; > : u2 þ gðh þ zÞ ¼ 6:386: 2
ð40Þ
The water surface and topography corresponding to this initial condition are shown in Fig. 6. We shall consider Grass model (2) with Ag = 0.005. As we see in Fig. 7, no major differences are observed between the two schemes. Both simulations have been computed
Fig. 9. Solution for (40), Ag = 0.07 with a coupled scheme.
with CFL = 0.95 and in the case of the splitting scheme, (27) has been used for the CFL condition. Table 1 shows the comparison between the solutions computed with both schemes and a reference solution. We remark that CPU time for the splitting technique is
. Fig. 8. Solution for (40) with Ag = 0.07.
.
986
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
Fig. 10. Comparison of splitting scheme using (29) and coupled scheme for (40) with Ag = 0.07 and CFL = 0.5.
Fig. 11. Eigenvalues for (40). True eigenvalues of the system (continuous line) and u
pffiffiffiffiffiffi gh approximation (dashed line).
Fig. 12. An example using Meyer-Peter and Müller bedload transport flux with a splitting technique and (27).
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
987
Fig. 13. An example using Meyer-Peter and Müller bedload transport flux.
Fig. 14. Transcritical steady-state for shallow water system.
smaller (16.73 s) than the CPU time for coupled scheme (19.06 s) which is due to the fact that the spectral decomposition of the
matrices Ai+1/2 is computed numerically in the case of the coupled scheme. Now, consider the same initial condition (40) but set Ag = 0.07. We remark in Fig. 8 that some instabilities arise when using the splitting scheme with (27) which disappear if the CFL is reduced. This is not the case when the coupled scheme is used as we see in Fig. 9. Doing the same computation with the splitting scheme and using now (29) gives similar results but now CFL = 0.95 can be used as it is shown in Fig. 10. Remark that still some oscillations may be seen in the topography with the splitting scheme. The main difference between the experiment with Ag = 0.005 and Ag = 0.07 is that, in the second case, the true eigenvalues pffiffiffiffiffiffi of the p full ffiffiffiffiffiffi system are far from the approximations u þ gh and u gh given by the splitting scheme. In particular, the pffiffiffiffiffimaximum ffi of the modulus of the eigenvalue is larger than juj þ gh. This implies that a CFL condition close to 1 should not be used with the splitting scheme (see Fig. 11) This the reason why the upper bound (29) works in general better than (27). In previous examples, we have used Grass model. Nevertheless, conclusions are not exclusive to such model and may be observed if we use any other formula for the bedload transport flux. The main difference observed between the different models is the strength of the interaction between the fluid and sediment.
Fig. 15. Solution for (41) with Ag = 0.0005 using splitting scheme.
988
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
Fig. 16. Solution for (41) with Ag = 0.0005 using full system.
state is reached, we let the sediment to evolve by using a Grass bedload flux with Ag = 0.0005. Fig. 15 shows that again we find instabilities by using the splitting technique with (27). These instabilities remain even for CFL = 0.05 (which requires almost 20 times the CPU time required for CFL = 0.95 with the coupled scheme). Using (29) the result obtained is similar. Note that theses instabilities do not appear when using the full coupled system even for CFL = 0.95 as it is shown in Fig. 16. This is due to the already known fact that one obtains always a negative eigenvalue even in supercritical regions. p Asffiffiffiffiffiwe see in ffi Fig. 17, in the supercritical region we obtain that u gh is always positive, while there exists a negative eigenvalue. As a consequence, the time-splitting method with upwind scheme does not take into account the information traveling backwards in general which explains the instabilities in the simulation shown before. 5. Conclusions Fig. p 17. ffiffiffiffiffiffi Eigenvalues for (41). True eigenvalues of the system (continuous line) and u gh approximation (dashed line).
The stronger the interaction is, the more likely to find oscillations. For instance, the Meyer-Peter and Müller formula models a weak interaction between fluid and sediment so that numerical instabilities due to splitting techniques are less frequent. Nevertheless, we may find that splitting technique fails. For example, consider the initial condition (40) but now impose q(x = 0, t) = 0.8. For the Meyer-Peter and Müller parameters we take sH cr ¼ 0:47, qs = 2612 and n = 0.095. As we see in Figs. 12 and 13, oscillations remain even for smaller CFL values. Remark that using too small CFL condition is not a good idea in general as this would increase the CPU time. For instance, using CFL = 0.125 we need 48.44 s of CPU time while the coupled scheme with CFL = 0.95 only needs 9.35 s. Remark that even using (29) we need to reduce CFL although in general this upper bound gives better results than (27). Now, let us consider an example with subcritical and supercritical regions. In order to do so, we begin with the initial condition
8 > < huðx; t ¼ 0Þ ¼ 0:6; 2 zb ðx; t ¼ 0Þ ¼ 0:1 þ 0:1eðx5Þ ; > : hðx; t ¼ 0Þ þ zb ðx; t ¼ 0Þ ¼ 0:4:
ð41Þ
We solve the shallow water system (Ag = 0.0) with this initial data until a steady state solution is reached (see Fig. 14) Once the steady
Exner equation coupled with shallow water equations results as a system that may lose hyperbolicity in some cases at least theoretically. A practical criterium has been introduced for the study of the hyperbolicity region. Nevertheless, the system remains hyperbolic in most of physical situations for classical definitions of the bedload transport flux. It could seem natural to do a splitting approach by solving first the shallow water system and then updating the topography. But even using a robust and well-balanced numerical scheme for shallow-water system the splitting technique may produce unphysical instabilities. In some cases the resulting instabilities may be avoided by reducing the CFL condition and using the upper bound (29) works better in general than (27). Nevertheless, one can observe that in certain cases these instabilities cannot be avoided with the splitting technique. Thus, in order to obtain a robust numerical scheme the full system must be considered and a coupled numerical scheme must be used. Acknowledgements This work has been partially supported by project METHODE ANR-07-BLAN-0232. One of the authors has been supported by the Spanish Government Research project MTM2009-11293. The authors also would like to acknowledge constructive comments of M. Castro, O. Cerdan, O. Delestre, C. Laguerre, C. Lucas, and P. Sochala.
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
References [1] Exner F. Über die wechselwirkung zwischen wasser und geschiebe in flüssen, Sitzungsber., Akad. Wissenschaften pt. IIa; 1925. Bd. 134. [2] Castro Díaz M, Fernández-Nieto E, Ferreiro AM. Sediment transport models in shallow water equations and numerical approach by high order finite volume methods. Comput Fluids 2008;37(3):299–316. doi:10.1016/j.compfluid.2007. 07.017. Available from:. [3] Grass A. Sediment transport by waves and currents, SERC London Cent. Mar. Technol Report No. FL29. [4] Hervouet J. Les équations de Navier–Stokes à surface libre et leurs formes simplifiées en eau peu profonde, Support de cours: Ecoulements peu profonds, INRIA; 2002. [5] Shields A. Application of similarity principles and turbulence research to bedload movement. Mitteilungen der Preußischen Versuchsanstalt für Wasserbau, Berlin; 1936. [6] Meyer-Peter E, Müller R. Formulas for bed-load transport. In: 2nd meeting IAHSR, Stockholm, Sweden; 1948. p. 1–26. [7] Luque RF, van Beek R. Erosion and transport of bedload sediment. J Hydraul Res 1976;14(2):127–44. [8] Nielsen P. Coastal bottom boundary layers and sediment transport. World Scientific Pub. Co. Inc.; 1992. [9] Ribberink J. Bedload transport for steady flows and unsteady oscillatory flows. Coastal Eng 1998;34:52–82. [10] Camenen B, Larson M. A general formula for non-cohesive bed load sediment transport. Estuar Coast Shelf Sci 2005;63:249–60. [11] Morales de Luna T, Castro Díaz MJ, Parés Madroñal C. A duality method for sediment transport based on a modified Meyer-Peter & Müller model, J Sci Comput; in press. doi:10.1007/s10915-010-9447-1. [12] Parés C. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J Numer Anal 2006;44(1):300–21. [13] Godlewski E, Raviart P. Numerical approximation of hyperbolic systems of conservation laws. Appl Math Sci 1996;118. [14] Bouchut F. Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources. Frontiers in mathematics. Birkhauser; 2004. [15] Lyn DA. Unsteady sediment-transport modeling. J Hydraul Eng 1987;113(1):1–15. doi:10.1061/(ASCE)0733-9429. Available from: . [16] Morris PH, Williams DJ. Relative celerities of mobile bed flows with finite solids concentrations. J Hydraul Eng 1996;122(6):311–5. Available from: . doi:10.1061(ASCE)0733-9429. [17] Lyn DA, Altinakar M. St. Venant–Exner equations for near-critical and transcritical flows. J Hydraul Eng 2002;128(6):579–87. Available from: . doi:10.1061(ASCE)0733-9429.
989
[18] Cordier S. Hyperbolicity of grad’s extension of hydrodynamic models of ionospheric plasma. part one: the single species case. Math Models Methods Appl Sci (M3AS) 1994;4(5):625–45. [19] Cordier S. Hyperbolicity of grad’s extension of hydrodynamic models of ionospheric plasma. part two: the two species case. Math Models Methods Appl Sci (M3AS) 1994;4(5):647–67. [20] Cordier S. Hyperbolicity of the hydrodynamical model of plasmas under the quasineutrality hypothesis. Math Models Appl Sci (M2AS) 1995;18:627–47. [21] Kubatko E, Westerink J, Dawson C. An unstructured grid morphodynamic model with a discontinuous Galerkin method for bed evolution. Ocean Model 2006;15:71–89. [22] Nord G, Esteves M. Psem_2d: a physically based model of erosion processes at the plot scale. Water Resour Res 41; 2005. 14 pp. doi:200510.1029/2004WR003690.. [23] Yang J. Assimilation de données variationnelle pour les problèmes de transport des sédiments en rivière, Ph.D. Thesis, Joseph-Fourier – Grenoble I University; 2004. [24] Toro EF. Shock-capturing methods for free-surface shallow flows. Wiley and Sons Ltd.; 2001. [25] Delestre O. Rain water overland flow on agricultural fields simulation. Ph.D. Thesis, Orleans University; 2010. . [29] Toumi I. A weak formulation of Roe’s approximate Riemann solver. J Comput Phys 1992;102(2):360–73. [30] Parés C, Castro M. On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems. M2AN Math Model Numer Anal 2004;38(5):821–52. [31] Weisstein EW. CRC concise encyclopedia of mathematics. 2nd ed. Chapman & Hall; 2002. [32] Morales de Luna T, Castro Díaz MJ, Parés Madroñal C, Fernández Nieto ED. On a shallow water model for the simulation of turbidity currents. Commun Comput Phys 2009;6(4):848–82. Available from: .
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help S. Cordier a,⇑, M.H. Le a,b, T. Morales de Luna c a
Laboratoire de Mathématiques - Analyse, Probabilités, Modélisation - Orléans (MAPMO), UMR CNRS 6628, Université d’Orléans , France Bureau de Recherches Géologiques et Minières (BRGM), France c Departamento de Matemáticas, Universidad de Córdoba, Spain b
a r t i c l e
i n f o
Article history: Received 23 December 2010 Received in revised form 8 May 2011 Accepted 8 May 2011 Available online 23 May 2011 Keywords: Shallow water system Exner equation Sediment transport Hyperbolicity Stability Splitting methods
a b s t r a c t In this paper, we are concerned with sediment transport models consisting of a shallow water system coupled with the so called Exner equation to describe the evolution of the topography. We show that, for some bedload transport models like the well-known Meyer-Peter and Müller model, the system is hyperbolic and, thus, linearly stable, only under some constraint on the velocity. In practical situations, this condition is hopefully fulfilled. Numerical approximations of such system are often based on a splitting method, solving first shallow water equation on a time step and, updating afterwards the topography. It is shown that this strategy can create spurious/unphysical oscillations which are related to the study of hyperbolicity. Using an upper bound of the largest eigenvector may improve the results although the instabilities cannot be always avoided, e.g. in supercritical regions. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction Soil erosion is caused by the movement of sediment due to mechanical actions of flow. In the context of bedload transport, a mass conservation law also called Exner equation [1] is used to update the bed elevation. This equation is often coupled with the shallow water equations describing the overland flows (see [2] and references therein). The complete system of PDE may be written in the form
8 > < @ t h þ @ x ðhuÞ ¼ 0; 2 2 @ t ðhuÞ þ @ x ðhu þ gh =2Þ ¼ gh@ x zb þ gh@ x H; > : @ t zb þ n@ x qb ¼ 0;
ð1Þ
where h is the water depth, u is the flow velocity and zb is the thickness of sediment layer which can be modified by the fluid. The sediment layer is assumed to be located on a fixed bedrock layer at depth H which is not modified by the fluid (see Fig. 1). n = (1 c)1 with c the porosity of the sediment layer. qb is the volumetric bedload sediment transport rate per unit time and width and g is the acceleration due to gravity. The conservative variable hu is also called water discharge and noted by q. In what follows and for the sake of simplicity, we shall assume that H is constant ⇑ Corresponding author. E-mail addresses: [email protected] (S. Cordier), [email protected] (M.H. Le), [email protected] (T. Morales de Luna). 0309-1708/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2011.05.002
so that the term @ xH can be neglected. We shall also assume that n = 1. Let us first quote some of the expressions proposed for the bedload transport rate in the literature. Many researches have developed different formulae for predicting and estimating bedload transport rate. Their expressions are usually proposed for granular non-cohesive sediments and quantified empirically. One of the simplest expression is the one proposed by Grass [3] which considers qb as a function of the flow velocity and a constant coefficient Ag which depends on soil properties
Grass ð1981Þ :
qb ¼ Ag ujujmg 1 ;
1 6 mg 6 4:
ð2Þ
In practice, estimations of bedload transport rate are mainly based on the bottom shear stress sb, i.e. the force of water acting on the bed during its routing. In the context of laminar flows, the bottom shear stress is given as
sb ¼ qghSf ;
ð3Þ
where q is the water density and Sf is the friction term that can be quantified by different empirical laws such as the Darcy–Weisbach or Manning formulae (see [4]) Darcy–Weisbach:
Sf ¼
fujuj ; 8gh
where f is the Darcy–Weisbach’s coefficient.
ð4Þ
981
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
0.07
Bedload rate (m^2/s)
0.06
Meyer-Peter & Muller Fernandez Luque & Van Beek Nielsen Ribberink Camenen & Larson
0.05 0.04 0.03 0.02 0.01 0
0
0.05
0.1
0.15
n2 ujuj h
4=3
ð5Þ
;
where n is the Manning’s coefficient. The bottom shear stress is usually used in dimensionless form, noted sH b , which is also called Shields parameter. It is defined in terms of the ratio between drag forces and the submerged weight by
sHb ¼
0.25
0.3
0.35
0.4
0.45
0.5
Fig. 2. Bedload rate as function of Shields parameter with ds = 32 mm, q = 1000 kg/ m3, qs = 2600 kg/m3 and sH cr ¼ 0:05.
Manning:
Sf ¼
0.2
Shields parameter
Fig. 1. Sketch of shallow water over an erodible bed.
sb ; ðqs qÞgds
ð6Þ
where qs is the sediment density and ds is the diameter of sediment. H The main hypothesis is that sH b must exceed a threshold value scr in order to initiate motion. The threshold value sH depends on the cr physical properties of sediment and is usually computed experimentally. One of the first works on this topic was done by Shields [5] in which sH cr is determined in relation with the boundary Reynolds number. Bedload transport rate may be represented as a function of sH b via a non-dimensional parameter U by
s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qs 3 qb ¼ U 1 gds :
ð7Þ
q
The followings expressions, illustrated by Fig. 2, have been often applied [6–10]:
€ ller ð1948Þ : U ¼ 8 Meyer-Peter & Mu
sHb sHcr
ndez Luque&Van Beek ð1976Þ : U ¼ 5:7 Ferna qffiffiffiffiffiffi H H Nielsen ð1992Þ : U ¼ 12 sH b sb scr þ ; H 1:65 Ribberink ð1998Þ : U ¼ 11 sH b scr þ ;
3=2
þ
ð8Þ
;
sHb sHcr
3=2 þ
or one could consider a modification of the classical formulae by including gravity effects as it was done in [11]. As we have said and for the sake of simplicity we shall assume here that H(x) is constant. Remark that the general case where H(x) is non-constant can be considered by introducing the trivially satisfied equation
@t H ¼ 0
ð13Þ
and following the ideas presented in [12]. PDE system (1) can be written in vectorial form as
@ t W þ @ x FðWÞ ¼ BðWÞ@ x W;
ð14Þ
where W = (h, q, zb)T is the state vector in conservative form, t 2 2 FðWÞ ¼ hu; hu þ 2g h ; qb ðh; huÞ and B(W) = ghE23 where Eij denotes the canonical base of the 3 3 matrices. System (1) can also be written in non conservative form as
@ t W þ AðWÞ@ x W ¼ 0;
ð15Þ
where A(W) = DWF B(W) is the matrix of transport coefficients. An important property of such systems is hyperbolicity [13] which requires that the matrix A(W) is R diagonalizable (or strictly hyperbolic when eigenvalues are distinct). Let us quickly give an interpretation of this property as a stability condition for the linearized system. Assume that W is a small wave perturbation of a constant state, i.e. of the form
Wðx; tÞ ¼ W 0 þ W 1 eiðkxxtÞ ;
1:
ð16Þ
This is a solution of the linearized problem
;
ð9Þ ð10Þ ð11Þ
1:5 H Camenen and Larson ð2005Þ : U ¼ 12 sH exp 4:5sH b cr =sb : ð12Þ Remark 1. Classical formulae for bedload transport do not take into account gravity effects. As a consequence, particles located at the advancing front of a dune do not fall due to gravity as the motion of particles only depends on the hydrodynamical variables. This means that vertical profiles may be observed in numerical simulations which are not found in physical situations, for example at the front of an advancing dune. In order to obtain more realistic profiles, a diffusion term may be added to the third equation of (1)
ixW 1 þ AðW 0 ÞðikÞW 1 ¼ 0
ð17Þ
if and only if W1 is an eigenvector associated to the eigenvalue xk . In other words, if one can write any arbitrary perturbation W1 as a linear combination of such eigenvector associated with real velocity, the solution will propagate without amplification of the perturbation. On the contrary, if any of the eigenvalues is complex, this will lead to instability. Let us remark that these eigenvalues are also important when using explicit upwind schemes to ensure stability. Indeed, the time step, for a constant mesh size Dx, has to satisfy the so called Courant–Friedichs–Lewy condition (CFL condition in what follows) (see [13,14] and references therein)
jkmax jDt < Dx;
ð18Þ
where jkmaxj is the maximum of the modulus of the eigenvalues or in other words, kmax is the maximum velocity for propagation of
982
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
information. We recall that when using finite volume schemes technique for (14), this condition can be seen as the definition of a time step sufficiently small so that the different Riemann problems at each intercell do not interact between each other, that is, information of each Riemann problem does not cross more than one cell. This paper is organized as follow: in next section, we study the hyperbolicity of system (1) and establish the main result. Next, we compare two different strategies for numerically solving (1): a time-splitting strategy and a coupled scheme. While some works have already been done comparing coupled and uncoupled schemes, we show how information on the eigenvalues may give better results. Finally, we present some numerical tests to justify the arguments mentioned.
Let us find relations between a and b in the different models proposed before: (2) and (8)–(12). Grass: qb = Agujujm1
u ¼ hq ) @qb @q
2
0
6 AðWÞ ¼ 4 gh u2 a
1
0
3
7 2u gh 5; b
ð19Þ
0
b b with a ¼ @q and b ¼ @q . The characteristic polynomial of A(W) can @h @q be written as
pA ðkÞ ¼ k
k 1 k 1 gh a b u þ gh 2u k
@u ; @q
@u @h @qb @h
)
¼ hq2
¼
ðqb Þ0u
)
@u @h
@qb q @qb ¼ : h @q @h
@qb @q
sb ¼ ðqb Þ0sb @@q ¼ ðqb Þ0sb ðsb Þ0u @u @q
@qb @h
sb ¼ ðqb Þ0sb @@h ¼ ðqb Þ0sb ðsb Þ0u @u @h
)
)
@qb q @qb ¼ h @q @h
Manning : 2 ujuj qjqj sb ¼ qgh nh4=3 ¼ a hujuj 1=3 ¼ a 7=3 h @ sb @q
jqj ¼ 2ha7=3
@ sb @h
¼ 76
9 =
q 2ajqj h h7=3
¼ 76
q @ sb h @q
;
)
@qb 7 q @qb ¼ : 6 h @q @h
Thus, we can summarize that we have two types of relation: a = ub or a ¼ 76 ub. We shall assume that b > 0 since sediment rate increases with that of the flow. The line d(k) can be rewritten as d(k) = ghb(k ku) with k = 1 or k = 7/6 and we shall assume the slope ghb > 0. Proposition 2.1. Consider system (1) with qb = qb(h, q) such that
@qb q @qb ¼ k : h @q @h
ð21Þ
For a given state (h, q), the system is strictly hyperbolic if and only if
a < ku < aþ ;
ð22Þ
where a± is defined by expression (26). More explicitly, In the case k = 1 (or a = ub), the system is always strictly hyperbolic. In the case k = 7/6 (or a ¼ 76 ub), a sufficient condition for system (1) to be strictly hyperbolic is
pffiffiffiffiffiffi juj < 6 gh:
2
¼ k½ðu kÞ2 gh þ ghðbk þ aÞ:
ðqb Þ0u
¼ 1h ;
Meyer-Peter and Müller, Fernández Luque and Van Beek, Nielsen, Ribberink, Camenen and Larson: qb = qb(sb) Darcy–Weisbach : sb ¼ qgh fujuj ¼ aujuj ða ¼ cstÞ 8gh
2. Domain of hyperbolicity We are concerned here by the hyperbolicity of system (1) for the different models proposed (2)–(11) for the bedload transport rate qb. It is known that, for Grass model, the system is always hyperbolic [2]. Nevertheless, to our knowledge, the study the hyperbolicity of system (1) for more general bedload transport fluxes has not been done. The usual approach in such cases relies b n on the typically small magnitude of @q which permits a perturba@u h tion solution of the cubic characteristic equation. This approach allows to give approximations of the eigenvalues of the system in that framework, see for instance [15–17]. Nevertheless this does not answer the fact of whether the model is always hyperbolic or not. In fact, Castro et al. [2] have stated numerically that hyperbolicity may be lost in some cases for the Meyer-Peter and Müller bedload transport flux. Moreover, it would be interesting to have an easier way to know if one would find complex eigenvalues or not rather than computing directly the roots of characteristic polynomial. Proposition 2.1 will give an answer to this fact and states a necessary and sufficient condition for system (1) to be hyperbolic for the different bedload transport fluxes proposed before. Note that qb = qb(h, q) and thus the matrix of transport coefficients is given by
¼
@u @q
ð23Þ
ð20Þ
System (1) is thus strictly hyperbolic if and only if pA(k) has three different solutions noted by k1 < k2 < k3. In other words if the curve f(k) = k[(u k)2 gh] and the line d(k) = gh(bk + a) have three points of intersection. This is illustrated in Fig. 3 for the case of a subcritical flow.
Proof. For simplicity, let us consider the case u > 0 so a = kub < 0. The case u < 0 can be treated with the same arguments. We define the two tangents of the curve f(k) which are parallel to d(k). Their intersections with f(k) are characterized by f0 (k) = ghb which yields to two values of k of the form def
k ¼
2u
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ 3ðgh þ ghbÞ : 3
ð24Þ
The two tangents are such that d±(k±) = f(k±). This implies the equations for the tangents are given by
d ðkÞ ¼ ghbðk a Þ;
ð25Þ
with def
a ¼ k Fig. 3. Eigenvalues of the transport coefficients matrix.
f ðk Þ : ghb
ð26Þ
The roots of pA(k) which correspond to the eigenvalues of A(W) are given as the intersection of f(k) and d(k) (see Fig. 3). Recall that
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
n pffiffiffiffiffiffio f(k) is a third order polynomial with roots 0; u gh . The equation pA(k) = d(k) f(k) will have 3 distinct solutions if and only if the line d(k) lies in between d-(k) and d+(k). This can be equivalently written as a- < ku < a+. pffiffiffiffiffiffi Itpcan ffiffiffiffiffiffi be checked that we always have a < u gh < u < u þ gh < aþ so the systempis hyperbolic =ffiffiffiffiffi 1.ffi ffiffiffiffiffialways ffi pffiffiffiffiffiffiin the case kp In the case k = 7/6, if juj < 6 gh, we have u gh < 76 u < u þ gh and thus the hyperbolicity condition is verified which concludes the proof. Note that similar arguments have been used by one of the author in [18–20] to characterize the hyperbolicity domains for systems of moments equations arising in plasma physics. h
Remark 2. Condition (22) is interesting because it gives a necessary and sufficient condition to check if the model remains hyperbolic without computing the eigenvalues. Moreover, condition (22) shows that the model is indeed hyperbolic in most physical situations. For instance, Fig. 4 shows the region where condition (22) is satisfied for the Meyer-Peter and Müller formula with values sHcr ¼ 0:47, qs = 2612 and n = 0.0196. 3. Time-splitting method The strategy of time-splitting consists in solving separately the shallow water system during a first step for a fixed topography zb, and updating afterward the topography on a second step using the third equation. See, for example, the references [21,22] or Chapter 3 of [23]. It seems natural such approach as the variation of topography is slow compared to the characteristic time associated to hydrodynamical variables. Using a time-splitting strategy allows to minimize computational costs and makes easier to implement the codes for different bedload transport fluxes in a modular way. Moreover, time-splitting techniques allows to use robust and efficient numerical schemes already developed for shallow water system like, for example, those presented in the publication [24]. In particular this would be an interesting approach for the code FullSWOF_2D (licence DL 03434-01) developed in the framework of the multidisciplinary project METHODE ANR-07-BLAN-0232 (see [25] and http://www.univ-orleans.fr/mapmo/methode/). However, it is a known fact that such an approach is not a good one as this may re-
sult in numerical instabilities, see for instance [26]. The main reason is related to the bad approximation of the true eigenvalues of the system. We shall propose in next section different numerical tests in order to study how reducing CFL condion or the new upper bound (29) may help. More explicitly: 1. As explained in the introduction, the stability condition for shallow water equations only requires the CFL limitation associated with
juj þ
pffiffiffiffiffiffi gh Dt < Dx:
ð27Þ
Remark that the CFL condition (27) is not enough to ensure the stability condition (18) for the coupled system (1) because the maximum eigenvalue k3 of the matrix pffiffiffiffiffiffi of transport coefficients A(W) is always larger than juj þ gh (see Fig. 3). This may lead to numerical instabilities in some numerical simulations when a splitting technique is used. These instabilities are less frequent or do not exist if the numerical scheme used is sufficiently diffusive. For instance, when Lax–Friedrichs scheme is used it is very difficult to find instabilities while for Roe type schemes they are easily found. We recall that something similar has been studied in [27,28] where the splitting technique for the two-layer shallow water system was studied. In [26] a comparison between coupled and decoupled schemes have also been done. Instead of computing the exact eigenvalues of A(W), an upper bound k03 of k3 may be used, defined asffi the intersection of d(k) pffiffiffiffiffi and the tangent dT(k) of f(k) at u þ gh (see Fig. 5)
pffiffiffiffiffiffi pffiffiffiffiffiffi dT ðkÞ ¼ f 0 u þ gh k u gh :
ð28Þ
k03 is thus of the form
def k03 ¼
uþ
pffiffiffiffiffiffi pffiffiffiffiffiffi gh f 0 u þ gh þ gha : pffiffiffiffiffiffi f 0 u þ gh ghb
ð29Þ
This value can be used to impose a CFL condition associated to the full coupled system (1). With this upper-bound the true CFL condition is now granted. Nevertheless we remark that now the true maximum propagation speed is overestimated. It will be interesting to study the splitting technique when this estimation is used. 2. Another problem is related to the sign of eigenvalues of the matrix A(W). It is known that there exist always a negative eigenvalue even in the supercritical case (see [15,16,26]). Indeed, rewriting the characteristic polynomial pA(k) in form
pA ðkÞ ¼ ðk k1 Þðk k2 Þðk k3 Þ; the product of three eigenvalues, in the case u > 0, satisfies
k1 k2 k3 ¼ pA ð0Þ ¼ gha < 0:
Fig. 4. Hyperbolicity region (in red) for Meyer-Peter and Müller. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
983
Fig. 5. Upper bound first order of the maximum eigenvalue.
984
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
pffiffiffiffiffiffi This means that in a supercritical flow, i.e. juj > gh, k2 and k3 are positive and consequently k1 < 0. This can be interpreted as a wave propagating upstream. The using of splitting strategy in such situation cannot account this information as we would have pffiffiffiffiffiffi take into pffiffiffiffiffi ffi 0 < u gh < u þ gh. This will lead to numerical instabilities when using splitting method with upwind schemes. 4. Numerical results
where
is assumed to be a Roe matrix for (1) in the particular case that
@qb q @qb ¼ k : h @q @h
4.1. Description of the numerical schemes Let us describe briefly the numerical schemes used hereafter. As usual, we consider a set of computing cells Ii ¼ ½xi1=2 ; xiþ1=2 ; i 2 Z. We shall assume that these cells have a constant size Dx and that xi+1/2 = iDx. xi = (i 1/2)Dx is the center of the cell Ii. Let Dt be the time step and tn = nDt. We denote by W ni the approximation of the cell averages of the exact solution
W ni ffi
1 Dx
Z
xiþ1=2
Wðx; tn Þdx:
ð30Þ
xi1=2
n
biþ1=2
pffiffiffiffi pffiffiffiffiffiffiffiffiffi g 1 hi þ hiþ1 mX Ag ¼ n pffiffiffiffi pffiffiffiffiffiffiffiffiffi ðuiþ1 Þk ðui Þmg 1k : hi hiþ1 þ hiþ1 hi k¼0
i Dt h þ n Ai1=2 W i W ni1 þ Aiþ1=2 W niþ1 W ni ; Dx ð31Þ
where
ð38Þ
For other models considered in this work, it is not always possible to obtain an explicit formula and its implementation is very costly. So, in practice the following approximation is used
biþ1=2 ¼
Splitting scheme: We use a Roe scheme for the shallow water system plus an upwinding technique based on the sign of velocity for the bedload transport equation, that is,
W inþ1=2 ¼ W ni
Remark 3. As it was stated in [2] it is possible to write a Roe matrix for to Grass model with (35) and
n
We consider two different schemes for (1):
ð37Þ
@qb ðhiþ1=2 ; hiþ1=2 uiþ1=2 Þ: @ðhuÞ
ð39Þ
As a consequence, we are using an approximated Roe matrix and we cannot use directly schemes in the form (36) but rather equivalent expressions based on the physical flux function F(W) as it is done in [2]. Remark 4. In order to compute the eigenstructure of matrix Ai+1/2, one could use the Cardano–Vieta’s formula (see [31]). Once the eigenvalues are known, one could express the eigenvectors as functions of these eigenvalues (see [32]). Nevertheless, numerical simulations in next section have been carried out by computing numerically the eigenstructure of Ai+1/2. 4.2. Numerical tests
is a Roe matrix for the shallow water system with non-flat topography, and then nþ1
hi
nþ1=2
¼ hi
;
nþ1=2
unþ1 ¼ ui i
ðzb Þnþ1 ¼ ðzb Þni i
;
ð32Þ
i Dt h nþ1 n ðqb Þnþ1 iþ1=2 ðqb Þi1=2 ; Dx
ð33Þ
with
ðqb Þnþ1 iþ1=2
8 nþ1 nþ1 > ; if uiþ1=2 P 0; < qb hi ; qi ¼ > : q hnþ1 ; qnþ1 ; if u iþ1=2 < 0 b iþ1 iþ1
In what follows and unless said otherwise, the computational domain is discretized using 400 cells. Boundary conditions have been set by defining a ghost cell duplicated from the first and last cell, respectively for each time step. We shall take the usual value mg = 3 for Grass model. First, we shall show a case where the splitting technique gives the same result as the fully coupled scheme. Let us consider as
ð34Þ
and ui+1/2 the velocity at the interface of the cells i and i + 1, for instance the value used in the Roe’s matrix. In particular, the usual definitions
hiþ1=2 ¼
hi þ hiþ1=2 ; 2
uiþ1=2 ¼
pffiffiffiffi pffiffiffiffiffiffiffiffiffi hi ui þ hiþ1 uiþ1 pffiffiffiffi pffiffiffiffiffiffiffiffiffi hi þ hiþ1
ð35Þ
is taken, which corresponds to the choice of segments as the path connecting two different states. See [27,29,30] for further details. Coupled Scheme: We use a Roe scheme for the full system which can be written in the form
W nþ1 ¼ W ni i
Dt þ n Ai1=2 W i W ni1 þ Aiþ1=2 ðW niþ1 W ni Þ ; Dx ð36Þ
Fig. 6. Initial condition for (40).
985
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
Fig. 7. Solution for (40) with Ag = 0.005.
Table 1 Comparison at time t = 10 between numerical solutions computed with a Roe scheme for the full system and a splitting scheme and a reference solution computed on a grid with 3200 points. # Points
50
100
200
400
Coupled L1 error L1 error L1 error
0.1056 0.1026 0.0385
0.0534 0.0503 0.0204
0.0264 0.0244 0.0107
0.0126 0.0114 0.0051
0.1021 0.1067 0.0330
0.0503 0.0534 0.0167
0.0240 0.0269 0.0084
0.0111 0.0135 0.0038
scheme in h in hu in z
Splitting L1 error in h L1 error in hu L1 error in z
initial condition a subcritical steady state for the classical shallow water system given by
8 huðx; t ¼ 0Þ ¼ 0:5; > < 2 zb ðx; t ¼ 0Þ ¼ 0:1 þ 0:1eðx5Þ ; > : u2 þ gðh þ zÞ ¼ 6:386: 2
ð40Þ
The water surface and topography corresponding to this initial condition are shown in Fig. 6. We shall consider Grass model (2) with Ag = 0.005. As we see in Fig. 7, no major differences are observed between the two schemes. Both simulations have been computed
Fig. 9. Solution for (40), Ag = 0.07 with a coupled scheme.
with CFL = 0.95 and in the case of the splitting scheme, (27) has been used for the CFL condition. Table 1 shows the comparison between the solutions computed with both schemes and a reference solution. We remark that CPU time for the splitting technique is
. Fig. 8. Solution for (40) with Ag = 0.07.
.
986
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
Fig. 10. Comparison of splitting scheme using (29) and coupled scheme for (40) with Ag = 0.07 and CFL = 0.5.
Fig. 11. Eigenvalues for (40). True eigenvalues of the system (continuous line) and u
pffiffiffiffiffiffi gh approximation (dashed line).
Fig. 12. An example using Meyer-Peter and Müller bedload transport flux with a splitting technique and (27).
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
987
Fig. 13. An example using Meyer-Peter and Müller bedload transport flux.
Fig. 14. Transcritical steady-state for shallow water system.
smaller (16.73 s) than the CPU time for coupled scheme (19.06 s) which is due to the fact that the spectral decomposition of the
matrices Ai+1/2 is computed numerically in the case of the coupled scheme. Now, consider the same initial condition (40) but set Ag = 0.07. We remark in Fig. 8 that some instabilities arise when using the splitting scheme with (27) which disappear if the CFL is reduced. This is not the case when the coupled scheme is used as we see in Fig. 9. Doing the same computation with the splitting scheme and using now (29) gives similar results but now CFL = 0.95 can be used as it is shown in Fig. 10. Remark that still some oscillations may be seen in the topography with the splitting scheme. The main difference between the experiment with Ag = 0.005 and Ag = 0.07 is that, in the second case, the true eigenvalues pffiffiffiffiffiffi of the p full ffiffiffiffiffiffi system are far from the approximations u þ gh and u gh given by the splitting scheme. In particular, the pffiffiffiffiffimaximum ffi of the modulus of the eigenvalue is larger than juj þ gh. This implies that a CFL condition close to 1 should not be used with the splitting scheme (see Fig. 11) This the reason why the upper bound (29) works in general better than (27). In previous examples, we have used Grass model. Nevertheless, conclusions are not exclusive to such model and may be observed if we use any other formula for the bedload transport flux. The main difference observed between the different models is the strength of the interaction between the fluid and sediment.
Fig. 15. Solution for (41) with Ag = 0.0005 using splitting scheme.
988
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
Fig. 16. Solution for (41) with Ag = 0.0005 using full system.
state is reached, we let the sediment to evolve by using a Grass bedload flux with Ag = 0.0005. Fig. 15 shows that again we find instabilities by using the splitting technique with (27). These instabilities remain even for CFL = 0.05 (which requires almost 20 times the CPU time required for CFL = 0.95 with the coupled scheme). Using (29) the result obtained is similar. Note that theses instabilities do not appear when using the full coupled system even for CFL = 0.95 as it is shown in Fig. 16. This is due to the already known fact that one obtains always a negative eigenvalue even in supercritical regions. p Asffiffiffiffiffiwe see in ffi Fig. 17, in the supercritical region we obtain that u gh is always positive, while there exists a negative eigenvalue. As a consequence, the time-splitting method with upwind scheme does not take into account the information traveling backwards in general which explains the instabilities in the simulation shown before. 5. Conclusions Fig. p 17. ffiffiffiffiffiffi Eigenvalues for (41). True eigenvalues of the system (continuous line) and u gh approximation (dashed line).
The stronger the interaction is, the more likely to find oscillations. For instance, the Meyer-Peter and Müller formula models a weak interaction between fluid and sediment so that numerical instabilities due to splitting techniques are less frequent. Nevertheless, we may find that splitting technique fails. For example, consider the initial condition (40) but now impose q(x = 0, t) = 0.8. For the Meyer-Peter and Müller parameters we take sH cr ¼ 0:47, qs = 2612 and n = 0.095. As we see in Figs. 12 and 13, oscillations remain even for smaller CFL values. Remark that using too small CFL condition is not a good idea in general as this would increase the CPU time. For instance, using CFL = 0.125 we need 48.44 s of CPU time while the coupled scheme with CFL = 0.95 only needs 9.35 s. Remark that even using (29) we need to reduce CFL although in general this upper bound gives better results than (27). Now, let us consider an example with subcritical and supercritical regions. In order to do so, we begin with the initial condition
8 > < huðx; t ¼ 0Þ ¼ 0:6; 2 zb ðx; t ¼ 0Þ ¼ 0:1 þ 0:1eðx5Þ ; > : hðx; t ¼ 0Þ þ zb ðx; t ¼ 0Þ ¼ 0:4:
ð41Þ
We solve the shallow water system (Ag = 0.0) with this initial data until a steady state solution is reached (see Fig. 14) Once the steady
Exner equation coupled with shallow water equations results as a system that may lose hyperbolicity in some cases at least theoretically. A practical criterium has been introduced for the study of the hyperbolicity region. Nevertheless, the system remains hyperbolic in most of physical situations for classical definitions of the bedload transport flux. It could seem natural to do a splitting approach by solving first the shallow water system and then updating the topography. But even using a robust and well-balanced numerical scheme for shallow-water system the splitting technique may produce unphysical instabilities. In some cases the resulting instabilities may be avoided by reducing the CFL condition and using the upper bound (29) works better in general than (27). Nevertheless, one can observe that in certain cases these instabilities cannot be avoided with the splitting technique. Thus, in order to obtain a robust numerical scheme the full system must be considered and a coupled numerical scheme must be used. Acknowledgements This work has been partially supported by project METHODE ANR-07-BLAN-0232. One of the authors has been supported by the Spanish Government Research project MTM2009-11293. The authors also would like to acknowledge constructive comments of M. Castro, O. Cerdan, O. Delestre, C. Laguerre, C. Lucas, and P. Sochala.
S. Cordier et al. / Advances in Water Resources 34 (2011) 980–989
References [1] Exner F. Über die wechselwirkung zwischen wasser und geschiebe in flüssen, Sitzungsber., Akad. Wissenschaften pt. IIa; 1925. Bd. 134. [2] Castro Díaz M, Fernández-Nieto E, Ferreiro AM. Sediment transport models in shallow water equations and numerical approach by high order finite volume methods. Comput Fluids 2008;37(3):299–316. doi:10.1016/j.compfluid.2007. 07.017. Available from:
989
[18] Cordier S. Hyperbolicity of grad’s extension of hydrodynamic models of ionospheric plasma. part one: the single species case. Math Models Methods Appl Sci (M3AS) 1994;4(5):625–45. [19] Cordier S. Hyperbolicity of grad’s extension of hydrodynamic models of ionospheric plasma. part two: the two species case. Math Models Methods Appl Sci (M3AS) 1994;4(5):647–67. [20] Cordier S. Hyperbolicity of the hydrodynamical model of plasmas under the quasineutrality hypothesis. Math Models Appl Sci (M2AS) 1995;18:627–47. [21] Kubatko E, Westerink J, Dawson C. An unstructured grid morphodynamic model with a discontinuous Galerkin method for bed evolution. Ocean Model 2006;15:71–89. [22] Nord G, Esteves M. Psem_2d: a physically based model of erosion processes at the plot scale. Water Resour Res 41; 2005. 14 pp. doi:200510.1029/2004WR003690.