**Volume of fluid methods** **for**

**compressible multiphase flow**

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

**Objective**

Discuss a class of volume of fluid (vs. level set, MAC,

particles) methods for interface problems with application to compressible multiphase flow

1** ^{.}** Adaptive moving grid approach (this lecture)
Cartesian grid embedded volume tracking

Useful also for turbulent flow (Glimm & Abarzhi’s talk)

Moving mapped grid interface capturing

Simple method if mesh handling is robust

2** ^{.}** Eulerian interface sharpening approach (next lecture)
Artificial interface compression method

Anti-diffusion method

**Incompressible 2-phase flow**

One model for free-surface (or 2-phase) flow can be described by incompressible Navier-Stokes equations:

∇ · ~u = 0 (Continuity)

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = ∇ · τ + ρ~g + ~f_{σ} (Momentum)

∂_{t}α + ~u · ∇α = 0 (Volume fraction transport )
Material quantities in 2-phase coexistent region are often

computed by α-weighted average (α ∈ [0, 1]) as

ρ = αρ_{1} + (1 − α) ρ_{2}, ǫ = αǫ_{1} + (1 − α) ǫ_{2},
where

τ = ǫ ∇~u + ∇~u^{T}

, f~_{σ} = −σκ∇α with κ = ∇ ·

∇α

|∇α|

**Incompressible 2-phase flow**

Accurate advection scheme for volume fraction α is of fundamental importance. Typical methods employed are

Simple line interface calculation (SLIC) Donor acceptor scheme

Flux corrected transport method Youngs (PLIC) method

Level set method

CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes)

THINC (Tangent of Hyperbola for INterface Capturing) or its variant

**Compressible flow in porous medium**

One model for gas flow in variable porous medium can be described by compressible Euler equations:

∂_{t} (αρ) + ∇ · (αρ~u) = 0

∂_{t} (αρ~u) + ∇ · (αρ~u ⊗ ~u) + ∇ (αp) = p∇α

∂_{t} (αρE) + ∇ · (αρE~u + αp~u) = 0

Embeded discontinuous porosity α ∈ [0, 1]

Well-balanced & entropy solution due to non-conservative term

**Baer-Nunziato 2-phase model**

Known 2-phase mixture model for deflagration-to-detonation transition (DDT) in reactive granular materials

(α_{1}ρ_{1})_{t} + ∇ · (α_{1}ρ_{1}~u_{1}) = 0

(α_{1}ρ_{1}~u_{1})_{t} + ∇ · (α_{1}ρ_{1}~u_{1} ⊗ ~u_{1}) + ∇(α_{1}p_{1}) = p_{0}∇α_{1} + λ(~u_{2} − ~u_{1})
(α_{1}ρ_{1}E_{1})_{t} + ∇ · (α_{1}ρ_{1}E_{1}~u_{1} + α_{1}p_{1}~u_{1}) = p_{0} (α_{2})_{t} + λ~u_{0} · (~u_{2} − ~u_{1})
(α_{2}ρ_{2})_{t} + ∇ · (α_{2}ρ_{2}~u_{2}) = 0

(α_{2}ρ_{2}~u_{2})_{t} + ∇ · (α_{2}ρ_{2}~u_{2} ⊗ ~u_{2}) + ∇(α_{2}p_{2}) = p_{0}∇α_{2} − λ(~u_{2} − ~u_{1})

(α_{2}ρ_{2}E_{2}) _{t} + ∇ · (α_{2}ρ_{2}E_{2}~u_{2} + α_{2}p_{2}~u_{2}) = −p_{0} (α_{2})_{t} − λ~u_{0} · (~u_{2} − ~u_{1})
(α_{2})_{t} + ~u_{0} · ∇α_{2} = µ (p_{2} − p_{1})

α_{k} = V_{k}/V : volume fraction (α_{1} + α_{2} = 1), ρ_{k}: density,

~u_{k}: velocity, p_{k}: pressure, E_{k} = e_{k} + ~u^{2}_{k}/2: specific total
energy, e_{k}: specific internal energy, k = 1, 2

**Baer-Nunziato 2-phase model**

