Christophe Chalons

Comptes Rendus Mathématique

C. R. Acad. Sci. Paris, Ser. I 342 (2006) 623–626http://france.elsevier.com/direct/CRASS1/Numerical AnalysisTransport-equilibrium schemes for computing nonclassical shocksChristophe ChalonsLaboratoire J.-L Lions, U.M.R. 7598, université P. et M. Curie, boîte courrier 187, 75252 Paris cedex 05, FranceReceived 23 June 2005 ; accepted after revision 7 February 2006Presented by Olivier PironneauAbstractThis Note presents a very efficient numerical strategy for computing weak solutions of a scalar conservation law which fails tobe genuinely nonlinear. In such a situation, the dynamics of shock solutions turns out to be mainly driven by a prescribed kineticfunction that imposes the speed of propagation of the discontinuities. We show how to enforce the validity of the kinetic criterion atthe discrete level. The resulting scheme provides, in addition, sharp profiles. Numerical evidence are included. To cite this article:C. Chalons, C. R. Acad. Sci. Paris, Ser. I 342 (2006). 2006 Académie des sciences. Published by Elsevier SAS. All rights reserved.RésuméSchémas transport-équilibre pour l’approximation des chocs nonclassiques. Cette Note présente un algorithme très efficacepour le calcul des solutions faibles d’une loi de conservation scalaire non vraiment nonlinéaire. Dans ce contexte, la dynamiquedes solutions choc repose principalement sur la donnée d’une fonction cinétique qui fixe la vitesse de propagation des disconti-nuités. Nous montrons comment forcer la validité du critère cinétique au niveau discret. Le schéma obtenu fournit par ailleurs desdiscontinuités sans diffusion numérique. Des résultats numériques sont présentés. Pour citer cet article : C. Chalons, C. R. Acad.Sci. Paris, Ser. I 342 (2006). 2006 Académie des sciences. Published by Elsevier SAS. All rights reserved.1. IntroductionWe are interested in computing nonclassical weak solutions of an initial-value problem for a scalar conservationlaw of the form{∂tu + ∂xf (u) = 0, u(x, t) ∈ R, (x, t) ∈ R × R+∗,u(x,0) = u0(x), (1)where f : R → R is a (smooth) nonconvex flux-function. Generally speaking, solutions of problem (1) may be discon-tinuous and are not uniquely determined by initial data u0. According to a general regularization principle, we thusask solutions of (1) to satisfy a single entropy inequality of the form∂tU(u) + ∂xF (u)  0, (2)E-mail address: chalons@math.jussieu.fr (C. Chalons).1631-073X/$ – see front matter  2006 Académie des sciences. Published by Elsevier SAS. All rights reserved.doi:10.1016/j.crma.2006.02.003624 C. Chalons / C. R. Acad. Sci. Paris, Ser. I 342 (2006) 623–626where U : R → R and F : R → R are specified functions such that U is strictly convex and F ′ = U ′f ′. When f isconvex, entropy condition (2) actually selects a unique classical solution of (1). When f fails to be convex, it is nec-essary to supplement (1)–(2) with an additional selection criterion called kinetic relation from [4]. More precisely, theRiemann problem associated with (1)–(2) still admits a one-parameter family of solutions, which may contain shockwaves violating Lax shock inequalities. Such discontinuities are referred as to undercompressive shocks or nonclas-sical shocks. In order for the uniqueness to be ensured, a kinetic relation needs to be added along each nonclassicaldiscontinuity connecting a left state u− to a right state u+. It takes the form u+ = ϕ(u−) or u− = ϕ−(u+) whereϕ is the so-called kinetic function and ϕ− its inverse. We refer to [4] for a general theory of nonclassical entropysolutions.The numerical approximation of nonclassical solutions is known to be very challenging. The main difficulty is therespect of the kinetic relation at the discrete level. At once, such nonclassical solutions are actually present in real andcomplex problems coming from the physics. Let us mention for instance the study of compressible fluids undergoingvapor-liquid phase changes, which is our main motivation.In this Note, we present a new scheme for capturing discontinuities whose dynamics is driven by a kinetic function.Our strategy deals directly with the kinetic function ϕ to tackle the nonclassical solutions. The resulting algorithmprovides numerical results in full agreement with exact ones, whatever the strength of the shocks are. In particular,our scheme leaves sharp isolated nonclassical shocks.2. The case of cubic flux and nonclassical Riemann solverWithout loss of generality, we take f (u) = u3 which is to some extent the simplest example of a nonconvex func-tion, and refer to [1] for more general flux functions. We consider weak solutions of (1) satisfying entropy inequality(2) with U(u) = u2 and F(u) = 34u4, and choose (again without restriction) ϕ(u) = −βu as a kinetic function, withβ ∈ [1/2,1) so that each nonclassical shock obeys (2).Given two constant states ul, ur such that ul > 0, we consider a Riemann initial data u0 defined by u0(x) = ul ifx < 0 and u0(x) = ur if x > 0. Following [4] and defining ϕ(u) = −u−ϕ(u), the solution of (1)–(2), supplementedwith the kinetic criterion associated with ϕ, is given as follows:(1) If ur  ul , the solution is a rarefaction wave connecting ul to ur .(2) If ur ∈ [ϕ(ul), ul), the solution is a classical shock wave connecting ul to ur .(3) If ur ∈ (ϕ(ul), ϕ(ul)), the solution contains a nonclassical shock connecting ul to ϕ(ul), followed by a classicalshock connecting ϕ(ul) to ur .(4) If ur  ϕ(ul), the solution contains a nonclassical shock connecting ul to ϕ(ul), followed by a rarefactionconnecting ϕ(ul) to ur .3. Numerical approximationWe now present a suitable algorithm for approximating these Riemann solutions. Let be given a time step t and aspace step x. Introducing xj+1/2 = jx for j ∈ Z and tn = nt for n ∈ N, we seek at each time tn an approximationunj of u on each interval [xj−1/2;xj+1/2). In this context, we choose without restriction a two-point numerical fluxfunction g : (u, v) → g(u, v) consistent with the flux function f .Let us first motivate our algorithm. Actually, it is observed in [1,2] that a classical conservative scheme associatedwith g is not able to propagate any isolated nonclassical shock. Instead, some spurious values are created by thescheme, eventually leading to a damaged numerical solution. Our algorithm aims at removing these spurious values.This is achieved in the first step of the method by systematically making stationary the nonclassical discontinuities.Which means in particular that the conservation property is lost, so that the use of a nonconservative update formulais completely natural in this first step. Then, the dynamics of the nonclassical discontinuities is taken into account inthe second step of the strategy, thus restoring (in a weak sense) the conservation property.First step (tn → tn+1−) As motivated above, this step aims at making stationary the nonclassical discontinuitiesof our problem. For that, we introduce the following nonconservative update formula:un+1− = un − λ(gL − gR ), j ∈ Z, with λ = t/x, (3)j j j+1/2 j−1/2C. Chalons / C. R. Acad. Sci. Paris, Ser. I 342 (2006) 623–626 625where the numerical fluxes gLj+1/2 and gRj+1/2 have to be suitably defined. In this note, we will focus ourselves on thenumerical approximation of the nonclassical discontinuities only (the most difficult ones to capture numerically), thatis those separating two states u− and u+ such that u+ = ϕ(u−) < ϕ(u−) when u− > 0. Then, it is easily checkedthat defining gLj+1/2 and gRj+1/2 as follows when unj > 0:gLj+1/2 ={g(unj , ϕ−(unj+1)) if unj+1 < ϕ(unj ),g(unj , unj+1) otherwise,(4)gRj+1/2 ={g(ϕ(unj ), unj+1) if unj+1 < ϕ(unj ),g(unj , unj+1) otherwise,(5)and, in a first approach at least, equal to gj+1/2 when unj  0 (see [1] for details) is sufficient to keep at stationaryequilibrium all the discontinuities separating two states u− and u+ such that u+ = ϕ(u−).Second step (tn+1− → tn+1) This step deals with the dynamics of the nonclassical discontinuities left stationaryduring the first step. We first recall that the speed of propagation σ(u−, u+) of a discontinuity between u− and u+is given by Rankine–Hugoniot conditions, that is σ(u−, u+) = [f (u+) − f (u−)]/[u+ − u−]. We then define at eachinterface xj+1/2 a speed of propagation σj+1/2 byσj+1/2 ={σ(un+1−j , un+1−j+1 ) if unj+1 < ϕ(unj ),0 otherwise,(6)and solve at each discontinuity xj+1/2 a transport equation with speed σj+1/2. In order to avoid the appearance ofnew spurious values and to get a new approximation un+1j at time tn+1 = tn + t , we propose to pick up randomlyon interval [xj−1/2, xj+1/2[ a value in the juxtaposition of these Riemann solutions at time t chosen small enoughto avoid wave interactions. Given a well distributed random sequence (an) within interval (0,1), it amounts to set:un+1j =un+1−j−1 if an+1 ∈ [0, λσ+j−1/2[,un+1−j if an+1 ∈ [λσ+j−1/2,1 + λσ−j+1/2[,un+1−j+1 if an+1 ∈ [1 + λσ−j+1/2,1[,with{σ−j+1/2 = min(σj+1/2,0),σ+j−1/2 = max(σj−1/2,0).(7)4. Numerical experimentsWe give two numerical evidences to validate our scheme. We take a Roe scheme as a basic numerical flux g, andthe van der Corput random sequence for (an). For the kinetic function ϕ, we set β = 34 . We consider the typicalnonclassical behaviors of the Riemann solution given in Section 2, when taking ul = 4 and ur equal to −2 (test 1)or −5 (test 2). Numerical solutions are plotted on Fig. 1 for x = 10−2.Fig. 1. Nonclassical solutions: test 1 (left) and test 2 (right).Fig. 1. Solutions nonclassiques : test 1 (gauche) et test 2 (droite).626 C. Chalons / C. R. Acad. Sci. Paris, Ser. I 342 (2006) 623–626We observe that the numerical solutions fully agree with the exact ones. In particular, the right states of the non-classical waves are exactly captured, which is remarkable, while there are not any points in their profile (see [1] andthe references therein). For test 1, we note that the classical shock contains some numerical diffusion induced by theRoe scheme. In [1], we show how to slightly modify the definitions of the numerical fluxes gLj+1/2 and gRj+1/2 in(4)–(5) in order to make sharp the classical shocks, too.To conclude, an efficient numerical strategy has been presented for computing nonclassical solutions of a particularscalar conservation law. It provides sharp nonclassical interfaces propagating at the right speed whatever the strengthof interfaces are. The method turns out to be nonconservative but measures in [1] have shown that the loss of conserva-tion is extremely low, while numerical solutions fully agree with exact ones. We emphasize that our approach is builtwithout explicitly using the knowledge of the underlying nonclassical Riemann solver, contrary to Glimm’s method forinstance. Which means that it is not an expensive method and that it may be used for more complex applications. Anapplication to pedestrian flows has been addressed in [2]. In [3], we consider the case of a hyperbolic-elliptic modelwith phase changes that has initially motivated the present study.References[1] C. Chalons, Transport-Equilibrium schemes for computing nonclassical shocks. I. Scalar conservation laws, Laboratoire J.-L. Lions, PreprintR05031, 2005. Available at http://www.ann.jussieu.fr.[2] C. Chalons, Numerical approximation of a macroscopic model of pedestrian flows, Laboratoire J.-L. Lions, Preprint R05032, 2005. Availableat http://www.ann.jussieu.fr.[3] C. Chalons, Transport-equilibrium schemes for computing nonclassical shocks, II. A hyperbolic-elliptic model of phase transitions, in prepara-tion.[4] P.G. LeFloch, Hyperbolic Systems of Conservation Laws: The Theory of Classical and Nonclassical Shock Waves, E.T.H. Lecture Notes Series,Birkhäuser, 2002.

Transport-equilibrium schemes for computing nonclassical shocks

https://comptes-rendus.academie-sciences.fr/mathematique/item/10.1016/j.crma.2006.02.003.pdf

Transport-equilibrium schemes for computing nonclassical shocks

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Mathématique