A new technique for the numerical solution of the compressible Euler - PowerPoint PPT Presentation

A new technique for the numerical solution of the compressible Euler equations with arbitrary Mach numbers ∗ Miloslav Feistauer and V´ aclav Kuˇ cera Charles University Prague Faculty of Mathematics and Physics ∗ Presented at conference HYP06, Lyon 17–21 July 2006.

Our goal: to develop a sufficiently accurate and robust method allowing the solution of problems with a wide range of Mach numbers Inviscid flow 2D flow: Basic equations: continuity eq., Euler equtions of motion, en- ergy equation = nonlinear hyperbolic system of conservation laws 2 ∂ w ∂ f s ( w ) � ∂t + = 0 in Q T = Ω × (0 , T ) , (1) ∂x s s =1 w = ( w 1 , . . . , w 4 ) T = ( ρ, ρv 1 , ρv 2 , E ) T (2) ( f i 1 , . . . , f im ) T f i ( w ) = (3) ( ρv i , ρv 1 v i + δ 1 i p, ρv 2 v i + δ 2 i p, ( E + p ) v i ) T = p = ( γ − 1) ( E − ρ | v | 2 / 2) . (4) Notation : ρ – density, p – pressure, E – total energy, v = ( v 1 , v 2 ) – velocity, δ sk – Kronecker symbol, γ > 1 – Pois- son adiabatic constant

Initial condition: w ( x , 0) = w 0 ( x ) , x ∈ Ω , (5) Boundary conditions chosen in such a way that the problem is linearly well–posed Important property of fluxes: A s ( w ) = D f s ( w ) , (6) D w – the Jacobi matrices of the mappings f s f s , – homogeneous mappings of order one = ⇒ f s ( w ) = A s ( w ) w (7)

Discontinuous Galerkin (DG) space discretization for N = 2 Mesh : Ω h = polygonal approximation of Ω T h = mesh in Ω h consisting of triangles or quadrilaterals K i ∈ T h , i ∈ I (= suitable index set) Γ ij = common edge between two neighbouring elements K i and K j or face ⊂ Ω h Introduce index sets s ( i ) so that j ∈ s ( i ) = ⇒ Γ ij ⊂ Ω h γ ( i ) so that i ∈ γ ( i ) = ⇒ Γ ij ⊂ ∂ Ω h S ( i ) = s ( i ) ∪ γ ( i ) . Then � � ∂K i = Γ ij , ∂K i ∩ ∂ Ω h = Γ ij . (8) j ∈ S ( i ) j ∈ γ ( i ) n ij = (( n ij ) 1 , ( n ij ) 2 ) = unit outer normal to ∂K i on the side Γ ij

Γ ij n K i ij K j Neighbouring elements K i , K j – mesh can be nonconforming

Space of approximate solutions Discontinuous piecewise polynomial functions: S h = [ S h ] 4 , (9) S h ≡ S r, − 1 (Ω , T h ) = { v ; v | K ∈ P r ( K ) ∀ K ∈ T h } , r ≥ 0 – integer and P r ( K ) denotes the space of all polyno- mials on K of degree ≤ r . Derivation of the discrete problem : w – exact sufficiently regular solution multiply (1) by any ϕ ∈ S h , integrate over any K i , i ∈ I , apply Green’s theorem, sum over all K i , i ∈ I

� ⇒ ∂ � = w · ϕ d x (10) ∂t K i K i ∈T h 2 � f s ( w ) · ∂ ϕ � � d x − − ∂x s K i s =1 K i ∈T h 2 � � � � + f s ( w ) · ϕ n s d S = 0 . Γ ij K i ∈T h s =1 j ∈ S ( i ) Numerical flux Approximation of fluxes through interfaces Γ ij – with the aid of a numerical flux H ( α, β, n ) : let w h , ϕ h ∈ S h 2 � � f s ( w h ) n s · ϕ h d S (11) Γ ij s =1 � ≈ H ( w h | Γ ij , w h | Γ ji , n ij ) · ϕ h d S Γ ij w h | Γ ij = the value of w h on Γ ij considered from the interior of K i , w h | Γ ji = the value of w h on Γ ij considered from the exterior of K i

