Late-time/stiff relaxation asymptotic-preserving approximations of hyperbolic equations
aa r X i v : . [ m a t h . A P ] S e p LATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVINGAPPROXIMATIONS OF HYPERBOLIC EQUATIONS
CHRISTOPHE BERTHON, PHILIPPE G. L E FLOCH, AND RODOLPHE TURPAULT
Abstract.
We investigate the late-time asymptotic behavior of solutions tononlinear hyperbolic systems of conservation laws containing stiff relaxationterms. First, we introduce a Chapman-Enskog-type asymptotic expansion andderive an effective system of equations describing the late-time/stiff relaxationsingular limit. The structure of this new system is discussed and the roleof a mathematical entropy is emphasized. Second, we propose a new finitevolume discretization which, in late-time asymptotics, allows us to recover adiscrete version of the same effective asymptotic system. This is achieved pro-vided we suitably discretize the relaxation term in a way that depends on amatrix-valued free-parameter, chosen so that the desired asymptotic behavioris obtained. Our results are illustrated with several models of interest in con-tinuum physics, and numerical experiments demonstrate the relevance of theproposed theory and numerical strategy. Introduction
Motivations.
We are interested in the numerical approximation of solutions tononlinear hyperbolic systems of conservation laws containing stiff relaxation terms.An extensive literature is available on such systems since they arise in many phys-ical problems of interest, for instance, in the modeling of complex multiphaseflows involving phase transitions or kinetic-type phenomena. (See, for instance,[12, 13, 14, 22, 30].) Stiff relaxation source terms are essential to model phenomenainvolving distinct physical time-scales. The derivation and the analysis (existence,stability) of an effective system of equations (also of hyperbolic type as the origi-nal system) as the relaxation time goes to zero is required for understanding thebehavior of general solutions. (See, for instance, [25, 31, 34, 35].)In the present paper, we go beyond these classical works and investigate the late-time behavior of solutions to systems with stiff relaxation and, specifically,we consider the class of systems ( ε > ε ∂ t U + ∂ x F ( U ) = − R ( U ) ε , t > , x ∈ R , where the main unknown U : R × R + → Ω takes its values in a convex set Ω ⊂ R N . The associated homogeneous first-order system —obtained by neglecting therelaxation term R : Ω → R N in the right-hand side of (1.1)— is assumed to behyperbolic in the sense that the N × N matrix A := D U F admits real eigenvaluesand a full basis of eigenvectors. Mathematics Subject Classification.
Key words and phrases.
Nonlinear hyperbolic system, stiff source term, late-time behavior,diffusive regime, finite volume scheme, asymptotic preserving. E FLOCH, AND R. TURPAULT
In comparison with classical works on this subject [11] as well as also [1, 4, 10,25], the main novelty in the present work lies in the fact that the term ε ∂ t U isrescaled to be proportional to ε . When ε tends to zero, this will lead us to limitsolutions whose behavior is quite distinct from the ones of standard relaxationlimits. Indeed, as we establish below, the effective problem associated with (1.1)turns out to be a diffusion problem, in which the diffusion operator is determinedby the (nonlinear) relaxation term R , in a way explicitly determined in the presentwork . The relaxation map R : Ω → R N is assumed to be sufficiently regularand satisfy the conditions introduced by Chen, Levermore, and Liu [11] for the“hyperbolic to hyperbolic” relaxation problem.First of all, we assume the existence of a (constant) n × N matrix Q with (max-imal) rank n < N such that(1.2) QR ( U ) = 0 , U ∈ Ω . From this, it follows that, if U is a solution to (1.1), then QU : R × R + → R N satisfies the system of n conservation laws(1.3) ε ∂ t ( QU ) + ∂ x ( QF ( U )) = 0 , in which the unknown takes values in the (convex) set ω := Q Ω ⊂ R n . We also assume that a map E : ω → Ω uniquely determines local equilibria, asdefined by the relations(1.4) Q E ( u ) = u, R ( E ( u )) = 0 , u ∈ ω. This suggests to introduce the equilibrium manifold M := (cid:8) U = E ( u ) (cid:9) . The dimension of the null space of the N × N matrix B := DR with the equilibriumsubmanifold is assumed to be “maximal” in the sense thatdim (cid:16) ker( B ( E ( u ))) (cid:17) = n, (1.5) ker (cid:0) B ( E ( u )) (cid:1) ∩ Im (cid:0) B ( E ( u )) (cid:1) = { } . (1.6)Finally, since we are interested in the late-time behavior of solutions it is necessaryto impose that QF ( E ( u )) = c, u ∈ ω, for some constant vector c ∈ R n . Indeed, in view of (1.3), we can formally write ε ∂ t u + ∂ x QF ( E ( u )) → , so that QF ( E ( u )) must be a constant, normalized so that(1.7) QF = 0 on M . After this paper was completed, P. Marcati pointed out to the authors the relevant references[15, 16] in which several convergence results are established.
ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 3
Main results.
Our objective in this paper is two-fold. First, we determine aneffective system of equations driving the late-time asymptotic behavior of generalsolutions to (1.1) in the singular regime ε →
0. Second, we introduce a novel nu-merical strategy allowing us to precisely recover the expected asymptotic behavior.The paper is organized as follows. Section 2 is devoted to the asymptotic analysis:we derive an effective system and, under the hypotheses stated in this introductionas well as the assumption of the existence of a mathematical entropy compatiblewith the relaxation source term (in a sense specified below), we study the structureof this system which is found to be of diffusion-type and we show that the associatedtotal entropy is non-increasing in time.Then, in Section 3 below, we generalize our analysis to a nonlinear version of(1.1), since this is required to encompass certain models arising in the aplicationswhen the relaxation has strong nonlinearities. The class of systems under consid-eration in this section reads(1.8) ε ∂ t U + ∂ x F ( U ) = − R ( U ) ε m , where the parameter m ≥ additional scale in the problem. Byimposing that R (cid:0) E ( u ) + ε U (cid:1) = ε m R (cid:0) E ( u ) + M U (cid:1) , U ∈ Ω , u ∈ ω, plus certain conditions on the N × N matrix M (see Section 3), we derive again aneffective system which, now, consists of nonlinear diffusion equations. Next, in Section 4, we discuss several specific examples of practical interest andshow that they fit into our general framework: (1) the Euler equations of compress-ible fluids with friction term, (2) the so-called M Late-time/stiff-relaxation framework
Derivation of the effective system.
Throughout this section, we consider thenonlinear system of balance laws (1.1), and we impose the conditions stated in theintroduction.Our objective is to exhibit the system of effective equations satisfied by localequilibria u = u ( t, x ) ∈ ω . In the spirit of Chapman-Enskog expansions of the C. BERTHON, P.G. L E FLOCH, AND R. TURPAULT kinetic theory [8], we consider a formal expansion of solutions U to (1.1) in theform:(2.1) U ε = E ( u ) + ε U + ε U + . . . , where U i is referred to as the i th -order corrector. Such an expansion is natural inview of the assumptions (1.4) and (1.7). We consider the system(2.2) ε ∂ t U ε + ∂ x F ( U ε ) = − ε R ( U ε ) ,QU ε = u, in which we plug the formal expansion (2.1) and match together terms of the sameorder of magnitude in ε .From QU ε = u we deduce that Q E ( u ) = u , and QU i = 0 for all i th -ordercorrectors. Next, for sufficiently regular flux F , we obtain the following expansion:(2.3) F ( U ε ) = F ( E ( u )) + εA ( E ( u )) U + ε D U F ( E ( u )) ( U , U )+ ε A ( E ( u )) U + O ( ε ) . Analogously, it is useful to write down the formal expansion of the relaxation term.Taking (1.4) into account, we obtain(2.4) 1 ε R ( U ε ) = B ( E ( u )) U + ε D U R ( E ( u )) . ( U , U ) + εB ( E ( u )) U + O ( ε ) . We then plug (2.3) and (2.4) into (2.2) and obtain(2.5) ε∂ t E ( u ) + ∂ x F ( E ( u )) + ε∂ x A ( E ( u )) U = − B ( E ( u )) U − ε D U R ( E ( u )) . ( U , U ) − εB ( E ( u )) U + O ( ε ) . Let us consider (2.5) and expand in ε . The zeroth-order terms give ∂ x F ( E ( u )) = − B ( E ( u )) U . At this stage, we can solve for U by recalling that QU = 0. Indeed, since we haveimposed the properties (1.5) and (1.6) on the null space of B ( E ( u )) then, for allfixed J ∈ R N , the system(2.6) B ( E ( u )) .V = J,QV = 0 , admits a unique solution V ∈ R n if and only if QJ = 0. (See Lemma 2.4 at the endof this section.) Here, we have J = − ∂ x F ( E ( u )) which satisfies Q∂ x F ( E ( u )) = 0by (1.7). As a consequence, for all u ∈ ω , we can uniquely determine U such that(2.7) B ( E ( u )) .U = − ∂ x F ( E ( u )) ,QU = 0 . Next, we consider the first-order terms in (2.5), and obtain ∂ t E ( u ) + ∂ x A ( E ( u )) U = − D U R ( E ( u )) . ( U , U ) − B ( E ( u )) U . Multiplying by Q and using the assumption (1.2), we obtain QD U R ( E ( u )) . ( U , U ) ≡ , QB ( E ( u )) U ≡ . ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 5
Since Q E ( u ) = u , we arrive at an effective system for the limit u , that is,(2.8) ∂ t u = − ∂ x ( QA ( E ( u )) U ) , where U is the unique solution to (2.7). The role of a mathematical entropy.
Now, we assume the existence of a suffi-ciently regular, mathematical entropy Φ : Ω → R so that the matrix D U Φ( U ) A ( U )is symmetric for all U in Ω and there exists an entropy-flux map Ψ : Ω → R suchthat(2.9) D U Φ( U ) A ( U ) = D U Ψ( U ) , U ∈ Ω . Hence, all smooth solutions to (2.2) satisfy(2.10) ε∂ t Φ( U ε ) + ∂ x Ψ( U ε ) = − ε D U Φ( U ε ) R ( U ε ) . As usual, we impose Φ to be convex by requiring the N × N matrix D U Φ( U ) to bepositive. In addition, we assume that the entropy is compatible with the relaxationin the sense that, for some map ν : M → R n , D U Φ( U ) R ( U ) ≥ , U ∈ Ω , ‘(2.11) D U Φ = ν Q on M . (2.12)We now analyze the nature of the limiting system (2.8). Theorem 2.1.
Consider the nonlinear system of balance laws (1.1) , under theassumptions (1.2) – (1.7) and (2.11) – (2.12) . Then, the associated limiting system(2.8) takes the form (2.13) ∂ t u = ∂ x (cid:16) S L − ( u ) S T ( ∂ x D u Φ( E ( u ))) T (cid:17) , where (2.14) S := QA ( E ( u )) and (2.15) L ( u ) := D U Φ( E ( u )) B ( E ( u )) , and, for all b ∈ R N with Qb = 0 , L − ( u ) .b denotes the unique solution to the system L ( u ) .V = b, QV = 0 . In addition, this system is dissipative with respect to the entropy Φ in the sense thatthe following positivity condition holds for all u in ω : (cid:0) ∂ x D u Φ( E ( u )) (cid:1) S L − ( u ) S T (cid:0) ∂ x D u Φ( E ( u )) (cid:1) T ≥ . Proof.
We follow the strategy in [11] which we adapt to our problem. We firstestablish (2.13). To simplify the notation, we set D ( u ) = − QA ( E ( u )) U , to restate (2.8) as follows: ∂ t u = ∂ x D ( u ) . Here, the vector U ∈ R N is the unique solution of (2.7). Since D U Φ( U ) is a positive N × N matrix for all U ∈ Ω, we can rewrite (2.7) as follows:(2.16) D U Φ( E ( u )) B ( E ( u )) .U = − D U Φ( E ( u )) ∂ x F ( E ( u )) ,QU = 0 . C. BERTHON, P.G. L E FLOCH, AND R. TURPAULT
Invoking (2.15), the state vector U is defined as the unique solution of L ( u ) U = − D U Φ( E ( u )) ∂ x F ( E ( u )) ,QU = 0 . As a consequence, by definition of L − ( u ) we have U = −L − ( u ) D U Φ( E ( u )) ∂ x F ( E ( u )) . From (2.14), we find D ( u ) = S L − ( u ) D U Φ( E ( u )) ∂ x F ( E ( u )) . To obtain (2.13), we then establish(2.17) D U Φ( E ( u )) ∂ x F ( E ( u )) = S T v, where we define v ∈ R n as follows:(2.18) v = ( ∂ x D u Φ( E ( u ))) T . To this end we differentiate (2.12) (which involves the map ν ) and we obtain( D u E ( u )) T D U Φ( E ( u )) = D u Φ( E ( u )) Q. By transposition, we thus get(2.19) D U Φ( E ( u )) D u E ( u ) = Q T D U Φ( E ( u )) . Now, we have D U Φ( E ( u )) ∂ x F ( E ( u )) = D U Φ( E ( u )) A ( E ( u )) D u E ( u ) ∂ x u Since the matrix D U Φ( E ( u )) A ( E ( u )) is symmetric, we can write D U Φ( E ( u )) ∂ x F ( E ( u )) = ( A ( E ( u ))) T D U Φ( E ( u )) D u E ( u ) ∂ x u. We use (2.19) and obtain(2.20) D U Φ( E ( u )) ∂ x F ( E ( u )) = ( A ( E ( u ))) T Q T D u Φ( E ( u )) ∂ x u = S T v, and the identity (2.13) follows.The proof will be completed as soon as we establish v T S L − ( u ) S T v ≥ , which is equivalent to v T D ( u ) ≥ D ( u ) = − SU . We thus have v T D ( u ) = − v T SU . But from (2.16) and (2.17), we deduce v T D ( u ) = (cid:0) D U Φ( E ( u )) B ( E ( u )) U (cid:1) T U . As a consequence, the expected inequality v T D ( u ) ≥ D U Φ( E ( u )) D U R ( E ( u )) is non-negative.Recalling the entropy assumption (2.11) and the equilibrium property R ( E ( u )) =0, we have D U Φ( U ) R ( U ) ≥ , U ∈ Ω , ( D U Φ( U ) R ( U )) | U = E ( u ) = 0 , u ∈ ω. ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 7
As a consequence, the matrix D U ( D U Φ( U ) R ( U )) | U = E ( u ) is non-negative. Next, acalculation using the chain rule and the fact that R vanishes on the equilibriumsubmanifold, that is, R ( E ( u )) = 0 , leads us easily to D U ( D U Φ( U ) R ( U )) | U = E ( u ) = D U Φ( E ( u )) B ( E ( u )) + (cid:0) D U Φ( E ( u )) B ( E ( u )) (cid:1) T , in which no third-order derivative terms arise since R precisely vanishes on theequilibrium manifold. From the above identity, it then follows that(2.21) (cid:0) ( D U Φ( E ( u )) B ( E ( u ))) U (cid:1) T U ≥ , U ∈ Ω . We thus obtain the expected inequality v T D ( u ) ≥ (cid:3) Monotonicity of the entropy.
We then study the asymptotic behavior of theentropy inequality. In order to exhibit the entropy law satisfied by the equilibriumsolution E ( u ), we need the following technical result. Lemma 2.2.
Under the assumptions of Theorem 2.1, the entropy-flux map re-stricted to the equilibrium Ψ( E ( u )) remains constant for all u in ω .Proof. We consider the map u Ψ( E ( u )) and, after differentiation, obtain D u Ψ( E ( u )) = D U Ψ( E ( u )) D u E ( u ) . By definition of Ψ given by (2.9), we have D u Ψ( E ( u )) = D U Φ( E ( u )) A ( E ( u )) D u E ( u ) . Then, the assumption (2.12) made on Φ yields the following relation:(2.22) D u Ψ( E ( u )) = D u Φ( E ( u )) QA ( E ( u )) D u E ( u ) , = D u Φ( E ( u )) D u QF ( E ( u )) . Since we have QF ( E ( u )) = 0 over ω , then D u QF ( E ( u )) = 0. As a consequence, D u Ψ( E ( u )) = 0 for all u in ω and the proof is completed. (cid:3) Equipped with this result, we can exhibit the asymptotic equation satisfied bythe equilibrium entropy Φ( E ( u )). Arguing the formal asymptotic expansion (2.1)satisfied by U : U ε = E ( u ) + εU + . . . , where U is defined as the unique solution to (2.7), we consider the formal expansionof each term in (2.10). First, since the entropy and entropy-flux are regular maps,we have(2.23) Φ( U ε ) = Φ( E ( u )) + εD U Φ( E ( u )) U + O ( ε ) , and(2.24) Ψ( U ε ) = Ψ( E ( u )) + εD U Ψ( E ( u )) U + O ( ε ) . But, by applying Lemma 2.2, we deduce the following relation:(2.25) ∂ x Ψ( U ε ) = ε∂ x D U Ψ( E ( u )) U + O ( ε ) . Similarly, concerning the entropy relaxation source term, we easily have:(2.26) D U Φ( U ε ) R ( U ε ) = ε D U Φ( E ( u )) D U R ( E ( u )) U + O ( ε ) . C. BERTHON, P.G. L E FLOCH, AND R. TURPAULT
Now, we plug the expressions (2.23), (2.25) and (2.26) into (2.10), and consider thefirst-order terms only:(2.27) ∂ t Φ( E ( u )) = − ∂ x ( D U Ψ( E ( u )) U ) − U T (cid:0) D U Φ( E ( u )) B ( E ( u )) (cid:1) U , where, once again, U is the unique solution to (2.7). From this entropy evolutionlaw, we can state an entropy decreasing principle. Indeed, we have established thatthe matrix D U Φ( E ( u )) D U R ( E ( u )) is positive (cf. (2.21) and, as a consequence, U T (cid:0) D U Φ( E ( u )) B ( E ( u )) (cid:1) U ≥ , U ∈ Ω . We have thus proven the following statement.
Proposition 2.3.
Under the assumptions of Theorem 2.1, the entropy is non-increasing in the following sense: ∂ t Φ( E ( u )) ≤ − ∂ x ( D U Ψ( E ( u )) U ) . To conclude this presentation of the asymptotic system of diffusion equationssatisfied by stiff relaxation term for late times, let us emphasize the role playedby the entropy. As recognized by Chen, Levermore and Liu [11], the existence ofa convex mathematical entropy provides an important structure to investigate theasymptotic regime satisfied by the model. However, the main discrepancy with[11] lies in the nature of the singular limit system. Indeed, in [11], the singularlimit system turns out to be an hyperbolic system supplemented by an ε first-orderdiffusive term. Here, the obtained asymptotic system defines a system of diffusionequations. Even if the limiting solution is smooth (due to the diffusive nature ofthe limiting system), the mathematical entropy is essential in order to establish thestability of the asymptotic regime. A technical lemma.
We end this section with a technical lemma that was usefulin the above derivation of the asymptotic system.
Lemma 2.4.
Let A be a real N × N matrix such that dim(ker( A )) = n and ker( A ) ∩ Im ( A ) = { } . Let also Q be a real n × N matrix such that rank ( Q ) = n and QA = 0 Then, for any b ∈ R N , the linear system (2.28) Ax = b,Qx = 0 , admits a unique solution if and only if Qb = 0 .Proof. If the system (2.28) has a solution then left-multiplying Ax = b by Q leadsto QAx = Qb . Since QA = 0 then Qb = 0.Now we suppose that Qb = 0. Since QA = 0 then the columns of Q T are elementsof the left null-space of A . Furthermore, dim(ker( A T )) = n = rank( Q ) thereforethe columns of Q T are a basis of the left null-space of A . This implies that if z ∈ ker( A T ) then there exists y ∈ R n such that z = Q T y .Let us consider z ∈ R N such that Qz = 0. Using the fundamental theorem oflinear algebra we have z = z + z where z ∈ ker( A T ) and z ∈ Im( A ).We then characterize z and z using their respective definitions: there exists y ∈ R N such that z = Ay , and there exists y ∈ R n such that z = Q T y . With thesecharacterizations we have:0 = Qz = Qz + Qz = QQ T y + QAy = QQ T y + 0 , ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 9 since QA = 0. Now, rank( Q ) = n implies that QQ T is a symmetric positive-definitematrix and therefore that y = 0. Therefore z ∈ Im( A ) if and only if Qz = 0.By considering first the existence issue, the above property on b allows us to saythat b ∈ Im( A ) and hence there exists x ∈ R N such that Ax = b . If x is such asolution and since ker( A ) ⊕ Im( A ) = R N , we have x = x + x where x ∈ ker( A )and x ∈ Im( A ). With these notation x is a solution to (2.28).Considering next the uniqueness issue and suppose that x and x are two solu-tions to (2.28). Then ( x − x ) is the solution to A ( x − x ) = 0 ,Q ( x − x ) = 0 , and therefore ( x − x ) ∈ ker( A ) ∩ Im( A ) = { } . Hence, the solution to (2.28) isunique. (cid:3) Nonlinear diffusive regime
Derivation of a nonlinear asymptotic system.
Some physical models involveseveral relaxation time-scales. More precisely, we suppose now that the ratio of therelaxation time and the late time under consideration is no longer constant. Theextended model thus reads(3.1) ε∂ t U + ∂ x F ( U ) = − ε m R ( U ) , t > , x ∈ R , with U ∈ Ω ⊂ R N . Here, m ≥ m = 1 has beendiscussed in the previous section, and we now assume m > m = 1 are kept here. Precisey, we assumethe existence of an n × N matrix Q with rank n < N satisfying (1.2), as wellas the existence of a map E : ω → Ω satisfying (1.4). We also impose that theflux F satisfies (1.7). Concerning the nonlinear relaxation map R , an additionalassumption must be imposed: there exists a N × N matrix, denoted M ( ε ), suchthat(3.2) R ( E ( u ) + εU ) = ε m R (cid:0) E ( u ) + M ( ε ) U (cid:1) , U ∈ Ω , u ∈ ω. The matrix M ( ε ) is assumed to be sufficiently smooth in ε ∈ [0 , B ( E ( u )) is irrelevantif m >
1. Indeed, this kernel assumption was imposed to ensure the existence anduniqueness of the solution to (2.6). We are going to see that this linear system isno longer relevant and must be replaced by a nonlinear problem. In this sense, theproposed extension is called nonlinear since the diffusive asymptotic regime willinvolve a nonlinear differential operator.First, to derive the effective system of equations satisfied by the local equilibrium u ∈ ω , we introduce again a Chapman-Enskog-type expansion: U ε = E ( u ) + εU + ε U + ... We plug this expansion into (3.1) and match terms of the same order in ε .Enforcing QU ε = u , the condition (1.4) on the local equilibrium implies QU i = 0for each corrector term. Note that the expansion (2.3) for the flux remains valid, E FLOCH, AND R. TURPAULT and we only have to evaluate the expansion of the relaxation term by recalling (3.2).Indeed, this assumption gives R ( U ε ) = ε m R (cid:0) E ( u ) + M ( ε ) U + εM ( ε ) U + ... (cid:1) , and so(3.3) 1 ε m R ( U ε ) = R ( E ( u ) + M (0) U ) + εB ( M (0) U ) . ( M (0) U )+ ε D ε ( M ( ε ) U ) | ε =0 .R ( E ( u ) + M (0) U ) + O ( ε ) . Setting (2.3) and (3.3) into (3.1), we obtain(3.4) ε∂ t E ( u ) + ∂ x F ( E ( u )) + εA ( E ( u )) U = − R ( E ( u ) + M (0) U ) − εB ( M (0) U ) . ( M (0) U ) − εD ε ( M ( ε ) U ) | ε =0 .R ( E ( u ) + M (0) U ) + O ( ε ) . Considering the zeroth-order terms, we get ∂ x F ( E ( u )) = − R (cid:0) E ( u ) + M (0) U (cid:1) , which turns out to be a nonlinear system of equations with unknown U , supple-mented by the condition QU = 0.At this level, we see that the assumption on the kernel of B ( E ( u )) is no longerrelevant (or sufficient) in order to solve (2.6), since we now have to consider(3.5) R ( E ( u ) + M (0) U ) = − ∂ x F ( E ( u )) ,QU = 0 . In view of the strong nonlinearities involved in this equation, we tacitly assume theexistence and uniqueness of the solution, denoted ¯ U , of (3.5). In the applications,this property will be checked directly.Then, we match first-order terms issuing from (3.4) to get ∂ t E ( u ) + ∂ x A ( E ( u )) . ¯ U = − B ( M (0) ¯ U ) . ( M (0) U ) − D ε ( M ( ε ) ¯ U ) | ε =0 .R ( E ( u ) + M (0) ¯ U ) . Since Q E ( u ) = u for all u ∈ ω , and since QR ( U ) = 0 and QB ( U ) = 0 for all U ∈ Ω,by multiplying the above relation by Q , we obtain(3.6) ∂ t u = − ∂ x (cid:0) QA ( E ( u )) ¯ U (cid:1) . Once again, this equilibrium equation involves a nonlinear differential operatorin the right-hand side since ¯ U is a nonlinear map applied to first-order spacederivatives.Similarly to the “linear” case governed by (2.8), we want interpret (3.6) as asystem of diffusion equations. We thus assume the existence of a convex entropyΦ : Ω → R which satisfies all the compatibility conditions imposed in the previoussection. The matrix D U Φ( U ) D U F ( U ) to be symmetric for all U ∈ Ω and, in addi-tion, we assume that the compatibility conditions (2.11) and (2.12) hold. Smoothsolutions to (3.1) satisfy the additional balance law(3.7) ε∂ t Φ( U ε ) + ∂ x Ψ( U ε ) = − ε m D U Φ( U ε ) R ( U ε ) . ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 11
Equipped with this convex entropy, we observe that ¯ U is, equivalently, a solutionto the nonlinear algebraic system D U Φ( E ( u )) R ( E ( u ) + M (0) U ) = − D U Φ( E ( u )) ∂ x F ( E ( u )) ,QU = 0 , for any convex entropy compatible with the relaxation term.Next, consider the matrix S defined by (2.14) and the vector v given by (2.18).In view of (2.20), the above system is equivalent to(3.8) D U Φ( E ( u )) R ( E ( u ) + M (0) U ) = − S T v,QU = 0 . With some abuse of notation and for the sake of clarity, we set N u ( U ) = D U Φ( E ( u )) R ( E ( u ) + M (0) U )and introduce the notation(3.9) ¯ U = N − u ( − S T v ) . Hence, the equilibrium equation (3.6) reads as follows:(3.10) ∂ t u = ∂ x (cid:0) − S N − u ( − S T v ) (cid:1) . Once again, we note the crucial role played by the convex entropies. Indeed, wenow exhibit the limiting system of equations satisfied by the equilibrium entropyΦ( E ( u )). We skip here the details of the computation which are similar to the linearcase. Lemma 2.2 still holds, as well as the entropy expansion (2.23) and (2.25). Infact, only the entropy relaxation source term expansion changes. We now obtain D U Φ( U ε ) R ( U ε ) = ε m +1 ¯ U T D U Φ( E ( u )) R ( E ( u ) + M (0) ¯ U ) . Plugging these expansions into (3.7) and considering first-order terms, we find(3.11) ∂ t Φ( E ( u )) = − ∂ x (cid:0) D U Ψ( E ( u )) ¯ U (cid:1) − ¯ U T D U Φ( E ( u )) R ( E ( u ) + M (0) ¯ U ) . To conclude this section, we show that the asymptotic system of equations (3.10) isof diffusive-type, and that an associated mathematical entropy is non-decreasing.
Lemma 3.1.
Let ¯ U be given by (3.5). Assume the existence of a non-negativemap c ( u ) ≥ such that (3.12) R ( E ( u ) + M (0) ¯ U ) = c ( u ) ¯ U , u ∈ ω. Then the limiting equation (3.6) is nonlinearly dissipative with respect to the entropyin the following sense: (3.13) − v T S N − u ( − S T v ) ≥ , u ∈ ω, where S is given by (2.14) and v by (2.18).Moreover, the entropy is decreasing as follows: (3.14) ∂ t Φ( E ( u )) ≤ − ∂ x (cid:0) D U Ψ( E ( u )) . ¯ U (cid:1) . Proof.
Arguing (3.8) and the definition (3.9), we have − v T S N − u ( − S T v ) = ¯ U T D U Φ( E ( u )) R ( E ( u ) + M (0) ¯ U ) . By recalling (3.12), we obtain − v T S N − u ( − S T v ) = c ( u ) ¯ U T D U Φ( E ( u )) ¯ U . E FLOCH, AND R. TURPAULT
The inequality (3.13) follows from the convexity of the entropy. Recalling (3.11)we obtain (3.14), and the proof is completed. (cid:3) Physical examples
Euler equations with friction term.
Many models involving distinct physicalscales enter the framework proposed in the present paper and, specifically, we willnow illustrate the interest of our framework with four examples. We begin withthe isentropic Euler equations supplemented with a friction term (cf. [27, 33] forfurther details). The asymptotics for this model has been already considered in theliterarure and, more recently, relevant numerical techniques have been proposed[5, 9]. (See also [26, 3].) Importantly, this model satisfies all of the conditionsrequired in Section 3, above.The Euler model with friction reads(4.1) ε∂ t ρ + ∂ x ρv = 0 ,ε∂ t ρv + ∂ x ( ρv + p ( ρ )) = − ε ρv, where ρ > v ∈ R the velocity of a compressible fluid.The pressure function p : R + → R + is assumed to be sufficiently regular andsatisfy p ′ ( ρ ) >
0, so that the first-order homogeneous system associated with (4.1)is strictly hyperbolic.In view of (1.1), we should set(4.2) U = (cid:18) ρρv (cid:19) , F ( U ) = (cid:18) ρvρv + p ( ρ ) (cid:19) , R ( U ) = (cid:18) ρv (cid:19) , which corresponds to the matrix Q = (1 0) , and scalar local equilibria u = ρ , with E ( u ) = (cid:18) ρ (cid:19) . As required, we also have QF ( E ( u )) = 0, and it is easily checked that all ourassumptions of the previous sections hold.Considering now the asymptotic diffusive regime in the limit ε →
0, we note thatthe equilibrium solution must satisfy (2.8), i.e. ∂ t ρ = − ∂ x ( QA ( E ( u )) U ) , where A ( E ( u )) = (cid:18) p ′ ( ρ ) 0 (cid:19) , and U is the unique solution to (2.7). Since B ( E ( u )) = (cid:18) (cid:19) , ∂ x F ( E ( u )) = (cid:18) ∂ x p ( ρ ) (cid:19) , the diffusive regime associated with this Euler model with friction is described bythe equation(4.3) ∂ t ρ = ∂ x (cid:0) p ( ρ ) (cid:1) . Based on Theorem 2.1, we observe that the diffusive nature of (4.3) follows fromthe existence of a convex entropy which is compatible with the relaxation source
ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 13 term in the sense (2.11)-(2.12). Indeed, by introducing the internal energy e ( ρ ) > e ′ ( ρ ) = p ( ρ ) ρ , we see that smooth solutions to (4.1) satisfy ε∂ t (cid:18) ρ v ρe ( ρ ) (cid:19) + ∂ x (cid:18) ρ v ρe ( ρ ) + p ( ρ ) (cid:19) v = − ε ρv . The function Φ( U ) = ρ v + ρe ( ρ ) is a convex entropy satisfying the compatibilityconditions with the relaxation terms. The M model for radiative transfer. The second example of interest relieson a more complex physical set-up, relevant in radiative transfer and referred to asthe M ε∂ t e + ∂ x f = 1 ε ( τ − e ) ,ε∂ t f + ∂ x (cid:16) χ ( f /e ) e (cid:17) = − ε f,ε∂ t τ = 1 ε ( e − τ ) , where e > f the radiative flux, restricted by the “fluxlimitation” condition (cid:12)(cid:12)(cid:12)(cid:12) fe (cid:12)(cid:12)(cid:12)(cid:12) ≤ , and τ > χ : [ − , → R + stands for theEddington factor defined by χ ( ξ ) = 3 + 4 ξ p − ξ . Once again, we rely on the notation introduced in the previous sections (cf. sys-tem (1.1)) and we write U = efτ , F ( U ) = fχ ( fe ) e , R ( U ) = e − τ fτ − e . The local equilibria are described by the map E ( u ) = τ τ , where u = τ + τ is now a scalar. We set Q = (1 0 1) and, in agreement with (1.7),we have QF ( E ( u )) = 0.The asymptotic regime is governed by (2.8) and to exhibit its explicit formulationwe need the expression of A ( E ( u )) and U . A straightforward calculation gives A ( E ( u )) = χ (0) χ ′ (0) 00 0 0 = , E FLOCH, AND R. TURPAULT while U is the solution to (2.7) which here reads − τ − τ U = ∂ x ( τ )0 , (1 0 1) U = 0 . We thus obtain U = τ ∂ x τ , and, according with (2.8), the asymptotic regime is governed by the following dif-fusion equation: ∂ t ( τ + τ ) = ∂ x (cid:18) τ ∂ x τ (cid:19) . A coupled Euler/ M model. We propose here an example of system that de-generates into a system of diffusion equations of dimension n >
1. To do so, wecouple the Euler model (4.1) with the M ε∂ t ρ + ∂ x ρv = 0 ,ε∂ t ρv + ∂ x ( ρv + p ( ρ )) = − κε ρv + σε f,ε∂ t e + ∂ x f = 0 ,ε∂ t f + ∂ x χ (cid:18) fe (cid:19) e = − σε f, in the notation previously introduced. Here, κ and σ denote positive constants.Even though this is a toy model, the pressure has to be sufficiently small in orderto represent an application of physical interest. We will therefore consider thefollowing pressure law: p ( ρ ) = C p ρ η , C p ≪ , η > . In the formalism (1.1) we need to set U = ρρvef , F ( U ) = ρvρv + p ( ρ ) fχ ( fe ) e , R ( U ) = κρv − σf σf . The local equilibrium is given by E ( u ) = ρ e , u = QU = (cid:18) ρe (cid:19) , Q = (cid:18) (cid:19) . Once again, one has QF ( E ( u )) = 0. To derive the asymptotic regime, we exhibit A ( E ( u )) and U : A ( E ( u )) = p ′ ( ρ ) 0 0 00 0 0 10 0 , U = κ (cid:2) − ∂ x p ( ρ ) − ∂ x e (cid:3) − σ ∂ x e . ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 15
The asymptotic diffusive regime of the system (4.5) is therefore given by(4.6) ∂ t ρ − κ ∂ x p ( ρ ) − κ ∂ x e = 0 ,∂ t e − σ ∂ x e = 0 . Shallow water with strong friction effect.
The last suggested example is de-voted to the well-known shallow-water model supplemented by a strong frictionterm in a late-time regime where the friction is assumed to dominate over theconvection.This model is given as follows:(4.7) ε∂ t h + ∂ x hv = 0 ,ε∂ t hv + ∂ x (cid:0) h v + p ( h ) (cid:1) = − κ ( h ) ε g hv | hv | , where h > v ∈ R the velocity and p ( h ) = g h the pressure law.Here, g > κ : R + → R + is a given and positive function. In [28], several examples of friction κ are proposed.A standard choice is κ ( h ) = κ h where κ > m = 2. According to the notation involved in (3.1), we haveset U = (cid:18) hhv (cid:19) , F ( U ) = (cid:18) hvhv + p ( h ) (cid:19) , R ( U ) = (cid:18) κ ( h ) ghv | hv | (cid:19) . Concerning the equilibrium, we have E ( u ) = (cid:18) h (cid:19) , where u = h is a scalar and Q = (1 0). The assumption (3.2) easily holds since R ( E ( u ) + εU ) = ε R (cid:0) E ( U ) + M ( ε ) U (cid:1) , where M ( ε ) = (cid:18) ε
00 1 (cid:19) . We turn our attention now to the asymptotic regime which is governed by thenonlinear diffusion equation (3.6). To get its explicit from, we have to exhibit U = ( α β ) T the solution to (3.5) which reads as follows: α = 0 ,κ ( h ) gβ | β | = − ∂ x p ( h ) . We easily find β = − √ h∂ x hκ ( h ) p | ∂ x h | , and the effective nonlinear diffusion equation is thus(4.8) ∂ t h = ∂ x √ hκ ( h ) ∂ x h p | ∂ x h | ! . This is a nonlinear Laplacian equation (for instance, see [23] and references therein). E FLOCH, AND R. TURPAULT
Lemma 3.1 applies here with the following entropy. By introducing the internalenergy e ( h ) = gh/
2, we see that smooth solutions to (4.7) satisfy ε∂ t (cid:18) h v g h (cid:19) + ∂ x (cid:18) h v gh (cid:19) v = − κ ( h ) ε ghv | hv | . The entropy Φ( U ) = h v + g h satisfies all the required properties, and the condi-tion (3.12) holds since R ( E ( u ) + M (0) ¯ U ) = (cid:18) ∂ x p ( h ) (cid:19) , where ¯ U = (0 β ) T . As a consequence, we obtain R ( E ( u ) + M (0) ¯ U ) = c ( u ) ¯ U with c ( u ) = gκ ( h ) p h | ∂ x h | ≥ . Hence, the limit equation (4.8) is of diffusive-type in the sense of Lemma 3.1.5.
Asymptotic-preserving schemes
Objective.
In this section, we consider the numerical approximation of solutionsto (1.1). Our goal is to derive a class of numerical schemes that restore the relevantasymptotic regime, given by (2.8), in the limit ε →
0. One of the main difficultieswhen deriving asymptotic-preserving schemes lies in the independent role playedby each ε and the mesh size. More precisely, the limit discrete diffusion equation(as ε tends to zero) must be obtained independently of the space mesh-size.Such a numerical problem was investigated during the last decade on severalspecific examples. For instance, in [5], the Euler equations with friction term areconsidered. Relating works to radiative transfer and the M M ε :(5.1) ∂ t U + ∂ x F ( U ) = 0 . Next, we derive a suitable correction to obtain a finite volume discretization ofthe source term. Hence, the corrected finite volume method gives approximatesolutions of the following system:(5.2) ∂ t U + ∂ x F ( U ) = − γR ( U ) , where γ > γ = ε .The asymptotic scheme is thus obtained in the limit of ε to zero. Modulo a suitablecorrection, it is thus proved to be asymptotic preserving. A discretization of (5.2).
As a first step, we suggest to consider the well-knownGodunov-type scheme introduced by Harten, Lax and van Leer [21] with a singleconstant intermediate state, to approximate the weak solutions to (5.1).We consider a uniform mesh made of cells [ x i − / , x i +1 / ) where x i +1 / = x i + ∆ x for all i ∈ Z with a constant cell size ∆ x . The time discretization is defined by ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 17 t n +1 = t n + ∆ t where the time increment will be restricted later on by a CFL likecondition.We define the discrete initial data as follows: U ( x ) = 1∆ x Z x i +1 / x i − / U ( x, dx, x ∈ [ x i − / , x i +1 / ) . We seek for a piecewise constant approximation of the exact solution of (5.1) at thetime t n , U n ( x ) = U ni , x ∈ [ x i − / , x i +1 / ) , with U ni ∈ Ω for all i ∈ Z .By considering a suitable sequence of approximate Riemann solvers, we canevolve this approximation and get a piecewise constant function ˜ U n +1 ( x ) which isan approximation of the solution to (5.1) at the time t n +∆ t . Following Harten, Laxand van Leer [21], at each cell interface we use the following approximate Riemannsolver:(5.3) ˜ U R ( xt ; U L , U R ) = U L , xt < − b, ˜ U ⋆ , − b < xt < b,U R , xt > b, where b > U ⋆ = 12 ( U L + U R ) − b ( F ( U R ) − F ( U L )) . As a consequence, as soon as the following CFL restriction holds:(5.5) b ∆ t ∆ x ≤ , we are considering a juxtaposition of non-interacting approximate Riemann solver(cf. Figure 1) denoted ˜ U n ∆ x ( x, t n + t ) for t ∈ [0 , ∆ t ). Figure 1.
HLL scheme: Juxtaposition of approximated Riemann problems. U ni − U ni U ni +1 ˜ U ⋆i − / ˜ U ⋆i +1 / − b b − b bx i − / x i +1 / xt n t The updated approximated solution at time t n +1 is thus defined as follows:(5.6) ˜ U n +1 i = 1∆ x Z x i +1 / x i − / ˜ U n ∆ x ( x, t n + ∆ t ) dx. Since we have(5.7) ˜ U ⋆i +1 / = 12 ( U ni + U ni +1 ) − b ( F ( U ni +1 ) − F ( U ni )) , E FLOCH, AND R. TURPAULT an easy computation gives the following standard conservation form:(5.8) ˜ U n +1 i = U ni − ∆ t ∆ x ( F HLLi +1 / − F HLLi − / ) , where we have set(5.9) F HLLi +1 / = 12 ( F ( U ni ) + F ( U ni +1 )) − b U ni +1 − U ni ) . For simplicity in the presentation, we have adopted here a constant numericalcone of dependence (cf. Figure 1) characterized by a single speed parameter b > b − i +1 / , b + i +1 / ) with b − i +1 / < b + i +1 / . However, for the sakeof simplicity and without genuine loss of generality, we present our strategy for thesimpler case b + i +1 / = − b − i +1 / = b > U ⋆i +1 / belongs to Ω for all i ∈ Z to deduce that U n +1 i in Ω for all i ∈ Z . Indeed, from (5.6), we have˜ U n +1 i = b ∆ t ∆ x ˜ U ⋆i − / + (cid:18) − b ∆ t ∆ x (cid:19) U ni + b ∆ t ∆ x ˜ U ⋆i +1 / and, based on the CFL restriction (5.5), the above relation is a convex combinationof states in Ω. Since Ω is a convex set, we deduce that ˜ U n +1 i belongs to Ω.Our main necessary condition for the above argument is that ˜ U ⋆i +1 / ∈ Ω for all i in Z . But, once again, ˜ U ⋆i +1 / can be seen as a convex combination, as follows:˜ U ⋆i +1 / = 12 (cid:18) U ni + 1 b F ( U ni ) (cid:19) + 12 (cid:18) U ni +1 − b F ( U ni +1 ) (cid:19) . Since Ω is an open convex set, we can choose b to be large enough so that to enforcethe condition ˜ U ⋆i +1 / ∈ Ω.Next, we modify this approximate Riemann solver and introduce a discretiza-tion of the source-term in order to approximate the solutions of (5.2). Similarmodifications were made for specific problems in [3] and [2, 6], while we proposehere a general approach based on matrix-valued free-parameters. We modify theapproximate Riemann solver (5.3) as follows:(5.10) U R ( xt ; U L , U R ) = U L , xt < − b,U ⋆L , − b < xt < ,U ⋆R , < xt < b,U R , xt > b, where we have set(5.11) U ⋆L = α ˜ U ⋆ + ( I − α )( U L − ¯ R ( U L )) ,U ⋆R = α ˜ U ⋆ + ( I − α )( U R − ¯ R ( U R )) . Here, α denotes a N × N matrix defined as follows:(5.12) α = (cid:18) I + γ ∆ x b ( I + σ ) (cid:19) − , ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 19 and(5.13) ¯ R ( U ) = ( I + σ ) − R ( U ) . The N × N matrices I and σ respectively denote the identity matrix and a parametermatrix to be defined.A choice for the matrix σ will be made later, which will turn out to govern theasymptotic regime. At this point, the matrix σ is assumed to be such that theinverse matrices in (5.12) and (5.13) are well-defined.We adopt the modified approximate Riemann solver (5.10) to derive a modifiedGodunov type scheme in the spirit of [3]. At each cell interface x i +1 / , we set theapproximate Riemann solver U R ( x − x i +1 / t − t n ; U ni , U ni +1 ) to define a juxtaposition ofmodified approximate Riemann solver, denoted by U n ∆ x ( x, t n + t ) for t ∈ [0 , ∆ t ).(See Figure 2.) Figure 2.
Juxtaposition of modified approximated Riemann problems. U ni − U ni U ni +1 ˜ U ⋆Li − / ˜ U ⋆Ri − / ˜ U ⋆Li +1 / ˜ U ⋆Ri +1 / − b b − b b x i − / x i +1 / xt n t Thanks to the CFL condition (5.5), such a juxtaposition is non-interacting.At time t n +1 , the updated approximated solution is given as follows for all i ∈ Z :(5.14) U n +1 i = Z x i +1 / x i − / U n ∆ x ( x, t n + ∆ t ) dx. We compute this integral form and obtain(5.15) 1∆ t ( U n +1 i − U ni ) + 1∆ x ( α i +1 / F HLLi +1 / − α i − / F HLLi − / )= 1∆ x ( α i +1 / − α i − / ) F ( U ni ) − b ∆ x ( I − α i − / ) ¯ R i − / ( U ni ) − b ∆ x ( I − α i +1 / ) ¯ R i +1 / ( U ni ) , where the numerical flux F HLLi +1 / is given by (5.9). The discretized source-term canbe rewritten in a more relevant form: b ∆ x ( I − α i +1 / ) ¯ R i +1 / ( U ni ) = b ∆ x α i +1 / ( α − i +1 / − I ) ¯ R i +1 / ( U ni ) . In view of the definition of α i +1 / (see (5.12)), we deduce that b ∆ x ( I − α i +1 / ) ¯ R i +1 / ( U ni ) = γ α i +1 / R ( U ni )and, similarly, b ∆ x ( I − α i − / ) ¯ R i − / ( U ni ) = γ α i − / R ( U ni ) . E FLOCH, AND R. TURPAULT
As a consequence, the scheme (5.15) reads(5.16) 1∆ t ( U n +1 i − U ni ) + 1∆ x ( α i +1 / F HLLi +1 / − α i − / F HLLi − / )= 1∆ x ( α i +1 / − α i − / ) F ( U ni ) − γ α i +1 / + α i − / ) R ( U ni ) . This scheme satisfies the following statement.
Theorem 5.1.
Assume that the matrix σ defines a spatially continuous map. Thenumerical scheme (5.16) is consistant with the equation (5.2).At time t n , assume U ni ∈ Ω for all i ∈ Z and, in addition, that the state vectors U ⋆Li +1 / and U ⋆Ri +1 / , defined by U ⋆Li +1 / = α i +1 / ˜ U ⋆i +1 / + ( I − α i +1 / )( U ni − ¯ R ( U ni )) ,U ⋆Ri +1 / = α i +1 / ˜ U ⋆i +1 / + ( I − α i +1 / )( U ni +1 − ¯ R ( U ni +1 )) , belong to Ω . Then, the updated state vector U n +1 i , defined by (5.16), belongs theset Ω for all i in Z .Proof. The consistency property follows from the definition of α i +1 / , given by(5.12). Indeed, we easily obtain α i +1 / = I + O (∆ x ) , to deduce the expected consistency of the flux and the relaxation source term. Onlythe term x ( α i +1 / − α i − / ) F ( U ni ) is left over. Recalling (5.12), we have1∆ x ( α i +1 / − α i − / ) F ( U ni ) = − γ b α i +1 / ( σ i +1 / − σ i − / ) α i − / F ( U ni ) , to obtain 1∆ x ( α i +1 / − α i − / ) F ( U ni ) = O (∆ x ) , as soon as σ i +1 / − σ i − / = O (∆ x ). The expected equation consistency is thereforeobtained.Concerning the robustness of the method, from (5.14), we have U n +1 i = b ∆ t ∆ x U ⋆Li − / + (cid:18) − b ∆ t ∆ x (cid:19) U ni + b ∆ t ∆ x U ⋆Ri +1 / , which is nothing but a convex combination in Ω. Then U n +1 i is in Ω and the proofis completed. (cid:3) To conclude this derivation, observe that the term x ( α i +1 / − α i − / ) F ( U ni )may seem to be a discrepancy in the method. In fact, this term is standard toderive asymptotic preserving schemes and it can be found in several works. (See,for instance, [2, 3, 6].) Our approach allows us recover a scheme proposed earlierin [18]. ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 21
The linear asymptotic regime.
The scheme (5.16) is now considered to suggesta discretization of our initial model (1.1). The expected scheme is thus easilyobtained when substituting ∆ t by ∆ tε and fixing γ = ε . The resulting scheme readsas follows:(5.17) ε ∆ t ( U n +1 i − U ni ) + 1∆ x ( α i +1 / F HLLi +1 / − α i − / F HLLi − / )= 1∆ x ( α i +1 / − α i − / ) F ( U ni ) − ε ( α i +1 / + α i − / ) R ( U ni ) , where(5.18) α i +1 / = (cid:18) I + ∆ x εb ( I + σ i +1 / ) (cid:19) − . For the sake of simplicity in the forthcoming asymptotic derivation, we propose tointroduce the N × N matrix α εi +1 / defined by(5.19) α εi +1 / = (cid:18) εI + ∆ x b ( I + σ i +1 / ) (cid:19) − , so that we have α i +1 / = εα εi +1 / . Recalling the definition, the scheme (5.16) takesthe form:(5.20) ε ∆ t ( U n +1 i − U ni ) + ε ∆ x ( α εi +1 / F HLLi +1 / − α εi − / F HLLi − / )= ε ∆ x ( α εi +1 / − α εi − / ) F ( U ni ) −
12 ( α εi +1 / + α εi − / ) R ( U ni ) . We observe that U ni remains close to the equilibrium state E ( u ni ) for ε small andwe thus consider the following expansion: U ni = E ( u ni ) + ε ( U ) ni + O ( ε ) , which we now plug into (5.20). We easily have1 ε R ( U ni ) = B ( E ( u ni )) . ( U ) ni + O ( ε ) ,F HLLi +1 / | E ( u )+ O ( ε ) = F HLLi +1 / | E ( u ) + O ( ε ) , where F HLLi +1 / | E ( u ) = 12 (cid:0) F ( E ( u ni )) + F ( E ( u ni +1 )) (cid:1) − b (cid:0) E ( u ni +1 ) − E ( u ni ) (cid:1) . In addition we have α εi +1 / = 2 b ∆ x ( I + σ i +1 / ) − + O ( ε ) . By considering the first-order terms issuing from (5.20), we obtain1∆ t ( E ( u n +1 i ) − E ( u ni ))+ 2 b ∆ x (cid:16) ( I + σ i +1 / ) − F HLLi +1 / | E ( u ) − ( I + σ i − / ) − F HLLi − / | E ( u ) (cid:17) = 2 b ∆ x (cid:16) ( I + σ i +1 / ) − − ( I + σ i − / ) − (cid:17) F ( E ( u ni )) − b ∆ x (cid:16) ( I + σ i +1 / ) − + ( I + σ i − / ) − (cid:17) B ( E ( u ni )) . ( U ) ni . E FLOCH, AND R. TURPAULT
Next, assume the existence of a n × n squared matrix, denoted by M i +1 / , suchthat(5.21) Q ( I + σ i +1 / ) − = 1 b M i +1 / Q. Recalling the assumptions (1.2), (1.4), and (1.7), we write1∆ t ( u n +1 i − u ni ) = − b ∆ x (cid:16) M i +1 / QF HLLi +1 / | E ( u ) − M i − / QF HLLi − / | E ( u ) (cid:17) , where QF HLLi +1 / | E ( u ) = 12 Q (cid:0) F ( E ( u ni )) + F ( E ( u ni +1 )) (cid:1) − b Q (cid:0) E ( u ni +1 ) − E ( u ni ) (cid:1) = − b u ni +1 − u ni ) . As a consequence, the asymptotic discrete regime is given by(5.22) 1∆ t ( u n +1 i − u ni ) = 1∆ x (cid:0) M i +1 / ( u ni +1 − u ni ) + M i − / ( u ni − − u ni ) (cid:1) . We thus have established the following result.
Theorem 5.2.
Consider any N × N matrix σ i +1 / such that the matrices I + σ i +1 / and (1 + ∆ x εb ) I + σ i +1 / are nonsingular for all ε > . Assume the existence of a n × n matrix M i +1 / such that (5.21) holds, and introduce the n × n matrix M ( u ) defined by M ( u ) = QA ( E ( u )) L − ( u ) D U Φ( E ( u )) A ( E ( u )) D u E ( u ) . In addition, assume that the matrix M i +1 / is a discrete form of M ( u ) at eachcell interface x i +1 / . Then, the asymptotic behavior of the scheme (5.20) coincideswith a discrete form for the limit diffusion equation (2.13).Proof. We directly deduce from (2.13) that the diffusive limit equation reads ∂ t u = ∂ x ( M ( u ) ∂ x u ) . Since the asymptotic regime satisfied by the scheme (5.20) is governed by (5.22),the proposed choice of the matrix M i +1 / leads us to the correct behavior of thescheme as ε goes to zero. (cid:3) The nonlinear asymptotic regime.
We propose to extend the above numericalscheme to consider the nonlinear asymptotic regime governed by the system (3.1).To address such an issue, once again we consider the scheme (5.16) where ∆ t issubstituted by ∆ tε and we set γ = ε . To be consistent, we substitute R ( U ni ) by ε m − R ( U ni ).Adopting such a strategy, the same arguments used to obtain (5.20) now give(5.23) 1∆ t ( U n +1 i − U ni ) + 1∆ x ( α εi +1 / F HLLi +1 / − α εi − / F HLLi − / )= 1∆ x ( α εi +1 / − α εi − / ) F ( U ni ) − ε m ( α εi +1 / + α εi − / ) R ( U ni ) . where α εi +1 / is defined by (5.19). The above scheme exactly coincides with (5.20)as soon as we fix m = 1. ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 23
Once again, as soon as ε tends to zero, the state vector U ni remains in a neigh-borhood of the equilibrium state E ( u ni ). As a consequence, we adopt a formalexpansion given by U ni = E ( u ni ) + ε ( U ) ni + O ( ε ) . Arguing the same calculations as used in the linear expansion case, the asymptoticdiscrete equation is given by (5.22). Hence, we have to choose M i +1 / in order toget a discretization of (3.6). Theorem 5.3.
Consider an N × N matrix σ i +1 / such that the matrices I + σ i +1 / and (1 + ∆ x εb ) I + σ i +1 / are non-singular for all ε > . Assume the existence of a n × n matrix M i +1 / such that (5.21) holds. Consider also the n × n matrix M ( u ) defined by (5.24) M ( u ) = QA ( E ( u )) R − ( u ) , where R − : ω → Ω defines the unique solution of (3.5). Assume that M i +1 / isa discrete form of M ( u ) at each cell interface x i +1 / . The asymptotic behavior ofthe scheme (5.23) defines a discrete form of the nonlinear diffusion equation (3.6).Proof. By definition of R − ( u ), we have ¯ U = R − ( u ) the unique solution to (3.5).As a consequence, we deduce the following rewriting of (3.6): ∂ t u = ∂ x ( M ( u ) ∂ x u ) , where M ( u ) is given by (5.24). The proposed definition of M i +1 / concludes theproof. (cid:3) Numerical experiments
Euler equations with friction.
To illustrate the interest of the asymptotic pre-serving scheme (5.17), we apply it now to the Euler equations with a friction term(4.1). Here, we suggest to fix the matrix parameter σ i +1 / as follows: σ i +1 / = σ i +1 / I, where σ i +1 / stands for a scalar parameter to be defined. As a consequence, thematrix α i +1 / , defined by (5.18), is now given by(6.1) α i +1 / = α i +1 / I, α i +1 / = 11 + ∆ x εb (1 + σ i +1 / ) . The scheme (5.17) thus reads ε ∆ t ( ρ n +1 i − ρ ni ) + 1∆ x (cid:16) α i +1 / F ρ,HLLi +1 / − α i − / F ρ,HLLi − / (cid:17) = − α i +1 / σ i +1 / − σ i − / εb α i − / ( ρv ) ni , (6.2) ε ∆ t (( ρv ) n +1 i − ( ρv ) ni ) + 1∆ x (cid:16) α i +1 / F ρv,HLLi +1 / − α i − / F ρv,HLLi − / (cid:17) = − α i +1 / σ i +1 / − σ i − / εb α i − / (cid:0) ρ ni ( v ni ) + p ( ρ ni ) (cid:1) − ε α i +1 / + α i − / ρv ) ni , (6.3)where the numerical flux vector ( F ρ,HLLi +1 / , F ρv,HLLi +1 / ) is defined by (5.9). E FLOCH, AND R. TURPAULT
First, applying Theorem 5.1, we observe that the proposed scheme is consistentwith the system (4.1) and preserves the admissible states. To address such an issue,in view of Theorem 5.1, we have to establish the positivity property ( i ∈ Z ) ρ ⋆Li +1 / = α i +1 / ˜ ρ ⋆i +1 / + (1 − α i +1 / ) ρ ni ,ρ ⋆Ri +1 / = α i +1 / ˜ ρ ⋆i +1 / + (1 − α i +1 / ) ρ ni +1 . Since α i +1 / , defined by (6.1), belongs to (0 , ρ ⋆Li +1 / > ρ ⋆Ri +1 / > ρ ⋆i +1 / > b (cf. relation(5.4)).Now, we study for the asymptotic behavior of the scheme (6.2)-(6.3). Put inother words, we have to fix the free parameter σ i +1 / to recover the expectedasymptotic regime in the limit of ε to zero. This required asymptotic behaviormust be governed by a discrete form of (4.3). First, observe thatlim ε → α i +1 / = 0 , lim ε → α i +1 / ε = 0From the momentum ( ρv ) n +1 i governing equation (6.2), we easily deduce the fol-lowing momentum behavior in the limit of ε to zero:( ρv ) ni = 0 , i ∈ Z . The density approximation thus admits the following asymptotic regime:1∆ t ( ρ n +1 i − ρ ni ) + 2 bδx (cid:18)
11 + σ i +1 / F ρ,HLLi +1 / | ρv =0 −
11 + σ i − / F ρ,HLLi − / | ρv =0 (cid:19) = 0 . But we have F ρ,HLLi − / | ρv =0 = − b (cid:0) ρ ni +1 − ρ ni (cid:1) , and therefore1∆ t ( ρ n +1 i − ρ ni ) = b δx (cid:18)
11 + σ i +1 / ( ρ ni +1 − ρ ni ) + 11 + σ i − / ( ρ ni − − ρ ni ) (cid:19) . We choose(6.4) σ i +1 / = b ρ ni +1 − ρ ni p ( ρ ni +1 ) − p ( ρ ni ) − , if ρ ni +1 = ρ ni ,b p ′ ( ρ ni ) − , otherwise , and arrive at the following discretization of the diffusion equation (4.3):1∆ t ( ρ n +1 i − ρ ni ) = 1 δx (cid:0) p ( ρ ni +1 ) − p ( ρ ni ) + p ( ρ ni − ) (cid:1) . To conclude, observe that Theorem 5.2 applies if the matrix M i +1 / is defined by M i +1 / = m i +1 / I n where I n denotes the n × n identity matrix and m i +1 / = p ( ρ ni +1 ) − p ( ρ ni ) ρ ni +1 − ρ ni , if ρ ni +1 = ρ ni ,p ′ ( ρ ni ) , otherwise . In [3, 9], the scheme (6.2)-(6.3) was derived by a completely different approach.Similarly, in [3, 2, 6], the application of our general asymptotic preserving scheme(5.17) is found in the framework of the radiative transfer (4.4).
ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 25
Figure 3.
Reference (full line) and proposed scheme (dashed line)solution comparison.
Figure 4.
Decrease of the entropy.As an illustration, the scheme (6.2)-(6.3) supplemented by the asymptotic pre-serving correction (6.4) is used to approximate the solution when the initial datais given by ( ρ, ρv ) T = ( (2 , T , x ∈ [1 . , . , (1 , T , otherwise . Furthermore, we choose the simple pressure law p ( ρ ) = ρ and ε = 10 − . InFigure 3, the solution on a 300 cells grid (∆ x = 10 − ) is compared at time t = E FLOCH, AND R. TURPAULT . with a numerical approximation of (4.3) computed with 600 cells. We notea fairly good agreement between the two approximate results. In agreement withLemma 2.3, Figure 4 shows that entropy is decreasing in time. For two choices ofspace grids, we plot P i (cid:16) ρ v + ρe ( ρ ) (cid:17) i versus time. Coupled Euler/ M equations. We now propose to approximate the solutionof the system (4.5) by adopting the asymptotic-preserving scheme (5.17) with thefollowing matrix parameter σ i +1 / : σ i +1 / = σ ,i +1 / − σ ,i +1 /
00 0 0 00 0 σ ,i +1 /
00 0 0 0 , where σ j,i +1 / are parameters to be defined later, in order to reach the requiredasymptotic regime (4.6). With such a definition α εi +1 / , defined by (5.19), becomes: α εk = b k εθ ,k b k εγ ∆ xσ ,k θ ,k θ ,k
00 2 b k ε b k ε + γ ∆ x b k εθ ,k
00 0 0 2 b k ε b k ε + γ ∆ x , where θ j,k = 2 b k ε + γ ∆ x (1 + σ j,k ) and γ = max( κ,σ ) ε . The asymptotic regimeassociated with this scheme is ρ n +1 i = ρ ni − ∆ t ∆ x b γ ∆ x ρ ni +1 − ρ ni σ ,i +1 / − ρ ni − ρ ni − σ ,i − / + σ ,i +1 / ( e ni +1 − e ni )(1 + σ ,i +1 / )(1 + σ ,i +1 / ) − σ ,i − / ( e ni − e ni − )(1 + σ ,i − / )(1 + σ ,i − / ) ! ,e n +1 i = e ni − ∆ t ∆ x b γ ∆ x e ni +1 − e ni σ ,i +1 / − e ni − e ni − σ ,i − / ! . To recover a discrete form of (4.6), we propose to set σ ,i +1 / = κγ ρ ni +1 − ρ ni p ni +1 − p ni − ,σ ,i +1 / = σγ ρ ni +1 − ρ ni p ni +1 − p ni ,σ ,i +1 / = 3 σγ − . ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 27
In conclusion, setting ∆ p ni +1 / = p ni +1 − p ni ρ ni +1 − ρ ni we obtain σ i +1 / = κγ ∆ p ni +1 / − − σγ ∆ p ni +1 /
00 0 0 00 0 3 σγ − , This scheme is now applied to the following numerical experiment. We choosean initial data given by ρ = 0 . , v = 0 , f = 0 , e = ( . , x ∈ [0 . , . , , otherwise.The parameters of the model are κ = 2 , σ = 1 , ε = 10 − , C p = 10 − and η = 2.The numerical solution is computed with ∆ x = 10 − and compared to Figure 5 witha reference solution obtained by solving (4.6). Once again, the solution is perfectlycaptured even on a coarse grid. We emphasize that this case is very challengingbecause of the very specific form of diffusion involved on the density.7. Concluding remarks
We have presented a general framerwork to investigate the late-time/stiff-rela-xation limit of a large class of hyperbolic systems. This framework was shown tocover the examples of interest arising in the modeling of complex fluid flows whenseveral time-scales are involved. A new class of schemes was proposed for theirapproximation, and we demonstrated that the proposed modification was crucialin order to ensure the correct asymptotic behavior for late times. It would beinteresting to make comparisons between the numerical results obtained here andconcrete experimental results, especially on radiative transfer problems. In futurework, we plan to generalize our continuous and discrete frameworks to problemswith several space dimensions, and to numerically investigate the robustness andaccuracy of such multi-dimensional finite volume discretizations in realistic physicalsituations.
Acknowledgments
The authors gratefully thank Eric Paturel and Friedrich Wagemann for fruitfuldiscussions. The second author (PLF) was partially supported by the Agence Na-tionale de la Recherche (ANR) through the grant 06-2-134423, and by the CentreNational de la Recherche Scientifique (CNRS).
References [1] D. Aregba-Driollet, R. Natalini, Convergence of relaxation schemes for conservation laws,Appl. Anal., 61, no. 1-2, 163–193 (1996).[2] C. Berthon, P. Charrier, B. Dubroca, An HLLC scheme to solve the M1 Model of radiativetransfer in two space dimensions, J. Sci. Comput., 31, no. 3, 347–389 (2007).[3] C. Berthon, R. Turpault, Asymptotic preserving HLL schemes, Numer. Methods PartialDifferential Equations, DOI: 10.1002/num.20586.[4] S. Bianchini, B. Hanouzet, R. Natalini, Asymptotic behavior of smooth solutions for partiallydissipative hyperbolic systems with a convex entropy, Comm. Pure Appl. Math., 60, 1559–1622 (2007). E FLOCH, AND R. TURPAULT
Density
Radiative energy
Figure 5.
Reference (full line) and proposed scheme (dashedlines) solution with and without asymptotic preserving correction.Density ρ is shown above and Radiative energy e is shown below. [5] F. Bouchut, H. Ounaissa, B. Perthame, Upwinding of the source term at interfaces for Eulerequations with high friction, J. Comput. Math. Appl. 53, No. 3-4, 361–375 (2007).[6] C. Buet, S. Cordier, An asymptotic preserving scheme for hydrodynamics radiative transfermodels: numerics for radiative transfer, Numer. Math. 108, 199–221 (2007).[7] C. Buet, B. Despr´es, Asymptotic preserving and positive schemes for radiation hydrodynam-ics, J. Comput. Phys. 215, 717–740 (2006). ATE-TIME/STIFF RELAXATION ASYMPTOTIC-PRESERVING APPROXIMATIONS 29 [8] C. Cercignani, The Boltzmann equation and its applications, Applied Math. Sciences, Vol. 67,Springer-Verlag, New York, 1988.[9] C. Chalons, F. Coquel, E. Godlewski, P.-A. Raviart, N. Seguin Godunov-type schemes forhyperbolic systems with parameter dependent source. The case of Euler system with friction,M2AN Math. Model. Numer. Anal., to appear.[10] C. Chalons, J.-F. Coulombel, Relaxation approximation of the Euler equations, J. Math.Anal. Appl., vol 348, no. 2, 872–893 (2008).[11] G.Q. Chen, C.D. Levermore, T.P. Liu, Hyperbolic conservation laws with stiff relaxationterms and entropy, Comm. Pure Appl. Math. 47, 787–830 (1995).[12] F. Coquel, B. Perthame, Relaxation of energy and approximate Riemann solvers for generalpressure laws in fluid dynamics, SIAM J. Numer. Anal. 35, 2223–2249 (1998).[13] P. Crispel, P. Degond, M.-H. Vignal, An asymptotic preserving scheme for the two-fluidEuler-Poisson model in the quasi-neutral limit, J. Comput. Phys. 223, 208–234 (2007).[14] P. Degond, F. Deluzet, A. Sangam, M.-H. Vignal, An asymptotic preserving scheme for theEuler equations in a strong magnetic field, J. Comput. Phys. 228, 3540–3558 (2009).[15] D. Donatelli and P. Marcati, Convergence of singular limits for multi-D semilinear hyperbolicsystems to parabolic systems, Trans. of AMS 356, 2093–2121 (2004).[16] P. Marcati and B. Rubino, Hyperbolic to parabolic relaxation theory for quasilinear firstorder systems, J. Differential Equations 162, 359–399 (2000).[17] B. Dubroca, J.-L. Feugeas, Entropic moment closure hierarchy for the radiative transferEquation, C. R. Acad. Sci. Paris, Ser. I, 329, 915–920 (1999).[18] L. Gosse, G. Toscani, Asymptotic-preserving well-balanced scheme for the hyperbolic heatequations, C. R. Acad. Sci. Paris, 334, 337–342 (2002).[19] L. Gosse, G. Toscani, Space localization and well-balanced schemes for discrete kinetic modelsin diffusive regimes, SIAM J. Numer. Anal., 41, 641–658 (2003).[20] T. Goudon, J.-E. Morel, B. Despr`es, C. Buet, C. Berthon, J. Dubois, R. Turpault, J.-F.Coulombel,
Mathematical models and numerical methods for radiative transfer,
Panoramaset synth`eses, Vol. 28, Soc. Fran¸caise Math., 2009.[21] A. Harten, P.D. Lax, B. van Leer, On upstream differencing and Godunov-type schemes forhyperbolic conservation laws, SIAM Review, 25, 35–61 (1983).[22] S. Jin, Z. Xin, The relaxation scheme for systems of conservation laws in arbitrary spacedimension, Comm. Pure Appl. Math., 45, 235–276 (1995).[23] N. Ju, Numerical analysis of parabolic p -Laplacian: approximation of trajectories. SIAM J.Numer. Anal. 37, 1861–1884 (2000).[24] R. J. LeVeque, Finite volume methods for hyperbolic problems,
Cambridge Texts in AppliedMathematics, Cambridge Univ. Press, 2002.[25] T.P. Liu, Hyperbolic conservation laws with relaxation, Comm. Math. Phys., 108, 153–175(1987).[26] P. Marcati, Approximate solutions to conservation laws via convective parabolic equations,Comm. Partial Differential Equations, 13, 321–344 (1988).[27] P. Marcati, A. Milani, The one-dimensional Darcy’s law as the limit of a compressible Eulerflow, Journal of Differential Equations, 84, no. 1, 129–146 (1990).[28] F. Marche, Derivation of a new two-dimensional viscous shallow water model with varyingtopography, bottom friction and capillary effects, Eur. J. Mech./B-Fluid, 26, 49–63 (2007).[29] G.N. Minerbo, Maximum entropy Eddington factors, J. Quant. Spectrosc. Radiat. Transfer,20, 541–545 (1978).[30] I. Suliciu, On modelling phase transitions by means of rate-type constitutive equations, shockwave structure, Int. J. Engrg. Sci. 28, 829–841 (1990).[31] I. Suliciu, On the thermodynamics of fluids with relaxation and phase transitions. Fluids withrelaxation, Int. J. Engag. Sci. 36, 921–947 (1998).[32] E.F. Toro,
Riemann solvers and numerical methods for fluid dynamics. A practical introduc-tion,
Second edition, Springer Verlag, Berlin, 1999.[33] J.L. Vazquez, Hyperbolic aspects in the theory of the porous medium equation, in Proceedingsof the Workshop on Metastability and PDE’s, Academic Press, Minneapolis, 1985.[34] W.-A. Yong, Singular perturbations of first-order hyperbolic systems with stiff source terms,J. Differential Equations, 155, 89–132 (1999).[35] W.-A. Yong, K. Zumbrun, Existence of relaxation shock profiles for hyperbolic conservationlaws, SIAM J. Appl. Math., 60, 1565–1575 (2000). E FLOCH, AND R. TURPAULT
Universit´e de Nantes, Laboratoire de Math´ematiques Jean Leray, CNRS UMR 6629,2 rue de la Houssini`ere, BP 92208, 44322 Nantes, France.
E-mail address : [email protected]
Laboratoire Jacques-Louis Lions, Centre National de la Recherche Scientifique,Universit´e Pierre et Marie Curie (Paris 6), 4 Place Jussieu, 75252 Paris, France.
E-mail address : [email protected] Universit´e de Nantes, Laboratoire de Math´ematiques Jean Leray, CNRS UMR 6629,2 rue de la Houssini`ere, BP 92208, 44322 Nantes, France.
E-mail address ::