compressible multiphase flow

(1)

Volume of fluid methods for

compressible multiphase flow

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

(2)

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

(3)

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 − α) ǫ₂, where

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

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

∇α

|∇α|

(4)

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

(5)

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

(6)

Baer-Nunziato 2-phase model

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

(α₁ρ₁)_t + ∇ · (α₁ρ₁~u₁) = 0

(α₁ρ₁~u₁)_t + ∇ · (α₁ρ₁~u₁ ⊗ ~u₁) + ∇(α₁p₁) = p₀∇α₁ + λ(~u₂ − ~u₁) (α₁ρ₁E₁)_t + ∇ · (α₁ρ₁E₁~u₁ + α₁p₁~u₁) = p₀ (α₂)_t + λ~u₀ · (~u₂ − ~u₁) (α₂ρ₂)_t + ∇ · (α₂ρ₂~u₂) = 0

(α₂ρ₂~u₂)_t + ∇ · (α₂ρ₂~u₂ ⊗ ~u₂) + ∇(α₂p₂) = p₀∇α₂ − λ(~u₂ − ~u₁)

(α₂ρ₂E₂) _t + ∇ · (α₂ρ₂E₂~u₂ + α₂p₂~u₂) = −p₀ (α₂)_t − λ~u₀ · (~u₂ − ~u₁) (α₂)_t + ~u₀ · ∇α₂ = µ (p₂ − p₁)

α_k = V_k/V : volume fraction (α₁ + α₂ = 1), ρ_k: density,

~u_k: velocity, p_k: pressure, E_k = e_k + ~u²_k/2: specific total energy, e_k: specific internal energy, k = 1, 2

(7)

Baer-Nunziato 2-phase model

Quantities p₀ & ~u₀ denote interfacial pressure & velocity Baer & Nunziato (1986)

p₀ = p₂, ~u₀ = ~u₁ Saurel & Abgrall (1999)

p₀ = P₂

k=1 α_kp_k, ~u₀ = P₂

k=1 α_kρ_k~u_k

P₂

k=1 α_kρ_k

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

pressure, respectively

(8)

Reduced 2-phase model

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

∂_t (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0

∂_t (α₂ρ₂) + ∇ · (α₂ρ₂~u) = 0

∂_t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

∂_t (ρE) + ∇ · (ρE~u + p~u) = 0

∂_tα₁ + ~u · ∇α₁ = ρ₂c²₂ − ρ₁c²₁ P2

k=1 ρ_kc²_k/α_k ∇ · ~u

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

Positivity & sharpness of α is important for accurate solution

(9)

Cartesian embedded boundary method

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

∂_tq + ∇ · ~f (q) = 0

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

q ∈ lR^m: vector of m conserved quantities

f = (f~ ₁, f₂, . . . , f_d) ∈ lR^m×d: flux matrix; f_j ∈ lR^m

Γ: embedded boundary

Γ →

(10)

Finite-volume discretization

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

Qⁿ_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 discretization ∆x = ∆x , j = 1, . . . , d

A B

C D

(11)

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_Sf (q)~ ^k · ~n^k = 1 κ_S∆x

X

k∈∂S

h^k_SF˘^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

(12)

Finite-volume discretization

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

Qⁿ⁺¹_S = Qⁿ_S − ∆tD_S · ~f (q) = Qⁿ_S − ∆t κ_S∆x

X

k∈∂S

h^k_SF˘^k

which satisfies discrete conservation identity X

∪S

κ_SQⁿ⁺¹ = X

∪S

κ_SQⁿ − ∆t

∆x

X

∪S

X

k∈∂S

h^k_SF˘^k

Method is subject to stringent CFL stability constraint

∆t = O

(κ_min)^1/d∆x

|λ|_max

≪ 1 if κ_min ≪ 1 (Small-cell problem)

(13)

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

(14)

Flux calculation: Wave propagation

On 2D uniform grid, high-resolution in fluctutation form Qⁿ⁺¹_ij = Qⁿ_ij − ∆t

∆x₁

A⁺₁ ∆Q_i−¹

2,j + A⁻₁ ∆Q_i+¹

2,j

− ∆t

∆x₂

A⁺₂ ∆Q_i,j−¹

2 + A⁻₂ ∆Q_i,j+¹

2

− ∆t

∆x₁

F˜_i+¹ 1

2,j − ˜F_i−¹ 1 2,j

− ∆t

∆x₂

F˜_i,j+² 1 2

− ˜F_i,j−² 1 2

A⁺₁ ∆Q

A⁻₁ ∆Q

A⁻₂ A⁺₁ ∆Q

A⁺₂ A⁻₁ ∆Q

A⁺₂ A⁺₁ ∆Q

A⁻₂ A⁻₁ ∆Q

(15)

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^±₁ ∆Q =

mw

X

k=1

λ^1,k^±

W^1,k, A^±₂ ∆Q =

mw

X

k=1

λ^2,k^±

W^2,k

A^±₂ A^±₁ ∆Q =

mw

X

p=1

λ^1,p^±

A^±₂ W^1,p =

mw

X

p=1

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

q=1

λ^2,q^±

W^2,q

F˜¹ := F˜¹ + 1 2

mw

X

k=1

_λ^1,k

1 − ∆t

∆x₁

_λ^1,k

Wf^1,k

To aviod oscillations near discontinuities, a wave limiter is applied leading to limited waves _Wf

(16)

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ⁿ⁺¹_S := Qⁿ⁺¹_S − M (Z^p ∩ S)

M(S) W^p

Z↓^p

A B

C

D Z↓^p

A B

C D

(17)

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

(18)

Wave propagation with boundary

High resolution with piecewise linear wave update

wave before propagation after propagation

a) b)