Properties of the numerical flux: – (locally) Lipschitz continuous – consistent: H ( α, α, n ) = � N s =1 f s ( α ) n s – conservative: H ( α, β, n ) = − H ( β, α, − n )

Definition of forms For w h , ϕ h ∈ S h we define � ( w h , ϕ h ) h = w h · ϕ h d x , (12) Ω h B h ( w h , ϕ h ) (13) 2 � f s ( w h ) · ∂ ϕ h � � = − d x ∂x s K s =1 K ∈T h � � � + H ( w h | Γ ij , w h | Γ ji , n ij ) · ϕ h d S. Γ ij K i ∈T h j ∈ S ( i ) Approximate DGFE semidiscrete solution: w h ∈ C 1 ([0 , T ] , S h ) , a) d b) d t ( w h ( t ) , ϕ h ) + B h ( w h ( t ) , ϕ h ) = 0 (14) ∀ ϕ h ∈ S h ∀ t ∈ (0 , T ) , w h (0) = Π h w 0 = L 2 − projection of w 0 on S h c)

(14) ≡ large system of ordinary differential equations – usually solved by explicit Runge-Kutta methods – work well for fast flow - conditionally stable!!! Slow flows – with very low Mach numbers – represent big obstacle for all existing methods: FD, FV, FE − → general opinion that this type of flow can be solved only on the basis of modified governing equations using physical variables and derived on the basis of a multiscale analysis R. Klein, C.-D. Munz, A. Meister,... Question: Is it possible to develop a method using conser- vation form of the Euler equations and robust with respect to the Mach number?

Successful technique allowing the solution of all Mach numbers flows: DGFEM & Ingredients: 1) Semi-implicit time discretization (M.F. & V. Dolejˇ s ´ ı, JCP 2004): For each k ≥ 0 find w k +1 such that h w k +1 a) ∈ S h , (15) h    w k +1 − w k h , w k +1 h h + b h ( w k b) , ϕ h , ϕ h ) = 0 , ∀ ϕ h ∈ S h , k = 0 , 1 , . . . ,  h τ h w 0 h = Π h w 0 . c) h , w k +1 , ϕ h ) - linear with respect to ϕ h and w k +1 b h ( w k h h

Definition of b h : Partial linearization of B h We have B h ( w k +1 , ϕ h ) (16) h 2 � ( x )) · ∂ ϕ h ( x ) � � f s ( w k +1 = − d x h ∂x s K s =1 K ∈T h � �� � =:˜ σ 1 � � � H ( w k +1 | Γ ij , w k +1 + | Γ ji , n ij ) · ϕ h d S . h h Γ ij K i ∈T h j ∈ S ( i ) � �� � =:˜ σ 2

Linearization of ˜ σ 1 : Homogeneity of f s ⇒ 2 2 � � f s ( w ) n s = A s ( w ) w n s . (17) s =1 s =1 − → σ 1 ≈ σ 1 ˜ (18) 2 � ( x ) · ∂ ϕ h ( x ) � � h ( x )) w k +1 A s ( w k = d x . h ∂x s K s =1 K ∈T h

Linearization of ˜ σ 2 = Γ ij H ( w k +1 | Γ ij , w k +1 � � � | Γ ji , n ij ) · ϕ h d S : K i ∈T h j ∈ S ( i ) h h Assume that H = Vijayasundaram numerical flux: The matrix P ( w , n ) = � 2 s =1 A s ( w ) n s is diagonalizable: P ( w , n ) = T DT − 1 , D = diag ( λ 1 , . . . , λ 4 ) , (19) λ 1 , . . . , λ 4 = the eigenvalues of P . “Positive” and “negative” part of P : P ± ( w , n ) = T D ± T − 1 , D ± = diag ( λ ± 1 , . . . , λ ± 4 ) . (20) Vijayasundaram numerical flux: H V ( w 1 , w 2 , n ) (21) � w 1 + w 2 � � w 1 + w 2 � = P + w 1 + P − , n , n w 2 . 2 2

= ⇒ � � P + � � � � w k +1 � w k ˜ σ 2 ≈ h � ij , n ij | Γ ij (22) h Γ ij K i ∈T h j ∈ S ( i ) + P − � � � w k +1 � w k h � ij , n ij | Γ ji · ϕ h d S, h where 1 � � � w k w k h | Γ ij + w k h � ij ≡ h | Γ ji . (23) 2

⇒ Linearized convection form : = h , w k +1 b h ( w k , ϕ h ) (24) h 2 � ( x ) · ∂ ϕ h ( x ) � � h ( x )) w k +1 A s ( w k = − d x h ∂x s K K ∈T h s =1 � � P + � � � � w k +1 � w k + h � ij , n ij | Γ ij h Γ ij K i ∈T h j ∈ S ( i ) + P − � � � w k +1 � w k · ϕ h d S, h � ij , n ij | Γ ji h which is linear with respect to w k +1 and ϕ h . h Semi-implicit linearized numerical scheme: for each k ≥ 0 find w k +1 such that h w k +1 a) ∈ S h , (25) h � � � � w k +1 h , w k +1 h + τ k b h ( w k w k b) , ϕ h , ϕ h ) = h , ϕ h h h h ∀ ϕ h ∈ S h , k = 0 , 1 , . . . , w 0 h = Π h w 0 . c)

