A unified moving grid method for

(1)

A unified moving grid method for

hyperbolic systems of

partial differential equations

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

Department of Applied Mathematics, Ta-Tung University, April 23, 2009 – p. 1/73

(2)

Talk objective

Describe a unified-coordinate moving grid approach for numerical approximation of first-order hyperbolic system

∂

∂tq (~x, t) +

N

X

j=1

A_j ∂q

∂x_j = 0 or ∂

∂tq (~x, t) +

N

X

j=1

∂

∂x_j f_j (q, ~x) = 0

with discontinuous initial data in general N ≥ 1 geometry

~x = (x₁, x₂, . . . , x_N): spatial vector, t: time q ∈ lR^m: vector of m state quantities

A_j ∈ lR^m×m: m × m matrix, f_j ∈ lR^m: flux vector

Model is assumed to be hyperbolic, where ^P^N_j=1 α_jA_j or

PN

j=1 α_j (∂f_j/∂q) is diagonalizable with real e-values, α_j ∈ lR

(3)

Talk outline

Preliminary

Sample models in Cartesian coordinates Cartesian cut-cell method & results

Mathematical model in unified coordinates Basic physical equations

Moving grid condition & geometric conservation law Finite volume approximation

Riemann problem & approximate solution Godunov-type method

Numerical examples Future work

(4)

Preliminary: model equations

Acoustics in heterogeneous media

∂

∂t





 p u₁ u₂





 +







0 K 0

1/ρ 0 0

0 0 0







∂

∂x₁





 p u₁ u₂





 +







0 0 K

0 0 0

1/ρ 0 0







∂

∂x₂





 p u₁ u₂







= 0

Shallow water equations with bottom topography

∂

∂t



 h hu_i



 +

N

X

j=1

∂

∂x_j





hu_j

hu_iu_j + ¹₂gh²δ_ij



 =





0

−gh_∂x^∂B

i



, i = 1, . . . , N

p: pressure, ρ: density, K: bulk modulus, u_i: x_i-velocity h: water height, δ_ij: Kronecker delta, B: bottom topo.

g: gravitational constant

(5)

Model equations (Cont.)

Compressible Euler equations

∂

∂t





 ρ ρu_i

E





 +

N

X

j=1

∂

∂x_j







ρu_j

ρu_iu_j + pδ_ij Eu_j + pu_j







= 0, i = 1, . . . , N

E = ρe + ρ P_N

j=1 u²_j/2: total energy, e(ρ, p): internal energy Note constitutive law for p is required to complete the

model, for example,

Polytropic gas: p = (γ − 1)ρe

Stiffened gas: p = (γ − 1)ρe − γB

van der Waals gas: p = _1−bρ^γ−1 ρe + aρ² − aρ²

(6)

Model equations (Cont.)

Compressible reduced 2-phase flow model

Proposed by Murrone & Guillard (JCP 2005)

Derive from Baer & Nunziato’s model by assuming 1-pressure & 1-velocity across interfaces

∂

∂t







α₁ρ₁ α₂ρ₂ ρu_i

E





 +

N

X

j=1

∂

∂x_j







α₁ρ₁u_j α₂ρ₂u_j ρu_iu_j + pδ_ij

Eu_j + pu_j







= 0, i = 1, . . . , N

∂α₁

∂t +

N

X

j=1

u_j ∂α₁

∂x_j = α₁α₂ ρ₁c²₁ − ρ₂c²₂ P2

k=1 α_kρ_kc²_k

! _N X

j=1

∂u_j

∂x_j

α_k: volume fraction for phase k, α1 + α2 = 1, c_k:

sound speed, ρ = α1ρ1 + α2ρ2: mixture (total) density

(7)

Reduced 2-phase model (Cont.)

Mixture equation of state: p = p(α2, α1ρ1, α2ρ2, ρe) Isobaric closure: p1 = p2 = p

For a class of EOS, explicit formula for p is available For some complex EOS, from (α2, ρ1, ρ2, ρe) in model equations we recover p by solving

p₁(ρ₁, ρ₁e₁) = p₂(ρ₂, ρ₂e₂) &

2

X

k=1

α_kρ_ke_k = ρe

Shyue (JCP 1998) & Allaire ^{et al.} (JCP 2002) proposed

∂α₁

∂t +

N

X

