moving-mesh relaxation scheme for

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 (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0 (Continuity α₁ρ₁)

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

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

α₁ + α₂ = 1 =⇒ α₁ρ₁

ρ₁(p) + α₂ρ₂

ρ₂(p) = 1

yielding scalar nonlinear equation for equilibrium p

Equilibrium speed of sound

Non-monotonic sound speed ceq defined as 1

ρc²eq

= α₁

ρ₁c²₁ + α₂

ρ₂c²₂ (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²

f ⁼

PY_kc²_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 (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0

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

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

∂_tα₁ + ~u · ∇α₁ = µ (p₁ (ρ₁) − p₂ (ρ₂)) (Transport α₁) When taking infinite pressure relaxation µ → ∞, we have

p₁(ρ₁) = p₂(ρ₂) =⇒ p₁

α₁ρ₁ α₁



− p₂

 α₂ρ₂ 1 − α₁



= 0 yielding scalar nonlinear equation for volume fraction α₁

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²_f = P

Y_kc²_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 (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0

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

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

∂_t (α₁ρ₁e₁) + ∇ · (α₁ρ₁e₁~u) + α₁p₁∇ · ~u = −pµ (p¯ ₁ − p₂)

∂_t (α₂ρ₂e₂) + ∇ · (α₂ρ₂e₂~u) + α₂p₂∇ · ~u = pµ (p¯ ₁ − p₂)

∂_tα₁ + ~u · ∇α₁ = µ (p₁ − p₂)

¯

p = (p₁Z₂ + p₂Z₁)/(Z₁ + Z₂); Z_k = ρ_kc_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²_gρlc²_l (ρl+α(ρg−ρl)) ρl(^ρgc2g−α(^{ρgc2g−ρlc2}l ))



if ρ < ρ_sat

with p_gl = ρ_gc²_gρ_lc²_l (ρ_g − ρ_l)

ρ²_gc²_g − ρ²_l c²_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

∂_tq + ∇ · 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₁ = x1(ξ1, ξ₂) x2 = x2(ξ1, ξ2) x1

x2

Mapped grid methods

On a curvilinear grid, a finite volume method takes Qⁿ⁺¹_ij = Qⁿ_ij− ∆t

κ_ij∆ξ₁

F_i+¹ 1

2,j − F_i−¹ 1 2,j



− ∆t κ_ij∆ξ₂

F_i,j+² 1 2

− F_i,j−² 1 2



∆ξ₁, ∆ξ₂ denote mesh sizes in ξ₁- & ξ₂-directions

κ_ij = M(C_ij)/∆ξ₁∆ξ₂ is area ratio between areas of grid cell in physical & logical spaces

F_i−¹ ₁

2,j = γ_i−¹

2,jF˘_i−¹

2,j, F_i,j−² ₁

2

= γ_i,j−¹

2

F˘_i,j−¹

2 are normal fluxes per unit length in logical space with γ_i−¹

2,j = h_i−¹

2,j/∆ξ₁ &

γ_i,j−¹

2 = h_i,j−¹

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ⁿ⁺¹_ij = Qⁿ_ij− ∆t κ_ij∆ξ₁



A⁺₁ ∆Q_i−¹

2,j + A⁻₁ ∆Q_i+¹

2,j



−

∆t κ_ij∆ξ₂



A⁺₂ ∆Q_i,j−¹

2 + A⁻₂ ∆Q_i,j+¹

2



with right-, left-, up-, & down-moving fluctuations ^A+₁ ∆Q_i−¹

2,j, A⁻₁ ∆Q_i+¹

2,j, A⁺₂ ∆Q_i,j−¹

2, & A⁻₂ ∆Q_i,j+¹

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ⁿ⁺¹_ij := Qⁿ⁺¹_ij − 1 κ_ij

∆t

∆ξ1

Fe_i+¹ ₁

2,j − Fe¹

i−¹₂,j



− 1 κ_ij

∆t

∆ξ2

Fe_i,j+² ₁

2

− Fe²

i,j−¹₂



For example, at cell edge (i − ¹₂, j) correction flux takes e

F¹

i−¹₂,j = 1 2

Nw

X

k=1

^λ1,k_i−¹

2,j

 ^{1 −} _κ ^∆t

i−¹₂,j∆ξ₁ 

^λ1,k_i−¹

2,j

! f W^1,k

i−₂¹,j

κ_i−¹

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±₁ ∆Q_i−¹

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

down-going ^A±₂ ^A+₁ ∆Q_i−¹

2,j & ^A±₂ ^A−₁ ∆Q_i−¹

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

2,j

  κ_ip,j∆ξ1

, 

λ^2,k_i,j−¹

2

  κ_i,jp∆ξ2



 ≤ 1,

i_p = i if λ^1,k

i−¹₂,j > 0 & i − 1 if λ^1,k

i−¹₂,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_jdξ

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_jG˘_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

∂_tq˜ + ∇_ξ · ˜f = G

˜

q = Jq, f˜_j = J (q ∂_tξ_j + ∇ξ_j · f) , J = det(∂ξ/∂x)⁻¹ G = q 

∂_tJ + ∇_ξ · (J∂_tξ_j) +

XN j=1

f_j∇_ξ · J∂_xjξ_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⁻⁶ 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⁻³

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, ρ₀, p₀)₁ = 

7, 3000 bar, 10³ kg/m³, 1 bar while parameters for gas phase as

(γ, B, ρ₀, p₀)₂ = 

1.4, 0, 1 kg/m³, 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⁸

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 (α₁ρ₁) + ∇ · (α₁ρ₁~u) = 0

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

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

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

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

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

Baer-Nunziato Two-Phase Flow Model

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

Baer-Nunziato Model (Cont.)

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

equilibrium