**Adaptive**

**moving-mesh relaxation scheme** **for**

**compressible two-phase barotropic** **flow with cavitation**

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

**Outline**

Introduce simple relaxation model for compressible 2-phase barotropic flow with & without cavitation

Describe finite-volume method for solving proposed model numerically on fixed mapped (body-fitted) grids Discuss adaptive moving mesh approach for efficient improvement of numerical resolution

Show sample results & discuss future extensions

**Conventional barotropic 2-phase model**

Assume homogeneous equilibrium barotropic flow with single pressure p & velocity ~u in 2-phase mixture region

Equations of motion

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

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

∂_{t} (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (Momentum ρ~u)
Saturation condition for volume fraction α_{k}, k = 1, 2

α_{1} + α_{2} = 1 =⇒ α_{1}ρ_{1}

ρ_{1}(p) + α_{2}ρ_{2}

ρ_{2}(p) = 1

yielding scalar nonlinear equation for equilibrium p

**Equilibrium speed of sound**

Non-monotonic sound speed ceq defined as 1

ρc^{2}eq

= α_{1}

ρ_{1}c^{2}_{1} + α_{2}

ρ_{2}c^{2}_{2} (Wood’s formula)

exists in 2-phase region, air-water case shown below;

yielding numerical difficulty such as inaccurate wave transmission across diffused interface

0 0.2 0.4 0.6 0.8 1

0 500 1000 1500

α

c eq

**Equilibrium vs. frozen sound speed**

Equilibrium (solid) & frozen (dashed) sound speeds,
c^{2}

f ^{=}

PY_{k}c^{2}_{k}, in case of passive advection of air-water
interface

**Equilibrium vs. frozen sound speed**

When acoustic wave interacts with numerical diffusion zone, sound speeds difference leads to time delay τ of transmitted waves through interface

**Relaxation barotropic 2-phase model**

To overcome these difficulties, we consider a relaxation

model proposed by Caro, Coquel, Jamet, Kokh, Saurel ^{et al.}

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

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

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

∂_{t}α_{1} + ~u · ∇α_{1} = µ (p_{1} (ρ_{1}) − p_{2} (ρ_{2})) (Transport α_{1})
When taking infinite pressure relaxation µ → ∞, we have

p_{1}(ρ_{1}) = p_{2}(ρ_{2}) =⇒ p_{1}

α_{1}ρ_{1}
α_{1}

− p_{2}

α_{2}ρ_{2}
1 − α_{1}

= 0
yielding scalar nonlinear equation for volume fraction α_{1}

**Relaxation barotropic model: Remarks**

This model can be viewed as isentropic version of

relaxation model proposed by Saurel, Petitpas, Berry (JCP 2009, see below) for 2-phase flow

This model is hyperbolic & has monotonic sound speed
c^{2}_{f} = P

Y_{k}c^{2}_{k}

Cavitation is modeled as a simplified mechanical

relaxation process, occurring at infinite rate & not as a mass transfer process

*i.e.*, cavitation pockets appear as volume fraction

increases for a small amount of gas present initially

**Relaxation 2-phase model**

Saurel, Petitpas, Berry (JCP 2009)

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

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

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

∂_{t} (α_{1}ρ_{1}e_{1}) + ∇ · (α_{1}ρ_{1}e_{1}~u) + α_{1}p_{1}∇ · ~u = −pµ (p¯ _{1} − p_{2})

∂_{t} (α_{2}ρ_{2}e_{2}) + ∇ · (α_{2}ρ_{2}e_{2}~u) + α_{2}p_{2}∇ · ~u = pµ (p¯ _{1} − p_{2})

∂_{t}α_{1} + ~u · ∇α_{1} = µ (p_{1} − p_{2})

¯

p = (p_{1}Z_{2} + p_{2}Z_{1})/(Z_{1} + Z_{2}); Z_{k} = ρ_{k}c_{k}

Model agrees with reduced 2-phase model of Kapila,

Menikoff, Bdzil, Son, Stewart (Phys. Fluid 2001) formally

**1-phase barotropic cavitation models**

Cutoff model

p =

(p(ρ) if ρ ≥ ρ_{sat}
p_{sat} if ρ < ρ_{sat}
Schmidt model

p =

p(ρ) if ρ ≥ ρ_{sat}

p_{sat} + p_{gl} ln

ρgc^{2}_{g}ρlc^{2}_{l} (ρl+α(ρg−ρl))
ρl(^{ρ}^{g}^{c}^{2}^{g}^{−}^{α}(^{ρ}^{g}^{c}^{2}^{g}^{−}^{ρ}^{l}^{c}^{2}l ))

if ρ < ρ_{sat}

with p_{gl} = ρ_{g}c^{2}_{g}ρ_{l}c^{2}_{l} (ρ_{g} − ρ_{l})

ρ^{2}_{g}c^{2}_{g} − ρ^{2}_{l} c^{2}_{l} & ρ = αρ_{g} + (1 − α)ρ_{l}
Modified Schmidt model & its variant

**Mapped grid method**

We want to use finite-volume mapped grid approach to solve proposed relaxation model in complex geometry Assume mapped grids are logically rectangular & will

review method for hyperbolic system of conservation laws

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

i − 1 i − 1

i

i j

j

j + 1 j + 1

Cˆ_{ij}
C_{ij}

ξ1

ξ2

mapping

∆ξ1

∆ξ2

logical grid physical grid

←−

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

x2

**Mapped grid methods**

On a curvilinear grid, a finite volume method takes
Q^{n+1}_{ij} = Q^{n}_{ij}− ∆t

κ_{ij}∆ξ_{1}

F_{i+}^{1} 1

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

− ∆t
κ_{ij}∆ξ_{2}

F_{i,j+}^{2} 1
2

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

∆ξ_{1}, ∆ξ_{2} denote mesh sizes in ξ_{1}- & ξ_{2}-directions

κ_{ij} = M(C_{ij})/∆ξ_{1}∆ξ_{2} is area ratio between areas of grid
cell in physical & logical spaces

F_{i−}^{1} _{1}

2,j = γ_{i−}^{1}

2,jF˘_{i−}^{1}

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

2

= γ_{i,j−}^{1}

2

F˘_{i,j−}^{1}

2 are normal fluxes
per unit length in logical space with γ_{i−}^{1}

2,j = h_{i−}^{1}

2,j/∆ξ_{1} &

γ_{i,j−}^{1}

2 = h_{i,j−}^{1}

2 /∆ξ_{2} representing length ratios

**Wave propagation method**

First order wave propagation method devised by LeVeque on mapped grid is a Godunov-type finite volume method

Q^{n+1}_{ij} = Q^{n}_{ij}− ∆t
κ_{ij}∆ξ_{1}

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

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

2,j

−

∆t
κ_{ij}∆ξ_{2}

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

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

2

with right-, left-, up-, & down-moving fluctuations ^{A}^{+}_{1} ∆Q_{i−}^{1}

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

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

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

2 that are entering into grid cell

To determine these fluctuations, one-dimensional Riemann problems in direction normal to cell edges are solved

**Wave propagation method**

Speeds & limited versions of waves are used to calculate second order correction terms as

Q^{n+1}_{ij} := Q^{n+1}_{ij} − 1
κ_{ij}

∆t

∆ξ1

Fe_{i+}^{1} _{1}

2,j − Fe^{1}

i−^{1}_{2},j

− 1
κ_{ij}

∆t

∆ξ2

Fe_{i,j+}^{2} _{1}

2

− Fe^{2}

i,j−^{1}_{2}

For example, at cell edge (i − ^{1}_{2}, j) correction flux takes
e

F^{1}

i−^{1}_{2},j = 1
2

Nw

X

k=1

^{λ}^{1,k}_{i−}^{1}

2,j

^{1 −} _{κ} ^{∆t}

i−^{1}_{2},j∆ξ_{1}

^{λ}^{1,k}_{i−}^{1}

2,j

!
f
W^{1,k}

i−_{2}^{1},j

κ_{i−}^{1}

2,j = (κ_{i−1,j} + κ_{ij})/2. To aviod oscillations near

discontinuities, a wave limiter is applied leading to limited waves Wf

**Wave propagation method**

Transverse wave propagation is included to ensure second
order accuracy & also improve stability that ^{A}^{±}_{1} ∆Q_{i−}^{1}

2,j are each split into two transverse fluctuations: up- &

down-going ^{A}^{±}_{2} ^{A}^{+}_{1} ∆Q_{i−}^{1}

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

2,j, at each cell edge

This method can be shown to be conservative & stable

under a variant of CFL (Courant-Friedrichs-Lewy) condition of form

ν = ∆t max

i,j,k

λ^{1,k}_{i−}^{1}

2,j

κ_{i}_{p}_{,j}∆ξ1

,

λ^{2,k}_{i,j−}^{1}

2

κ_{i,j}_{p}∆ξ2

≤ 1,

i_{p} = i if λ^{1,k}

i−^{1}_{2},j > 0 & i − 1 if λ^{1,k}

i−^{1}_{2},j < 0

**Extension to moving mesh**

To extend mapped grid method to solution adaptive moving grid method one simple way is to take approach proposed by

H. Tang & T. Tang, Adaptive mesh methods for one- and two-dimensional hyperbolic conservation laws, SIAM J.

Numer. Anal., 2003

In each time step, this moving mesh method consists of three basic steps:

(1) Mesh redistribution

(2) Conservative interpolation of solution states (3) Solution update on a fixed mapped grid

**Mesh redistribution scheme**

Winslow’s approach (1981)

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

for ξ(x). Coefficient D is a positive definite matrix which may depend on solution gradient

Variational approach (Tang & many others)

Solve ^{∇}_{ξ} ^{·} (D∇_{ξ}x_{j}) = 0, j = 1, . . . , N_{d}
for x(ξ) that minimizes “energy” functional

E(x(ξ)) = 1 2

Z

Ω
N_{d}

X

j=1

∇^{T}_{ξ} D∇x_{j}dξ

Lagrangian (ALE)-type approach (* ^{e.g.}*, CAVEAT code)

**Conservative interpolation**

Numerical solutions need to be updated conservatively, ^{i.e.}

XM C^{k+1}

Q^{k+1} = X

M C^{k}
Q^{k}

after each redistribution iterate k. This can be done by Finite-volume approach

M C^{k+1}

Q^{k+1} = M C^{k}

Q^{k} −

Ns

X

j=1

h_{j}G˘_{j}, G˘ = ( ˙x · n)Q

Geometric approach

"

X

S

M C_{p}^{k+1} ∩ S_{p}^{k}#

Q^{k+1}_{C} = X

S

M C_{p}^{k+1} ∩ S_{p}^{k}
Q^{k}_{S}

C_{p}, S_{p} are polygonal regions occupied by cells C & S

**Interpolation-free moving mesh**

If we want to derive an interpolation-free moving mesh method, one may first consider coordinate change of equations via (x, t) 7→ (ξ, t), yielding transformed

conservation law as

∂_{t}q˜ + ∇_{ξ} · ˜f = G

˜

q = Jq, f˜_{j} = J (q ∂_{t}ξ_{j} + ∇ξ_{j} · f) , J = det(∂ξ/∂x)^{−}^{1}
G = q

∂_{t}J + ∇_{ξ} · (J∂_{t}ξ_{j})
+

XN j=1

f_{j}∇_{ξ} · J∂_{x}_{j}ξ_{k}

= 0 (if GCL & SCL are satisfied)

Numerical method can be devised easily to solve these equations

**Relaxation solver on moving meshes**

In each time step, our numerical method for solving

barotropic 2-phase flows on a moving mesh consists of following steps:

(1) Moving mesh step

Determine cell-interface velocity & cell-interface location in physical space over a time step

(2) Frozen step µ → 0

Solve homogeneous part of relaxation model on a

moving mapped grid over same time step as in step 1 (3) Relaxation step µ → ∞

Solve model system with only source terms in infinite relaxation limit

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

**Water-vapor cavitation**

Physical grid in x-t plane

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1x 10^{−3}

time

x

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

Physical grid in x-t plane

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

time

x

**Oscillating water column**

**Oscillating water column**

Time evolution of pressure at left & right boundaries

0 1 2 3 4 5 6 7 8 9 10

0.8 1 1.2 1.4 1.6 1.8

p L

0 1 2 3 4 5 6 7 8 9 10

0.8 1 1.2 1.4 1.6 1.8

p R

mmesh mesh0

mmesh mesh0

time

**Tait equation of state & parameters**

Each fluid phase (liquid & gas) satisfies Tait equation of state

p_{k}(ρ) = (p_{0k} + B_{k})

ρ
ρ_{0k}

γk

− B_{k} for k = 1, 2.
with parameters for liquid phase as

(γ, B, ρ_{0}, p_{0})_{1} =

7, 3000 bar, 10^{3} kg/m^{3}, 1 bar
while parameters for gas phase as

(γ, B, ρ_{0}, p_{0})_{2} =

1.4, 0, 1 kg/m^{3}, 1 bar

**Supersonic flow over a bluntbody**

Formation of cavitation zone

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Volume fraction

0
0.5
1
1.5
2
2.5
x 10^{8}

Pressure

**Underwater explosion**

Contours of density and pressure at selected times

Density t=0.24ms

air

water

Pressure

t=0.4ms

**Underwater explosion**

Contours of density and pressure at selected times

t=0.8ms

t=1.2ms

**Underwater explosion**

Physical grid at selected times

−1 0 1

−1

−0.5 0 0.5

−1 0 1

−1

−0.5 0 0.5

t = 0.2ms t = 0.4ms

−1

−0.5 0 0.5

−1

−0.5 0 0.5

t = 0.8ms t = 1.2ms

**Final remarks**

Show preliminary results obtained using relaxation moving mesh methods for barotropic 2-phase flow problem

Cavitation problems are often occurred in low Mach

scenario & so suitable fixed up such as preconditioning

& others are necessary for solution accuracy improvement

Extension to non-barotropc cavitation with phase transition

. . .

## Thank you

**Reduced 2-phase flow model**

Reduced 2-phase flow model of Kapila * ^{et al.}* is zero-order
approximation of Baer-Nunziato equations with stiff

mechanical relaxation that takes

∂_{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}
P_{2}

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

**Baer-Nunziato Two-Phase Flow Model**

Baer & Nunziato (J. Multiphase Flow 1986)
(α_{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 : specific internal energy, k = 1, 2

**Baer-Nunziato Model (Cont.)**

p_{0} & ~u_{0}: 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}
λ & µ (> 0): relaxation parameters that determine rates at
which velocities and pressures of two phases reach

equilibrium