Quantities p_{0} & ~u_{0} denote interfacial pressure & velocity
Baer & Nunziato (1986)

p_{0} = p_{2}, ~u_{0} = ~u_{1}
Saurel & Abgrall (1999)

p_{0} = P_{2}

k=1 α_{k}p_{k}, ~u_{0} = P_{2}

k=1 α_{k}ρ_{k}~u_{k}

P_{2}

k=1 α_{k}ρ_{k}

Quantities λ & µ (> 0) are velocity & pressure relaxation parameters that give rates to equilibrium for velocitiy &

pressure, respectively

**Reduced 2-phase model**

For some problems, it is simple to use reduced 2-phase model of Kapila, Menikoff, Bdzil, Son, & Stewart (2001)

∂_{t} (α_{1}ρ_{1}) + ∇ · (α_{1}ρ_{1}~u) = 0

∂_{t} (α_{2}ρ_{2}) + ∇ · (α_{2}ρ_{2}~u) = 0

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

∂_{t} (ρE) + ∇ · (ρE~u + p~u) = 0

∂_{t}α_{1} + ~u · ∇α_{1} = ρ_{2}c^{2}_{2} − ρ_{1}c^{2}_{1}
P2

k=1 ρ_{k}c^{2}_{k}/α_{k} ∇ · ~u

Model is zero-order approximation of Baer-Nunziato model with stiff mechanical relaxation

Positivity & sharpness of α is important for accurate solution

**Cartesian embedded boundary method**

Towards volume tracking for sharp interface resolution, we consider conservation laws

∂_{t}q + ∇ · ~f (q) = 0

with embedded boundary (* ^{i.e.}*, interface with rigid BCs)

q ∈ lR^{m}: vector of m
conserved quantities

f = (f~ _{1}, f_{2}, . . . , f_{d}) ∈ lR^{m×d}:
flux matrix; f_{j} ∈ lR^{m}

Γ: embedded boundary

Γ →

**Finite-volume discretization**

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

Q^{n}_{S} ≈ 1
M(S)

Z

S

q(x, t_{n}) dV = 1

κ_{S} (∆x)^{d}
Z

S

q(x, t_{n}) dV
M(S): measure (area in 2D

or volume in 3D) of cell S
κ_{S} = M(S)/ (∆x)^{d}: volume
fraction

h^{k}_{S} = M(∂S_{k})/∆x^{(d−1)}: area
fraction

Assume uniform discretiza- tion ∆x = ∆x , j = 1, . . . , d

A B

C D

**Finite-volume discretization**

Divergence theorem over each control volume leads to finite volume approximation for

∇ · ~f (q)

S ≈ 1

M(S) Z

S

∇ · ~f (q) dV = 1 M(S)

Z

∂S

f (q) · ~n dA~

= 1

κ_{S}∆x

X

k∈∂S

h^{k}_{S}f (q)~ ^{k} · ~n^{k} = 1
κ_{S}∆x

X

k∈∂S

h^{k}_{S}F˘^{k}

≡ D_{S} · ~f (q)

Away from boundaries,

method reduces to standard conservative finite difference discretization

τ^{∆x} =

∇ · ~f

S − D_{S} · ~f =
O(∆x^{p−1}/κ ) if F^{˘}^{k} = O(∆^{p})

A B

C D

**Finite-volume discretization**

For conservation laws, On each control volume S, a fully discretized finite volume method can be written as

Q^{n+1}_{S} = Q^{n}_{S} − ∆tD_{S} · ~f (q) = Q^{n}_{S} − ∆t
κ_{S}∆x

X

k∈∂S

h^{k}_{S}F˘^{k}

which satisfies discrete conservation identity X

∪S

κ_{S}Q^{n+1} = X

∪S

κ_{S}Q^{n} − ∆t

∆x

X

∪S

X

k∈∂S

h^{k}_{S}F˘^{k}

Method is subject to stringent CFL stability constraint

∆t = O

(κ_{min})^{1/d}∆x

|λ|_{max}

≪ 1 if κ_{min} ≪ 1 (Small-cell problem)

**Small-cell problem**

Methods proposed for small-cell problem include Non-conservation approach: Glimm & coworker