j=1

u_j ∂α₁

∂x_j = 0

(8)

Preliminary: model problem

Moving cylindrical vessel with ~u_b = (1, 0)

initial condition interface

media 1 media 2

(9)

Model problem: grid system

Typical discrete grid systems for cylindrical vessel

−0.5 0 0.5

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8

−1 −0.5 0 0.5 1

−1

−0.5 0 0.5

1 Cartesian grid Quadrilateral grid

x1

x 2

x2

(10)

Model problem: Cartesian results

Compressible flow case with air-helium interface

initial condition interface

air helium

(11)

Model problem: Cartesian results

Compressible flow case with air-helium interface Solution at time t = 0.25

(12)

Model problem: Cartesian results

(13)

Model problem: Cartesian results

(14)

Model problem: Cartesian results

Compressible flow case with air-helium interface Solution at time t = 1

(15)

Cartesian cut-cell method

Finite volume formulation of wave propagation method, Qⁿ_S gives approximate value of cell average of solution q over cell S at time t_n

Qⁿ_S ≈ 1 M(S)

Z

S

q(X, t_n) dV

M(S): measure (area in 2D or volume in 3D) of cell S

C E

D F

G H

(16)

Cartesian cut-cell method (Cont.)

First order version: Piecewise constant wave update Godunov-type method: Solve Riemann problem at each cell interface in normal direction & use resulting waves to update cell averages

Qⁿ⁺¹_S := Qⁿ⁺¹_S −M (W_p ∩ S)

M(S) R_p, R_p being jump from RP

↓

W_p W_p

(17)

Cartesian cut-cell method (Cont.)

First order version: Transverse-wave included

Use transverse portion of equation, solve Riemann problem in transverse direction, & use resulting

waves to update cell averages as usual

Stability of method is typically improved, while conservation of method is maintained

↓ ↓

W_pq W_pq

(18)

Cartesian cut-cell method (Cont.)

High resolution version: Piecewise linear wave update wave before propagation after propagation

a) b)

c) d)

α_pr

p/2

α_pr

p/2

λ_p∆ t

(19)

Embedded boundary conditions

For tracked segments representing rigid (solid wall)

boundary (stationary or moving), reflection principle is used to assign states for fictitious subcells in each time step:

z_C := z_E (z = ρ, p, α)

~u_C := ~u_E − 2(~u_E · ~n)~n + 2(~u_b · ~n)

~u_b: moving boundary velocity

C

E F

G H

~n

(20)

Cylinder lift-off problem

Moving speed of cylinder is governed by Newton’s law Pressure contours are shown with a 1000 × 200 grid

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15 0.2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15 0.2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15 0.2

t = 0

t = 0.1641s

t = 0.30085s

(21)

Shock-Bubble Interaction Problem

Cartesian grid results

(22)

Shock-Bubble Interaction Problem

(23)

Shock-Bubble Interaction Problem

(24)

Shock-Bubble Interaction Problem

(25)

Shock-Bubble Interaction Problem

(26)

Shock-Bubble Interaction Problem

(27)

Shock-Bubble Interaction Problem

(28)

Shock-Bubble Interaction Problem

(29)

Shock-Bubble Interaction Problem

(30)

Shock-Bubble Interaction Problem

(31)

Shock-Bubble Interaction (cont.)

Approximate locations of interfaces

time=55µs

air R22

time=115µs time=135µs

time=187µs time=247µs time=200µs

time=342µs time=417µs time=1020µs

(32)

Cartesian cut-cell: Remarks

Small cell problems Stability

Accuracy

Numerical implementation

Challenging task for embedded 3D geometry

Challenging task for interface tracking in general geometry (even in 2D)

(33)

Cartesian cut-cell: Remarks

Small cell problems Stability

Accuracy

Numerical implementation

Challenging task for embedded 3D geometry

Challenging task for interface tracking in general geometry (even in 2D)

This work is aimed at devising a more robust moving grid method than the aforementioned Cartesian cut-cell method

To begin with, take unified coordinate method proposed by Hui & coworkers (JCP 1999, 2001)

(34)

Model system in unified coord.

To begin with, consider a general non-rectangular domain Ω (N = 2 shown below) & introduce coordinate change

(~x, t) 7→ (~ξ, τ ) via

