Divergence-free discontinuous Galerkin method for ideal compressible MHD equations Praveen Chandrashekar praveen@math.tifrbng.res.in http://cpraveen.github.io Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore-560065, India http://math.tifrbng.res.in Seminar, Dept. of Mathematics, IIT Delhi 24 October 2019 1 / 37
Maxwell Equations Linear hyperbolic system ∂ B ∂ D ∂t + ∇ × E = 0 , ∂t − ∇ × H = − J B = magnetic flux density D = electric flux density E = electric field H = magnetic field J = electric current density µ, ε ∈ R 3 × 3 symmetric B = µ H , D = ε E , J = σ E ε = permittivity tensor µ = magnetic permeability tensor σ = conductivity ∂ρ ∇ · B = 0 , ∇ · D = ρ (electric charge density) , ∂t + ∇ · J = 0 2 / 37
Two fluid MHD Non-linear hyperbolic system Conservation laws for each species ∂ρ α ∂t + ∇ · ( ρ α v α ) =0 ∂ ( ρ α v α ) + ∇ · ( ρ α v α ⊗ v α + p α I ) = 1 ρ α q α ( E + v α × B ) , α = i, e ∂t m α ∂ E α ∂t + ∇ · [( E α + p α ) v α ] = 1 ρ α q α E · v α m α γ α − 1 + 1 p α 2 ρ α | v α | 2 Total energy: E α = Coupled with Maxwell’s equations ∂ B 1 ∂ E ∂t + ∇ × E = 0 , ∂t − ∇ × B = − µ 0 ( ρ i q i v i + ρ e q e v e ) c 2 together with the constraints ∇ · E = 1 ∇ · B = 0 , ( ρ i q i + ρ e q e ) ǫ 0 3 / 37
Ideal compressible MHD equations Nonlinear hyperbolic system Compressible Euler equations with Lorentz force ∂ρ ∂t + ∇ · ( ρ v ) = 0 ∂ ( ρ v ) + ∇ · ( P I + ρ v ⊗ v − B ⊗ B ) = 0 ∂t ∂E ∂t + ∇ · (( E + P ) v + ( v · B ) B ) = 0 ∂ B ∂t − ∇ × ( v × B ) = 0 P = p + 1 γ − 1 + 1 2 ρ | v | 2 + 1 p 2 | B | 2 , 2 | B | 2 E = Magnetic monopoles do not exist: = ⇒ ∇ · B = 0 4 / 37
Divergence constraint ∂ B Rotated shock tube ∂t + ∇ × E = 0 Powell 0 . 1 0 . 0 ∇ · ∂ B ∂t + ∇ · ∇× E = 0 − 0 . 1 � �� � Relative error on B � cleaning =0 0 . 1 ∂ 0 . 0 ∂t ∇ · B = 0 − 0 . 1 Athena If 0 . 1 ∇ · B = 0 at t = 0 0 . 0 − 0 . 1 then − 0 . 4 − 0 . 2 0 . 0 0 . 2 0 . 4 x ∇ · B = 0 for t > 0 Guillet et al., MNRAS 2019 Discrete div-free = ⇒ positivity Intrinsic property, not a dynamical (Kailiang Wu) equation 5 / 37
Objectives • Based on conservation form of the equations • Upwind-type schemes using Riemann solvers • Divergence-free schemes for Maxwell’s and compressible MHD ◮ Cartesian grids at present ◮ Divergence preserving schemes (RT) ◮ Divergence-free reconstruction (BDM) 1 • High order accurate ◮ discontinuous-Galerkin FEM • Non-oscillatory schemes for MHD ◮ using limiters ◮ div-free reconstruction using BDM • Explicit time stepping ◮ local mass matrices • Based on previous work for induction equation ◮ J. Sci. Comp., Vol. 79, pp, 79-102, 2019 1 Hazra et al., JCP, Vol. 394, 2019 6 / 37
Some existing methods Exactly divergence-free methods • Yee scheme (Yee (1966)) • Projection methods (Brackbill & Barnes (1980)) • Constrained transport (Evans & Hawley (1989)) • Divergence-free reconstruction (Balsara (2001)) • Globally divergence-free scheme (Li et al. (2011), Fu et al, (2018)) Approximate methods • Locally divergence-free schemes (Cockburn, Li & Shu (2005)) • Godunov’s symmetrized version of MHD (Powell, Gassner et al., C/K) • Divergence cleaning methods (Dedner et al.) 7 / 37
MHD equations in 2-D ∂ U ∂t + ∂ F x ∂x + ∂ F y ∂y = 0 ρ ρv x ρv y P + ρv 2 x − B 2 ρv x v y − B x B y ρv x x P + ρv 2 y − B 2 ρv x v y − B x B y ρv y y ρv x v z − B x B z ρv y v z − B y B z ρv z U = , F x = , F y = E ( E + P ) v x − B x ( v · B ) ( E + P ) v y − B y ( v · B ) B x 0 E z B y − E z 0 B z v x B z − v z B x v y B z − v z B y where P = p + 1 γ − 1 + 1 p 2 ρ | v | 2 + 1 2 | B | 2 , 2 | B | 2 B = ( B x , B y , B z ) , E = E z is the electric field in the z direction E z = − ( v × B ) z = v y B x − v x B y 8 / 37
MHD equations in 2-D Split into two parts U = [ ρ, ρ v , E , B z ] ⊤ , B = ( B x , B y ) ∂ U ∂B x ∂t + ∂E z ∂B y ∂t − ∂E z ∂t + ∇ · F ( U , B ) = 0 , ∂y = 0 , ∂x = 0 The fluxes F = ( F x , F y ) are of the form ρv x ρv y P + ρv 2 x − B 2 ρv x v y − B x B y x P + ρv 2 y − B 2 ρv x v y − B x B y y F x = , F y = ρv x v z − B x B z ρv y v z − B y B z ( E + P ) v x − B x ( v · B ) ( E + P ) v y − B y ( v · B ) v x B z − v z B x v y B z − v z B y 9 / 37
Constraint preserving finite difference Store magnetic field on the faces: ( B x ) i + 1 2 ,j , ( B y ) i,j + 1 2 ( E z ) i + 1 2 − ( E z ) i + 1 ∂B x ∂t + ∂E z d 2 ,j + 1 2 ,j − 1 ∂y = 0 = ⇒ d t ( B x ) i + 1 2 ,j + 2 = 0 ∆ y ( E z ) i + 1 2 − ( E z ) i − 1 ∂B y ∂t − ∂E z d 2 ,j + 1 2 ,j + 1 ∂x = 0 = ⇒ d t ( B y ) i,j + 1 2 − 2 = 0 ∆ x Measure divergence at cell center ( B x ) i + 1 2 ,j − ( B x ) i − 1 ( B y ) i,j + 1 2 − ( B y ) i,j − 1 2 ,j 2 ∇ h · B i,j = + ∆ x ∆ y Then d d t ∇ h · B i,j = 0 The corner fluxes cancel one another !!! 10 / 37
Approximation of magnetic field B h ∈ V h = some finite element space If we want ∇ · B h = 0 , it is natural to look for approximations in B h ∈ V h ⊂ H ( div, Ω) = { B ∈ L 2 (Ω) : div( B ) ∈ L 2 (Ω) } For conformal approximate of functions in H ( div, Ω) on a mesh T h with piecewise polynomials, we need B · n continuous across element faces 11 / 37
Approximation spaces: Degree k ≥ 0 Map cell K to reference cell ˆ K = [ − 1 2 , + 1 2 ] × [ − 1 2 , + 1 2 ] P r ( ξ ) = span { 1 , ξ, ξ 2 , . . . , ξ r } , Q r,s ( ξ, η ) = P r ( ξ ) ⊗ P s ( η ) Hydrodynamic variables in each cell k k � � U ( ξ, η ) = U ij φ i ( ξ ) φ j ( η ) ∈ Q k,k i =0 j =0 Normal component of B on faces k � on vertical faces : b x ( η ) = a j φ j ( η ) ∈ P k ( η ) j =0 k � on horizontal faces : b y ( ξ ) = b j φ j ( ξ ) ∈ P k ( ξ ) j =0 { φ i ( ξ ) } are orthogonal polynomials on [ − 1 2 , + 1 2 ] , with degree φ i = i . 12 / 37
Approximation spaces: Degree k ≥ 0 For k ≥ 1 ,define certain cell moments � + 1 � + 1 1 2 2 α ij = α ij ( B x ) := B x ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η, 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k m ij − 1 − 1 � �� � 2 2 Q k − 1 ,k � + 1 � + 1 1 2 2 β ij = β ij ( B y ) := B y ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η, 0 ≤ i ≤ k, 0 ≤ j ≤ k − 1 m ij − 1 − 1 � �� � 2 2 Q k,k − 1 � + 1 � + 1 � + 1 2 2 2 [ φ i ( ξ ) φ j ( η )] 2 d ξ d η = m i m j , [ φ i ( ξ )] 2 d ξ m ij = m i = − 1 − 1 − 1 2 2 2 α 00 , β 00 are cell averages of B x , B y Solution variables { U ( ξ, η ) } , { b x ( η ) } , { b y ( ξ ) } , { α, β } The set { b x , b y , α, β } are the dofs for the Raviart-Thomas space. 13 / 37
RT reconstruction: b ± x ( η ) , b ± y ( ξ ) , α, β → B ( ξ, η ) Given b ± x ( η ) ∈ P k and b ± 2 3 y ( ξ ) ∈ P k , b + y ( ξ ) and set of cell moments x ( η ) x ( η ) α { α ij , 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k } b − b + β { β ij , 0 ≤ i ≤ k, 0 ≤ j ≤ k − 1 } b − y ( ξ ) 0 1 Find B x ∈ Q k +1 ,k and B y ∈ Q k,k +1 such that 2 , η ) = b ± 2 ) = b ± B x ( ± 1 η ∈ [ − 1 2 , 1 B y ( ξ, ± 1 ξ ∈ [ − 1 2 , 1 x ( η ) , 2 ] , y ( ξ ) , 2 ] � + 1 � + 1 1 2 2 B x ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η = α ij , 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k m ij − 1 − 1 2 2 � + 1 � + 1 1 2 2 B y ( ξ, η ) φ i ( ξ ) φ j ( η ) d ξ d η = β ij , 0 ≤ i ≤ k, 0 ≤ j ≤ k − 1 m ij − 1 − 1 2 2 (1) ∃ unique solution. (2) B · n continuous. (3) Data div-free = ⇒ reconstructed B is div-free. 14 / 37
DG scheme for B on faces On every vertical face of the mesh: ∂B x ∂t + ∂E z ∂y = 0 � + 1 � + 1 ∂t φ i d η − 1 d φ i d η d η + 1 ∂b x 2 2 ˆ ∆ y [ ˜ E z E z φ i ] = 0 , 0 ≤ i ≤ k ∆ y − 1 − 1 2 2 On every horizontal face of the mesh: ∂B y ∂t − ∂E z ∂x = 0 � + 1 � + 1 ∂t φ i d ξ + 1 d φ i d ξ d ξ − 1 ∂b y 2 2 ˆ ∆ x [ ˜ E z E z φ i ] = 0 , 0 ≤ i ≤ k ∆ x − 1 − 1 2 2 Numerical fluxes ˆ E z : on face, 1-D Riemann solver ˜ E z : at vertex, 2-D Riemann solver 15 / 37
DG scheme for B on cells � + 1 � + 1 m ij d α ij ∂B x 2 2 = ∂t φ i ( ξ ) φ j ( η ) d ξ d η, 0 ≤ i ≤ k − 1 , 0 ≤ j ≤ k d t − 1 − 1 2 2 � + 1 � + 1 − 1 ∂E z 2 2 = ∂η φ i ( ξ ) φ j ( η ) d ξ d η ∆ y − 1 − 1 2 2 � + 1 − 1 2 [ ˆ 2 ) − ˆ E z ( ξ, 1 2 ) φ i ( ξ ) φ j ( 1 E z ( ξ, − 1 2 ) φ i ( ξ ) φ j ( − 1 = 2 )] d ξ ∆ y − 1 2 � + 1 � + 1 + 1 2 2 E z ( ξ, η ) φ i ( ξ ) φ ′ j ( η ) d ξ d η ∆ y − 1 − 1 2 2 Numerical fluxes ˆ E z : on face, 1-D Riemann solver Not a Galerkin method, test functions ( Q k − 1 ,k ) different from trial functions ( Q k +1 ,k ) 16 / 37
Recommend
More recommend
