Construction of semi-implicit methods for magnetohydrodynamics based on multiple scale asymptotics

magnetohydrodynamics based on multiple scale asymptotics

Friedemann Kemm IAG Universit¨at Stuttgart kemm@iag.uni-stuttgart.de

Claus-Dieter Munz IAG Universit¨at Stuttgart munz@iag.uni-stuttgart.de

The non dimensional equations in primitive variables

ρ_t + uρ_x + ρu_x = 0 , u_t + uu_x + 1

M² p_x

ρ + 1 Av²

1

2ρ(B²)_x = 0 , v_t + uv_x − 1

Av² B₁

ρ B_2x = 0 , w_t + uw_x − 1

Av² B₁

ρ B_3x = 0 , B_2t + B₂u_x − B₁v_x + uB_2x = 0 , B_3t + B₃u_x − B₁w_x + uB_3x = 0 ,

p_t + γpu_x + up_x = 0

M := v₀

a₀ , Av := v₀ c_A,0 .

The eigenvalues

u , u ± c_s , u ± c_A , u ± c_f with

c_A = 1

Avb₁ , c_{f /s} = 1

√2 s

1

M²a² + 1

Av²b² ±

r 1

M²a² + 1

Av²b²²

− 4 · 1

M²a² 1

Av²b²₁ , the abbreviations b₁ = p

B₁²/ρ, b² = B²/ρ and the acoustic speed a = γp/ρ.

The asymptotic Expansion

For small Mach numbers:

q(x, t; M) = q⁽⁰⁾(x, ξ, t) + M q⁽¹⁾(x, ξ, t) + M²q⁽²⁾(x, ξ, t) + . . . with

q⁽ⁱ⁾ : R × R × [0,∞) → R^m

(x, ξ, t) 7→ q⁽ⁱ⁾(x, ξ, t)

(i = 0,1,2, . . .)

and

ξ = M x .

For small Alfven numbers:

q(x, t; Av) = q⁽⁰⁾(x, ξ, t) + Av q⁽¹⁾(x, ξ, t) + Av²q⁽²⁾(x, ξ, t) + . . . with

q⁽ⁱ⁾ : R × R × [0,∞) → R^m

(x, ξ, t) 7→ q⁽ⁱ⁾(x, ξ, t)

(i = 0,1,2, . . .)

and

ξ = Av x . With the restriction

q⁽ⁱ⁾(x, ξ, t)

x → 0 for x → ∞ ∀i, t .

The velocity equation

Orders -2 and -1:

p⁽⁰⁾_x = 0 , p⁽⁰⁾_ξ + p⁽¹⁾_x = 0 .

The velocity equation

Orders -2 and -1:

p⁽⁰⁾_x = 0 , p⁽⁰⁾_ξ + p⁽¹⁾_x = 0 .

0 = 1

|I| Z

I

p⁽⁰⁾_ξ (ξ, t) dx + 1

|I| Z

I

p⁽¹⁾_x (x, ξ, t) dx

= p⁽⁰⁾_ξ + 1

|I| Z

I

p⁽¹⁾_x (x, ξ, t) dx

= p⁽⁰⁾_ξ + 1

|I|

p^(1)b

a , and therefore

p⁽⁰⁾_ξ = 0 f¨ur |I| → ∞ .

The behavior and the meaning of the pressure variables

The meaning:

p⁽⁰⁾= p⁽⁰⁾(t) thermodynamic , p⁽¹⁾= p⁽¹⁾(ξ, t) acoustic ,

p⁽²⁾= p⁽²⁾(x, ξ, t) incompressible .

The behavior and the meaning of the pressure variables

The meaning:

p⁽⁰⁾= p⁽⁰⁾(t) thermodynamic , p⁽¹⁾= p⁽¹⁾(ξ, t) acoustic ,

p⁽²⁾= p⁽²⁾(x, ξ, t) incompressible .

The time behavior of p⁽⁰⁾:

p⁽⁰⁾_t + 1

|I|γp⁽⁰⁾ u^b

a = 0 .

The long wave length acoustic:

(ρu)_t + p⁽¹⁾_ξ = 0 , p⁽¹⁾_t + γp⁽⁰⁾u_ξ = 0 .

The zero Mach number limit

The pressure is the only variable, whose higher order terms we must take into consideration. So we skip the upper indices of the others.

ρ_t + ρu_x + ρ_xu = 0 , u_t + uu_x + 1

ρp⁽²⁾_x + 1

2Av²ρ(B²)_x = −1

ρp⁽¹⁾_ξ , v_t + uv_x − 1