ξ = (ξ~ ₁, ξ₂, . . . , ξ_N) , ξ_j = ξ_j(~x, t), τ = t,

that maps a physical domain Ω to a logical one Ω^ˆ

−1 0 1

−1.5

−1

−0.5 0

−1 −0.5 0 0.5

−2

−1.5

−1

−0.5

x₁ x2

ξ₁

ξ2

Ω −→ Ωˆ

mapping

ξ₁ = ξ₁(x₁, x₂) ξ₂ = ξ₂(x₁, x₂)

logical domain physical domain

(35)

Unified coord. model (Cont.)

To derive hyperbolic conservation laws, for example, in this generalized coordinate (~ξ, τ ), using chain rule of partial

differentiation, derivatives in physical space become

∂

∂t = ∂

∂τ +

N

X

i=1

∂ξ_i

∂t

∂

∂ξ_i, ∂

∂x_j =

N

X

i=1

∂ξ_i

∂x_j

∂

∂ξ_i for j = 1, 2, . . . , N ,

yielding the equation

∂q

∂τ +

N

X

j=1

∂ξ_j

∂t

∂q

∂ξ_j +

N

X

i=1

∂ξ_i

∂x_j

∂f_j

∂ξ_i

!

= 0

Note this is not in divergence form, and hence is not conservative.

(36)

Unified coord. model (Cont.)

To obtain a strong conservation-law form as

∂ ˜q

∂τ +

N

X

j=1

∂ ˜f_j

∂ξ_j = ˜ψ

for some q˜, f^˜_j, & ψ^˜, we first multiply J = det

∂~ξ/∂~x⁻1

to the aforementioned non-conservative equations, and have

J ∂q

∂τ +

N

X

j=1

J ∂ξ_j

∂t

∂q

∂ξ_j +

N

X

i=1

∂ξ_i

∂x_j

∂f_j

∂ξ_i

!

= Jψ(q)

Then use differentiation by parts, u dv = d(uv) − v du, yielding

∂ ˜q

∂τ +

N

X

j=1

∂ ˜f_j

∂ξ_j = ˜ψ + G with q = Jq˜ , f˜_j = J

q^∂ξ_∂t^j + PN

k=1 f_k _∂x^∂ξ^j

k

, ψ = Jψ,˜ & G (see next)

(37)

Unified coord. model (Cont.)

Here we have

G = q





∂J

∂τ +

N

X

j=1

∂

∂ξ_j

J ∂ξ_j

∂t



 +

N

X

j=1

f_j

" _N X

k=1

∂

∂ξ_k

J ∂ξ_k

∂x_j

#

With the use of basic grid-metric relations, it is known that

∂J

∂τ +

N

X

j=1

∂

∂ξ_j

J ∂ξ_j

∂t

= 0 (geometric conservation law)

N

X

k=1

∂

∂ξ_k

J ∂ξ_k

∂x_j

= 0 ∀ j = 1, 2, . . . , N (compatibility condition)

and hence G = 0

(38)

Unified coord. model (Cont.)

Shallow water equations

∂

∂τ





hJ hJu_i



 +

N

X

j=1

∂

∂ξ_j J





hU_j

hu_iU_j + ¹₂gh²δ_ij ^∂ξ_∂x^j

i



 =





0

−ghJ _∂x^∂B

i





Compressible Euler equations

∂

∂τ







ρJ ρJu_i

JE





 +

N

X

j=1

∂

∂ξ_j J







ρU_j

ρu_iU_j + p^∂ξ_∂x^j

i

EU_j + pU_j − p^∂ξ_∂t^j







=







0

−ρJ _∂x^∂φ

i

−ρJ~u · ∇φ







U_j = ∂_tξ_j + PN

i=1 u_i∂_x_iξ_j: contravariant velocity in ξ_j-direction

φ: gravitational potential

(39)

Unified coord.: Geometric claw

With non-trivial ∂_τ~x, we should impose conditions on grid metrics ∂_t~ξ & ∇_~_xξ~ to have the fulfillment of geometrical conservation law

∂J

∂τ +

N

X

j=1

∂

∂ξ_j

J ∂ξ_j

∂t

= 0

Here we are interested in an approach that is based on the compatibility condition of ∂_τ∂_ξ_jx_i & ∂_ξ_j∂_τx_i, ^i.e.,

∂

∂τ

