Asymptotically balanced schemes for non-homogeneous hyperbolic systems - application to the Shallow Water equations

Numerical Analysis Asymptotically balanced schemes for non-homogeneous hyperbolic systems – application to the Shallow Water equations Tomás Chacón Rebolloa, Antonio Domínguez Delgadob, Enrique D. Fernández Nietob aDpto Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, c/Tarfia s/n, 41012 Sevilla, Spain bDpto Matemática Aplicada I, Universidad de Sevilla, Avda. Reina Mercedes N. 2, 41012 Sevilla, Spain Received 30 June 2003; accepted after revision 4 November 2003 Presented by Olivier Pironneau Abstract In this work we introduce a class of balanced numerical schemes, up to second order, for the solution of general nonhomogeneous hyperbolic systems of conservation laws. We give a general technique to build such schemes. We also prove that they balance up to second order a large class of steady solutions in the whole domain but some subset whose measure tends to zero as the grid size decreases to zero. We finally present an application to Shallow Water equations that exhibit the good performances of some of the schemes introduced. To cite this article: T. Chacón Rebollo et al., C. R. Acad. Sci. Paris, Ser. I 338 (2004). Schémas asymptotiquement equilibrés pour systèmes hyperboliques non-homogénes – application aux équations de Saint-Venant. Dans ce travail nous introduisons une classe de schémas numériques equilibrés au second ordre pour la solution de systèmes hyperboliques de lois de conservation. Nous donnons une technique générale pour construire ce type de schémas. Nous prouvons que ces schémas équilibrent au second ordre une grande classe de solutions stationaires, dans tout le domaine, excepté un sous-ensemble de mesure qui tend vers zéro lorsque la taille de la maille tend vers zéro. Nous présentons finalement une application aux équations de Saint-Venant qui montre les bonnes performances de quelques-uns des schémas présentés. Pour citer cet article:T. Chacón Rebollo et al., C. R. Acad. Sci. Paris, Ser. I 338 (2004). Version française abrégée Dans ce travail nous introduisons quelques schémas numériques equilibrés pour la résolution de systèmes hyperboliques de lois de conservation non homogènes. C’est un fait bien connu, et largement traité dans la littérature à cet égard, que pour qu’un schéma numérique fournisse des solutions avec un niveau de précision acceptable, il doit approcher les équilibres (solutions E-mail addresses: [email protected] (T. Chacón Rebollo), [email protected] (A. Domínguez Delgado), [email protected] (E.D. Fernández Nieto). 1631-073X/$ – see front matter doi:10.1016/j.crma.2003.11.008 86 T. Chacón Rebollo et al. / C. R. Acad. Sci. Paris, Ser. I 338 (2004) 85–90 stationnaires) du système, au moins au deuxième ordre. Autrement, les solutions obtenues présentent de fortes erreurs, étant souvent physiquement inacceptables. Le point-clé pour construire des schémas equilibrés est d’introduire un décentrage du terme source compatible avec le terme de diffusion numérique du schéma (cf. [1,5,6]). Dans ce travail nous présentons d’abord une technique générale pour construire des termes sources décentrés, de façon à garantir que les équilibres sont calculés au moins au second ordre (schéma (2)–(7)). L’idée consiste à rajouter un terme correcteur au système équivalent. Ceci s’applique a des systèmes hyperboliques de lois de conservationsnon-homogènesquelconques,à des schémas de type Flux-differenceet Flux-splitting et à une grande classe de solutions stationnaires (Théorème3.3). Essentiellement, nous démontronsque notre schéma équilibre les états stationnaires à l’ordre 2 dans les zones du domaine où les vitesses caractéristiques gardent un signe constant (même zéro). Asymptotiquement (lorsque la taille de la maille tends vers zéro), on équilibre ces états partout sauf sur un ensemble de mesure nulle. L’intérêt de ce résultat est qu’on n’a pas besoin de connaître les solutions stationnaires auparavant. Nousdonnonségalementdesconditionssuffisantespourquenotreschémaéquilibredefaçonexacteune solution stationnaire connue. Finalement, la technique proposée est appliquée aux équations de Saint-Venant : nous proposons une discrétisation concrète du terme source, et montrons quelques résultats numériques assez satisfaisants pour des tests significatifs (Section 4). 1. Introduction In this paper we study the numerical discretization of non-homogeneous hyperbolic systems of conservation laws. It is well known that in order that the scheme is accurate enough, it should solve some equilibria of the system up to second order. There is a wide literature which treats this problem. Bermúdez–Vázquez in [1] define the C-property; which asks that the scheme calculates a given stationary solution for a given system at grid nodes, exactly or up to second order. In [3] an extension of a flux-splitting scheme to non-homogeneous 1D SWE is derived. In [2], based on the ideas of [3], a systematic technique to build numerical schemes which verify the C-property is introduced. In [5] Greenberg–Leroux introduce the concept of a well-balanced scheme. In this paper, for non-homogeneous scalar equations, the authors introduce a numerical scheme based on the resolution of regularized Riemann problems, which preserve the balance between internal forces and the source term of the hyperbolic conservation law. The numerical scheme balances the equilibria of non-homogeneousequations. In the case of such a system, in [6], Jin introduces the interface-method. He proves that the proposed scheme balances up to second order all stationary solutions for the system of Shallow Water equations, at cell interfaces. In [7] Perthame–Simeoni prove an extension of the theorem of Lax–Wendroff to non-homogeneous scalar conservation laws. The remarkablefact, in our context, is that the theorem holds,providedthe scheme balances certain stationary solutions. In the present work we introduce, for general non-homogeneoushyperbolic systems, numerical schemes which balance virtually all stationary solutions up to second order. We consider a general formulation of numerical schemes, under which can be written a wide class of known and new schemes. We present an application to Shallow Water equations that exhibits the good performances of the schemes introduced. 2. Motivation In this section, we give the general structure of the schemes that we propose. We consider a non-homogeneoushyperbolic system, ∂W ∂t +∂F ∂x =G. (1) T. Chacón Rebollo et al. / C. R. Acad. Sci. Paris, Ser. I 338 (2004) 85–90 87 The structure of schemes in conservative form for the hyperbolic system is obtained by integrating the differential system in the control volume, in the 1D case, (xi−1/2,xi+1/2),wherexi=ix. Then, for non-homogeneous hyperbolic systems with an explicit discretization in time, if t is the time step, we have Wn+1 i−Wn i t + φn i+1/2−φn i−1/2 x =Gn C,i,(2) where Gn C,i is a centered (second order)approximationof 1 x xi+1/2 xi−1/2G(x, Wn h)dx,Wn hbeing a function built from {Wn i}i. The numerical flux function, φi+1/2is an approximation of Fat x=xi+1/2. We consider upwind schemes, that treat the convection dominance by introducing some numerical diffusion. For these schemes, the numerical flux function has the structure φi+1/2=FC(Wi,W i+1)−1 2D(Wi,W i+1)(Wi+1−Wi), (3) where by FC(Wi,W i+1), we denote a centered (second order) approximation of F(W) in x=xi+1/2,and D(Wi,W i+1)is a first-order approximation of the diffusion matrix D( Wi+1/2),where Wi+1/2is an intermediate state between Wiand Wi+1. These schemes can be interpreted as second order approximationsin the space of a parabolic equivalent system. If we also take a centered approximation of the source term, then this equivalent system is ∂W ∂t +∂ ∂xF(W)−ν∂ ∂xD(W) ∂ ∂xW=G(x, W), (4) where νis equal to half of the space step (ν=x/2). Notice that if Wis a stationary solution of (1) then, it verifies ∂ ∂xF(W)=G(x, W).ThenWis not a stationary solution of (4), as thereis a first order error dueto the upwindingof the flux. It is well knownthat this is a source of large errors in the numerical solution. These large errors disappear when the scheme has a high precision in space (cf. [1,5]). We shall develop here schemes with second order accuracy in space. For this purpose, we propose to add a correcting term (that we denote by VST, Viscosity Source Term) to the equivalent system: ∂W ∂t +∂ ∂xF(x,W)−ν∂ ∂xD(x, W )∂W ∂x =G(x, W ) +VST.(5) The term VST, for stationary solutions of (1) must verify −ν∂ ∂x(D(W)∂W ∂x )=VST. To give the definition of VST we consider a smooth stationary solution Wof (1). By A(W), we denote the Jacobian matrix of F.IfA(W) is non singular, then A(W)∂W ∂x =G(x, W ) ⇒−ν∂ ∂xD(W)∂W ∂x =−ν∂ ∂xD(W)A−1(W)G(x, W). With this definition of VST, we can write the system (5) (with a conservative structure) as ∂W ∂t +∂ ∂xF(W)−νD(W)∂W ∂x −A−1G(x, W )=G(x, W). (6) We can conclude then that to build balanced numerical schemes, up to second order, for the hyperbolic system (1), it is enough to use approximationsof order two in space of (6). This motivates the definition of φi+1/2, as a second order approximation of F(W)−νD(W)(∂W ∂x −A−1G(x, W )) in x=xi+1/2. Concretely, we define φi+1/2=FC(Wi,W i+1)−νD(Wi,W i+1)Wi+1−Wi x − A−1(Wi,W i+1)GD(xi,xi+1,W i,W i+1),(7) where GDand  A−1are first-order approximations of Gand A−1( Wi+1/2), respectively. 88 T. Chacón Rebollo et al. / C. R. Acad. Sci. Paris, Ser. I 338 (2004) 85–90 There is still the difficulty of defining VST when Ais singular. If by λj,j=1,...,N, we denote the eigenvalues of A, we shall assume that there exists Ncontinuous surfaces γj⊂RN,j=1,...,N, such that λjis invertible and C2in RN\γj. We propose to define the matrix  A−1(U, V ) as  A−1(U, V ) =X( W(U,V)) Λ−1(U, V )X−1( W(U,V)) where  W(U,V) is an intermediate state between Uand V,byXwe denote the matrix defined by the eigenvectors of A and  Λ−1(U, V ) =Diag λ−1 1(U, V ), . . . ,  λ−1 N(U, V )with  λ−1 j=λ−1 j Wif L[U,V]⊂RN\γj, 0otherwise, where L[U,V]is the segment in RNthat connects Uand V. This definition of  A−1states that when some eigenvalueof Avanishes, then no upwindingof the source term in characteristics variables must be performed. For scalar equations, we may prove that this definition of  A−1ensures the stability of the scheme. This choice allows us to balance the scheme, even in some cases when Ais singular. 3. Balance properties In this section we prove that the scheme (2)–(7) balances up to second order the system (1) for a wide class of stationary solutions. We also prove that under some natural hypotheses upon matrix D, the defined scheme balances system (1) for all stationary solutions in all [0,L]but in a set whose measure tends to zero as x →0. We will call them ‘asymptotically balanced schemes’. Firstly, we observe that matrix D(W) is only a Lipschitz function, with non-continuous derivative at the points where its eigenvaluesvanish. A possible choice of the matrix is D(W) =|A(W)|; for example this happens for the scheme of Roe. Thus, we cannot hope that the introduced scheme balances all stationary solutions, with order two, in all the domain [0,L]. We introduce the following definitions, Definition 3.1 (Balanced schemes).We say that the numerical scheme (2)–(7) is balanced on a smooth stationary solution W(x) of the hyperbolic system (1) if the associated consistency error has second order accuracy. Definition 3.2 (Asimptotically balanced schemes).We say that the scheme (2)–(7) is asymptotically balanced on a smooth stationary solution W(x) of the hyperbolic system (1) if there exists an increasing sequence {Kn}nof compact subsets of [0,L]such that (1) µ([0,L]nKn)=0, where by µwe denote Lebesgue measure in R. (2) For all nthere exists a δn>0 such that if 0 <x<δ n, then the scheme balances the system in Kn. Assuming that A−1and Dpresent singularities only on curves in RN, we obtain the following result: Theorem 3.3. We consider W:[0,L]→RNa stationary solution of class C2of the system (1). Then, (a) If A(W) does not have singular points, the scheme (2)–(7) balances system (1) on Wup to second order. (b) If the set of points xwhere A(W(x)) is singular has measure zero, then the scheme (2)–(7) asymptotically balances system (1) on W. (c) If Dhas the same eigenvectors as Aand their eigenvalues vanish at the same points, then the scheme (2)–(7) asymptotically balances system (1) on W. The proof of this theorem is based on the fact that in regions of [0,L]where Ais inversible, a smooth solution of the original system (1) is also a solution of the equivalentsystem (6). Point (c) applies to steady solutions Wfor T. Chacón Rebollo et al. / C. R. Acad. Sci. Paris, Ser. I 338 (2004) 85–90 89 which the set of point where A(W) is singular has non-zero measure. The complete proof of the theorem can be found in [4]. The interest of point (c) is that it provides a wide class of schemes that balance all (smooth) steady solutions of system (1). We also havethe following resultthat yields sufficient conditionsfor the exactcalculation of stationary solutions. Theorem 3.4. Let Wbe a stationary solution of (1). Assume that FC,GC, and GDevaluated on Wverify FC(Wi,W i+1)−FC(Wi−1,W i) x =GC,i,Wi+1−Wi x = A−1(Wi,W i+1)GD(xi,xi+1,W i,W i+1). Then, the scheme defined by (2)–(7) balances in an exact way the system (1) for the stationary solution at grid points. This happens independently of the choice of the upwinding matrix D(U,V ). 4. Application to Shallow Water equations In this section we apply the proposed scheme to Shallow Water equations, which model the behaviour of water flows in channels. The unknowns are h, the height of the water, and qthe discharge. The physical flux functionand the source term for variable depth, are defined by F(W)=q q2 h+1 2gh2,G(x,W)=0 ghH(x) where W=h q, gis the gravity constant and H(x)=¯ h−zb(x) is the depth of the channel, about a fixed reference level. The different elements that determine scheme (2)–(7) that we use are defined as follows: FC(Wi,W i+1)=F(W i+α)+F(W i+1−α)/2,where Wi+α=(1−α)Wi+αWi+1with α∈[0,1]. The first component of GDand GCis equal to zero, in both cases, and the second components are GD(xi,xi+1,W i,W i+1)2=ghi+1/2Hi+1−Hi x , [GC,i]2=g 2αhi+α/2+(1−α)hi+(1−α)/2Hi+1−Hi x +αhi−α/2+(1−α)hi+(α−1)/2Hi−Hi−1 x . Matrix Dhas the same eigenvectors as matrix Aand its eigenvalues are given by dj(Wi,W i+1)=|λj(Wi+1/2)|. With this choice, at least for the linear homogeneouscase, the scheme is L2stable. Using Theorem 3.3 we prove that the scheme is asymptotically balanced for all smooth stationary solutions. Moreover, the scheme exactly calculates water at rest, (h, q) =(H, 0). It is enough to observe that the conditions of the Theorem 3.4 are verified. Numerical test. We test our scheme for two stationary solutions, in a domain of length L=25 m and bottom function defined by the maximum between zero and the ‘bump’ zb(x) =0.2−0.05(x −10)2. For this test the discharge must be constant, so this test is interesting in order to observe the performance of the scheme near the bottom bump. We use a CFL condition equal to 0.8, a step discretization equal to x =0.25 and as initial condition water at rest with constant free surface at 0.5. As boundary conditions, we consider two different cases: for Test 1.a we impose h=0.66 m downstream when the flow is subcritical and q=1.53 m2/s upstream. Test 1.b is obtained by imposing h=0.33 downstream and q=0.18 m2/s upstream. In both cases we consider that the scheme converges to a stationary solution when the relative error between two consecutive approximations in time is less than 10−5. 90 T. Chacón Rebollo et al. / C. R. Acad. Sci. Paris, Ser. I 338 (2004) 85–90 Fig. 1. Test 1.a. Fig. 2. Test 1.b. In Figs. 1 and 2 we compare the calculated discharge for the Surface Gradient Method (SGM) (cf. [8]) and the proposedschemewith α=1/8.Inbothcasestheproposedschemeprovidessomeimprovement.Thisisparticularly relevant for Test 1.b, as in this case there is a hydraulic jump (discontinuity of the free surface) near x=12. Acknowledgements This work has been partially financed by the Spanish Government Research projects REN 2000-1162-C02-01, REN 2000-1168-C02-01and CYTMAR MAR97-1055-C02-02. References [1] A. Bermúdez, M.E. Vázquez Cendón, Upwind methods for hyperbolic conservation laws with source terms, Comput. & Fluids 23 (8) (1994) 1049–1071. [2] T. Chacón, A. Domínguez, E.D. Fernández-Nieto, A family of stable numerical solvers for the shallow water equations with source terms, Comput. Methods Appl. Mech. Engrg. 192 (2003) 203–225. [3] T. Chacón, E.D. Fernández-Nieto, M. Gómez, A flux-splitting solver for shallow water equations with source terms, Int. J. Numer. Methods Fluids 43 (2003) 23–55. [4] E.D. Fernández-Nieto, Aproximación numérica de leyes de conservación hiperbólicas no homogéneas. Aplicación a las ecuaciones de Aguas Someras, Ph.D. Thesis Universidad de Sevilla, 2003. [5] J.M. Greenberg, A.Y. Leroux, A well-balanced scheme for the numerical processing of source terms in hyperbolic equations, SIAM J. Numer. Anal. 33 (1) (1996) 1–16. [6] S. Jin, A steady-state capturing method for hyperbolic systems with geometrical source terms, Math. Modelling Numer. Anal. 35 (4) (2001) 631–645. [7] B. Perthame, C. Simeoni, Convergence of the Upwind Interface Source method for hyperbolic conservation laws, in: T. Hou, E. Tadmor (Eds.), Proceeding of Hyp 2002, Springer, 2003. [8] 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 (2001) 1–25.