**Wave-propagation based**

**generalized Lagrangian method** **for**

**hyperbolic balance laws**

**Application to inviscid compressible flow**

**Application to inviscid compressible flow**

*Keh-Ming Shyue*

Department of Mathematics National Taiwan University

Taiwan

**Objective**

Describe simple Lagrangian-like moving grid approach for numerical resolution of nonlinear hyperbolic balance laws of the form

∂

∂tq (~x, t) +

N

X

j=1

∂

∂x_{j} f_{j} (q, ~x) = ψ (q, ~x)

with discontinuous initial data in general N ≥ 1 rectangular or non-rectangular geometry

~x = (x_{1}, x_{2}, . . . , x_{N}): spatial vector, t: time
q ∈ lR^{m}: vector of m state quantities

f_{j} ∈ lR^{m}: flux vector, j = 1, 2, . . . , N, ψ ∈ lR^{m}: source term
Model equation is hyperbolic if ^{P}^{N}_{j=1}^{d} α_{j} (∂f_{j}/∂q) is

diagonalizable with real eigenvalues, α_{j} ∈ lR

**Outline**

Mathematical model for general balance laws Eulerian formulation

Generalized Lagrangian formulation

Example to single component compressible flow Wave-propagation based finite volume methods

Generalized Riemann problem & approximate solver Flux-based wave decomposition method

Sample numerical examples

Extension to compressible two-phase flow Future work

**Mathematical Model**

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

**Mathematical Model (Cont.)**

To derive hyperbolic balance laws 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}

!

= ψ(q)

Note this is not in divergence form, and hence is not conservative, in case the source term ψ is ignored.

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

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

#

Note with the use of basic grid-metric relations (see next), it is known that

∂J

∂τ +

N

X

j=1

∂

∂ξ_{j}

J ∂ξ_{j}

∂t

= 0 (geometrical conservation law)

N

X

k=1

∂

∂ξ_{k}

J ∂ξ_{k}

∂x_{j}

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

and hence G = 0

**Basic Grid-Metric Relations**

Assume existence of inverse transformation

t = τ, x_{j} = x_{j}(~ξ, t) for j = 1, 2, . . . , N ,

To find basic grid-metric relations between different coordinates, employ elementary differential rule

∂(τ, ~ξ)

∂(t, ~x) = ∂(t, ~x)

∂(τ, ~ξ)

−1

yielding in N = 3 case, for example, as

1 0 0 0

∂_{t}ξ_{1} ∂_{x}_{1}ξ_{1} ∂_{x}_{2}ξ_{1} ∂_{x}_{3}ξ_{1}

∂_{t}ξ_{2} ∂_{x}_{1}ξ_{2} ∂_{x}_{2}ξ_{2} ∂_{x}_{3}ξ_{2}

∂_{t}ξ_{3} ∂_{x} ξ_{3} ∂_{x} ξ_{3} ∂_{x} ξ_{3}

= 1 J

J 0 0 0

J_{01} J_{11} J_{21} J_{31}
J_{02} J_{12} J_{22} J_{32}
J_{03} J_{13} J_{23} J_{33}

**Grid-Metric Relations (Cont.)**

Here

J =

∂(x_{1}, x_{2}, x_{3})

∂(ξ_{1}, ξ_{2}, ξ_{3})

= det ∂(x_{1}, x_{2}, x_{3})

∂(ξ_{1}, ξ_{2}, ξ_{3})

,
J_{11} =

∂(x_{2}, x_{3})

∂(ξ_{2}, ξ_{3})

, J_{21} =

∂(x_{1}, x_{3})

∂(ξ_{3}, ξ_{2})

, J_{31} =

∂(x_{1}, x_{2})

∂(ξ_{2}, ξ_{3})

,
J_{12} =

∂(x_{2}, x_{3})

∂(ξ_{3}, ξ_{1})

, J_{22} =

∂(x_{1}, x_{3})

∂(ξ_{1}, ξ_{3})

, J_{32} =

∂(x_{1}, x_{2})

∂(ξ_{3}, ξ_{1})

,
J_{13} =

∂(x_{2}, x_{3})

∂(ξ_{1}, ξ_{2})

, J_{23} =

∂(x_{1}, x_{3})

∂(ξ_{2}, ξ_{1})

, J_{33} =

∂(x_{1}, x_{2})

∂(ξ_{1}, ξ_{2})

,

J_{0j} = −

Nd

X

i=1

J_{ij}∂_{τ}x_{i}, j = 1, 2, 3,

and so grid-metric relations between different coordinates

∇ξ = (∂ ξ , ∇ ξ ) = (∂ ξ , ∂ ξ , ∂ ξ , ∂ ξ ) = 1

(J , J , J , J )

**Grid-Metric Relations (Cont.)**

Note in two dimensions N = 2, we have

∂ξ_{1}

∂t , ∂ξ_{1}

∂x_{1} , ∂ξ_{1}

∂x_{2}

= 1 J

−∂x_{1}

∂τ

∂x_{2}

∂ξ_{2} + ∂x_{2}

∂τ