∂x_i

∂ξ_j

+ ∂

∂ξ_j

−∂x_i

∂τ

= 0 for i, j = 1, 2, . . . , N.

for unknowns ∂x_i/∂ξ_j, yielding easy computation of J & ∇ξ_j

(40)

Unified coord.: Grid movement

For fluid flow problems, to improve numerical resolution of interfaces (material or slip lines), it is popular to take ∂_τ~x as

Lagrangian case: ∂_τ~x = ~u (flow velocity)

Lagrangian-like case: ∂_τ~x = h0~u (pseudo velocity) h0 ∈ [0, 1] (fixed piecewise const.)

Unified coordinate case: ∂_τ~x = h~u

h ∈ [0, 1] but is determined from a PDE constraint arising from such as grid-angle or grid-Jacobian preserving condition

ALE-like case: ∂_τ~x = ~U (arbitrary velocity)

(41)

Unified coord.: Grid movement

For fluid flow problems, to improve numerical resolution of interfaces (material or slip lines), it is popular to take ∂_τ~x as

Lagrangian case: ∂_τ~x = ~u (flow velocity)

Lagrangian-like case: ∂_τ~x = h0~u (pseudo velocity) h0 ∈ [0, 1] (fixed piecewise const.)

Unified coordinate case: ∂_τ~x = h~u

h ∈ [0, 1] but is determined from a PDE constraint arising from such as grid-angle or grid-Jacobian preserving condition

ALE-like case: ∂_τ~x = ~U (arbitrary velocity)

For simplicity, we will focus on Lagrangian-like case

(42)

Unified coord. model: Summary

With ∂_τ~x = h0~u, unified coordinate model for single component compressible flow problems consists of

Physical balance laws

∂

∂τ







ρJ ρJu_i

JE





 +

N

X

j=1

∂

∂ξ_j J







ρU_j

ρu_iU_j + p^∂ξ_∂x^j

i

EU_j + pU_j − p^∂ξ_∂t^j







=







0

−ρJ _∂x^∂φ

i

−ρJ~u · ∇φ







Geometrical conservation laws

∂

∂τ

∂x_i

∂ξ_j

+ ∂

∂ξ_j

−∂x_i

∂τ

= 0 for i, j = 1, 2, . . . , N.

pressure law p(ρ, e)

(43)

Unified coord. model: Remarks

For unified coordinate models mentioned above, it is known that

when h0 = 0 (Eulerian case), the model is hyperbolic when h0 = 1 (Lagrangian case), the model is weakly hyperbolic (do not possess complete eigenvectors) when h0 ∈ (0, 1) (Lagrangian-like case), the model is hyperbolic

(44)

Unified coord. model: Remarks

For unified coordinate models mentioned above, it is known that

when h0 = 0 (Eulerian case), the model is hyperbolic when h0 = 1 (Lagrangian case), the model is weakly hyperbolic (do not possess complete eigenvectors) when h0 ∈ (0, 1) (Lagrangian-like case), the model is hyperbolic

If a prescribed velocity ~u_b for a rigid body motion is included in the formulation ^i.e., with ∂_τ~x = h⁰~u + ~u_b, we should be able to use the model to solve some moving body problems as well.

(45)

Unified coord. model: Review

The work presented here is related to

W.H. Hui ^{et al.} (JCP 1999, 2001): Unified coordinated system for Euler equations

W.H. Hui (Comm. Phys. Sci. 2007): Unified coordinate system in CFD

C. Jin & K. Xu (JCP 2007): Moving grid gas-kinetic method for viscous flow

P. Jia ^{et al.} (Computers and Fluids 2006) Unified

coordinated system for compressible milti-material flow Z. Chen ^{et al.} (Int J. Numer. Meth Fluids 2007): Wave speed based moving coordinates for compressible flow equations

(46)

Finite volume approximation

Employ finite volume formulation of numerical solution

Qⁿ_ijk ≈ 1

∆ξ₁∆ξ₂∆ξ₃ Z

Cijk

q(ξ₁, ξ₂, ξ₃, τ_n) dV

that gives approximate value of cell average of solution q over cell C_ijk at time τ_n (sample case in 2D shown below)

i − 1 i − 1

i

i j

j

j + 1 j + 1

C_ij Cˆ_ij

ξ₁ ξ₂