Flux-redistribution approach: Chern, Colella, & Modiano Cell-merging approach: Quirk, Powell, Causon

H-box approach: Berger, Helzel, LeVeque (see below)

Q¯

R( ¯Q) Q¯_{R}

Q¯_{L}

**Flux calculation: Wave propagation**

On 2D uniform grid, high-resolution in fluctutation form
Q^{n+1}_{ij} = Q^{n}_{ij} − ∆t

∆x_{1}

A^{+}_{1} ∆Q_{i−}^{1}

2,j + A^{−}_{1} ∆Q_{i+}^{1}

2,j

− ∆t

∆x_{2}

A^{+}_{2} ∆Q_{i,j−}^{1}

2 + A^{−}_{2} ∆Q_{i,j+}^{1}

2

− ∆t

∆x_{1}

F˜_{i+}^{1} 1

2,j − ˜F_{i−}^{1} 1
2,j

− ∆t

∆x_{2}

F˜_{i,j+}^{2} 1
2

− ˜F_{i,j−}^{2} 1
2

A^{+}_{1} ∆Q

A^{−}_{1} ∆Q

A^{−}_{2} A^{+}_{1} ∆Q

A^{+}_{2} A^{−}_{1} ∆Q

A^{+}_{2} A^{+}_{1} ∆Q

A^{−}_{2} A^{−}_{1} ∆Q

**Flux calculation: Wave propagation**

Fluctuation terms in wave propagation form are obtained
from solving Riemann problems at each cell edge, yielding
speeds λ^{j,m} & waves W^{j,k}, j = 1, 2, k = 1, . . . , m_{w},

A^{±}_{1} ∆Q =

mw

X

k=1

λ^{1,k}^{±}

W^{1,k}, A^{±}_{2} ∆Q =

mw

X

k=1

λ^{2,k}^{±}

W^{2,k}

A^{±}_{2} A^{±}_{1} ∆Q =

mw

X

p=1

λ^{1,p}^{±}

A^{±}_{2} W^{1,p} =

mw

X

p=1

λ^{1,p}^{±} X^{m}^{w}

q=1

λ^{2,q}^{±}

W^{2,q}

F˜^{1} := F˜^{1} + 1
2

mw

X

k=1

_{λ}^{1,k}

1 − ∆t

∆x_{1}

_{λ}^{1,k}

Wf^{1,k}

To aviod oscillations near discontinuities, a wave limiter is
applied leading to limited waves _{W}f

**Wave propagation with boundary**

First order upwind with normal-wave update

Solve Riemann problem in direction normal at cell edge,

& use resulting waves to update nearby cell averages
Q^{n+1}_{S} := Q^{n+1}_{S} − M (Z^{p} ∩ S)

M(S) W^{p}

Z↓^{p}

A B

C

D Z↓^{p}

A B

C D

**Wave propagation with boundary**

First order upwind with transverse-wave update

Use transverse portion of equation, solve Riemann

problem in transverse direction, & use resulting waves to update neighboring cell averages as usual

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

Z↓^{pq}

A B

C

D Z↓^{pq}

A B

C D

**Wave propagation with boundary**

High resolution with 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

**Moving boundary tracking**

With a prescribed moving-object velocity ~u_{b} at current
time, tracked boundaries are advanced in a Lagrangian
manner over time step ∆t

Final moving boundary locations are set by where new segments intersect grid lines; interpolation is required with care to initialize cut-cell values

DC BA new object↓

↑

old object

**Embedded boundary condition**

Rigid embedded boundary (stationary or moving), reflection
principle is used to assign states for fictitious subcells, * ^{i.e.}*,

z := z_{E} (z = ρ, p, α)

~u := ~u_{E} − 2(~u_{E} · ~n)~n + 2(~u_{b} · ~n),
where ~u_{b} is velocity of object

If motion of object is governed by Newton’s second law,
velocity of object ~u^{n+1}_{b} at time step t_{n+1} is computed by

~u^{n+1}_{b} = ~u^{n}_{b} + ∆t

M_{b}F~_{b}^{n},

where M_{b} is mass of object & F^{~}_{b}^{n} is force on object

