Simulation of Field-Aligned Simulation of Field-Aligned Ideal MHD - - PowerPoint PPT Presentation

von Karman Institute for Fluid Dynamics

11th INTERNATIONAL CONFERENCE ON HYPERBOLIC PROBLEMS: THEORY, NUMERICS, APPLICATIONS Mehmet Sarp YALIM David Vanden ABEELE, Andrea LANI, Herman DECONINCK

Department of Aeronautics and Aerospace von Karman Institute for Fluid Dynamics (VKI) 72, Chaussée de Waterloo B-1640 Rhode-Saint-Genèse BELGIUM

Simulation of Field-Aligned Simulation of Field-Aligned Ideal MHD Flows Around deal MHD Flows Around Perfectly Conducting Cylinders Using An Artificial Perfectly Conducting Cylinders Using An Artificial Compressibility Approach Compressibility Approach

von Karman Institute for Fluid Dynamics

OUTLINE

• Ideal Magnetohydrodynamics (MHD) Equations
• An Overview of Solenoidal Constraint Satisfying Techniques
• Artificial Compressibility Approach (ACA)
• The COOLFluiD Framework
• Upwind Finite Volume Method (FVM) MHD Solver
• Results
• Conclusion

von Karman Institute for Fluid Dynamics

IDEAL MHD EQUATIONS

Assumptions

• quasi-neutral plasma
• infinite conductivity
• non-relativistic plasma
• negligible Hall term

( )

2 2 = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ + + − − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ + + ⋅ ∇ + ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ B u B u B B Bu uB BB B B I uu u B u p E p E t ρ ρ ρ ρ

von Karman Institute for Fluid Dynamics

DIMENSIONLESS NUMBERS IN MHD

2

2 pressure magnetic pressure acoustic μ β B p = =

η Lu Rm =

Plasma Beta Magnetic Reynolds Number

so plasmas cs astrophysi for ≈ ⎭ ⎬ ⎫ ∞ → ∞ → η

m

R L

von Karman Institute for Fluid Dynamics

MHD WAVES

Euler equations 3 waves (shock, contact surface, rarefaction) Ideal MHD equations 7 waves MORE COMPLEX

• 1 entropy wave
• 2 Alfven waves (incompressible wave, strongly anisotropic)
• 2 fast magnetoacoustic waves (compressible wave, strongly

anisotropic)

• 2 slow magnetoacoustic waves (compressible wave, strongly

anisotropic)

von Karman Institute for Fluid Dynamics

von Karman Institute for Fluid Dynamics

divB CONSTRAINT ?

• Numerically divB can never be exactly zero.
• Discretization errors accumulate at each iteration.
• After a significant amount of error accumulation, the

algorithm breaks down. Precaution: Errors in divB should be damped. For this purpose:

• Powell source term
• Balsara’s scheme
• Projection scheme

von Karman Institute for Fluid Dynamics

POWELL SOURCE TERM (1)

Sources of Instabilities for Ideal MHD Equations

• Jacobian matrix is singular (one of the modes remain

undamped)

• Terms proportional to divB act as extra unphysical forces

Solution Proposed by Powell A source term is added to the RHS

• to remove the singularity in a mathematically appropriate way
• to cancel the terms proportional to divB

von Karman Institute for Fluid Dynamics

POWELL SOURCE TERM (2)

= ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ ∇ ∇ ⋅ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ ∇ ∂ ∂ ρ ρ B u B t

• no change at PDE level (i.e. mathematically divB = 0 is a constraint)
• non-conservative source term introduces non-conservativity into the

PDE system

• new divergence wave (i.e. 8. wave for ideal MHD equations)
• errors in divB will be swept away by the flow

( )

⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ ⋅ −∇ = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ + + − − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ + + ⋅ ∇ + ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ B u u B B B u B u B B Bu uB BB B B I uu u B u 2 2 p E p E t ρ ρ ρ ρ

von Karman Institute for Fluid Dynamics

POWELL SOURCE TERM (3)

Problem with Non-Conservative Nature of Powell Source Term