mapping

∆ξ₁

∆ξ₂ logical domain physical domain

←−

x₁ = x₁(ξ₁, ξ₂) x₂ = x₂(ξ₁, ξ₂) x₁

x₂

(47)

Finite volume (Cont.)

In three dimensions N = 3, equations to be solved take

∂

∂τ q ~ξ, τ + X^N

j=1

∂

∂ξ_j f_j

q, ~ξ

= ψ

q, ~ξ

A simple dimensional-splitting method based on f-wave approach of LeVeque ^{et al.} is used for approximation, ^i.e.,

Solve one-dimensional Riemann problem normal at each cell interfaces

Use resulting jumps of fluxes (decomposed into each wave family) of Riemann solution to update cell

averages

Introduce limited jumps of fluxes to achieve high resolution

(48)

Finite volume (Cont.)

Basic steps of a dimensional-splitting scheme ξ1-sweeps: solve

∂q

∂τ + f₁ ∂

∂ξ, q, ∇~ξ

= 0 updating Qⁿ_ijk to Q^∗_ijk

ξ2-sweeps: solve

∂q

∂τ + f₂

∂

∂ξ₂, q, ∇~ξ

= 0 updating Q^∗_ijk to Q^∗∗_ijk

ξ3-sweeps: solve

∂q

∂τ + f₃

∂

∂ξ₃ , q, ∇~ξ

= 0 updating Q^∗∗_ijk to Qⁿ⁺¹_ijk

(49)

Finite volume (Cont.)

Consider ξ₁-sweeps, for example, First order update is

Q^∗_ijk = Qⁿ_ijk − ∆τ

∆ξ₁

h A⁺₁ ∆Qn

i−1/2,jk + A⁻₁ ∆Qn

i+1/2,jk

i

with the fluctuations

(A⁺₁ ∆Q)ⁿ_i−1/2,jk = X

m:(λ1,m)ⁿ_i−1/2,jk>0

(Z_1,m)ⁿ_i−1/2,jk

and

(A⁻₁ ∆Q)ⁿ_i+1/2,jk = X

m:(λ1,m)ⁿ_i+1/2,jk<0

(Z_1,m)ⁿ_i+1/2,jk

(λ_1,m)ⁿ_ι−1/2,jk & (Z_1,m)ⁿ_ι−1/2,jk are in turn wave speed and f-waves for the mth family of the 1D Riemann problem solutions

(50)

Finite volume (Cont.)

High resolution correction is

Q^∗_ijk := Q^∗_ijk − ∆τ

∆ξ₁

˜F₁n

i+1/2,jk − ˜F₁n

i−1/2,jk

with ( ˜F₁)ⁿ_i−1/2,jk = 1 2

mw

X

m=1

sign(λ_1,m)

1 − ∆τ

∆ξ₁ |λ_1,m|

Z˜_1,m

n

i−1/2,jk

Z˜_ι,m is a limited value of Z_ι,m

It is clear that this method belongs to a class of upwind schemes, and is stable when the typical CFL (Courant-Friedrichs-Lewy) condition:

ν = ∆τ max_m (λ_1,m, λ_2,m, λ_3,m)

min (∆ξ₁, ∆ξ₂, ∆ξ₃) ≤ 1,

(51)

Finite volume: Riemann problem

Riemann problem for our model equations at cell interface ξ_i−1/2 consists of the equation











∂q_i−1,jk

∂τ + f₁

∂

∂ξ₁ , q_i−1,jk

= 0 if ξ₁ < (ξ₁)_i−1/2,

∂q_ijk

∂τ + f₁

∂

∂ξ₁ , q_ijk

= 0 if ξ₁ > (ξ₁)_i−1/2,

together with piecewise constant initial data

q(ξ₁, 0) =







Qⁿ_i−1,jk for ξ < ξ_i−1/2 Qⁿ_ijk for ξ > ξ_i−1/2

q_ijk = q|_(∂_ξ2_~_{x, ∂}_ξ3_~_x)_ijk & f₁(∂_ξ₁, q_ijk) = f₁(∂_ξ₁, q)|_(∂_ξ2_~_{x, ∂}_ξ3_~_x)_ijk

(52)

Riemann problem

Riemann problem at time τ = 0