**Moving piston problem**

Numerical results agree well with exact solution
Piston moves with speed ~u_{b} = (0, 300)m/s

0 0.5 1

0.5 1 1.5 2

x (m) ρ (kg/m3 )

0 0.5 1

0 0.05 0.1 0.15 0.2 0.25 0.3

x (m)

u (km/s)

0 0.5 1

0.5 1 1.5 2 2.5

x (m)

p (Bar)

piston

**Cylinder lift off**

Cylinder motion 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.05 0.1 0.15 0.2

t = 0

t = 0.1641s

t = 0.30085s

**Cylinder lift off**

A convergence study of center of cylinder & relative mass loss at final stopping time t = 0.30085s

Mesh size Center of cylinder Relative mass loss 250 × 50 (0.618181, 0.134456) −0.257528

500 × 100 (0.620266, 0.136807) −0.131474 1000 × 200 (0.623075, 0.138929) −0.066984

Results are comparable with numerical appeared in literature

**Falling rigid object in water tank**

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2

Density Pressure 3 Volume fraction air

water t = 1ms

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

air

water t = 3ms

**VOF interface tracking**

Given volume fractions, VOF tracking consists of:

1** ^{.}** Volume fraction update

Take a time step on current grid using wave propagation method to update volume fractions α governed by

∂_{t}α + ~u · ∇α = 0

at next time step; ~u underlying velocity field

2** ^{.}** Interface reconstruction (by PLIC or MOF, see Z. Jia’s
talk at 17:40-18:00, 9/15)

Standard PLIC method does:

Compute interface normal, using method such as Gradient or least squares of Youngs or Puckett Determine interface location by iterative bisection

**Interface reconstruction: Example**

Cell-averaged volume fractions (left) & reconstructed interface (right); a set of disconnected line segments

0 0

0 0

0 0

0 0

0 0

0

0 0

0 0

0 0

0 0

0

1 0.29

0.68 0.09

0.51 0.51

0.09 0.68

0.29 interface↓

**Interface reconstruction: Example**

With α in previous slide & ~u = (1, 1) updated volume
fraction over ∆t = 0.06 shown here ( κ_{1} = 5.7 × 10^{−}^{3} &

κ_{2} = 1.3 × 10^{−}^{3})

New reconstructed interface location

0

0 0 0 0

0 0 0 0

0 0

0 0

0 0

0 0

0

1 0.01

0.38 0.11

0.25 0.72

0.06 0.85

0.74 κ^{1}
κ_{2}

↑

↓

old interface

new interface

**Solution update**

Solution is double-valued on cells subdivided by VOF
tracked interface; Again use finite volume formulation of
solution Q^{n}_{S} for approximate value of solution q

Q^{n}_{S} ≈ 1
M(S)

Z

S

q(x, t_{n}) dV

C E

D F

G H

**Wave propagation method**

Normal- & transverse-edge wave propagation

↓

↓

Z_{p}
Z_{p}

↓ ↓

Z_{pq}
Z_{pq}

**Stability issues**

Choose time step ∆t based on uniform mesh size ∆x as
ν = |λ|_{max}∆t

∆x ≤ 1

Use wave propagation method described above which expands range of influence of solution to maintain

stability of method even in presence of small Cartesian cut cells

Use smoothing operator (such as, h-box approach of
Berger * ^{et al.}* ) locally for cell averages in irregular cells
for accuracy

**Shock in air & R22 bubble**

**Shock in air & R22 bubble**

**Shock in air & R22 bubble**

**Shock in air & R22 bubble**

Quantitative assessment of prominent flow velocities:

Velocity (m/s) V_{s} V_{R} V_{T} V_{ui} V_{uf} V_{di} V_{df}
Haas & Sturtevant 415 240 540 73 90 78 78
Quirk & Karni 420 254 560 74 90 116 82
Our result (tracking) 411 243 538 64 87 82 60
Our result (capturing) 411 244 534 65 86 98 76

V_{s} (V_{R}, V_{T}) Incident (refracted, transmitted) shock
speed t ∈ [0, 250]µs (t ∈ [0, 202]µs, t ∈ [202, 250]µs )
V_{ui} (V_{uf}) Initial (final) upstream bubble wall speed
t ∈ [0, 400]µs (t ∈ [400, 1000]µs)