c) d)

α_pr_p/2

λ_p∆ t

(19)

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

(20)

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ⁿ⁺¹_b at time step t_n+1 is computed by

~uⁿ⁺¹_b = ~uⁿ_b + ∆t

M_bF~_bⁿ,

where M_b is mass of object & F^~_bⁿ is force on object

(21)

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

(22)

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

t = 0

t = 0.1641s

t = 0.30085s

(23)

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

(24)

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

(25)

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

(26)

Interface reconstruction: Example

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

0 0

0

0 0

0

1 0.29

0.68 0.09

0.51 0.51

0.09 0.68

0.29 interface↓

(27)

Interface reconstruction: Example

With α in previous slide & ~u = (1, 1) updated volume fraction over ∆t = 0.06 shown here ( κ₁ = 5.7 × 10⁻³ &

κ₂ = 1.3 × 10⁻³)

New reconstructed interface location

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 κ¹ κ₂

↑

↓

old interface

new interface

(28)

Solution update

Solution is double-valued on cells subdivided by VOF tracked interface; Again use finite volume formulation of solution Qⁿ_S for approximate value of solution q

Qⁿ_S ≈ 1 M(S)

Z

S

q(x, t_n) dV

C E

D F

G H

(29)

Wave propagation method

Normal- & transverse-edge wave propagation

↓

Z_p Z_p

↓ ↓

Z_pq Z_pq

(30)

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

(31)

Shock in air & R22 bubble

(32)

Shock in air & R22 bubble

(33)

Shock in air & R22 bubble

(34)

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)

(35)

Underwater explosion

(36)

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

(37)

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)⁻¹ should be set

(38)

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)⁻¹ 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

(39)

Reduced 2-phase model (revisit)

Physical balance laws

∂

∂τ







α₁ρ₁J α₂ρ₂J ρu_iJ

ρEJ





 ⁺

Xd j=1

∂

∂ξ_j J







α₁ρ₁U_j α₂ρ₂U_j ρu_iU_j + p_∂x^∂ξ^j

i

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





 ^{= 0}

∂α₁

∂τ +

Xd j=1

U_j ∂α₁

∂ξ_j = ρ₂c²₂ − ρ₁c²₁ P₂

k=1 ρ_kc²_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(α₂, α₁ρ₁, α₂ρ₂, ρe) with isobaric closure: p₁ = p₂ = p

(40)

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

(41)

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

∂τ = ε∇²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

(42)

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

(43)

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 →

(44)

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

(45)

Water-vapor cavitation