∂x_{1}

∂ξ_{2} , ∂x_{2}

∂ξ_{2} , −∂x_{1}

∂ξ_{2}

∂ξ_{2}

∂t , ∂ξ_{2}

∂x_{1} , ∂ξ_{2}

∂x_{2}

= 1 J

∂x_{1}

∂τ

∂x_{2}

∂ξ_{1} − ∂x_{2}

∂τ

∂x_{1}

∂ξ_{1} , −∂x_{2}

∂ξ_{1} , ∂x_{1}

∂ξ_{1}

J = ∂x_{1}

∂ξ_{1}

∂x_{2}

∂ξ_{2} − ∂x_{1}

∂ξ_{2}

∂x_{2}

∂ξ_{1}

Thus to have G = 0 fulfilled, grid-metrics should obey

∂J

∂τ + ∂

∂ξ_{1}

J ∂ξ_{1}

∂t

+ ∂

∂ξ_{2}

J ∂ξ_{2}

∂t

= 0

∂

∂ξ_{1}

J ∂ξ_{1}

∂x_{1}

+ ∂

∂ξ_{2}

J ∂ξ_{2}

∂x_{1}

= ∂

∂ξ_{1}

∂x_{2}

∂ξ_{2}

+ ∂

∂ξ_{2}

−∂x_{2}

∂ξ_{1}

= 0

∂

∂ξ

J ∂ξ_{1}

∂x

+ ∂

∂ξ

J ∂ξ_{2}

∂x

= ∂

∂ξ

−∂x_{1}

∂ξ

+ ∂

∂ξ

∂x_{1}

∂ξ

= 0

**Mathematical Model (Cont.)**

As an example, with gravity effect included, Euler equations for single component compressible gas flow take

Cartesian coordinate case

∂

∂t

ρ
ρu_{i}

E

+

N

X

j=1

∂

∂x_{j}

ρu_{j}

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

=

0

−ρ_{∂x}^{∂φ}

i

−ρ~u · ∇φ

, i = 1, . . . , N

Generalized coordinate case

∂

∂τ

ρ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 · ∇φ

ρ: density, p = p(ρ, e): pressure , e: internal energy

E = ρe + ρPN

j=1 u^{2}_{j}/2: total energy, φ: gravitational potential

**Mathematical Model (Cont.)**

Note that to complete the model, we must

1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially

Depending on how complex the geometry is, this can be done by various numerical means

**Mathematical Model (Cont.)**

Note that to complete the model, we must

1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially

Depending on how complex the geometry is, this can be done by various numerical means

2. choose a moving grid strategy for ∂_{τ}~x

**Mathematical Model (Cont.)**

Note that to complete the model, we must

1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially

Depending on how complex the geometry is, this can be done by various numerical means

2. choose a moving grid strategy for ∂_{τ}~x

When ∂_{τ}~x = 0 or ∂_{τ}~x = ~u_{b}(t) (rigid-body motion)

∂_{t}ξ~ & ∇_{~}_{x}ξ~ are time-independent; no need to have
more additional condition

**Mathematical Model (Cont.)**

Note that to complete the model, we must

1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially

Depending on how complex the geometry is, this can be done by various numerical means

2. choose a moving grid strategy for ∂_{τ}~x

When ∂_{τ}~x = 0 or ∂_{τ}~x = ~u_{b}(t) (rigid-body motion)

∂_{t}ξ~ & ∇_{~}_{x}ξ~ are time-independent; no need to have
more additional condition

While ∂_{τ}~x 6= 0 (flow-dependent motion) (see next)

∂_{t}ξ~ & ∇_{~}_{x}ξ~ would be time-dependent; require

additional conditions to determine ∇_{ξ}_{~}~x (N^{2} of them
in total) over time (see below)

**Lagrangian-Like Moving Grid**

For compressible flow, 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 = h^{0}~u (pseudo velocity)
h^{0} ∈ [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)

**Lagrangian-Like Moving Grid**

For compressible flow, 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 = h^{0}~u (pseudo velocity)
h^{0} ∈ [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)

Here we will focus on the simple Lagrangian-like case

**Unified Coordinate System**

Consider N = 2 case, for example, and use simplified

notation ~x = (x, y), ξ = (ξ, η)^{~} . At given time instance, free
parameter h can be chosen based on

Grid-angle preserving condition (Hui * ^{et al.}* JCP 1999)

∂

∂τ cos^{−1} ∇ξ

|∇ξ| · ∇η

|∇η|

= ∂

∂τ cos^{−1}

−y_{η}x_{η} − y_{ξ}x_{ξ}
qy_{ξ}^{2} + y_{η}^{2}q

x^{2}_{ξ} + x^{2}_{η}

= · · ·

= Ah_{ξ} + Bh_{η} + Ch = 0 (1st order PDE )

with

A = q

x^{2}_{η} + y_{η}^{2} (vx_{ξ} − uy_{ξ}) , B = q

x^{2}_{ξ} + y_{ξ}^{2} (uy_{η} − vx_{η})
C = q