V_{di} (V_{df}) Initial (final) downstream bubble wall speed
t ∈ [200, 400]µs (t ∈ [400, 1000]µs)

**Underwater explosion**

**Flying Aluminum plate**

Track vacuum-Al boundary Impact velocity 400m/s

Pressure contours are shown

Al target Al flyer

vacuum

time = 0.5µs time = 1µs time = 2µs time = 4µs

**Moving mapped grid method**

Towards sharp interface resolution, one may also consider moving mapped grid in that hyperbolic conservation laws

∂q

∂t +

Xd j=1

∂f_{j}(q)

∂x_{j} = 0

in transformed coordinates (x, t) 7→ (ξ, τ ) becomes conservation laws with patially varying fluxes:

∂

∂τ (Jq) +

Xd j=1

∂

∂ξ_{j} Jq ∂ξ_{j}

∂t +

Xd k=1

Jf_{k}(q) ∂ξ_{j}

∂x_{k}

!

= 0

To close system, conditions on grid motion ∂_{t}ξ & so grid
metrics ∂ξ_{j}/∂x_{k} with J = det (∂ξ/∂x)^{−}^{1} should be set

**Computing geometric quantities**

To set ∂_{t}ξ, ∂ξ_{j}/∂x_{k}, & J, one may use
Hui’s unified coordinate approach

∂x_{i}

∂τ = hu_{i} & ∂

∂τ

∂x_{i}

∂ξ_{j}

+ ∂

∂ξ_{j}

−∂x_{i}

∂τ

= 0

for i, j = 1, 2, . . . , d & compute J = det (∂ξ/∂x)^{−}^{1} directly.

Quantity h may be kept as a constant or varied with time by geometric measure

Variant ALE approach, see next talk by B. Chen or B.

Tian’s talk at 11:50-12:10, 9/15

∂J

∂τ +

Xd _{∂}

∂ξ_{j}

J ∂ξ_{j}

∂t

= 0

**Reduced 2-phase model (revisit)**

Physical balance laws

∂

∂τ

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

ρEJ

^{+}

Xd j=1

∂

∂ξ_{j} J

α_{1}ρ_{1}U_{j}
α_{2}ρ_{2}U_{j}
ρu_{i}U_{j} + p_{∂x}^{∂ξ}^{j}

i

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

^{= 0}

∂α_{1}

∂τ +

Xd j=1

U_{j} ∂α_{1}

∂ξ_{j} = ρ_{2}c^{2}_{2} − ρ_{1}c^{2}_{1}
P_{2}

k=1 ρ_{k}c^{2}_{k}/α_{k}

Xd j=1

∂U_{j}

∂ξ_{j}
U_{j} = ∂_{t}ξ_{j} + P_{d}

i=1 u_{i}∂_{x}_{i}ξ_{j}: ξ_{j}’s contravariant velocity

Mixture equation of state: p = p(α_{2}, α_{1}ρ_{1}, α_{2}ρ_{2}, ρe) with
isobaric closure: p_{1} = p_{2} = p

**Numerical approximation**

To ensure free-stream (constant states) preservation, which is one of basic properties that a moving mesh method

should have, a consistent discretization of

∂J

∂τ +

Xd j=1

∂

∂ξ_{j}

J ∂ξ_{j}

∂t

= 0

is required which may be difficult for 2D quadrilateral or 3D hexahedral mesh

In this regards, non-conservative form is often employed

∂q

∂τ + 1 J

Xd j=1

∂

∂ξ_{j}

Xd

k=1

Jf_{k}(q) ∂ξ_{j}

∂x_{k} + Jq∂ξ_{j}

∂t

!

= 0

**PLIC-VOF adaptive moving mesh**

We remark that Menciger & ˘Zun (JCP 2011) developed a PLIC-VOF adaptive moving grid method in that monitor function D employed in elliptic grid generator, say in form proposed by Ceniceros & Hou (JCP 2001),

∇_{ξ} · (D∇_{ξ}x_{j}) = 0, j = 1, . . . , d