• Jump conditions across the discontinuities can be predicted

wrongly in flows involving strong shocks

• Consequently, the scheme can produce incorrect shock

positions and speeds Problem with Divergence Wave Equation

• divB errors cannot be swept away in stagnant regions of the

flow

von Karman Institute for Fluid Dynamics

BALSARA’S SCHEME

• discrete satisfaction of the divB = 0 condition for the magnetic field
• hydrodynamic state variables at the cell center, magnetic field at the

cell faces and electric field at the cell edges (i.e. staggered grid approach)

• divB = 0 should be imposed as an initial condition in the whole

domain

• once imposed, divB = 0 will automatically be satisfied in time
• divB = 0 condition is satisfied up to the machine accuracy level

Implementation is exceedingly complex

von Karman Institute for Fluid Dynamics

PROJECTION SCHEME

• magnetic field is projected onto a discrete space in which divB = 0 is

satisfied in a numerical sense:

• divB = 0 should be imposed as an initial condition in the whole

domain

• once imposed, divB = 0 will automatically be satisfied in time
• divB = 0 condition is satisfied up to the machine accuracy level

Mixture of hyperbolic and elliptic numerical methods Computationally expensive (necessity of a Poisson solver)

( ) ( )

n n

φ φ ∇ − ← = ∇ − ⋅ ∇ = ⋅ ∇

+ ∗

B B B B

1

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (1)

• Hyperbolic system of 9 equations with eigenvalues: u, u±cf, u±cs,

u±cA, ±βref

• is a scalar potential function, βref is a reference velocity (e.g. inlet

flow velocity)

• Unlike the classical projection scheme approach, in this case, the

system is purely hyperbolic and can be solved with any conventional upwind finite volume numerical scheme suitable for hyperbolic PDE systems (e.g. Lax-Friedrichs scheme)

