**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

**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_{1}, x_{2}, . . . , 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} α_{j}A_{j} or

PN

j=1 α_{j} (∂f_{j}/∂q) is diagonalizable with real e-values, α_{j} ∈ lR

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

**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

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

**Preliminary: model equations**

Acoustics in heterogeneous media

∂

∂t

p
u_{1}
u_{2}

+

0 K 0

1/ρ 0 0

0 0 0

∂

∂x_{1}

p
u_{1}
u_{2}

+

0 0 K

0 0 0

1/ρ 0 0

∂

∂x_{2}

p
u_{1}
u_{2}

= 0

Shallow water equations with bottom topography

∂

∂t

h
hu_{i}

+

N

X

j=1

∂

∂x_{j}

hu_{j}

hu_{i}u_{j} + ^{1}_{2}gh^{2}δ_{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

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

**Model equations (Cont.)**

Compressible Euler equations

∂

∂t

ρ
ρu_{i}

E

+

N

X

j=1

∂

∂x_{j}

ρu_{j}

ρu_{i}u_{j} + pδ_{ij}
Eu_{j} + pu_{j}

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

E = ρe + ρ P_{N}

j=1 u^{2}_{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ρ^{2} − aρ^{2}

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

**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

α_{1}ρ_{1}
α_{2}ρ_{2}
ρu_{i}

E

+

N

X

j=1

∂

∂x_{j}

α_{1}ρ_{1}u_{j}
α_{2}ρ_{2}u_{j}
ρu_{i}u_{j} + pδ_{ij}

Eu_{j} + pu_{j}

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

∂α_{1}

∂t +

N

X

j=1

u_{j} ∂α_{1}

∂x_{j} = α_{1}α_{2} ρ_{1}c^{2}_{1} − ρ_{2}c^{2}_{2}
P2

k=1 α_{k}ρ_{k}c^{2}_{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

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

**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_{1}(ρ_{1}, ρ_{1}e_{1}) = p_{2}(ρ_{2}, ρ_{2}e_{2}) &

2

X

k=1

α_{k}ρ_{k}e_{k} = ρe

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

∂α_{1}

∂t +

N

X

j=1

u_{j} ∂α_{1}

∂x_{j} = 0

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

**Preliminary: model problem**

Moving cylindrical vessel with ~u_{b} = (1, 0)

initial condition interface

media 1 media 2

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

**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

x1

x 2

x2

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

**Model problem: Cartesian results**

Compressible flow case with air-helium interface

initial condition interface

air helium

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

**Model problem: Cartesian results**

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

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

**Model problem: Cartesian results**

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

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

**Model problem: Cartesian results**

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

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

**Model problem: Cartesian results**

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

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

**Cartesian cut-cell method**

Finite volume formulation of wave propagation method, Q^{n}_{S}
gives approximate value of cell average of solution q over
cell S at time t_{n}

Q^{n}_{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

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

**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^{n+1}_{S} := Q^{n+1}_{S} −M (W_{p} ∩ S)

M(S) R_{p}, R_{p} being jump from RP

↓

↓

W_{p}
W_{p}

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

**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}

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

**Cartesian cut-cell method (Cont.)**

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

a) b)

c) d)

α_{p}r

p/2

α_{p}r

p/2

λ_{p}∆ t

λ_{p}∆ t

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

**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

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

**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

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**Shock-Bubble Interaction Problem**

Cartesian grid results

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

**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

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

**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)

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

**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)

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

**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

