A consistent intermediate wave speed for a well-balanced HLLC solver.

Numerical Analysis /Analyse Num´erique. A consistent intermediate wave speed for a well-balanced HLLC solver. E.D. Fern´ andez-Nietoa, D. Breschb, J. Monnierc. aDpto Matem´atica Aplicada I, Universidad de Sevilla, Avda. Reina Mercedes N. 2, 41012 Sevilla (Spain). bLAMA, UMR5127 CNRS, Univ. Savoie, 73376 Le Bourget du Lac (France). cLJK - MOISE project-team, INP-Grenoble & INRIA, BP 53, F-38041 Grenoble (France) e-mails : [email protected],[email protected],[email protected] Abstract. In this work, we present a HLLC scheme modification for application to nonhomogeneous shallow-water equations with pollutant transport. This new version is related to the definition of a consistent approximation of the intermediate wave speed. Numerical results are presented to illustrate the importance of such approximation to get appropriate pollutant concentration profiles. Une vitesse interm´ediaire consistante pour un sch´ema HLLC bien ´equilibr´e. R´esum´e - Nous pr´esentons une modification du sch´ema HLLC pour application aux ´equations de Saint-Venant non homog`enes avec transport de polluants. Cette nouvelle version du sch´ema est reli´ee `a la d´efinition d’une vitesse d’onde interm´ediaire consistante. Un exemple est donn´e afin d’illuster comment une mauvaise approximation de la vitesse interm´ediaire peut entraˆıner de mauvais r´esultats sur l’approximation de la concentration en polluant mˆeme avec un sch´ema bien ´equilibr´e. Version fran¸caise abr´eg´ee - Nous proposons une extension du sch´ema HLLC pour les ´equations de Saint-Venant avec topographie et propagation polluants (voir [6]). Pour ce syst`eme le sch´ema HLLC peut s’´ecrire en fonction du flux num´erique provenant du sch´ema HLL. Pour le cas homog`ene, ce sch´ema est bas´e sur le fait que la troisi`eme composante du flux peut s’´ecrire en fonction de la premi`ere composante du flux et de la concentration en polluant. La d´efinition de la troisi`eme composante du flux est [Fhllc i+1/2]3= [Fhll i+1/2]1ϕ∗ La d´efinition de ϕ∗est simplement la concentration du polluant `a gauche ou `a droite de l’intervalle x=xi+1/2en fonction du signe de la vitesse interm´ediaire S∗. Pour le cas non homog`ene nous consid´erons un sch´ema de type HLL ´ecrit comme fonction du flux num´erique qui d´epend de la d´efinition du terme de front (voir [4], [8]). Ceci va nous permettre d’obtenir une extension naturelle du sch´ema HLLC au cas non homog`ene `a partir du cas homog`ene. Nous prouvons que le sch´ema propos´e calcule exactement la solution stationnaire de l’eau au repos et qu’il est asymptotiquement bien ´equilibr´e (voir [4]) ind´ependamment de la d´efinition de S∗. Nous proposons alors une modification de S∗, qui d´epend de la topographie, qui vaut z´ero au repos et qui ´evite des pics anormaux de polluants. En guise d’exemple, nous exhibons un test num´erique sur lequel nous observons l’importance de cette modification. L’expression de S∗est donn´ee par (5). 1 Introduction In this paper we study the extension of the HLLC solver to the Shallow Water Equations (SWE) with topography term, either in the 1D case with pollutant or equivalently in the 2D case. The HLLC solver is an improvement of HLL solvers. HLL solver is low cost-computing since no characteristic decomposition is required. But, HLL solver considers only the maximum and the minimum eigenvalues of the system, hence the scheme presents a large diffusion related to the 1 intermediate field. For 1D SWE with pollutant or 2D SWE, we have three characteristics fields. At contrary, HLLC solver takes into account the intermediate velocity, see e.g. [15]. In order to extend the HLLC solver to the non-homogeneous case, we study its well-balance properties. Also, HLLC solver needs to predict the intermediate wave speed at intercells, thus we propose a consistent definition. The well-balanced properties are related to the stationary solutions of the system. In our case, we seek numerical schemes that preserves at least the solution of water at rest. This is the so-called C-property introduced by Berm´udez and V´azquez in [2]. In [9], Greenberg and Leroux introduce the concept of well-balanced numerical schemes : they define a numerical scheme that balances the different terms for a non-homogeneous hyperbolic equation. The well-balanced property for hyperbolic systems with source terms has been studied for kinetic schemes (see Perthame et al., [13]), relaxation solvers (see Bouchut, [3]) and a family of Q-schemes (see Chac´on et al., [5]). An other approach consist to rewrite the hyperbolic system with source term as a purely nonconservative term, by including a new equation. In that approach, because of the presence of nonconservative products, a new definition of weak solution must be introduced, depending on the choice of a family paths, see DalMaso et al. [7], Castro et al. [12], Gosse et al [6], LeFloch-Tzavaras [10]. Also the well-balanced properties of a numerical scheme can be related to its convergence. In [14], the authors prove an extension of Lax-Wendroff theorem which is related to the well-balance properties. In [4], the authors introduce the concept of asymptotically well-balanced. A numerical scheme verifies this property if it preserves all regular stationary solutions of the non-homogeneous hyperbolic system over a set whose the measure of the complementary tends to zero when the space step ∆xtends to zero. We will consider this property in the present paper. In the present paper, we present a natural extension of the HLLC solver to the non-homogeneous SWE with pollutant. We show that, independently of the definition of the intermediate wave speed S∗, the scheme is asymptotically well-balanced. Anyway, we show that we need to define a consistent S∗depending on the source term that means an expression which vanish when water is at rest. To illustrate such phenomena, we present a numerical test which compare the same well-balanced scheme but with two different definitions of S∗, consistent or not. This paper is organized as follows. In Section 2, we present the 1D SWE equations we consider. In Section 3, we present the well-knowned HLL and HLLC solvers in the homogeneous case. In Section 4, we propose the extension of HLLC solvers, based on given HLL solvers. We show some properties of the proposed HLLC solver and propose a consistent definition of the intermediate wave speed at interfaces S∗. In Section 5, we present a numerical test which show the good results obtained when a well-balanced HLLC solver with a consistent definition of S∗are considered and bad results if S∗is not consistent. 2 The equations We denote the vectors of unknows by W:R×R+→Rm, and by F:Rm→Rmthe flux function. The conservation law we consider is the following : ∂tW(x, t) + ∂xF(W(x, t)) = G(x, W ), x ∈[0, L], t ∈[0, T ].(1) with : W= (h, q, r)T,h=h(x, t) is the height of the water column at instant tand position x, q=hu is the discharge, uis the averaged horizontal velocity, ris the third unknown, F(W) = (q, q2 h+1 2gh2,q r h)T,gis the gravity constant and G(x, W ) = (0,−ghz0 b(x),0)Tis the source term, where zbdenotes the topography. In the case where we consider the 1D SWE with pollutant, r=hϕ with ϕthe pollutant concentration. In the case where we consider the 2D SWE, we have : ∂tW+ div F(W) = 0, with W= (h, ~q)T,hthe water depth, ~q = (q1, q2) the discharge, ~q =h~u,~u the depth-averaged velocity, 2 F(W) = G(W), H(W)Tis the flux vector. Then, if we use the rotational invariance property of the equations, see [15], then the 2D SWE are reduced to a sum of 1D local Riemann problems, which are the same equations than the previous 1D SWE with pollutant. 3 The standard HLL and HLLC schemes for homogeneous case Using a conservative finite volume method, the general structure is the following. We define the points xi=i∆xas a partition of [0, L] with constant step ∆x. For a given time t=tn, we set tn+1 =tn+ ∆t. The time step ∆tis such that stability condition is satisfied. Then, if we denote Wn i≈1 ∆xRxi+1/2 xi−1/2W(x, tn)dx, a conservative three point finite volume scheme is : Wn+1 i−Wn i ∆t+Fn i+1/2−Fn i−1/2 ∆x= 0. where Fn i+1/2=Fi+1/2(Wn i, Wn i+1) is the numerical flux function. The eigenvalues of the jacobian matrix of the flux are : λ1= (u−c), λ2=u,λ3= (u+c) with c=√g h. The HLL solver is defined by : Fhll i+1/2=SRF(Wi)−SLF(Wi+1) + SLSR(Wi+1 −Wi) SR−SL where SLand SRare approximation of lower and upper bounds respectively, of the smallest and largest local speeds. We set : SL= min{λ1(Wi),0}and SR= max{λ3(Wi+1),0}, but other definitions are possible, see [15]. The HLLC solver takes into account the intermediate eigenvalue λ2=u. A possible version of HLLC solvers (for the homogeneous SWE) is such that the flux function relates the flux corresponding to the passive scalar to the flux of the mass conservation equation, see [15]. Namely, we set : [F]3= [F1]·ϕ(with ϕ=r/h). So, the first and second component of the flux are approximated by a given HLL solver. And the third component is approximated using the first component of the numerical flux of HLL solver as follows : [Fhllc i+1/2]3= [Fhll i+1/2]1·ϕ∗ where ϕ∗=ri/hiif S∗≥0 and ϕ∗=ri+1/hi+1 if S∗<0. The value of S∗must approximate the intermediate wave speed at intercell xi+1/2. In [15], different definitions of S∗can be found. The following definition is a common choice : S∗=SLhi+1(ui+1 −SR)−SRhi(ui−SL) hi+1(ui+1 −SR)−hi(ui−SL)(2) 4 Extension to the non-homogeneous case In this section we present the extension of the HLLC solvers to the non-homogeneous case. It is based on two points : i) the structure of the scheme in conservative form ; ii) a consistent definition of S∗. In Subsection 4.1, we present the extension of the HLLC solver and a consistent definition of S∗. In Subsection 4.2, we summarize some possible choices for HLL solvers. In Subsection 4.3, we show the well-balanced properties of the proposed HLLC scheme. 4.1 HLLC scheme with source term If we consider numerical schemes in conservative form, the structure is obtained by integrating (1) on the control volume (xi−1/2, xi+1/2). This gives : Wn+1 i−Wn i ∆t+Fn G,i+1/2−Fn G,i−1/2 ∆x=Gn i,(3) 3 where Gn iis a centered approximation of Gat x=xi:Gi= (Gi−1/2+Gi+1/2)/2. If we consider a Godunov solver, Fn G,i+1/2is the approximation of F(W(xi+1/2, tn)) where Wis the solution of the non-homogeneous Riemann problem. Let us point out that we do not follow the standard approach which consists to consider the homogeneous flux then upwind the source term, see [2]. Here, Fn G,i+1/2depends on the source term G. Moreover, we notice that under structure (3), the extension of HLLC solver is simple and natural. If we consider a given HLL numerical flux, we can apply the same idea as for the homogeneous system : let us denote by Fhll Gthe HLL numerical flux in the presence of source term G. Then, the first and second components of Fhllc Gare equal to Fhll G, and we set : [Fhllc G]3= [Fhll G]1·ϕ∗(4) with ϕ∗=ri/hiif S∗>0 and ϕ∗=ri+1/hi+1 if S∗<0. Concerning the definition of S∗, we can consider (2). Nevertheless, we can easily prove that for water at rest solution, this definition leads to S∗6= 0. Thus, we propose the following definition : S∗=SLqi+1 −SRqi−SLSR(hi+1 −hi−∆x[g A−1(Wi+1/2)Gi+1/2]1) hi+1(ui+1 −SR)−hi(ui−SL).(5) where Gi+1/2is approximation of Gat x=xi+1/2. Concretely we set [Gi+1/2]l= 0, for l=1,l= 3 and [Gi+1/2]2=−ghi+hi+1 2 zb(xi+1)−zb(xi) ∆x. By e A−1(Wi+1/2) we denote an approximation of the inverse matrix of Aat xi+1/2, and if an eigenvalue of Avanish then the corresponding eigenvalue of e A−1is setted to zero. 4.2 Well-balanced HLL solvers In this subsection, we present briefly the existing HLL scheme extended to the non-homogeneous case which respect the well-balance property. We rewrite them under the structure (3)-(4) we are interested in. We can classify these HLL solvers for SWE with topography term in two types : i) well-balanced schemes for the water at rest solution only ; ii) well-balanced schemes for all stationary solutions. i) The idea followed by many studies consists to use a scheme for the homogeneous case but evaluated at different states taking into account the topography term. That is, Fhll G,i+1/2=Fhll i+1/2(W− i+1/2, W+ i+1/2),verifying (W+ i+1/2−W− i+1/2) = 0 for water at rest. For example, Zhou et al. propose in [16] a linear approximation satisfying the previous property. LeVeque in [11] proposes a technique that build the states preserving a desired stationary solution. In the case of water at rest, LeVeque proposes to replace hby the water surface. Also, the hydrostatic reconstruction proposed by Audusse et al. in [1] can be applied. In this case, the water column is evaluated by adding and substituing two evaluations of the topography. Another technique consists to rewrite the SWE in function of the water surface and not in function of the water. ii) A well-balanced HLL solver for all stationary solutions is presented in [8]. It is defined as follows : Fhll G,i+1/2=F(Wi) + F(Wi+1) 2−1 2(SR+SL SR−SL (F(Wi+1)−F(Wi)−∆xG)+ +SRSL SR−SL (Wi+1 −Wi−∆xg A−1(Wi+1/2)Gi+1/2) (6) This scheme can be applied to any hyperbolic system with source term, hence to SWE with topography term. It is asymptotically well-balanced in the sense of [4]. That is, this scheme preserves any steady state up to second order. 4 4.3 Well-balanced and consistent HLLC solvers We can use any of the previous versions of HLL solvers. Then, following the structure (3)-(4), we obtain the corresponding HLLC solver. Let us study their well-balanced properties. We have : Theorem 4.1 If the HLL solver is well-balanced, then the HLLC solver defined by (4) is wellbalanced too, independently of the definition of S∗. The proof is straightforward since all stationary solutions verify q= 0 or ϕconstant. Corollary 4.2 a) Any HLLC solver builded from a HLL solver of type i) (see subsection 4.2) preserves the water at rest solution. b)Any HLLC solvers builded from HLL solver of type ii) is well-balanced for any stationary solutions (in the sense of [4]). This remains true independently of the definition of S∗. As we wrote above, the well-balanced properties are independent of the definition of S∗. Nevertheless, we proposed in (5) a modification of the classical definition of S∗in order to take into account the topography term. Then, we have the following result : Theorem 4.3 The intermediate wave speed S∗defined by (5) is an approximation of uat third order in space for any stationary solution. Furthermore, S∗vanishes when water is at rest. The proof is based on the fact that any stationary solution of SWE verifies qconstant. Moreover (hi+1 −hi−∆x[g A−1Gi+1/2]1) = O(∆x3) for any stationary solution and it vanishs for water at rest. 5 Numerical tests Let us consider the 1D SWE with pollutant : r=hϕ with ϕthe pollutant concentration. Then, we show an influence of the definition of S∗on the pollutant concentration. We consider a domain of 4 meters long discretized with 50 points, a CFL condition equal to 0.9, the initial conditions h(x, 0) = (18 −zb(x)), q(x, 0) = h(x, 0). We consider the topography zband rdefined by, Fig. 1 : zb(x) = 1 1.5<x<2.5 0 otherwise r(x) = h(x) 1.5<x<2.5 0 otherwise We set the final time to T= 0.5s. In Fig.1 a), we plot the computed pollutant concentration and the topography. We can observe the two different concentration values obtained with S∗defined by (2) (black dots) and (5) (squares) respectively. In the first case, S∗is not a consistent approximation of the velocity. In that case, we obtain a negative value near x= 2.5, see Fig. 1 b). At left of x= 2.5, S∗is positive and it is negative at x= 2.5. Thus, the pollutant is transported to the right until this point and transported to the left at x= 2.5 only. This produces the pollutant pick observed in Fig.1 a). At contrary, if we consider the consistent definition of S∗(5) then S∗is positive everywhere and the pollutant is transported smoothly, without any pick. Furthermore, we can notice in Fig.1 b) that S∗is defined at intercells x=xi+1/2, and the velocity is defined at points x=xi(line with cross). Then, we can observe that S∗defined by (5) is effectively a good approximation of the velocity at intermediate points. R´ef´erences [1] E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein, B. Perthame. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J.Sci. Comp 25, 6, 2050- 2065 (2004). 5 a) 01234 0 5 10 15 20 Water surface Discharge Conc. (S* non consistent) Conc. (S* consistent) Bottom b) 123 −1 0 1 2 3 S* non consistent S* consistent Velocity Fig. 1 – a) Pollutant concentration and topography; b) Value of S∗ [2] A. Berm´udez, M. E. V´azquez Cend´on. Upwind Methods for Hyperbolic Conservation Laws with Source Terms. Computers Fluids 23-8 1049-1071 (1994). [3] F. Bouchut. Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources. Birkh¨auser, (2004). [4] T. Chac´on Rebollo, A. Dom´ınguez Delgado, E.D. Fern´andez-Nieto, Asymptotically balanced schemes for non-homogeneous hyperbolic systemsapplication to Shallow Water equations, C.R. Acad. Sci. Paris. 338, 85-90(2004) . [5] T. Chac´on, A. Dom´ınguez, E.D. Fern´andez-Nieto. A family of stable numerical solvers for the shallow water equations with source terms. Comput. Methods Appl. Mech. Engrg. Vol. 192 pp. 203-225 (2003) [6] L. Gosse. A well-balanced scheme using non-conservative products designed for hyperbolic system of conservation laws with source terms. Mat. Mod. Meth. Appl. Sc. 11 : 339–365, (2001). [7] G. Dal Maso, P.G. LeFloch, F. Murat. Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74 : 483–548, (1995). [8] A. Dom´ınguez, E.D. Fern´andez-Nieto, S. Mart´ın, A well balanced third order HLL scheme for hyperbolic systems with source term. Application to Shallow Water Equations. Submitted (2006). [9] J.M. Greenberg, A. Y. Leroux. A Well-Balanced scheme for the numerical processing of source terms in hyperbolic equations. SIAM J. Numer. Anal. Vol. 33 No. 1, pp. 1-16 (1996). [10] P.G. LeFloch, A.E. Tzavaras. Representation and weak limits and definition of nonconservative products. SIAM J. Math. Anal. 30 : 1309–1342, (1999). [11] R.J. LeVeque Balancing Source Terms and Flux Gradients in High-Resolution Godunov Methods : The Quasi-Steady Wave-Propagation Algorithm. J. Comput. Phys. 146, 346-365 (1998). [12] C. Par´es, M.J. Castro.On the well-balanced property of Roe method for nonconservative hyperbolic systems. Applications to shallow water systems. ESAIM : M2AN, Vol 38(5) : 821–852, (2004). [13] B. Perthame, C. Simeoni. A kinetic scheme for the Saint-Venant system with a source term. Calcolo 38(4) : 201–231, (2001). [14] B. Perthame, C. Simeoni. Convergence of the Upwind Interface Source method for hyperbolic conservation laws. Proceeding of Hyp 2002, T Hou and E. Tadmor editor, Springer (2003). [15] E F. Toro. Shock-Capturing Methods for Free-Surface Shalow Flows. Springer (1999). [16] J.G. Zhou, D.M. Causon, C.G. Mingham, and D. M. Ingram. The Surface Gradient Method for the Treatment of Source Terms in the Sallow-Water Equations. J. Comput. Phys., Vol. 168 pp. 1-25 (2001). 6