A family of stable numerical solvers for the shallow water equations with source terms

A family of stable numerical solvers for the shallow water equations with source terms q Tom aas Chac oon Rebollo a,* , Antonio Dom ıınguez Delgado b , Enrique D. Fern aandez Nieto b a Departamento de Ecuaciones Diferenciales y Anaalisis Numeerico, Universidad de Sevilla, C/Tarfia, s/n. 41080 Sevilla, Spain b Departamento de Matem aatica Aplicada I, Universidad de Sevilla, ETS Arquitectura Avda, Reina Mercedes, N. 2, 41012 Sevilla, Spain Received 23 November 2001; received in revised form 3 June 2002 Abstract In this work we introduce a multiparametric family of stable and accurate numerical schemes for 1D shallow water equations.These schemes are based upon the splitting of the discretization of the source term into centered and decentered parts.These schemes are specifically designed to fulfill the enhanced consistency condition of Bermuudez and Vaazquez, necessary to obtain accurate solutions when source terms arise.Our general family of schemes contains as particular cases the extensions already known of Roe and Van Leer schemes, and as new contributions, extensions of Steger–Warming, Vijayasundaram, Lax–Friedrichs and Lax–Wendroff schemes with and without flux-limiters. We include some meaningful numerical tests, which show the good stability and consistency properties of several of the new methods proposed.We also include a linear stability analysis that sets natural sufficient conditions of stability for our general methods. Keywords: Finite volume method; Upwinding; Shallow water; Source terms 1. Introduction and motivation This paper deals with the numerical solution of 1D shallow water equations for channels with variable depth and width. These equations are a couple of conservation laws linking the depth hand the discharge q, which in condensated form read as follows: q This research was partially supported by Spanish Government Research Projects REN2000-1162-C02-01 and REN2000-1168-C0201. * Corresponding author. Tel.: +34-54-557989; fax: +34-54-552898. E-mail addresses: [email protected] (T. Chaco ´n Rebollo), [email protected] (A. Dom ıınguez Delgado), [email protected] (E.D. Ferna ´ndez Nieto). 0045-7825/03/$ - see front matter . PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 5 5 1 - 0 oW otþo oxFðWÞ¼G1ðx;WÞþG2ðx;WÞ;in 0;L½0;T½:ð1Þ Here W¼h q is the unknown, while Fis the flux function, FðWÞ¼ q q2 hþ1 2gh2 0 @1 A:ð2Þ Also, G1,G2are the source terms that respectively arise due to variable depth and width of the channel. These are defined by G1ðx;WÞ¼ 0 ghH0ðxÞ  ;G2ðx;WÞ¼ qb0ðxÞ bðxÞ q2 h b0ðxÞ bðxÞ 0 B B @1 C C A;ð3Þ where HðxÞis a function which describes the bottom of the channel with respect to a reference height, and bðxÞis a function yielding the width of the channel. Both Hand bare considered known. The numerical solution of shallow water equations with source terms faces the problem that low-accuracy solvers yield quite inaccurate solutions, exhibiting in particular large errors in the computation of wave speed (cf. [2]). This difficulty is overcomed if the numerical scheme solves some steady solution at least with order two. This is the Berm uudez–V aazquez enhanced consistency condition (cf. [12]). The extension of usual solvers for homogeneous conservation laws to shallow water equations with source terms has been reached for several flux-difference (cf. [1,4]) and flux-splitting schemes (cf. [5]). However, the techniques employed in these papers are rather specific for the solvers considered, and it is not clear how to extend a given scheme in order to satisfy the enhanced consistency condition. In this paper we address the question of determining systematic techniques to build solvers which verify the enhanced consistency condition. The main contribution is to show that this is reached if the numerical source term is split into a centered and a decentered part, similarly to the numerical flux, in such a way that the decentered part balances up to second order the decentered part of the flux. We may motivate this by considering a linear scalar conservation law with source term, ov otþaov ox¼gðxÞ;in 0;L½0;T½;ð4Þ with constant velocity a6¼ 0. We consider a smooth steady solution vðxÞ. As usual, given a space step Dx, we denote by vjan approximation to vðxjÞ,withxj¼jDx. We consider the simplest upwind scheme, corresponding to the discretization aov oxx¼xj ’ avjvj1 Dxif a>0; avjþ1vj Dxif a<0; 8 > < > : which we write as the sum of a centered part plus a decentered part, in the form 204 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 aov oxx¼xj ’avjþ1vj1 2Dx1 2jajDxvjþ12vjþvj1 ðDxÞ2:ð5Þ For smooth v, the consistency errors are given by avjþ1vj1 2Dx¼aov oxðxjÞþOðDx2Þ;1 2jajDxvjþ12vjþvj1 ðDxÞ2¼1 2jajDxo2v ox2ðxjÞþOðDx2Þ: If we consider a centered approximation of the source term, for example gjx¼xj’gðxjþ1=2Þþgðxj1=2Þ 2; then the consistency error is of second order, gðxjþ1=2Þþgðxj1=2Þ 2¼gðxjÞþOðDx2Þ: Thus, there is a first order error, stemming from the upwinding of the flux, which is not compensated. However, if we perform an upwind discretization of the source term, gjx¼xj’gðxj1=2Þif a>0; gðxjþ1=2Þif a<0;  we obtain––for instance, when a>0, gðxj1=2Þ¼aov oxðxj1=2Þ¼aov oxðxjÞ1 2aDxo2v ox2ðxjÞþOðDx2Þ: Then, the overall consistency error of our fully upwind scheme is OðDx2Þ, so the enhanced consistency condition is satisfied. Notice that for homogeneous equations we automatically recover a scheme of second-order consistency error for the steady equation. Indeed, in this case ov=ox¼0 and the numerical viscosity term in (5) is proportional to o2v=o2x. The key point is that centered discretizations automatically yield second-order accuracy, while a specific upwind of the source term is needed to compensate the first order error introduced by the upwinding of the flux. In the following sections we shall introduce and test a multiparametric family of numerical schemes satisfying the enhanced consistency condition, constructed upon the above observation. The general structure of this family of schemes will be a generalization of the fully upwind scheme introduced above. To do this, we write this scheme––for the evolution equation (4)––as vnþ1 j¼vn jDt/ðvn j;vn jþ1Þ/ðvn j1;vn jÞ Dx Gðxj1=2;xjþ1=2Þ;ð6Þ where the numerical flux /and the numerical source term Gare split into centered and decentered parts as /ðu;vÞ¼/Cðu;vÞþ/Dðu;vÞ;with /Cðu;vÞ¼auþv 2;/Dðu;vÞ¼1 2jajðvuÞ; Gðx;yÞ¼GCðx;yÞþGDðx;yÞ; ð7Þ with GCðx;yÞ¼gðxÞþgðyÞ 2;GDðx;yÞ¼1 2sgnðaÞðgðyÞgðxÞÞ:ð8Þ T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 205 The structure of the paper is the following. In Section 2 we introduce a general formulation of flux splitting and flux difference schemes for homogeneous shallow water equations, as a generalization of (7). In Section 3 we introduce our general upwind discretization for shallow water equations with source terms, as a generalization of (8). We respectively prove in Sections 4–6 the enhanced consistency condition for variable depth and width of the channel and for friction terms. We adapt in Section 7 a standard linear stability analysis, which proves the stability of our general schemes for homogeneous linear equations under reasonable CFL conditions. Finally, we present in Section 8 some numerical tests of several of the new methods introduced, by either comparison with analytical solutions (for small Froude numbers), experimental measurements or highly performing already known methods. 2. Numerical flux Our strategy will be to start from general upwind numerical schemes for homogeneous hyperbolic systems and to extend them to systems with source term for shallow water equations in a way they verify the enhanced consistency property. We write system (1) as oW otþo oxFðWÞ¼Gðx;WÞ; where the flux function Fcan be written as FðWÞ¼AðWÞW, AðWÞ¼ 01 q2 h2þ1 2gh 2q h! : We consider the following numerical scheme generalization of (6), Wnþ1 i¼Wn iDt/þðWn i;Wn iþ1Þ/ðWn i1;Wn iÞ Dx Gðxi1;xi;xiþ1;Wi1;Wi;Wiþ1Þ;ð9Þ where the functions /þand /are defined by: /ðWj;Wjþ1Þ¼FðW1;jþ1=2;W2;jþ1=2ÞþFðW3;jþ1=2;W4;jþ1=2Þ 21 2DðW 5;jþ1=2ÞWjþ1þ1 2DðW 6;jþ1=2ÞWj;ð10Þ and the function FðU;VÞis defined as: FðU;VÞ¼AðUÞV; W1;jþ1=2,W2;jþ1=2,W3;jþ1=2,W4;jþ1=2,W 5;jþ1=2,W 6;jþ1=2are convex combinations of Wjand Wjþ1, and DðWÞis a matrix function, yielding the upwinding of the numerical flux. The first and second summand in (10) are respectively the centered and decentered parts of the numerical fluxes /þand /. This formulation includes flux-splitting schemes, taking D¼jAj(cf. [9,10]), besides schemes of fluxdifference, with D¼jAj(cf. [11]), where AðWÞdenotes the Jacobian matrix of FðWÞ, AðWÞ¼ 01 q2 h2þgh 2q h! : In particular, the methods of Steger–Warming, Vijayasundaram, Roe and Van Leer for specific choices of the Wk;jþ1=2. For instance, for D¼jAjand W1;jþ1=2¼W2;jþ1=2¼W3;jþ1=2¼W4;jþ1=2¼W 5;jþ1=2¼W 6;jþ1=2¼1 2ðWjþWjþ1Þ 206 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 we obtain the Vijayasundaram method, while for the choice D¼jAjand W1;jþ1=2¼W2;jþ1=2¼Wj;W3;jþ1=2¼W4;jþ1=2¼Wjþ1; W 5;jþ1=2¼W 6;jþ1=2¼e WWjþ1=2, where e WWjþ1=2is the intermediate value of Roe, we obtain RoeÕs method. 3. Structure of numerical source In this section we construct the numerical source function associated to the general scheme (9) and (10). Firstly, in order that the scheme which is built be consistent and due to the fact that we start from the numerical flux functions that provide consistent schemes for homogeneous hyperbolic systems, we have to ask function Gto verify: Gðx;x;x;W;W;WÞ¼Gðx;WÞ: The construction of Gmust reflect an upwinding of the source term according to the construction of the numerical flux function. Following this idea we start from the following expression: G¼AA1G¼A 2 þD 2þA 2D 2A1G¼1 2GþDA1G  |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} ðGLÞ þ1 2GDA1G  |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} ðGRÞ :ð11Þ We stress that here Dis the matrix responsible for the upwinding in the expression of the numerical flux term (10). In the same way, we may use matrix Ainstead A. In this last case, matrices of re-scaling must be introduced, with the aim of achieving the condition of enhanced consistency: G¼PAA1P1G¼PA 2 þD 2þA 2D 2A1P1G ¼1 2GþPDA1P1G  |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðGLÞ þ1 2GPDA1P1G  |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðGRÞ ;ð12Þ where P¼c10 0c2  and D¼PDP1with P¼c 10 0c 2  :ð13Þ Notice that the effect of matrix Pis to re-schale the eigenvector of the diffusion matrix D, while keeping the same eigenvalues. The key point is that these eigenvalues contain the information about the upwinding of the scheme. This technique was introduced in [5]. The matrix Pallows to increase the number of parameter in the scheme in order the condition of enhanced consistency to be verified. We can write (11) and (12) under the compact expression: G¼^ PP ^ AA ^ AA1^ PP1G¼1 2Gþ^ PP ^ DD ^ AA1^ PP1G  |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðc GLGLÞ þ1 2G^ PP ^ DD ^ AA1^ PP1G  |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðc GRGRÞ ;ð14Þ where ^ DD ¼^ PPD^ PP1;with ^ PP ¼^ cc10 0^ cc2  ;^ PP¼^ cc 10 0^ cc 2 ! T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 207 and ^ cc i¼^ cci¼1if ^ AA ¼A; ^ cc i¼c i;^ cci¼ciif ^ AA ¼A: ( In the next sections, we employ the structure of d GLGL and d GRGR for the construction of the left and right numerical source functions, respectively, for each actual source term corresponding to variable depth and variable width. 4. Enhanced consistency: variable depth In this section we apply the general scheme introduced in Sections 2 and 3 to discretize the source term G1appearing in Eq. (1) which, we recall, corresponds to variable depth. We give simple sufficient conditions on the nodes Wk;iþ1=2and the matrices Pand Pthat guarantee the enhanced consistency condition. We are going to use the following notation: As Wk;jþ1=2,k¼1;...;4 are convex combinations of Wj,Wjþ1, we will denote as xk;jþ1=2and Hk;jþ1=2,k¼1;...;4 the same convex combinations of xj,xjþ1and HðxjÞ, Hðxjþ1Þrespectively: if Wk;jþ1=2¼aWjþð1aÞWjþ1for some a2½0;1; then xk;jþ1=2¼axjþð1aÞxjþ1and Hk;jþ1=2¼aHðxjÞþð1aÞHðxjþ1Þ: We define the following numerical source function corresponding to G1: G1ðxi1;xi;xiþ1;Wi1;Wi;Wiþ1Þ¼1 2ððG1Þi;LþðG1Þi;DÞþ1 2ððG1Þi;RðG1Þi;DþÞ;ð15Þ where ðG1Þi;L¼ 0 gh1;i1=2þhi 2 HiH2;i1=2 xix2;i1=2 0 @1 Axix2;i1=2 Dxþ 0 gh3;i1=2þhi 2 HiH4;i1=2 xix4;i1=2 0 @1 Axix4;i1=2 Dx; ð16Þ ðG1Þi;R¼ 0 gh1;iþ1=2þhi 2 H2;iþ1=2Hi x2;iþ1=2xi 0 @1 Ax2;iþ1=2xi Dxþ 0 gh3;iþ1=2þhi 2 H4;iþ1=2Hi x4;iþ1=2xi 0 @1 Ax4;iþ1=2xi Dx; ð17Þ ðG1Þi;D¼1 Dx ^ PP ^ DDðW 5;i1=2Þ^ AA1Wi1þWi 2  ^ PP1 0 ghi1þhi 2HðxiÞ !" ^ PP ^ DDðW 6;i1=2Þ^ AA1Wi1þWi 2  ^ PP1 0 ghi1þhi 2Hðxi1Þ !# ;ð18Þ ðG1Þi;Dþ¼1 Dx ^ PP ^ DDðWþ 5;iþ1=2Þ^ AA1WiþWiþ1 2  ^ PP1 0 ghiþhiþ1 2Hðxiþ1Þ !" ^ PP ^ DDðWþ 6;iþ1=2Þ^ AA1WiþWiþ1 2  ^ PP1 0 ghiþhiþ1 2HðxiÞ !# :ð19Þ 208 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 Here, P¼P¼Iif ^ AA ¼A;and P¼10 02  ;P¼20 01  if ^ AA ¼A: 8 < :ð20Þ These values are given in such a way that the property of exact consistency is verified for the stationary solution h q H 0  :ð21Þ We have: Theorem 1. The scheme defined by (10), (15), (16), (17), (18), (19) and (20) calculates in an exact way the stationary solution (21) under the conditions x2;iþ1=2xi Dxþx4;iþ1=2xi Dx¼1;ð22Þ W1;iþ1=2þW3;iþ1=2¼WiþWiþ1:ð23Þ Remark 1. Condition (22) is required to have a consistent approximation of the source term G1. Proof of Theorem 1. We have to check the following equality: /þðWi;Wiþ1Þ/ðWi1;WiÞ Dx¼Gðxi1;xi;xiþ1;Wi1;Wi;Wiþ1Þ; when W¼H 0  . We begin calculating the expression of the numerical flux /defined by (10) for this stationary solution. FðW1;iþ1=2;W2;iþ1=2Þ¼AðW1;iþ1=2ÞW2;iþ1=2¼0 1 2gH1;iþ1=2H2;iþ1=2  :ð24Þ For the decentered part of /we have: DðW 5;jþ1=2ÞWjþ1¼D11ðW 5;jþ1=2ÞHjþ1 D21ðW 5;jþ1=2ÞHjþ1 ! ;ð25Þ then, /ðWj;Wjþ1Þ¼ 1 2D11ðW 5;jþ1=2ÞHjþ1þ1 2D11ðW 6;jþ1=2ÞHj 1 4gH 1;jþ1=2H2;jþ1=2þH3;jþ1=2H4;jþ1=2   1 2D21ðW 5;jþ1=2ÞHjþ1þ1 2D21ðW 6;jþ1=2ÞHj 0 B B @1 C C A: In order to simplify the notation, we perform the calculations by components. Using (24) and (25) the flux difference is /þðWi;Wiþ1Þ/ðWi1;WiÞ Dx¼Ui;1D Ui;2CþUi;2D  ; T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 209 where Ui;1D¼1 2DxðD11ðW 5;i1=2ÞHiD11ðW 6;i1=2ÞHi1D11ðWþ 5;iþ1=2ÞHiþ1þD11ðWþ 6;iþ1=2ÞHiÞ;ð26Þ Ui;2C¼1 Dx 1 4gH1;i1=2H2;i1=2H3;i1=2H4;i1=2þH1;iþ1=2H2;iþ1=2þH3;iþ1=2H4;iþ1=2;ð27Þ Ui;2D¼1 2DxðD21ðW 5;i1=2ÞHiD21ðW 6;i1=2ÞHi1D21ðWþ 5;iþ1=2ÞHiþ1þD21ðWþ 6;iþ1=2ÞHiÞ:ð28Þ We observe that the centered part of the first component of the flux difference is null. To calculate the numerical flux function, we distinguish with analogous notation, its two components: Gðxi1;xi;xiþ1;Wi1;Wi;Wiþ1Þ¼ di;1D di;2Cþdi;2D  : Using (16) and (17), we obtain the following expression of the centered component di;2C: di;2C¼1 Dx 1 4gðH1;i1=2 þHiÞðHiH2;i1=2ÞþðH3;i1=2þHiÞðHiH4;i1=2Þ þðH1;iþ1=2þHiÞðH2;iþ1=2HiÞþðH3;iþ1=2þHiÞðH4;iþ1=2HiÞ:ð29Þ To calculate the decentered part of the numerical source in an independent way from the possible choice of ^ AA, we notice that both Aand Acan be written under the same expression: ^ AAðWÞ¼ 01 ^ ggH 0  if W¼H 0  ; where ^ gg ¼gif ^ AA ¼A; g 2if ^ AA ¼A: ( We now perform several partial calculations: ^ AA1WiþWiþ1 2  ^ PP10 gHi1þHi 2Hi ! ¼02 ^ ggðHiþHiþ1Þ 10 0 @1 A 1 ^ cc1 0 01 ^ cc2 0 B B @1 C C A 0 gHi1þHi 2Hi ! ¼ 1 ^ cc2 g ^ ggHi 0 0 @1 A: To obtain the expression of ^ PP ^ DD, for the sake of simplicity, we do not specify the point at which matrix Dis evaluated. ^ PP ^ DD ¼ ^ cc1D11 ^ cc1^ cc 1 ^ cc 2 D12 ^ cc2^ cc 2 ^ cc 1 D21 ^ cc2D22 0 B B @1 C C A: From these last two identities, we obtain ^ PP ^ DD ^ AA1WiþWiþ1 2  ^ PP10 gHi1þHi 2Hi ! ¼ ^ cc1D11 ^ cc1^ cc 1 ^ cc 2 D12 ^ cc2^ cc 2 ^ cc 1 D21 ^ cc2D22 0 B B @1 C C A 1 ^ cc2 g ^ gg Hi 0 0 @1 A¼ ^ cc1 ^ cc2 g ^ gg HiD11 ^ cc 2 ^ cc 1 g ^ gg HiD21 0 B B @1 C C A: 210 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 Now, using (18) and (19) we achieve the expressions of di;1Dand di;2D: di;1D¼1 Dx 1 2 ^ cc1 ^ cc2 g ^ gg ðD11ðW 5;i1=2ÞHiD11ðW 6;i1=2ÞHi1D11ðWþ 5;iþ1=2ÞHiþ1þD11ðWþ 6;iþ1=2ÞHiÞ;ð30Þ di;2D¼1 Dx 1 2 ^ cc 2 ^ cc 1 g ^ gg ðD21ðW 5;i1=2ÞHiD21ðW 6;i1=2ÞHi1D21ðWþ 5;iþ1=2ÞHiþ1þD21ðWþ 6;iþ1=2ÞHiÞ:ð31Þ Let us analyze when the equality between the flux difference and the numerical source is reached. We start with the first component, that is, we have to see when Ui;1D¼di;1D. Comparing expressions (26) and (30), it is observed these are the same when ^ cc1 ^ cc2 g ^ gg ¼1:ð32Þ If ^ AA ¼Athen, by (20) c1¼c2¼1 and ^ gg ¼g, what yields (32). When ^ AA ¼Athen ^ gg ¼g=2 and, according to (20), c1¼1 and c2¼2, what also implies (32). With regard to the second component, we separately match the centered and the decentered parts. To obtain Ui;2D¼di;2D, we compare (28) and (31) and observe that it is enough to have ^ cc 2 ^ cc 1 g ^ gg ¼1:ð33Þ Using the definition of Pin (20), this identity is again fulfilled for ^ AA ¼Aand ^ AA ¼A. Finally, operating over the expression (29) we get: di;2C¼Ui;2Cþ1 Dx 1 4gHiH2;iþ1=2 þH4;iþ1=2ðH1;iþ1=2þH3;iþ1=2ÞþH1;i1=2þH3;i1=2 ðH2;i1=2þH4;i1=2Þ: From hypothesis (22) we deduce W2;iþ1=2þW4;iþ1=2¼WiþWiþ1. Then, di;2CUi;2C¼1 Dx 1 4gHiHi þHiþ1ðH1;iþ1=2þH3;iþ1=2ÞþH1;i1=2þH3;i1=2ðHi1þHiÞ; where the right-hand side vanishes due to hypothesis (23).  Remark 2. To obtain the enhanced consistency property we have not needed to impose any condition upon the matrix D, neither on the points in which W 5;iþ1=2,W 6;iþ1=2is evaluated. This is a consequence of the structure centered part plus decentered part that we have given to our numerical source. Remark 3. For Roe and Van LeerÕs schemes there have already been constructed numerical source functions in such a way that the given extension calculates in an accurate way the stationary solution (21) (see [11]). The construction we have suggested of the numerical source function coincides with the existing one for these two methods. This occurs because instead of using an evaluation of ^ AA1in ðWiþWiþ1Þ=2, it can also be used another point of evaluation, whenever the value of ^ AA1on that point coincides with ^ AA1ðWiþWiþ1Þ=2ðÞfor the stationary solution. For example, this happens with ~ WWi;iþ1, intermediate value of Roe (see [6,7]). Remark 4. As we can see in expressions (32) and (33), the relevant values of the matrices of re-scaling are c1=c2and c 1=c 2. Remark 5. When some eigenvalues of matrix ^ AA becomes null, the expression of ðG1Þi;Ddoes not make sense as matrix ^ AA1does not exist. T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 211 Now we present some particular cases of election of the matrix D: 1. D¼jAjupwind (flux difference). 2. D¼jAjupwind (flux splitting). 3. D¼Dt DxA2Lax–Wendroff (flux difference). 4. D¼Dt DxA2Lax–Wendroff (flux splitting). 5. D¼uðrÞjAjþð1uðrÞÞ Dt DxA2where uðrÞis a flux limiter (flux difference). 6. D¼uðrÞjAjþð1uðrÞÞ Dt DxA2(flux splitting). 7. D¼Dx DtI(I––identity matrix) Lax–Friedrich. 8. D¼jnAþð1nÞAjwith n2R. Theorem 4. For the above particular choices of the matrix Dthe schemes are L2stable if Dt DxqðDÞ61 where by qðDÞwe denote the spectral radium of the matrix D. Under the same condition, the schemes corresponding to the cases 1,2,7 and 8 are moreover L1stable. In the case of the full non-lineal equation (shallow water), numerical experiments show that stability is only attempted if ðG2Þi;Land ðG2Þi;Rare convex combinations of approximations of G2in ðxi1;xiÞand ðxi;xiþ1Þrespectively, it has to be fulfilled that 0 6rða;bÞ61. This occurs when ða;bÞlies in the set A1[A2 defined by (see Fig. 1): A1¼ða;bÞ2½0;1=2 f½0;1=2such that if bP1=4 then aP2b1=2g: A2¼ða;bÞ2½1=2;1 f½1=2;1such that if b63=4 then a62b1=2g: ð47Þ In the definition of ðG2Þi;Land ðG2Þi;R, since rða;bÞ¼rð1a;1bÞ, we can merely take ða;bÞ2A1: The most usual case is a¼b¼0. For example, this occurs with Roe and Van Leer methods. Fig. 1. Sets A1and A2. 218 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 8. Numerical tests We present in this section the results of three tests. The first one compares one of our new methods with an analytical approximated solution, accurate for small Froude number. Secondly we compare another of our methods with experimental measurements of a dam breaking experiment. Finally, we compare a third new method with two stationary transcritical solutions provided by the Van Leer method and the Surface Gradient Method (cf. [13]). All of them yield excellent performances, practically the same as the known extensions of Roe and Van LeerÕs methods. Test 1 (Analytical test). In [12] it is reported a limit solution of shallow water equations with source terms, for small Froude numbers and ‘‘short’’ domains. We compare this analytical solution with the method corresponding to a¼b¼1=2 and D¼jAj, which we call ‘‘modified Vijayasundaram method’’. If we take the initial and boundary condition as hðx;0Þ¼HðxÞ;qðx;0Þ¼0; hð0;tÞ¼uðtÞþHð0Þ;qðL;tÞ¼wðtÞ;ð48Þ with uðtÞ¼4þ4 sin p4t 86400  1 2;wðtÞ¼0; then, the analytical approximated solution is hðx;tÞ¼uðtÞþHðtÞ;qðx;tÞ¼wðtÞþu0ðtÞ bðxÞZL x bðsÞds: We observe that the surface water level is horizontal at each fixed time. We have taken L¼1500 (in m), T¼10,800 (in s), Dx¼7:5 and a CFL condition equal to 0.8. The Manning coefficients for bottom and sidewalls taken are Mb¼Mw¼0:1. Also, we have taken the profile bottom and width functions proposed in [12], these are represented in Fig. 2. This corresponds to a Froude number near 0.1. We present in Fig. 3 the computed velocity, compared to the analytical solution, at time t¼10800. Also we have obtained that the free surface is effectively horizontal. Fig. 2. Test 1: depth function and width function. T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 219 Test 2 (Dam break experiment). This test consists in comparing the results of our solver with the measurements of a dam break experiment, performed in the Laboratoire de Recherches Hydrauliques of the Universit ee Libre de Bruxelles, under the direction of Prof. Hiver. These measurements have been used in BrufauÕs Ph.D. thesis [3] to test the extension of RoeÕs solver developed by V aazquez Cend oon in [11]. The experiment is as follows. A reservoir of 15.5 m length is filled with water up to a height of 0.75 m. A floodgate separates the reservoir from a straight channel of 22.5 m length, with a triangular obstacle of 0.4 m height and 6 m length. The width is constant 1.75 m. Two kinds of outflow boundary are considered: free exit (Test 2.1) and large vertical wall (Test 2.2, see Fig. 5). The experimental data available are the free surface position along the time interval ½0;40(s), at 20 points, represented in Fig. 4. The Manning coefficients for bottom and sidewalls are Mb¼0:0125 and Mw¼0:011. The interest of this test is that it provides an analysis of accuracy of not only our numerical schemes, but also of the full modelling process. We recall that shallow water equations are based upon the hydrostatic pressure assumption and model the energy dissipation effects through ManningÕs law. This may lead to inaccuracies in the numerical results which are not due to the numerical model. However, this test is particularly well suited to test the accuracy in the computation of speeds of waves and shocks. In our experiments we have taken 152 points along the channel, corresponding to Dx¼25 cm. We have adapted a second-order four-step explicit Runge–Kutta scheme for the time discretization (cf. [5]). This provides a solver with an overall second-order accuracy. Thus, we may expect that the discrepancies between numerical results and experimental measurements are mainly due to the continuous model, rather than to the numerical solver. The increase of computational complexity involved for using this secondFig. 3. Test 1: computed versus analytical velocity at t¼10,800. Fig. 4. Test 2: measurement points. 220 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 order solver is not high, as it provides stable solutions with CFL numbers up to 2.5, while the explicit Euler scheme is stable for CFL numbers smaller than 0.8. We represent in Figs. 6 and 7 the time evolution of the computed free surface at points G10, G11, G13 and G20 for sub-Tests 2.1 and 2.2. These points likely are the more meaningful, as the first three are situated along the obstacle (in particular, the top of the obstacle is situated at G13), and the fourth one is between the obstacle and the outflow boundary. Test 2.1 (Free outflow condition). In this test we have used the method corresponding to a¼b¼1=4 and D¼jAj. We may observe a good accuracy at all points considered. The speed of propagation of discontinuities is well computed. Even the overall pattern of the free surface evolution at the top of the obstacle (point G13) is well reproduced. In general, depths are overestimated. We think that this is a consequence of the lack of energy dissipation mechanisms in our model. Test 2.2 (Large vertical wall at outflow boundary). We have used the method corresponding to a¼b¼1=2 and matrix Dcorresponding to case 8 presented in Section 7. In this case, the presence of the wall at the outflow boundary produces several reflections which travel back and forth along the channel and are in turn reflected by the obstacle. The overall pattern of the time behaviour of the free surface is again correctly represented by the numerical solution. The four fronts passing along the point G13 are recovered, so as their speed of propagation. The zero depth in the time interval ½29;33, is quite accurately reproduced. In point G20 these characteristics also are correctly simulated, with a quite satisfactory accuracy. Fig. 5. Test 2: initial condition and outflow boundary, Tests 2.1 and 2.2 respectively. Fig. 6. Test 2.1: free surface evolution at measurement points. Solid line: numerical results (a¼b¼1=4). Dotted line: experimental measurements. T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 221 However, for the measurement points situated G10 and G13, there are relatively large deviations with respect to the experimental measurements. This probably occurs because the modelling of friction effects by ManningÕs law that we include in our model does not produce an energy damping enough to correctly reproduce the turbulence effects in this flow. Very likely, this produces high levels of turbulence in this flow, especially close to these points. Moreover, the point G10 is situated in a corner where, as we mentioned, possibly the boundary layer detaches and a recirculation zone takes place. This could explain the larger discrepancies between experimental measurements and numerical results at point G10, with respect to the other points considered. Test 3 (Stationary solutions). In this test we calculate two stationary flows over a bump defined by zbðxÞ¼ 0:20:05ðx10Þ2if 8 <x<12; 0 otherwise:  The width is constant equal to 1 and the channel length is L¼25 m. We have taken Dx¼0:25 m and a CFL condition equal to 0.8. As initial condition we have taken h¼0:5zband q¼0. We compare the solution provided by the method corresponding to a¼b¼1=8 and D¼jAj, versus h and qprovided by the method proposed in [13] (the Surface Gradient Method), the analytical solution (see [8]) and the method proposed in [11] (RoeÕs method). We look for stationary flows. We consider that the scheme converges to a stationary solution if the relative error Rproposed in [13], R¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X i hn ihn1 i hn i  2 v u u t verifies R<5106. Transcritical flow without shock: This case is obtained by imposing h¼0:66 m downstream when the flow is subcritical and q¼1:53 m2/s upstream. In Figs. 8 and 9 we represent the height of the flow and discharge respectively. Fig. 7. Test 2.2: free surface evolution at measurement points. Solid line: numerical results (a¼b¼1=2). Dotted line: experimental measurements. 222 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 The method corresponding to a¼b¼1=8 has calculated a stationary solution after 1645 iterations, RoeÕs method after 1611 iterations and the surface gradient method after 1490 iterations. We observe a non-physical decrease of discharge near the obstacle reproduced by RoeÕs and SGM method. Method with a¼b¼1=8 does not suffer this effect. Transcritical flow with shock: This case is obtained by imposing h¼0:33 m downstream and q¼1:53 m2/s upstream. In Figs. 10 and 11 we represent the height of the flow and discharge respectively. The method corresponding to a¼b¼1=8 has calculated a stationary solution after 3356 iterations, RoeÕs method after 3689 iterations and the surface gradient method after 3748 iterations. Fig. 8. Test 3, h: transcritical flow without shock. Fig. 9. Test 3, q: transcritical flow without shock. T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 223 We observe in all methods some errors in discharge due to the interaction obstacle-shock. However, the errors due to method with a¼b¼1=8 are remarkably smaller than those due to RoeÕs and SGM methods. 9. Conclusion In this paper we have introduced a family of solvers for 1D shallow water equations with source terms, which satisfies a strengthened consistency condition for stationary solutions. This family includes as parFig. 10. Test 3, h: transcritical flow with shock. Fig. 11. Test 3, q: transcritical flow with shock. 224 T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 ticular cases all known extensions of usual solvers that satisfy this condition. Our main methodological innovation is the decentering of source terms which allows to compute the numerical viscosity of the flux up to second order. Several new methods are tested which present performances quite close to the best already known methods, which even improve in the case of transcritical flow on an obstacle. Acknowledgements The authors wish to thank Professors Pilar Brufau and Pilar Garcia Navarro for their technical help, and also Professors Manuel Castro D ııaz, Macarena G oomez M aarmol, Carlos Par ees Madro~ nnal and Maria Elena V aazquez Cend oon for their valuable remarks and in general their interest in the development of this work. References [1] A. Bermudez, A. Dervieux, J.A. Desideri, M.E. V aazquez Cend oon, Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes, Comput. Meth. Appl. Mech. Engrg. 155 (49) (1998). [2] A. Berm uudez, M.E. V aazquez Cend oon, Upwind methods for hyperbolic conservation laws with source terms, Comput. Fluids 23 (8) (1994) 1049–1071. [3] P. Brufau. Simulaci oon bidimensional de flujos hidrodin aamicos transitorios en gemotr ııas irregulares, Ph.D. thesis, Universidad de Zaragoza, 2000. [4] J. Burguete, P. Garc ııa-Navarro, Efficient construction of high-resolution TVD conservative schemes for equations with source terms, Application to shallow water flows, Int. J. Numer. Meth. Fluids 37 (2001) 209–248. [5] T. Chac oon Rebollo, E.D. Fer nnandez Nieto, M. G oomez M aarmol, A flux-splitting solver for shallow water equations with source terms, Int. J. Numer. Meth. Fluids, in press. [6] E. Godlewski, P.A. Raviart, Hyperbolic systems of conservation laws, Mathematiques et Applications, Ellipses, Paris, 1991. [7] E. Godlewski, P.A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer-Verlag, 1996. [8] N. Goutal, F. Maurel, in: Proceedings of the 2nd Workshop on Dam-Break Wave Simulation, Technical Report HE-43/97/016/A, Electricit ee de France, D eepartement Laboratoire National dÕHydraulique, Groude Hydraulique Fluviale, 1997. [9] LeVeque, H.C. Yee, A study of numerical methods for hyperbolic conservation laws with stiff source terms, J. Comput. Phys. 86 (1990) 187–210. [10] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, 1997. [11] M.E. V aazquez Cendon, Estudio de esquemas descentrados para su aplicacion a las leyes de conservaci oon hiperb oolicas con t eerminos fuente, Ph.D. thesis, Universidad de Santiago de Compostela, 1994. [12] M.E. V aazquez Cend oon, Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry, J. Comput. Phys. 148 (1999) 497–526. [13] J.G. Zhou, D.M. Causon, C.G. Mingham, D.M. Ingram, The surface gradient method for the treatment of source terms in the shallow-water equations, J. Comput. Phys. 168 (1–25) (2001). T. Chac oon Rebollo et al. / Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225 225