Regularization and relaxation tools for interface coupling

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Regularization and relaxation tools for interface coupling F. Coquel1, E. Godlewski1, N. Seguin1 1Universit´e Pierre et Marie Curie-Paris6, UMR 7598, Laboratoire Jacques-Louis Lions, Paris, F-75005; CNRS, UMR 7598, LJLL, Paris, F-75005 France. E-mails: [email protected], [email protected], [email protected]. Palabras clave: Hyperbolic systems, coupling, relaxation, Riemann problems, finite volume methods, fluid model Resumen We analyze a relaxation method for approximating the coupling of two Euler systems at a fixed interface and more generally for approximating fluid systems. 1. Introduction We have been considering, in a series of papers, the coupling of systems, both from a theoretical and from a numerical point of views. The actual problem arises from the coupling of codes which model liquid-vapor flows and the systems under consideration are systems of conservation laws of hyperbolic nature. The coupling problem might be interpreted as solving conservation laws with discontinuous coefficients, in which case the flux is assumed to be continuous at the interface. Many works have been recently devoted to the study of conservation laws with discontinuous coefficients, see [1] [4] and the references therein. These conservative approaches involve rather naturally some interface entropy condition along the line of discontinuity. However, we have followed an approach where the coupling is a priori non conservative, which we have named state coupling since (in general) it ensures the continuity or transmission of the state variables, as opposed to flux coupling. In some cases the two may coincide as we will see. In state coupling, we may chose to transmit the conservative variables say u, or any set of state variables vsuch that u→vis an admissible change of variables. We will then speak of v−state coupling (we refer to [2] for detailed definitions). The precise coupling conditions, introduced in [10], impose that two boundary value problems be well-posed. Since boundary value problems for hyperbolic systems are a difficult subject (see [11]), these coupling conditions cannot always be explicited, and 1 F. Coquel, E. Godlewski, N. Seguin moreover may lead to ill-posed problems. However, we can justify our approach from a theoretical point of view in some cases. In fact, the state-coupling approach can handle different models of coupling including flux coupling (via relaxation). First, the coupling conditions can indeed be explicited in the case of Lagrangian systems whose particular structure has enabled us to exhibit a special set of transmitted variables for rather general fluid systems in [2]. In particular, for the Euler system in Lagrangian coordinates, transmitting the primitive variables (v= (%, u, p), %density, u velocity, ppressure) yields the continuity of u, p at the interface and thus results in a conservative coupling [2]. The case of gas dynamics in Eulerian coordinates is more complex since the eigenvalues may change sign (another way of expressing the difficulty is that the flux does not make an admissible change of variables). However, if we choose to transmit the conservative variables, or any set of variables which does not yield continuity of the flux, the coupling is non conservative and moreover there is no natural criteria such as an entropy condition at the interface to select a unique admissible solution. We may think of some regularization procedure to treat the problem. Dafermos regularization was introduced in [5] in the scalar case. This regularization allows the approximation of coupled Riemann problems by smooth profiles and enlights the possible discontinuous behavior of the limit solutions at the interface, but does not give uniqueness as expected in the nonconservative case. In another direction, we have a now classical numerical procedure for approximating the problem which in some cases is justified by convergence results. Indeed, the two flux finite volume method introduced in [10] provides a natural tool for the numerical coupling associated to our state coupling approach. It is based on finite volume methods characterized by their numerical flux and a key point lies in the definition of the two fluxes which express the transmission at the interface of coupling. We have also used approximation by relaxation systems in a global numerical procedure. In [3] we introduce a larger but simpler system on which a state coupling procedure is applied. It ensures a conservative numerical coupling of two Euler systems (in the isentropic case) which avoids resonance (though we do not have uniqueness). This very powerful tool can be used for rather general fluid systems as we will see below. Note also that in the scalar case, we can link the results to those of [4] for conservation laws with discontinuous coefficients. This link should exist for systems which enter into the framework of [12]. We first adress the coupling of two Euler systems (with energy) and then investigate the extension of our relaxation approach to general fluid systems [9]. We present the Lagrangian coordinate case and forget the coupling for this part. However, the results extend to systems in Eulerian coordinates thanks to Wagner’s equivalence result. The approach can also be done directly in Eulerian coordinates but the computations are obviously simpler with the Lagrangian form. For details we refer to [8]. 2. Relaxation solver for flux-coupling of Euler systems We want to ensure a conservative coupling of two Euler systems    ∂tρ+∂x(ρu) = 0, ∂t(ρu) + ∂x(ρu2+p) = 0, x ∈R, t > 0 ∂t(%e) + ∂x((%e +p)u) = 0, (1) 2 Normas para la edici´on (CEDYA 2007) differing by their pressure law at a fixed interface x= 0 p=p(x) = pα(%, ε) = ˜pα(τ, s), α =Lin x < 0, α =Rin x > 0,(2) or equivalently p(x)≡(1 −H(x))pL+H(x)pR, and for some given initial condition, say u0(x). This problem is a simple model for coupling two codes which simulate the same kind of multiphase flow but take into account different closure laws. In (2), pαis a given function expressed either in terms of density and internal energy (%, ε) or equivalently (with a tilde) in terms of specific volume and entropy (τ, s) (with τ= 1/% and the physical specific entropy is −s), satisfying usual assumptions derived from laws of thermodynamics, in particular we can write ε= ˜ε(τ, s) with ˜ετ=−p < 0, ˜εs=−T < 0, and ˜ε(τ, s) strictly convex (see [11]). The internal energy satisfies ε=e−u2/2, the set of states is {u= (%, %u, %e); % > 0, u ∈R, e −u2/2>0}. Discontinuous solutions of (1) are required to satisfy an inequality ∂t%s +∂x(%su)≤0, which becomes an equality for smooth solutions ∂t%s +∂x(%su) = 0.(3) 2.1. The relaxation system We approximate the solutions of (1) by way of the relaxation system        ∂t%+∂x(%u) = 0 ∂t(%u) + ∂x(%u2+ Π) = 0 ∂t%s +∂x(%su) = 0, ∂t(ρT) + ∂x(%Tu) = λρ(τ− T ), (4) with a singular perturbation, where λ > 0 stands for the relaxation parameter. The set of states for (4) is {U= (%, %u, %s, %T)T;% > 0, u ∈R, s > 0,T>0}. The closure relation for Π is Π = Πα(U) = ˜ Πα(τ, s, T)≡˜pα(T, s) + a2(T − τ), α =Lin x < 0, α =Rin x > 0,(5) i.e., the exact pressure law ˜pof (2) modified by a correction term. We will note Fα(U) = (%u, %u2+Π, %su, %Tu)T. In (5), ais a positive constant which plays the role of a Lagrangian sound speed and is required to upper bound the exact sound speed a2>−∂τ˜pα(T, s),(6) for all the states Uunder consideration, this is Whitham’s or subcharacteristic condition. Formally as the relaxation parameter λ→ ∞,T − τ→0, so that Π →pand we indeed recover system (1) at equilibrium where the energy equation has been replaced by the entropy (3). In order to recover a consistent approximation of (1), we follow the approach developped in [7]. The system (4) without source term (λ= 0) when the law pα(pLor pR) is used on the whole line (i.e. without coupling) has four linearly degenerate fields u−aτ, u, u, u +aτ, and explicit solutions for Riemann problems which we note Wα(x/t;U`,Ur). A coupled Riemann problem (CRP) is then a Cauchy problem for (4 with λ= 0) (5), corresponding to a piecewise constant initial data U(x, 0) = U`, x < 0, Ur, x > 0,(7) 3 F. Coquel, E. Godlewski, N. Seguin and a coupling condition associated to a transmitted set of variables say V. The solution of a CRP depends on the choice of V. We have proved existence and uniqueness of the solution of the CRP with transmission of V= (τ, u, Π,Σ), corresponding to data at equilibrium, i.e., such that Tα=τα,α=L, R. The proof assumes (6) is satisfied for all the states occuring in the solution of the Riemann problem. The law Tα(τ, Π, s) is obtained by inverting the relation ˜pα(T, s)+a2T= Π+a2τ, for asatisfying (6) and this yields that Vis an admissible change of variables obtained by composing (τ, u, s, T)T→(τ, u, s, Π)Tand (τ, u, s, Π)T→(τ, u, Σ,Π)T, and we may write U=ϕα(V). Theorem 1. Assuming (6), for any given equilibrium states U`,Ur, the coupled Riemann problem for (4), with transmission of Vadmits a unique solution WV c(x/t;U`,Ur)which coincides with the solution of the classical Riemann problem in each half-space WV c(x/t;U`,Ur) = WL(x/t;U`, ϕ`(Vr)), x < 0, WR(x/t;ϕR(V`),Ur), x > 0.(8) Moreover, %u,%u2+ Π and (%Σ + Π)u, are continuous at the interface. The solution we obtain for the above CRP expressed in terms of Vcoincides with the solution of the Riemann problem of the conservative system        ∂t%+∂x(%u) = 0, ∂t(%u) + ∂x(%u2+ Π) = 0, ∂t(%Π) + ∂x((%Π + a2)u) = 0, ∂t(%Σ) + ∂x((%Σ + Π)u) = 0, (9) where we have replaced the third equation in (4, λ= 0) by an equation on Π which is easy to verify. In the last one, the energy Σ is defined by Σ(τ, u, s, T) = ˜ε(T, s) + u2/2+(˜ Π2(τ, s, T)−˜p2(T, s))/2a2.(10) Since all the characteristic fields are linearly degenerate, the two systems are equivalent. For expressing the solution in terms of U, we use two changes of variables Tαand s(see [8] for details). System (4) is endowed with an entropy (the energy Σ). The relaxation procedure, outside coupling, satisfies some stability principles. First we have a minimization principle. Proposition 1. For a given u= (%, %u, %s)T, let us note Ueq = (%, %u, %s, 1) = (u,1). We have, noting U= (u, %T), %e =%Σ(Ueq) = minT ∈K%Σ(U).(11) Here, Kis a compact interval so that τ∈ ◦ K, and we take a2>maxT ∈K(−∂τ˜p(T, s)). Then a Chapman-Enskog type dissipativity principle is valid. Introducing a first order correction term for Tin inverse powers of λ:T(λ)=τ(λ)+λ−1τ(ν) 1+O(λ−2) (if we emphasize the dependence on λof the solution of (4) and then drop the superscript λ), yields, ∂t%u +∂x%u2+ ˜p(τ, s)=λ−1∂x(∂τ˜p(τ, s) + a2)τ∂xu+O(λ−2). This last equation is a dissipative approximation of the momentum equation of (1) if the stability or Whitham criterion (6) holds. 4 Normas para la edici´on (CEDYA 2007) 2.2. The global relaxation solver The aim of this section is to ensure the conservative coupling of two Euler systems (1) differing by their pressure law (2). In [3], we have introduced a global relaxation approximation of the coupling problem in the isentropic case and we follow the same approach. The scheme associates a relaxation procedure and finite volume numerical schemes. We first give the main lines of the relaxation part of the scheme which involves a fractional step method to advance the solution in time from tnto tn+1 =tn+∆twith three steps: reconstruction, evolution, projection. We describe the main lines of the procedure for n= 0, since they are similar at all other time. Let u0(x)=(%0, %0u0, %0e0)T(x) be an initial datum for system (1), (2): 1. Define the extended initial datum U0(x) = (%0, %0u0, %0s0, %0T0)T(x) for the relaxation system (4), where s0=s(u0) and defining T0≡1/%0:U0is at equilibrium. 2. Solve (approximately) the Cauchy problem (4), (5) with the initial data U0for t∈ (0,∆t] and with appropriate V−coupling condition, we obtain U− 1(x) = U(x, ∆t). 3. Project U− 1= (%1, %1u1, %1s1, %1T1)Ton the equilibrium set of system (4) (instantaneous relaxation) to get U1= (%1, %1u1, %1s1,1). 4. Define u1(x) = (%1, %1u1, %1e1)T(x). In steps 2 and 3 we have solved the coupling of the relaxation systems (4) by an operator splitting method. In the evolution step 2, we will use a Godunov solver, this results in our global relaxation solver (GRS). In fact step 4 is straightforward at the computational level but requires a careful analysis to justify the scheme for u= (%, %u, %e)Tfrom U= (%, %u, %s, %T)T. Moreover u0(x) is first discretized u0 j+1/2=Z(j+1)∆x j∆x u0(x)dx. Let us summarize the resulting global relaxation method for approximating the coupling of two Euler systems (1) differing by their pressure law (2). In each half space, we use a Godunov solver, and at the interface too. Starting from an initial condition u0(x) discretized by u0 j+1/2= (%, %u, %e)0 j+1/2, j ∈Z, setting µ= ∆t/∆x, - Define Un j+1/2= ((%, %u, %s)n j+1/2,1)Tthe extended equilibrium state. - Solve the Riemann problems Wα(0; Un j−1/2,Un j+1/2), α =Lfor j < 0, α =Rfor j > 0. Define the usual numerical fluxes of Godunov’s method GGod,n α,j =Fα(Wα(0; Un j−1/2,Un j+1/2)), j 6= 0.(12) Solve the CRP (8) with data Un −1/2,Un +1/2and prescribed coupling condition, we get an explicit formula for WV c(x/t;Un −1/2,Un +1/2), and define the corresponding numerical flux GGod,n α,0=Fα(WV c(0; Un −1/2,Un +1/2)).(13) 5 F. Coquel, E. Godlewski, N. Seguin - Update Un+1− j+1/2by Godunov’s scheme, we get what we call the Godunov scheme with V−state coupling    Un+1 j−1/2=Un j−1/2−µGGod,n L,j −GGod,n L,j−1, j ≤0, n ≥0, Un+1 j+1/2=Un j+1/2−µGGod,n R,j+1 −GGod,n R,j , j ≥0, n ≥0.(14) - Keep the two first components for the two first components of un+1 j+1/2. This results in a conservative discrete equation for (%, %u)n+1 j+1/2. Here, we skip the details of step 4 which follows [7] for defining the energy in the global solver so as to obtain a standard conservative scheme with an energy component of the flux. Moroever, using (11), an entropy inequality can be established. Proposition 2. The global relaxation solver satisfies a discrete entropy inequality %sn+1 j+1/2≤%sn j+1/2−µGn %s,j+1 − Gn %s,j,(15) where the discrete entropy flux Gn %s,j =Gn α,%s,j is defined by the third component of Godunov’s flux (12) (j6= 0) (resp. of (13), for j= 0). Finally, we have a formal Lax-Wendroff-type convergence result. Denoting by u∆the piecewise constant function associated in a classical way to the scheme, we prove Theorem 2. Assume that u∆is bounded in L∞(R+×R)and that the scheme converges in the sense that u∆→uin L1 loc(R+; L1 loc(R)) and a.e. with u∆(0, t)→u(0, t)in L1 loc(R) and a.e. Then the limit uis solution of system (1) (2) with initial condition u0. 3. Relaxation for Lagrangian systems 3.1. The form of general Lagrangian systems We consider systems of nconservation laws in Lagrangian coordinates which we write ∂tu+∂xf(u) = 0,(xstands for a mass variable) that meet some common properties (we refer to [9] for a detailed description): (i) they are endowed with a strictly convex entropy s(u), with null associated entropy flux, so that for smooth solutions ∂ts= 0. (ii) uis made of n−1−dstate variables wand dvelocity variables v. The last component of uis the total energy which we will denote e un≡e=1 2|v|2+(16) where the internal energy (w, s) is a state variable, then sis also a state variable. We will assume that s(u), satisfies ∂s ∂e(u)≡se(u)<0. The model is then called a fluid model. We also assume (iii) Galilean invariance, and (iv) reversibility in time for smooths solutions. Then, such a system can be written in a canonical form (see [2]), involving the entropy variables that symmetrize the system u∗≡s0(u)T= (su1,· · · , sun−1, se)T.Here we will only use the following form (equivalent for smooth solutions) of the equations    ∂tw−N∂xv=0, ∂tv−NT∂xφ(w, s) = 0, ∂ts= 0, (17) 6 Normas para la edici´on (CEDYA 2007) where Nis a constant d×(n−1−d) matrix, v∈Rd, is the velocity vector as defined above, w∈Rn−d−1represents the state variable vector, and φ(w, s) = ew(u), note that ev(u) = v. We use the notation ew≡ ∇wefor the vector of partial derivatives (ew1, ew2, ..., ewn−d−1), and ewi=∂e ∂wi. Finally the spectrum of f0(u) is symmetric and there are at least n−2dnull eigenvalues, and n−2d≥1 (there is at least one which is associated to the conservation of entropy). We may reduce even further the system in order to distinguish all the null eigenvalues (see [8]). Here we assume for simplicity that Nis a square invertible matrix and n= 2d+ 1. 3.2. The relaxation system Consider, again for simplicity, the isentropic case, s=s0, set (w) = (w, s0), 0=ew, then (17) reads ∂tw−N∂xv=0, ∂tv−NT∂x0(w) = 0.(18) We introduce a larger relaxation system    ∂tw−N∂xv=0, ∂tv−NT∂xX=0, ∂tW=λ(w− W), (19) with some initial condition U0(x) = (w0,v0,W0)(x), where X=X(w,W) = 0(W) + θ0(w)−θ0(W),(20) W ∈ Rd,θ:Rd→Rand we note either θwor θ0its derivative and (., .) the scalar product in Rd. For the system with entropy, all the computations can be done, only replacing the derivatives by partial derivatives, for instance θ=θ(w, s) will depend on s(see again [8] for details). In (20), θwill be chosen in order to assume entropy dissipation. By assumption w,w=00 is symmetric positive definite and we assume θw,w=θ00 is also symmetric positive definite. For example we may take θ(w)=(w,Λw) with Λ a positive definite matrix, even in the simplest case Λ may be diagonal with positive entries. Formally, as the relaxation parameter λ→ ∞,w− W → 0 and at equilibrium, X(w,w) = 0(w), we recover system (18). In order to justify the relaxation procedure, we introduce an energy Σ(w,v,W), Σ(w,v,W) = 1 2|v|2+(W) + θ(w)−θ(W) + ((0−θ0)(W),w− W),(21) which coincides with eat equilibrium, i.e. Σ(w,v,w) = 1 2|v|2+, and should be dissipated. More precisely, we have for fixed w,v e(w,v) = minW∈KΣ(w,v,W) i.e., the minimum of Σ is attained at equilibrium W=w. This holds under the following assumption (H) θ00(W)−00(W) is positive definite for all W ∈ K(where Kdenotes a compact set of the phase space such that w∈ ◦ K). 7 F. Coquel, E. Godlewski, N. Seguin 3.3. Approximation results As in [6], following [13][14] we have convergence results for smooth solutions of (19) to smooth solutions of (18). As in these papers, we address the periodic case and denote by Hs(T) the Sobolev space of functions with period 1. If we apply the results of [13], we can prove a local in time result, without assuming that the initial data is close to equilibrium. Theorem 3. Let s≥2and consider an initial data U0= (w0,v0,W0)in Hs+2(T)that takes values in a compact subset. There exist θ, with θ00 constant, and T > 0such that -∀λ > 0,∃U(λ)= (w(λ),v(λ),W(λ))∈ C([0, T], Hs(T)) unique solution of (19), - system (18) admits a unique solution (we,ve)∈ C([0, T], Hs+2(T)) with initial data (w0,v0), -(v(λ),u(λ))converges towards (we,ve)in C([0, T], Hs(T)) as λ→ ∞. Note that θ00 constant, which corresponds to the constant ain Whitham’s condition (6), is a little restrictive but important for the applications. We can also follow [14], and prove the existence of a (unique global smooth) solution to (19), under an initial data close to an equilibrium state Ue. The proof of these results (see [8]) relies on checking that our system satisfies the structural properties stated in [13][14]. Referencias [1] Adimurthi ; Mishra, Siddhartha ; Gowda, G. D. Veerappa, Optimal entropy solutions for conservation laws with discontinuous flux-functions. J. Hyperbolic Differ. Equ. 2no. 4, 783–837 (2005) [2] A. Ambroso, C. Chalons, F. Coquel, E. Godlewski, F. Lagouti`ere, P.-A. Raviart, N. Seguin, Coupling of general Lagrangian systems. Math. of Computation (to appear, 2007) [3] A. Ambroso, C. Chalons, F. Coquel, E. Godlewski, F. Lagouti`ere, P.-A. Raviart, N. Seguin, Relaxation methods and coupling procedures. Proceedings of ICFD Conference on numerical methods for fluid dynamics, University of Reading, March 2007 [4] F. Bachmann, J. Vovelle, Existence and uniqueness of entropy solution of scalar conservation laws with a flux function involving discontinuous coefficients. Comm. Partial Differential Equations, 31 no. 1-3 (2006), 371–395 [5] B. Boutin, F. Coquel, E. Godlewski, Dafermos regularization for interface coupling of conservation laws. In Eleventh International Conference on Hyperbolic Problems, Theory, Numerics, Applications. Hyp2006 proceedings, Springer (2007) [6] C. Chalons, J.-F. Coulombel, Relaxation approximation of the Euler equations (submitted) [7] F. Coquel, E. Godlewski, B. Perthame, A. In, P. Rascle, Some new Godunov and relaxation methods for two-phase flow problems. In Godunov methods (Oxford, 1999), Kluwer/Plenum, 2001;179–188 [8] F. Coquel, E. Godlewski, N. Seguin, work in preparation. [9] B. Despr´es, Lagrangian systems of conservation laws. Invariance properties of Lagrangian systems of conservation laws, approximate Riemann solvers and the entropy condition. Numer. Math., 89, Vol. 1, 99–134 (2001) [10] E. Godlewski, K.-C. Le Thanh, P.-A. Raviart, The numerical interface coupling of nonlinear hyperbolic systems of conservation laws: II. The case of systems. ESAIM:M2AN, 39 49–692 (2005) [11] E. Godlewski, P.-A. Raviart. Numerical approximation of hyperbolic systems of conservation laws. Applied Math. Science 118, Springer, New York, 1996 [12] E. Isaacson and B. Temple, Nonlinear resonance in systems of conservation laws, Siam J. Appl. Math., 52, no.5, 1260–1278 (1992) [13] Yong W.-A., Singular perturbations of first-order hyperbolic systems with stiff source terms. Journal of Differential Equations, 155, 89–132 (1999) [14] Yong W.-A., Entropy and global existence for hyperbolic balance laws. Arch. Rational. Mech. Anal., 172, 247–266 (2004) 8