Well-balanced schemes for Euler equations with gravity Praveen. C - PowerPoint PPT Presentation

Well-balanced schemes for Euler equations with gravity Praveen. C praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065 Conference on Partial Differential Equations LNMIT Jaipur, 10-11 December, 2014 This work supported by Airbus Foundation Chair at TIFR-CAM/ICTS 1 / 61

Acknowledgements • Christian Klingenberg Dept. of Mathematics Univ. of W¨ urzburg, Germany • Fritz R¨ opke and Phillip Edelmann Dept. of Physics Univ. of W¨ urzburg, Germany • Airbus Foundation Chair at TIFR-CAM 2 / 61

Euler equations with gravity Flow properties ρ = density , u = velocity p = pressure , E = total energy Gravitational potential φ ; force per unit volume of fluid − ρ ∇ φ System of conservation laws ∂ρ ∂t + ∂ ∂x ( ρu ) = 0 ∂t ( ρu ) + ∂ ∂ − ρ∂φ ∂x ( p + ρu 2 ) = ∂x ∂E ∂t + ∂ − ρu∂φ ∂x ( E + p ) u = ∂x 3 / 61

Euler equations with gravity Perfect gas assumption � E − 1 � γ = c p 2 ρu 2 p = ( γ − 1) , > 1 c v In compact notation   0 ∂ q ∂t + ∂ f  ∂φ ∂x = − ρ  ∂x ρu where     ρ ρu  , p + ρu 2 q = ρu f =    E ( E + p ) u 4 / 61

Hydrostatic solutions • Fluid at rest u e = 0 • Mass and energy equation satisfied • Momentum equation d p e d φ d x = − ρ e (1) d x • Need additional assumptions to solve this equation • Assume ideal gas and some temperature profile T e ( x ) p e ( x ) = ρ e ( x ) RT e ( x ) , R = gas constant integrate (1) to obtain � x φ ′ ( s ) � � p e ( x ) = p 0 exp − RT e ( s )d s x 0 5 / 61

Hydrostatic solutions • If the hydrostatic state is isothermal , i.e., T e ( x ) = T e = const , then � φ ( x ) � p e ( x ) exp = const (2) RT e Density ρ e ( x ) = p e ( x ) RT e • If the hydrostatic solution is polytropic then we have following relations ν 1 − − p e ρ − ν ν − 1 ν − 1 = const , p e T = const , ρ e T = const (3) e e e where ν > 1 is some constant. From (1) and (3), we obtain νRT e ( x ) + φ ( x ) = const (4) ν − 1 E.g., pressure is ν − 1 p e ( x ) = C 1 [ C 2 − φ ( x )] ν 6 / 61

Existing schemes • Isothermal case: Xing and Shu [2], well-balanced WENO scheme • If ν = γ we are in isentropic case h ( x ) + φ ( x ) = const has been considered by Kappeli and Mishra [1]. • Desveaux et al: Relaxation schemes, general hydrostatic states 7 / 61

Well-balanced scheme • Scheme is well-balanced if it exactly preserves hydrostatic solution. • General evolutionary PDE ∂ q ∂t = R ( q ) • Stationary solution q e R ( q e ) = 0 • We are interested in computing small perturbations q ( x, 0) = q e ( x ) + ε ˜ q ( x, 0) , ε ≪ 1 • Perturbations are governed by linear equation ∂ ˜ q ∂t = R ′ ( q e )˜ q 8 / 61

Well-balanced scheme • Some numerical scheme ∂ q h ∂t = R h ( q h ) • q h,e = interpolation of q e onto the mesh • Scheme is well balanced if ∂ q h R h ( q h,e ) = 0 = ⇒ ∂t = 0 • Suppose scheme is not well-balanced R h ( q h,e ) � = 0 . Solution q h ( x, t ) = q h,e ( x ) + ε ˜ q h ( x, t ) 9 / 61