for x(ξ), make use of volume fraction α explicitly, * ^{i.e.}*,

∂D

∂τ = ε∇^{2}D with D(x, 0) = 1 + βM (∇α) & D|_{∂Ω} = 0
where M (z) = |z|R

dV / R

|z|dV

This can be used readily in moving-mesh method developed by T. Tang & many others

**Woodward-Colella’s problem**

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

0 0.5 1

0 5 10 15 20

ρ

Fine grid
h_{0}=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

0 2 4 6 8

ρ

−5 0 5 10 15

0 200 400 600

u p

t = 0.038 t = 0.038

t = 0.038

**Oscillating water column**

Initially, in closed shock tube, water column moves at u = 1 from left to right, yielding air compression at right

& air expansion at left

Subsequently, pressure difference built up across water column resulting deceleration of column of water to

right, makes a stop, & then acceleration to left; a reverse pressure difference built up across water column redirecting flow from left to right again

Eventually, water column starts to oscillate

air

air water

u →

**Oscillating water column**

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

time

oscilw grid

0 1 2 3 4 5 6 7 8 9 10

0.8 1 1.2 1.4 1.6 1.8

pL

0.8 1 1.2 1.4 1.6 1.8

pR

mmesh mesh0

mmesh mesh0

**Water-vapor cavitation**

Initially, in closed shock tube, flow is homogeneous
(contains α = 10^{−}^{6} gas in bulk liquid) at standard

atmospheric condition & exists interface separating flow with opposite motion (u = 100 m/s)

Result in pressure drop & formation of cavitation zone in middle; shocks form also from both ends

Eventually, shock-cavitation collision occurs

−u u

← Inteface

**Water-vapor cavitation**

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1x 10^{−3}

time

cavi grid

**Underwater explosions**

t=0.2ms

t=0.4ms

t=0.8ms

t=1.2ms

Moving mesh Eulerian mesh air

water

**Underwater explosions**

Moving mesh (coarsen by factor 5)

−1 0 1

−1

−0.5 0 0.5

−1 0 1

−1

−0.5 0 0.5

time = 0.2ms time = 0.4ms

−1

−0.5 0 0.5

−1

−0.5 0 0.5

time = 0.8ms time = 1.2ms

**Accuracy test in 2D**

Consider 2D compressible Euler equations with ideal gas law as governing equations

Take smooth vortex flow with initial condition

ρ =

1 − 25(γ − 1)

8γπ^{2} exp (1 − r^{2})

1/(γ−1)

,
p = ρ^{γ},

u_{1} = 1 − 5

2π exp ((1 − r^{2})/2) (x_{2} − 5) ,
u_{2} = 1 + 5

2π exp ((1 − r^{2})/2) (x_{1} − 5) ,

& periodic boundary conditions as an example, where r = p

(x_{1} − 5)^{2} + (x_{2} − 5)^{2}

**Accuracy test in 2D**

Grids used for this smooth vortex flow test

0 5 10

0 2 4 6 8 10

0 5 10

0 2 4 6 8 10

0 5 10

0 2 4 6 8

Grid 1 Grid 2 10 Grid 3

kE_{z}k_{1,∞} = kzcomput ^{− z}exact^{k}^{1,∞} denotes discrete 1- or
maximum-norm error for state variable z

Results shown below are at time t = 10 on N × N mesh

**Accuracy results in 2D: Grid 1**

N E_{1}(ρ) Order E_{1}(u_{1}) Order E_{1}(u_{2}) Order E_{1}(p) Order

40 0.6673 2.3443 1.7121 0.8143

80 0.1792 1.90 0.6194 1.92 0.4378 1.97 0.2128 1.94 160 0.0451 1.99 0.1537 2.01 0.1104 1.99 0.0536 1.99 320 0.0113 2.00 0.0384 2.00 0.0276 2.00 0.0134 2.00

N E^{∞}(ρ) Order E^{∞}(u_{1}) Order E^{∞}(u_{2}) Order E^{∞}(p) Order

40 0.1373 0.3929 0.1810 0.1742