ξ = (ξ~ _{1}, ξ_{2}, . . . , ξ_{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_{1}
x2

ξ_{1}

ξ2

Ω −→ Ωˆ

mapping

ξ_{1} = ξ_{1}(x_{1}, x_{2})
ξ_{2} = ξ_{2}(x_{1}, x_{2})

logical domain physical domain

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

**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.

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

**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)

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

**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

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

**Unified coord. model (Cont.)**

Shallow water equations

∂

∂τ

hJ
hJu_{i}

+

N

X

j=1

∂

∂ξ_{j} J

hU_{j}

hu_{i}U_{j} + ^{1}_{2}gh^{2}δ_{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_{i}U_{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

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

**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 ∂_{τ}∂_{ξ}_{j}x_{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}

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

**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)

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

**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

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

**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_{i}U_{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)

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

**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

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

**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

^{0}~u + ~u

_{b}, we should be able to use the model to solve some moving body problems as well.

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

**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

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

**Finite volume approximation**

Employ finite volume formulation of numerical solution

Q^{n}_{ijk} ≈ 1

∆ξ_{1}∆ξ_{2}∆ξ_{3}
Z

Cijk

q(ξ_{1}, ξ_{2}, ξ_{3}, τ_{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}

ξ_{1}
ξ_{2}

mapping

∆ξ_{1}

∆ξ_{2}
logical domain
physical domain

←−

x_{1} = x_{1}(ξ_{1}, ξ_{2})
x_{2} = x_{2}(ξ_{1}, ξ_{2})
x_{1}

x_{2}

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

**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

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

**Finite volume (Cont.)**

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

∂q

∂τ + f_{1} ∂

∂ξ, q, ∇~ξ

= 0 updating Q^{n}_{ijk} to Q^{∗}_{ijk}

ξ2-sweeps: solve

∂q

∂τ + f_{2}

∂

∂ξ_{2}, q, ∇~ξ

= 0 updating Q^{∗}_{ijk} to Q^{∗∗}_{ijk}

ξ3-sweeps: solve

∂q

∂τ + f_{3}

∂

∂ξ_{3} , q, ∇~ξ

= 0 updating Q^{∗∗}_{ijk} to Q^{n+1}_{ijk}

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

**Finite volume (Cont.)**

Consider ξ_{1}-sweeps, for example,
First order update is

Q^{∗}_{ijk} = Q^{n}_{ijk} − ∆τ

∆ξ_{1}

h A^{+}_{1} ∆Qn

i−1/2,jk + A^{−}_{1} ∆Qn

i+1/2,jk

i

with the fluctuations

(A^{+}_{1} ∆Q)^{n}_{i−1/2,jk} = X

m:(λ1,m)^{n}_{i−1/2,jk}>0

(Z_{1,m})^{n}_{i−1/2,jk}

and

(A^{−}_{1} ∆Q)^{n}_{i+1/2,jk} = X

m:(λ1,m)^{n}_{i+1/2,jk}<0

(Z_{1,m})^{n}_{i+1/2,jk}

(λ_{1,m})^{n}_{ι−1/2,jk} & (Z_{1,m})^{n}_{ι−1/2,jk} are in turn wave speed and f-waves
for the mth family of the 1D Riemann problem solutions

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

**Finite volume (Cont.)**

High resolution correction is

Q^{∗}_{ijk} := Q^{∗}_{ijk} − ∆τ

∆ξ_{1}

˜F_{1}n

i+1/2,jk − ˜F_{1}n

i−1/2,jk

with ( ˜F_{1})^{n}_{i−1/2,jk} = 1
2

mw

X

m=1

sign(λ_{1,m})

1 − ∆τ

∆ξ_{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}, ∆ξ_{2}, ∆ξ_{3}) ≤ 1,

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

**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_{1}

∂

∂ξ_{1} , q_{i−1,jk}

= 0 if ξ_{1} < (ξ_{1})_{i−1/2},

∂q_{ijk}

∂τ + f_{1}

∂

∂ξ_{1} , q_{ijk}

= 0 if ξ_{1} > (ξ_{1})_{i−1/2},

together with piecewise constant initial data

q(ξ_{1}, 0) =

Q^{n}_{i−1,jk} for ξ < ξ_{i−1/2}
Q^{n}_{ijk} for ξ > ξ_{i−1/2}

q_{ijk} = q|_{(∂}_{ξ2}_{~}_{x, ∂}_{ξ3}_{~}_{x)}_{ijk} & f_{1}(∂_{ξ}_{1}, q_{ijk}) = f_{1}(∂_{ξ}_{1}, q)|_{(∂}_{ξ2}_{~}_{x, ∂}_{ξ3}_{~}_{x)}_{ijk}

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

**Riemann problem**

Riemann problem at time τ = 0

∂_{τ}q_{L} + f_{1} (∂_{ξ}_{1}, q_{L}) = 0 ∂_{τ}q_{R} + f_{1} (∂_{ξ}_{1}, q_{R}) = 0
ξ_{1}
τ

Q^{n}_{i−1,jk} Q^{n}_{ijk}

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

**Riemann problem**

Exact Riemann solution: basic structure

∂_{τ}q_{L} + f_{1} (∂_{ξ}_{1}, q_{L}) = 0 ∂_{τ}q_{R} + f_{1} (∂_{ξ}_{1}, q_{R}) = 0
ξ_{1}
τ

Q^{n}_{i−1,jk} Q^{n}_{ijk}

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

**Riemann problem**

Shock-only approximate Riemann solution: basic structure

∂_{τ}q_{L} + f_{1} (∂_{ξ}_{1}, q_{L}) = 0 ∂_{τ}q_{R} + f_{1} (∂_{ξ}_{1}, q_{R}) = 0
ξ_{1}
τ

Q^{n}_{i−1,jk} Q^{n}_{ijk}

λ_{0}

λ_{1} λ_{2}

λ_{3}
q_{mL}^{−} q_{mL}^{+}

q_{mR}

Z^{1} = f_{L}(q_{mL}^{−} ) − f_{L}(Q^{n}_{i−1,jk}) Z^{2} = f_{R}(q_{mR}) − f_{R}(q_{mL}^{+} )

Z^{3} = f_{R}(Q^{n}_{ijk}) − f_{R}(q_{mR})

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

**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_{0})υ_{m} − M_{L}(p_{m})
ρ_{mL}(p_{m})

∇_{X}_{~} ξ_{1}

i−1/2,jk

(λ_{1,2})_{i−1/2,jk} = (1 − h_{0})υ_{m}

∇_{X}_{~}ξ_{1}

i−1/2,jk

(λ_{1,3})_{i−1/2,jk} =

(1 − h_{0})υ_{m} + M_{R}(p_{m})
ρ_{mR}(p_{m})

∇_{X}_{~} ξ_{1}

i−1/2,jk

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

**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

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

**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

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

**Woodward-Colella’s 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 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

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

**Woodward-Colella’s problem**

h0 = 0 Eulerian result

h0 = 0.99 Lagrangian-like result

sharper resolution for contact discontinuity

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

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

**Woodward-Colella’s problem**

h0 = 0 Eulerian result

h0 = 0.99 Lagrangian-like result

sharper resolution for contact discontinuity

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

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

**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

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

### 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

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

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

### 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

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

### 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

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