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

Well-balanced schemes for Euler equations with gravity Praveen Chandrashekar praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065 Dept. of Mathematics, IIT Delhi 26 October, 2016 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore 1 / 47

Acknowledgements • Dept. of Mathematics, Univ. of W¨ urzburg, Germany ◮ Markus Zenk (PhD student) ◮ Christian Klingenberg (professor) ◮ Dominik Zoar (master student) ◮ Jonas Berberich (master student) • Heidelberg Institute of Theoretical Science ◮ Fritz R¨ opke (professor) ◮ Phillip Edelmann (post-doc) • Airbus Foundation Chair at TIFR-CAM 2 / 47

Euler equations with gravity Flow properties ρ = density , u = velocity E = total energy = e + 1 2 ρu 2 p = pressure , 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 / 47

Euler equations with gravity In compact notation  0  ∂ q ∂t + ∂ f  ∂φ ∂x = − ρ  ∂x ρu where     ρ ρu  , p + ρu 2 q = ρu f =    E ( E + p ) u 4 / 47

Hydrostatic solutions • Fluid at rest u e = 0 • Momentum equation d p e d φ d x = − ρ e (1) d x • 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 / 47

Hydrostatic solutions • Isothermal hydrostatic state, i.e., T e ( x ) = T e = const , then � φ ( x ) � p e ( x ) exp = const (2) RT e ρ e ( x ) = p e ( x ) Density RT e • Polytropic hydrostatic state 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 / 47

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 7 / 47

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 = ⇒ q h (0) = q h,e , ∂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 ) 8 / 47

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 and/or r ≫ 1 9 / 47

Evolution of small perturbations The initial condition is taken to be the following φ = 1 2 x 2 , u = 0 , ρ ( x ) = exp( − φ ( x )) Add small perturbation to equilibrium pressure p ( x ) = exp( − φ ( x )) + ε exp( − 100( x − 1 / 2) 2 ) , 0 < ε ≪ 1 Non-well-balanced scheme ∂φ ∂x ( x i ) ≈ φ i +1 − φ i − 1 , reconstruct ρ, u, p 2∆ x Using exact derivative of potential does not improve results. In practice, φ is only available at grid points. 10 / 47

Evolution of small perturbations −4 −6 10 x 10 10 x 10 Initial Initial Well−balanced Well−balanced Non well−balanced Non well−balanced 8 8 Pressure perturbation Pressure perturbation 6 6 4 4 2 2 0 0 −2 −2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x ε = 10 − 3 , N = 100 cells ε = 10 − 5 , N = 100 cells 11 / 47

Evolution of small perturbations −6 −6 10 x 10 10 x 10 Initial Initial Well−balanced WB, cells=100 Non well−balanced WB, cells=500 8 8 Pressure perturbation Pressure perturbation 6 6 4 4 2 2 0 0 −2 −2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x ε = 10 − 5 , N = 500 cells ε = 10 − 5 12 / 47

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 • Ghosh and Constantinescu: CRWENO scheme 13 / 47

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

1-D finite volume scheme See: Chandrashekar & Klingenberg, SIAM J. Sci. Comput., Vol. 37, No. 3, 2015 • 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 d t + 2 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 f i + 1 2 ) i + 1 i + 1 15 / 47

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, etc. 16 / 47

1-D finite volume scheme • New set of independent variables ρ e − ψ , u, p e − ψ � ⊤ � w = w ( q , x ) := • We can compute q = q ( w , x ) • States for computing numerical flux q L 2 = q ( w L q R 2 = q ( w R 2 , x i + 1 2 ) , 2 , x i + 1 2 ) i + 1 i + 1 i + 1 i + 1 • First order scheme w L w R 2 = w i , 2 = w i +1 i + 1 i + 1 • Higher order scheme: Reconstruct w variables, e.g., 2 = w i + 1 w L 2 m ( θ ( w i − w i − 1 ) , ( w i +1 − w i − 1 ) / 2 , θ ( w i +1 − w i )) i + 1 17 / 47

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 18 / 47

Well-balanced property 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 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 19 / 47

Well-balanced property 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 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. 20 / 47

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 . 21 / 47

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 22 / 47

Approximation of source term • 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 • Gravitational potential required at faces φ i + 1 2 2 = 1 φ i + 1 2( φ i + φ i +1 ) 23 / 47