x^{2}_{ξ} + y_{ξ}^{2} (u_{η}y_{η} − v_{η}x_{η}) − q

x^{2}_{η} + y_{η}^{2} (u_{ξ}y_{ξ} − v_{ξ}x_{ξ})

**Unified Coordinate System**

Consider N = 2 case, for example, and use simplified

notation ~x = (x, y), ξ = (ξ, η)^{~} . Or alternatively, based on
Grid-Jacobian preserving condition

∂J

∂τ = ∂

∂τ (x_{ξ}y_{η} − x_{η}y_{ξ})

= x_{ξτ} y_{η} + x_{ξ} y_{ητ} − x_{ητ} y_{ξ} − x_{η} y_{ξτ}

= · · ·

= Ah_{ξ} + Bh_{η} + Ch = 0 (1st order PDE )

with

A = uy_{η} − vx_{η}, B = vx_{ξ} − uy_{ξ}, C = u_{ξ}y_{η} + v_{η}x_{ξ} − u_{η}y_{ξ} − v_{ξ}x_{η}

**Lagrangian-Like Grid (Cont.)**

Now with the temporal motion of the coordinate system
governed by ∂_{τ}~x = h^{0}~u. 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}

**Mathematical Model (Cont.)**

In summary, our Lagrangina-like model system 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.

Moving grid condition ∂_{τ}~x = h^{0}~u & pressure law p(ρ, e)

**Axisymmetric Compressible Flow**

Physical balance laws (ξ_{1}: axisymmetric direction)

∂

∂τ

ρJ
ρJu_{i}

JE

+

2

X

j=1

∂

∂ξ_{j} J

ρU_{j}

ρu_{i}U_{j} + p^{∂ξ}_{∂x}^{j}

i

EU_{j} + pU_{j} − p^{∂ξ}_{∂t}^{j}

=

−_{x}^{1}ρJu_{1}

−_{x}^{1}ρJu_{i}u_{j} − ρJ _{∂x}^{∂φ}

i

−_{x}^{1}J(E + p)u_{1} − ρJ~u · ∇φ

for i = 1, 2

Geometrical conservation laws

∂

∂τ

∂x_{i}

∂ξ_{j}

+ ∂

∂ξ_{j}

−∂x_{i}

∂τ

= 0 for i, j = 1, 2.

Moving grid condition ∂_{τ}~x = h_{0}~u & pressure law p(ρ, e)

**Mathematical Models: Remarks**

For single component compressible flow model mentioned above, it is known that under some thermodynamic stability conditions

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

**Mathematical Models: Remarks**

For single component compressible flow model mentioned above, it is known that under some thermodynamic stability conditions

when h^{0} = 0 (Eulerian case), the model is hyperbolic
when h^{0} = 1 (Lagrangian case), the model is weakly
hyperbolic (do not possess complete eigenvectors)
when h^{0} ∈ (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.

**Review of Previous Work**

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

**Numerical Methods**

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}

**Numerical Methods (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 generalized Riemann problem (defined below) 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

**Numerical Methods (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}

**Numerical Methods (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

**Numerical Methods (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,

**Generalized Riemann Problem**

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

**Generalized Riemann Problem**

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

**Generalized Riemann Problem**

Exact generalized Riemann solution: basic structure

∂ q + f (∂ , q ) = 0 ∂ q + f (∂ , q ) = 0
ξ_{1}
τ

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

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

**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 − h )υ + M_{R}(p_{m})

∇ ξ

**Lax’s Riemann Problem**

Ideal gas EOS with γ = 1.4
h^{0} = 0 Eulerian result

h^{0} = 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

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

**Woodward-Colella’s Problem**

Ideal gas EOS with γ = 1.4
h^{0} = 0 Eulerian result

h^{0} = 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

**Woodward-Colella’s Problem**

h^{0} = 0 Eulerian result

h^{0} = 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

**Woodward-Colella’s Problem**

h^{0} = 0 Eulerian result

h^{0} = 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

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

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

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

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

**Radially Symmetric Problem**

0 0.2 0.4

0 0.1 0.2 0.3 0.4 0.5

0 0.2 0.4

0 0.1 0.2 0.3 0.4 0.5

0 0.2 0.4

0 0.1 0.2 0.3 0.4

Density Pressure 0.5 Physical grid
a) h^{0} = 0.99

0.1 0.2 0.3 0.4 0.5

0.1 0.2 0.3 0.4 0.5

0.1 0.2 0.3 0.4

Density Pressure 0.5 Physical grid
b) h^{0} = 0

**Radially Symmetric Prob. (Cont.)**

0 0.1 0.2 0.3 0.4 0.5

0.2 0.4 0.6 0.8 1 1.2

one−d h0=0 h0=0.99

0 0.1 0.2 0.3 0.4 0.5

0 0.1 0.2 0.3 0.4

0 0.1 0.2 0.3 0.4 0.5

0 0.5 1 1.5

0 0.1 0.2 0.3 0.4 0.5

0.5 1 1.5 2

r(m) r(m)

ρ(Mg/m3 ) ¯u(km/s)

p(GPa) J