80 0.0377 1.87 0.1014 1.95 0.0502 1.85 0.0482 1.85 160 0.0093 2.02 0.0248 2.03 0.0123 2.03 0.0119 2.02 320 0.0022 2.07 0.0062 2.00 0.0030 2.04 0.0029 2.04

**Accuracy results in 2D: Grid 2**

N E_{1}(ρ) Order E_{1}(u_{1}) Order E_{1}(u_{2}) Order E_{1}(p) Order

40 0.9298 2.6248 2.1119 1.2104

80 0.2643 1.81 0.7258 1.85 0.5296 2.00 0.3277 1.89 160 0.0674 1.97 0.1833 1.99 0.1309 2.02 0.0845 1.96 320 0.0169 2.00 0.0458 2.00 0.0327 2.00 0.0212 1.99

N E^{∞}(ρ) Order E^{∞}(u_{1}) Order E^{∞}(u_{2}) Order E^{∞}(p) Order

40 0.1676 0.4112 0.2259 0.2111

80 0.0471 1.83 0.1242 1.73 0.0645 1.79 0.0586 1.85 160 0.0126 1.91 0.0333 1.90 0.0162 2.02 0.0149 1.97 320 0.0033 1.93 0.0085 1.97 0.0040 2.00 0.0038 1.98

**Accuracy results in 2D: Grid 3**

N E_{1}(ρ) Order E_{1}(u_{1}) Order E_{1}(u_{2}) Order E_{1}(p) Order

40 4.8272 4.7734 5.3367 5.4717

80 1.5740 1.62 1.5633 1.61 1.5660 1.77 1.5634 1.81 160 0.4536 1.79 0.4559 1.78 0.4537 1.79 0.4560 1.78 320 0.1215 1.90 0.1221 1.90 0.1222 1.89 0.1221 1.90

N E^{∞}(ρ) Order E^{∞}(u_{1}) Order E^{∞}(u_{2}) Order E^{∞}(p) Order

40 0.4481 0.4475 0.4765 0.4817

80 0.1170 1.94 0.1181 1.92 0.1196 1.99 0.1191 2.02 160 0.0434 1.43 0.0431 1.45 0.0442 1.43 0.0440 1.44 320 0.0117 1.89 0.0119 1.86 0.0119 1.89 0.0118 1.89

**Future perspectives**

Extension of Cartesian embedded boundary method in 3D will be based on either FronTier or Chombo

framework developed by Stony Brook, or ANAG, LBL Free-stream perserving PLIC-VOF moving mesh

method is under development

Low Mach solver for cavitation problems is ongoing

## Thank you

**Basic grid-metric relations**

Assume existence of inverse transformation

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

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

∂(τ, ~ξ)

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

∂(τ, ~ξ)

−1

yielding in d = 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}

= 1 J

J 0 0 0

J_{01} J_{11} J_{21} J_{31}
J_{02} J_{12} J_{22} J_{32}

**Basic grid-metric relations**

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} = −

N^{d}

X

i=1

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

and so grid-metric relations between different coordinates

∇ξ_{j} = (∂_{t}ξ_{j}, ∇_{~}_{x}ξ_{j}) = (∂_{t}ξ_{j}, ∂_{x}1ξ_{j}, ∂_{x}2ξ_{j}, ∂_{x}3ξ_{j}) = 1

J (J_{0j}, J_{1j}, J_{2j}, J_{3j})

**Basic grid-metric relations**

Note in two dimensions d = 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}

+ ∂

J ∂ξ_{2}

= ∂

−∂x_{1}

+ ∂

∂x_{1}

= 0

**Basic grid-metric relations**

Here we have

G = q

∂J

∂τ +

Xd j=1

∂

∂ξ_{j}

J ∂ξ_{j}

∂t

+

Xd j=1

f_{j}

" _{d}
X

k=1

∂

∂ξ_{k}

J ∂ξ_{k}

∂x_{j}

#

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

∂J

∂τ +

Xd j=1

∂

∂ξ_{j}

J ∂ξ_{j}

∂t

= 0 (geometrical conservation law) Xd

k=1

∂

∂ξ_{k}

J ∂ξ_{k}

∂x_{j}

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

and hence G = 0