( ) ( )

) 2 ( 2

2

= ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ − ⋅ + + + − − ⋅ + + ⋅ ∇ + ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ B B u B u B B I Bu uB BB B B I uu u B u

ref

p E p E t β φ ρ ρ φ ρ ρ

φ

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (2)

( )

1

2

= ⋅ ∇ + ∂ ∂ u ρ β t p

ref

Incompressible CFD Modified System of Ideal MHD Equations

1

2

= ⋅ ∇ + ∂ ∂ B t

ref

φ β

Incompressibility Condition Solenoidal Constraint

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (3)

( )

B u B × × ∇ = ∇ + ∂ ∂ φ t

φ

taking the divergence, differentiating w.r.t. time and substituting for divB from the evolution of , we get Physical Interpretation Wave equation for divB

( ) ( )

[ ]

2 2 2

= ⋅ ∇ ∇ ⋅ ∇ − ∂ ⋅ ∇ ∂ B B

ref

t β

Unlike the divergence wave of Powell’s approach, divB errors are radiated away with ±βref , wave speeds independent from the flow velocity, thereby guaranteeing the removal of divB errors even in stagnant regions of flow

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (4)

(1) Do the modified PDEs converge towards a correct solution of the

• riginal ideal MHD system ?

(2) Do the discretized equations converge to a correct weak solution of the original ideal MHD sysem ? Consistency Check

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (5)

( )

[ ]

2

= × × ∇ ⋅ ∇ = ∇ B u φ

boundaries

• ther

all at

• utlets

at = ∂ ∂ = n φ φ

= φ

For smooth flows Take the divergence of the modified induction equation, at steady state: Boundary conditions: By using separation of variables, the above Laplacian equation yields

• everywhere. Hence the artificial term in the induction

equation vanishes and we get the correct solution.

φ ∇

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (6)

R R L L

F F φ φ + = +

R L

φ φ =

For shocks Let’s consider only steady shocks for simplicity. The jump relations then become: The jump relations of the original ideal MHD system are satisfied by the hyperbolized projection form of the equations: put Couldn’t the modified form of the equations allow for entropy- satisfying unphysical shocks ?

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (7)

[ ] ( ) [ ]

B B u B × ∇ × ∇ − ∇ − × × ∇ = ∂ ∂ η φ t

η

φ ∇

Let’s consider the modified system with non-zero resistivity, . The modified induction equation will be: For smooth flows, taking the divergence of both sides of the above equation, for each value of the resistivity, with the previously defined boundary conditions, disappears and ACA converges towards the correct entropy-satisfying solution of the original system of ideal MHD equations.

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (8)

For shocks, the entropy equation for the modified system is derived from the momentum and energy equations, which are the same in the

• riginal ideal MHD system, as follows:

OR Entropy jump conditions across a shock for original ideal MHD system are derived without using the induction equation in which one of the two modifications are made in the original system to obtain the modified system; therefore correct entropy-satisfying weak solutions of the

• riginal ideal MHD system are preserved for the modified one. (Sutton

& Sherman, Engineering Magnetohydrodynamics, 1965, pp. 320-324)

( ) ( )

1

2 ≥

× ∇ − = − B η γ ρ ρ γ dt d p dt dp

Entropy term

von Karman Institute for Fluid Dynamics

ARTIFICIAL COMPRESSIBILITY APPROACH (9)

Finally, at the discrete level, we have a usual conservative cell-centered finite volume discretization which can be shown to convergence to a correct weak solution of the underlying PDEs on regular meshes. This is not the case with Powell’s method, for which the source term leads to consistency problems in shocks.

von Karman Institute for Fluid Dynamics

HYPERBOLIC DIVERGENCE CLEANING (1)

• Dedner et. al., JCP, 2002 proposed a solenoidal constraint satisfying

technique called Hyperbolic Divergence Cleaning (HDC).

• Mathematically, it is fully equivalent to ACA. (i.e. the system of

PDEs, the method of discretization are the same)

• ACA has a rigorous analytical consistency proof which clearly

indicates that the modified ideal MHD system boils down to the original system at convergence whereas such a consistency proof is not given for HDC.

• ACA gives the best performance with implicit schemes whereas such a

fact is not mentioned about HDC.

von Karman Institute for Fluid Dynamics

HYPERBOLIC DIVERGENCE CLEANING (2)

• Physical insight of the artificial compressibility and the inspiration

from the projection scheme constitute the basis of the modified system

• f PDEs of ACA and therefore of ACA itself whereas HDC is based on

a more mathematical basis.

• The idea of divB errors being removed with the two new characteristic

speeds is not clearly presented (i.e. a wave equation in terms of divB is not presented).

von Karman Institute for Fluid Dynamics

THE COOLFluiD FRAMEWORK

Objectives of COOLFluiD (Computational Object Oriented Library for Fluid Dynamics)

• Modular object oriented framework (under development at the

VKI)

• Efficient simulation of multi-physics (Euler, Navier-Stokes,

MHD, Reactive Flows ...)

• Unstructured grids
• Advanced numerical algorithms (RDS, FEM, FVM, ...)
• Life-Time of at least 8 years

von Karman Institute for Fluid Dynamics

THE STRUCTURE OF COOLFluiD

von Karman Institute for Fluid Dynamics

NUMERICAL SCHEMES (1)

( )

( )

( )⎟

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − − + + + =

+ + +

+ + +

i i i ref i i

φ φ λ β φ φ

1 2 1 max 2 1

2 1 2 2 1

1 i i 2 1 i 2 1 i

x x LF_Powell LF_ACA

B B F F

Lax-Friedrichs Scheme

( ) ( ) ( )

i 1 i i 1 i 2 1 i

Q Q Q F Q F F − − + =

+ + + + 2 1 max

2 1 2

i

λ

von Karman Institute for Fluid Dynamics

NUMERICAL SCHEMES (2)

( )

( )

( )

2 1 2 1

2 1 2 2 1

+ ± ± ± +

− ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + + + =

+ + +

i ref i i

ref ref ref 1 i i 2 1 i 2 1 i

β x x Roe_Powell Roe_ACA

K B B F F

β β

λ α β φ φ

Roe’s Scheme

( ) ( )

( )

∑

= + + +

− + =

8 1 2 1

2 1 2

l i l l l

K λ α

i 1 i 2 1 i

Q F Q F F

von Karman Institute for Fluid Dynamics

BOUNDARY CONDITIONS

In COOLFluiD, boundary conditions (BC) implemented for FVM MHD are

• Perfect Conducting Wall BC (normal velocity and magnetic field

components are cancelled)

• Superfast Inlet BC (Qghost = 2Q∞ - Qboundary)
• Superfast Outlet BC (Qghost = Qboundary)

Boundary cell Ghost cell Boundary face

von Karman Institute for Fluid Dynamics

SUPERSONIC FIELD-ALIGNED LOW-β FLOW OVER A PERFECTLY CONDUCTING QUARTER CYLINDER (2-D) (M∞ = 2.6, β = 0.4, # of Elements: 83008, Element Type: Triangular) Mesh in front of the quarter cylinder of radius 0.125 units and center located at the origin 2nd order TVD LF with Powell 2nd order TVD LF with ACA

von Karman Institute for Fluid Dynamics

SUPERSONIC FIELD-ALIGNED LOW-β FLOW OVER A PERFECTLY CONDUCTING QUARTER CYLINDER (2-D) (M∞ = 2.6, β = 0.4, # of Elements: 83008, Element Type: Triangular) 2nd order TVD LF with Powell (Zoomed) 2nd order TVD LF with ACA (Zoomed) Hans De Sterck’s Ph.D. Thesis, KU Leuven, 1999 2nd order TVD LF with Powell

von Karman Institute for Fluid Dynamics

SUPERSONIC FIELD-ALIGNED LOW-β FLOW OVER A PERFECTLY CONDUCTING QUARTER CYLINDER (2-D) (M∞ = 2.6, β = 0.4, # of Elements: 83008, Element Type: Triangular) Density variation along the stagnation line

von Karman Institute for Fluid Dynamics

SUPERSONIC FIELD-ALIGNED LOW-β FLOW OVER A PERFECTLY CONDUCTING QUARTER CYLINDER (2-D) (M∞ = 2.6, β = 0.4, # of Elements: 83008, Element Type: Triangular) 2nd order TVD LF with ACA Consistency Study with 1st order LF ACA

von Karman Institute for Fluid Dynamics

SUPERSONIC FIELD-ALIGNED LOW-β FLOW OVER A PERFECTLY CONDUCTING QUARTER CYLINDER (2-D) (M∞ = 2.6, β = 0.4, # of Elements: 17052, Element Type: Triangular) 1st order LF with Powell Backward Euler 1st order LF with ACA Backward Euler

von Karman Institute for Fluid Dynamics

SUPERSONIC FIELD-ALIGNED LOW-β FLOW OVER A PERFECTLY CONDUCTING QUARTER CYLINDER (2-D) (M∞ = 2.6, β = 0.4, # of Elements: 38208, Element Type: Triangular) 1st order LF with Powell & ACA CFL: 0.5 Forward Euler

von Karman Institute for Fluid Dynamics

CONCLUSION (1)

Properties of ACA

• Purely hyperbolic (unlike classical projection methods)
• Stable and consistent and works even in the vicinity of stagnant flow

regions (unlike Powell’s scheme)

• Easy to implement (unlike Balsara’s scheme)
• Poorer performance than Powell’s scheme in explicit schemes (due

to the reflecting superfast outlet BC) whereas comparable performance with Powell’s scheme in implicit schemes with the fully implicit treatment

von Karman Institute for Fluid Dynamics

CONCLUSION (2)

Properties of ACA (Cont.)

• Appears to be fully equivalent to Hyperbolic Divergence Cleaning

technique proposed by Dedner et. al., JCP, 2002. However, we believe that the derivation leading to our scheme and its practical implementation are more intuitive and straightforward.

von Karman Institute for Fluid Dynamics

END OF PRESENTATION

QUESTIONS SECTION