Initially, in closed shock tube, flow is homogeneous (contains α = 10⁻⁶ 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

(46)

Water-vapor cavitation

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1x 10⁻³

time

cavi grid

(47)

Underwater explosions

t=0.2ms

t=0.4ms

t=0.8ms

t=1.2ms

Moving mesh Eulerian mesh air

water

(48)

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

(49)

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γπ² exp (1 − r²)

1/(γ−1)

, p = ρ^γ,

u₁ = 1 − 5

2π exp ((1 − r²)/2) (x₂ − 5) , u₂ = 1 + 5

2π exp ((1 − r²)/2) (x₁ − 5) ,

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

(x₁ − 5)² + (x₂ − 5)²

(50)

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

(51)

Accuracy results in 2D: Grid 1

N E₁(ρ) Order E₁(u₁) Order E₁(u₂) Order E₁(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₁) Order E^∞(u₂) 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

(52)

Accuracy results in 2D: Grid 2

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

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

(53)

Accuracy results in 2D: Grid 3

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

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

(54)

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

(55)

Thank you

(56)

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ξ₁ ∂_x1ξ₁ ∂_x2ξ₁ ∂_x3ξ₁

∂_tξ₂ ∂_x1ξ₂ ∂_x2ξ₂ ∂_x3ξ₂







= 1 J







J 0 0 0

J₀₁ J₁₁ J₂₁ J₃₁ J₀₂ J₁₂ J₂₂ J₃₂







(57)

Basic grid-metric relations

Here

J =

∂(x₁, x₂, x₃)

∂(ξ₁, ξ₂, ξ₃)

= det

∂(x₁, x₂, x₃)

∂(ξ₁, ξ₂, ξ₃)

, J₁₁ =

∂(x₂, x₃)

∂(ξ₂, ξ₃)

, J₂₁ =

∂(x₁, x₃)

∂(ξ₃, ξ₂)

, J₃₁ =

∂(x₁, x₂)

∂(ξ₂, ξ₃) ,

J₁₂ =

∂(x₂, x₃)

∂(ξ₃, ξ₁)

, J₂₂ =

∂(x₁, x₃)

∂(ξ₁, ξ₃)

, J₃₂ =

∂(x₁, x₂)

∂(ξ₃, ξ₁) ,

J₁₃ =

∂(x₂, x₃)

∂(ξ₁, ξ₂)

, J₂₃ =

∂(x₁, x₃)

∂(ξ₂, ξ₁)

, J₃₃ =

∂(x₁, x₂)

∂(ξ₁, ξ₂) ,

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, ∂_x1ξ_j, ∂_x2ξ_j, ∂_x3ξ_j) = 1

J (J_0j, J_1j, J_2j, J_3j)

(58)

Basic grid-metric relations

Note in two dimensions d = 2, we have

∂ξ₁

∂t , ∂ξ₁

∂x₁ , ∂ξ₁

∂x₂

= 1 J

−∂x₁

∂τ

∂x₂

∂ξ₂ + ∂x₂

∂τ

∂x₁

∂ξ₂ , ∂x₂

∂ξ₂ , −∂x₁

∂ξ₂

∂t , ∂ξ₂

∂x₁ , ∂ξ₂

∂x₂

= 1 J

∂x₁

∂τ

∂x₂

∂ξ₁ − ∂x₂

∂τ

∂x₁

∂ξ₁ , −∂x₂

∂ξ₁ , ∂x₁

∂ξ₁

J = ∂x₁

∂ξ₁

∂x₂

∂ξ₂ − ∂x₁

∂ξ₂

∂x₂

∂ξ₁

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

∂J

∂τ + ∂

∂ξ₁

J ∂ξ₁

∂t

+ ∂

∂ξ₂

J ∂ξ₂

∂t

= 0

∂

∂ξ₁

J ∂ξ₁

∂x₁

+ ∂

∂ξ₂

J ∂ξ₂

∂x₁

= ∂

∂ξ₁

∂x₂

∂ξ₂

+ ∂

∂ξ₂

−∂x₂

∂ξ₁

= 0

∂

J ∂ξ₁

+ ∂

J ∂ξ₂

= ∂

−∂x₁

+ ∂

∂x₁

= 0

(59)

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