∂_τq_L + f₁ (∂_ξ₁, q_L) = 0 ∂_τq_R + f₁ (∂_ξ₁, q_R) = 0 ξ₁ τ

Qⁿ_i−1,jk Qⁿ_ijk

(53)

Riemann problem

Exact Riemann solution: basic structure

(54)

Riemann problem

Shock-only approximate Riemann solution: basic structure

λ₀

λ₁ λ₂

λ₃ q_mL⁻ q_mL⁺

q_mR

Z¹ = f_L(q_mL⁻ ) − f_L(Qⁿ_i−1,jk) Z² = f_R(q_mR) − f_R(q_mL⁺ )

Z³ = f_R(Qⁿ_ijk) − f_R(q_mR)

(55)

Shock-only Riemann solver

Rotate velocity vector in Riemann data normal to each cell interface Find midstate velocity υ_m and pressure p_m by solving

φ(p_m) = υ_mR(p_m) − υ_mL(p_m) = 0

derived from Rankine-Hugoniot relation iteratively, where

υ_mL(p) = υ_L − p − p_L

M_L(p), υ_mR(p) = υ_R + p − p_R M_R(p)

Propagation speed of each moving discontinuity is determined by

(λ_1,1)_i−1/2,jk =

(1 − h₀)υ_m − M_L(p_m) ρ_mL(p_m)

∇_X_~ ξ₁

i−1/2,jk

(λ_1,2)_i−1/2,jk = (1 − h₀)υ_m

∇_X_~ξ₁

i−1/2,jk

(λ_1,3)_i−1/2,jk =

(1 − h₀)υ_m + M_R(p_m) ρ_mR(p_m)

∇_X_~ ξ₁

i−1/2,jk

(56)

Lax’s Riemann problem

Ideal gas EOS with γ = 1.4 h0 = 0 Eulerian result

h0 = 0.99 Lagrangian-like result

sharper resolution for contact discontinuity

0 0.5 1

0 0.5 1 1.5 2

x

ρ

Exact h0=0 h0=0.99

0 0.5 1

0 0.5 1 1.5 2

x

u

0 0.5 1

0 1 2 3 4

x

p

(57)

Lax’s Riemann problem

Physical grid coordinates at selected times

Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15

x

time

(58)

Woodward-Colella’s problem

Ideal gas EOS with γ = 1.4 h0 = 0 Eulerian result

0 0.5 1

0 2 4 6 8

ρ

0 0.5 1

−10 0 10 20

0 0.5 1

0 100 200 300 400

u p

t = 0.016 t = 0.016

t = 0.016

(59)

Woodward-Colella’s problem

h0 = 0 Eulerian result

0 0.5 1

0 5 10 15 20

ρ

Fine grid h0=0 h0=0.99

0 0.5 1

−5 0 5 10 15

0 0.5 1

0 200 400 600

u p

t = 0.032 t = 0.032

t = 0.032

(60)

Woodward-Colella’s problem

h0 = 0 Eulerian result

0 0.5 1

0 2 4 6 8

ρ

0 0.5 1

−5 0 5 10 15

0 0.5 1

0 200 400 600

u p

x x

x

t = 0.038 t = 0.038

t = 0.038

(61)

Woodward-Colella’s problem

Physical grid coordinates at selected times

Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.01 0.02 0.03 0.04

x

time

(62)

2 D Riemann problem

With initial 4-shock wave pattern

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1





 ρ u v p







=





 1.5

0 0 1.5













0.532 1.206

0 0.3













0.138 1.206 1.206 0.029













0.532 0 1.206

0.3







(63)

2 D Riemann problem

With initial 4-shock wave pattern Lagrangian-like result

Occurrence of simple Mach reflection

0.2 0.4 0.6 0.8 0.2

0.4 0.6 0.8

0.2 0.4 0.6 0.8 0.2

0.4 0.6 0.8

0.2 0.4 0.6 0.8 0.2

0.4 0.6 0.8

Density Pressure Physical grid

(64)

2 D Riemann problem

With initial 4-shock wave pattern Eulerian result

Poor resolution around simple Mach reflection

0.2 0.4 0.6 0.8 0.2

0.4 0.6 0.8

0.2 0.4 0.6 0.8 0.2

0.4 0.6 0.8

0.2 0.4 0.6 0.8 0.2

0.4 0.6 0.8

Density Pressure Physical grid