Well-balanced scheme • Linearize the scheme around q h,e ∂ q h ) = R h ( q h,e ) + εR ′ ∂t ( q h,e + ε ˜ q h ) = R h ( q h,e + ε ˜ h ( q h,e )˜ q h or ∂ ˜ ∂t = 1 q h εR h ( q h,e ) + R ′ h ( q h,e )˜ q h • Scheme is consistent of order r : R h ( q h,e ) = Ch r � q h,e � ∂ ˜ ∂t = 1 q h εCh r � q h,e � + R ′ h ( q h,e )˜ q h • ε ≪ 1 then first term may dominate the second term; need h ≪ 1 • Canonical approach ∂ρu ∂t + ∂ ∂x ( p + ρu 2 ) = − ρ∂φ ∂x d t ( ρu ) i + 1 d φ i +1 − φ i − 1 ∆ x [ ˆ 2 − ˆ f i + 1 f i − 1 2 ] = − ρ i 2∆ x 10 / 61

Scope of present work • Second order finite volume scheme • Ideal gas model: well-balanced for both isothermal and polytropic solutions • Most numerical fluxes can be used 11 / 61

Source term [2] Define � x φ ′ ( s ) ψ ( x ) = − RT ( s )d s, x 0 is arbitrary x 0 Then � x RT ( s )d s = − φ ′ ( x ) φ ′ ( s ) ∂ψ ∂x = − ∂ ∂x RT ( x ) x 0 and ∂x = − exp( ψ ( x )) φ ′ ( x ) ∂x exp( ψ ( x )) = exp( ψ ( x )) ∂ψ ∂ RT ( x ) so that − ρ ( x ) ∂φ ∂x = p ( x ) exp( − ψ ( x )) ∂ ∂x exp( ψ ( x )) Euler equations  0  ∂ q ∂t + ∂ f  exp( − ψ ( x )) ∂ ∂x = p ∂x exp( ψ ( x ))  pu 12 / 61

1-D finite volume scheme • Divide domain into N finite volumes each of size ∆ x • i ’th cell = ( x i − 1 2 , x i + 1 2 ) • semi-discrete finite volume scheme for the i ’th cell ˆ 2 − ˆ �  0  ψ i + 1 ψ i − 1 � f i + 1 f i − 1 2 − e d q i e 2 = e − ψ i 2 d t + p i (5)   ∆ x ∆ x p i u i • ψ i , ψ i + 1 2 etc. are consistent approximations to the function ψ ( x ) • consistent numerical flux ˆ 2 = ˆ f ( q L 2 , q R 2 ) f i + 1 i + 1 i + 1 13 / 61

1-D finite volume scheme Def: Property C The numerical flux ˆ f is said to satisfy Property C if for any two states q L = [ ρ L , 0 , p/ ( γ − 1)] q R = [ ρ R , 0 , p/ ( γ − 1)] and we have ˆ f ( q L , q R ) = [0 , p, 0] ⊤ • states q L , q R in the above definition correspond to a stationary contact discontinuity. • Property C = ⇒ numerical flux exactly support a stationary contact discontinuity. • Examples of such numerical flux: Roe, HLLC 14 / 61

1-D finite volume scheme • First order scheme q L q R 2 = q i , 2 = q i +1 i + 1 i + 1 • Higher order scheme: To obtain the states q L 2 , q R 2 , reconstruct i + 1 i + 1 the following set of variables ρ e − ψ , u, p e − ψ � ⊤ � w = • Once w L 2 etc. are computed, the primitive variables are obtained as i + 1 ψ i + 1 ψ i + 1 ρ L 2 ( w 1 ) L u L 2 = ( w 2 ) L p L 2 ( w 3 ) L 2 = e 2 , 2 , 2 = e 2 , etc. i + 1 i + 1 i + 1 i + 1 i + 1 i + 1 15 / 61

Well-balanced property Theorem The finite volume scheme (5) together with a numerical flux which satisfies property C and reconstruction of w variables is well-balanced in the sense that the initial condition given by u i = 0 , p i exp( − ψ i ) = const , ∀ i (6) is preserved by the numerical scheme. Proof : Start computation with an initial condition that satisfies (6). Since we reconstruct the variables w , at any interface i + 1 2 we have ( w 2 ) L 2 = ( w 2 ) R ( w 3 ) L 2 = ( w 3 ) R 2 = 0 , i + 1 i + 1 i + 1 i + 1 2 Hence u L 2 = u R p L 2 = p R 2 = 0 , 2 = p i exp( ψ i + 1 2 − ψ i ) =: p i + 1 i + 1 i + 1 i + 1 i + 1 2 16 / 61

Well-balanced property and at i − 1 2 u L 2 = u R p L 2 = p R 2 = 0 , 2 = p i exp( ψ i − 1 2 − ψ i ) =: p i − 1 i − 1 i − 1 i − 1 i − 1 2 Since the numerical flux satisfies property C, we have ˆ ˆ 2 , 0] ⊤ , 2 , 0] ⊤ f i − 1 2 = [0 , p i − 1 f i + 1 2 = [0 , p i + 1 Mass and energy equations are already well balanced, i.e., d q (1) d q (3) i i = 0 , = 0 d t d t Momentum equation: on the left we have f (2) f (2) ˆ 2 − ˆ p i + 1 2 − p i − 1 i + 1 i − 1 2 2 = ∆ x ∆ x 17 / 61

Well-balanced property while on the right 2 − ψ i − p i e ψ i + 1 ψ i − 1 ψ i + 1 ψ i − 1 2 − ψ i p i + 1 2 − p i − 1 2 − e p i e − ψ i e = p i e 2 2 = ∆ x ∆ x ∆ x and hence d q (2) i = 0 d t This proves that the initial condition is preserved under any time integration scheme. 18 / 61

Approximation of source term • How to approximate ψ i , ψ i + 1 2 , etc. ? Need some quadrature • well-balanced property independent of quadrature rule to compute ψ . • To preserve isothermal/polytropic solutions exactly, the quadrature rule has to be exact for these cases. • To compute the source term in the i ’th cell, we define the function ψ ( x ) as follows � x φ ′ ( s ) ψ ( x ) = − RT ( s )d s x i where we chose the reference position as x i . 19 / 61

Approximation of source term • To approximate the integrals we define the piecewise constant temperature as follows T ( x ) = ˆ T i + 1 2 , x i < x < x i +1 (7) where ˆ T i + 1 2 is the logarithmic average given by T i +1 − T i ˆ T i + 1 2 = log T i +1 − log T i • The integrals are evaluated using the approximation of the temperature given in (7) leading to the following expressions for ψ . ψ i = 0 � x i − 1 φ i − φ i − 1 1 2 φ ′ ( s )d s = 2 ψ i − 1 = − R ˆ R ˆ T i − 1 T i − 1 2 x i 2 2 � x i + 1 φ i − φ i + 1 1 2 φ ′ ( s )d s = 2 ψ i + 1 = − R ˆ R ˆ T i + 1 T i + 1 2 x i 2 2 20 / 61

Approximation of source term • Gravitational potential required at faces φ i + 1 2 • φ is governed by Poisson equation and hence is a smooth function. We can interpolate 2 = 1 φ i + 1 2( φ i + φ i +1 ) Sufficient to obtain second order accuracy. Then 2 = φ i − φ i − 1 2 = φ i − φ i +1 ψ i − 1 , ψ i = 0 , ψ i + 1 (8) 2 R ˆ 2 R ˆ T i − 1 T i + 1 2 2 21 / 61

Approximation of source term Theorem The source term discretization given by (8) is second order accurate. Proof : The source term in (5) has the factor � � � � φ i − φ i +1 φ i − φ i − 1 exp − exp ψ i + 1 ψ i − 1 2 R ˆ 2 R ˆ 2 − e e − ψ i e T i + 1 T i − 1 2 2 2 = using (8) ∆ x ∆ x Using a Taylor expansion around x i we get 1 = 1 1 = 1 [1 + O (∆ x 2 )] , [1 + O (∆ x 2 )] ˆ ˆ T i T i T i − 1 T i + 1 2 2 22 / 61