Av²ρB₁B_2x = 0 , w_t + uw_x − 1

Av²ρB₁B_3x = 0 , B_2t + B₂u_x − B₁v_x + uB_2x = 0 , B_3t + B₃u_x − B₁w_x + uB_3x = 0 , p⁽⁰⁾_t + γp⁽⁰⁾u_x = 0 .

The behavior of the magnetic induction variables

The meaning:

B_i⁽⁰⁾ = B_i⁽⁰⁾(t) , B_i⁽¹⁾ = B_i⁽¹⁾(ξ, t) ,

B_i⁽²⁾ = B_i⁽²⁾(x, ξ, t) . (i = 2,3)

The behavior of the magnetic induction variables

The meaning:

B_i⁽⁰⁾ = B_i⁽⁰⁾(t) , B_i⁽¹⁾ = B_i⁽¹⁾(ξ, t) ,

B_i⁽²⁾ = B_i⁽²⁾(x, ξ, t) . (i = 2,3)

The time behavior of the leading order terms:

B_2,t⁽⁰⁾ + B₂⁽⁰⁾ 1

|I|

u^b

a − B₁ 1

|I|

v^b

a = 0 , B_3,t⁽⁰⁾ + B₃⁽⁰⁾ 1

|I|

u^b

a − B₁ 1

|I|

w^b

a = 0 .

Evolution equations for long wave magnetic phenomena:

(ρu)_t + B₂⁽⁰⁾B_2,ξ⁽¹⁾ + B₃⁽⁰⁾B_3,ξ⁽¹⁾

= 0 , (ρv)_t − B₁B_2,ξ⁽¹⁾ = 0 , (ρw)_t − B₁B_3,ξ⁽¹⁾ = 0 , B_2,t⁽¹⁾ + B₂⁽⁰⁾u_ξ − B₁v_ξ = 0 , B_3,t⁽¹⁾ + B₃⁽⁰⁾u_ξ − B₁v_ξ = 0 .

The zero Alfven number limit

ρ_t + ρu_x + ρ_xu = 0 , u_t + uu_x + 1

M² 1

ρp_x + 1

ρ(B²)⁽²⁾_x = −1

ρ(B²)⁽¹⁾_ξ , v_t + uv_x − 1

ρB₁B_2,x⁽²⁾ = 1

ρB₁B_2,ξ⁽¹⁾ , w_t + uw_x − 1

ρB₁B_3,x⁽²⁾ = 1

ρB₁B_3,ξ⁽¹⁾ , B_2,t⁽⁰⁾ + B₂⁽⁰⁾u_x − B₁v_x = 0 ,

B_3,t⁽⁰⁾ + B₃⁽⁰⁾u_x − B₁w_x = 0 , p_t + γpu_x + p_xu = 0 .

The hyperbolic-elliptic splitting

ρ_t + uρ_x + ρu_x = 0 , u_t + uu_x + 1

M² p_x

ρ + 1

2Av²ρ(B²)_x = 0 , v_t + uv_x − B₁

Av²ρB_2x = 0 , w_t + uw_x − B₁

Av²ρB_3x = 0 , B_2t + B₂u_x − B₁v_x + uB_2x = 0 , B_3t + B₃u_x − B₁w_x + uB_3x = 0 ,

p_t + γpu_x + up_x = 0

M := v₀

a₀ , Av := v₀ c_A,0 .

A semi-implicit method;

System with fast waves

For low Mach numbers: Discretization with backward differences in time:

ρⁿ⁺¹ − ρⁿ

∆t + ρⁿ⁺¹uⁿ⁺¹_x = 0 , uⁿ⁺¹ − uⁿ

∆t + 1 M²

1

ρⁿ⁺¹pⁿ⁺¹_x = 0 , pⁿ⁺¹ − pⁿ

∆t + γpⁿ⁺¹uⁿ⁺¹_x = 0 .

The scheme:

Step 1 Choose estimates ρ^∗ and p^∗ for density and pressure.

The scheme:

Step 1 Choose estimates ρ^∗ and p^∗ for density and pressure.

Step 2 Compute an estimate for the velocity via u^∗ = uⁿ − ∆t

M²ρ^∗p^∗_x .

Step 3 If the estimates are good enough → ready. If not then make correction Ansatz:

uⁿ⁺¹ = u^∗ + δu , pⁿ⁺¹ = p^∗ + δp , ρⁿ⁺¹ = ρ^∗ + δρ write δu in terms of δp, ρ^∗, M and ∆t:

δu = − ∆t