2) The use of isoparametric elements in order to obtain an accurate solution near curved boundaries 3) Special treatment of boundary conditions transparent for accoustic waves, based on the use of characteristic vari- ables 4) Shock capturing for avoiding the Gibbs phenomenon manifested by spurious overshoots and undershoots in com- puted quantities near discontinuities (shock waves, contact discontinuities) Artificial viscosity proposed by Jaffr´ e, Johnson and Szepessy - introduces some amount of artificial viscosity everywhere - nonphysical entropy production We introduce stabilization terms combining ideas of Jaffr´ e, Johnson and Szepessy and M.F., Dolejˇ s ´ ı and Schwab.

a)Define the discontinuity indicator g k ( i ) proposed by M.F., Dolejˇ s ´ ı and Schwab: � h ] 2 d S /( h K i | K i | 3 / 4 ) , g k ( i ) = [ ρ k K i ∈ T h . (26) ∂K i [ u ] | Γ ij = u Γ ij − u Γ ji = the jump on Γ ij of a function u ∈ S h . b)Define the discrete indicator G k ( i ) = 0 if g k ( i ) < 1 , G k ( i ) = 1 if g k ( i ) ≥ 1 , K i ∈ T h . (27) c)To the left-hand side of (15), b) we add the artificial viscosity form � � h , w k +1 ∇ w k +1 β h ( w k h K i G k ( i ) , ϕ ) = ν 1 · ∇ ϕ d x (28) h h K i i ∈ I with ν 1 ≈ 1 . d)Augment the left-hand side of (15), b) by adding the form � 1 � � h , w k +1 [ w k +1 J h ( w k 2( G k ( i )+ G k ( j )) , ϕ ) = ν 2 ] · [ ϕ ] d S , h h Γ ij i ∈ I j ∈ s ( i ) (29)

where ν 2 ≈ 1 . Resulting scheme: w k +1 a) ∈ S h , (30) h    w k +1 − w k h , w k +1 h + b h ( w k h b) , ϕ h , ϕ h )  h τ k h h , w k +1 h , w k +1 + β h ( w k , ϕ h ) + J h ( w k , ϕ h ) = 0 , ∀ ϕ h ∈ S h , k = 0 , 1 , . . . , h h w 0 h = Π h w 0 . c) This method successfully overcomes problems with the Gibbs phenomenon in the context of the semi-implicit scheme. Important: G k ( i ) vanishes in regions where the solution is regular. = ⇒ The scheme does not produce any nonphysical entropy in these regions.