M²ρ^∗δp_x and put this into the pressure equation

pⁿ⁺¹ − pⁿ

∆t + γp^∗uⁿ⁺¹_x = 0 .

This gives the “elliptic” condition

δp − ∆t²

M² γp^∗

1

ρ^∗δp_x

x

= pⁿ − p^∗ − ∆t γp^∗u^∗_x

for the pressure correction. Solve this equation numerically.

This gives the “elliptic” condition

δp − ∆t²

M² γp^∗

1

ρ^∗δp_x

x

= pⁿ − p^∗ − ∆t γp^∗u^∗_x

for the pressure correction. Solve this equation numerically.

Step 4 Compute corrected pressure, velocity and density. Use these values as the new estimates and return to Step 3.

For very small Mach numbers put the multiple pressure variables into the scheme.

The “elliptic” condition then gets the form

δp⁽²⁾ − ∆t²

M² γp^∗

1

ρ^∗δp⁽²⁾_x

x

= pⁿ − p^∗ − ∆t γp^∗u^∗_x . The corrected pressure is then computed by

pⁿ⁺¹ = p^∗ + M²δp⁽²⁾ .

For low Alfven numbers: Analogous precedure. Now the “elliptic” condition is δB₂ − ∆t²

Av²

B₂^∗(B₂^∗ + δB₂) + B₁²

ρ^∗ δB_2,x

x

− ∆t² Av²

B₂^∗(B₃^∗ + δB₃)

ρ δB_3,x

x

− ∆t²

Av²B₂^∗ δB₂B_2,x^∗ + δB₃B_3,x^∗

= B₂ⁿ − B₂^∗ − ∆t B₂^∗u^∗_x − B₁v_x^∗ , δB₃ − ∆t²

Av²

B₃^∗(B₃^∗ + δB₃) + B₁²

ρ^∗ δB_3,x

x

− ∆t² Av²

B₃^∗(B₂^∗ + δB₂)

ρ δB_2,x

x

− ∆t²

Av²B₃^∗ δB₃B_3,x^∗ + δB₂B_2,x^∗

= B₃ⁿ − B₃^∗ − ∆t B₃^∗u^∗_x − B₁v_x^∗ .

If we use multiple induction variables it is

δB₂⁽²⁾ − ∆t²

B₂^∗B₂^∗ + B₁²

ρ^∗ δB_2,x⁽²⁾

x

− ∆t²

B₂^∗B₃^∗

ρ δB_3,x⁽²⁾

x

− ∆t²B₂^∗ δB₂⁽²⁾B_2,x^∗ + δB₃⁽²⁾B_3,x^∗

= B₂ⁿ − B₂^∗ − ∆t B₂^∗u^∗_x − B₁v_x^∗ , δB₃⁽²⁾ − ∆t²

B₃^∗B₃^∗ + B₁²

ρ^∗ δB_3,x⁽²⁾

x

− ∆t²

B₃^∗B₂^∗

ρ δB_2,x⁽²⁾

x

− ∆t²B₃^∗ δB₃⁽²⁾B_3,x^∗ + δB₂⁽²⁾B_2,x^∗

= B₃ⁿ − B₃^∗ − ∆t B₃^∗u^∗_x − B₁v_x^∗ .

In the case of multiple induction variables the corrected magnetic field is computed by

B₂ⁿ⁺¹ = B₂^∗ + Av²δB₂⁽²⁾ , B₃ⁿ⁺¹ = B₃^∗ + Av²δB₃⁽²⁾ .

Conclusions

• Non dimensional equations

Conclusions

• Non dimensional equations

• Multiple scale asymptotics (one space dimension)

Conclusions

• Non dimensional equations

• Multiple scale asymptotics (one space dimension)

• Elliptic-hyperbolic splitting

Conclusions

• Non dimensional equations

• Multiple scale asymptotics (one space dimension)

• Elliptic-hyperbolic splitting

• Semi-implicit method:

Conclusions

• Non dimensional equations

• Multiple scale asymptotics (one space dimension)

• Elliptic-hyperbolic splitting

• Semi-implicit method:

? slow waves explicit

Conclusions

• Non dimensional equations

• Multiple scale asymptotics (one space dimension)

• Elliptic-hyperbolic splitting

• Semi-implicit method:

? slow waves explicit

? fast waves implicit

Conclusions

• Non dimensional equations

• Multiple scale asymptotics (one space dimension)

• Elliptic-hyperbolic splitting

• Semi-implicit method:

? slow waves explicit

? fast waves implicit

• correction of pressure or magnetic induction via elliptic equation