u 1 ρp u 2

Computing

Interface Motion in

Hyperbolic Problems

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

Problem Setup

Consider an interface I in 2D (shown below) that separate the domain into two parts with different states q_L and q_R Motion and dynamics of interface depends, clearly, on

Mathematical model Interface conditions

Initial & boundary cond.

Discretization scheme

q_L

q_R

~t ~n I

Mathematical Model

As an example, we consider problems governed by compressible Euler equations in N ≥ 1 dimension

∂

∂t





 ρ ρu_i

E





 +

N

X

j=1

∂

∂x_j







ρu_j

ρu_iu_j + pδ_ij Eu_j + pu_j





 =







0

−ρ∂_xiφ

−ρ~u · ∇φ





 , i = 1, . . . , N

ρ p = p(ρ, e) e E = ρe + ρ PN

j=1 u²_j/2 φ

density pressure int. energy total energy gr. potential

Mathematical Model (Cont.)

For the ease of the latter reference, we write the above compressible flow model as

∂q

∂t +

N

X

j=1

∂f_j(q)

∂x_j = ψ(q),

where

q = [ρ, ρu1, . . . , ρu_N, E]^T ,

f_j = ρu_j, ρu1u_j + pδ_j1, . . . , ρu_Nu_j + pδ_jN, Eu_j + pu_jT

, ψ = [0, −ρ∂_x1φ, . . . , −ρ∂_xNφ, −ρ~u · ∇φ]^T

Interface Conditions

Interface conditions of interest are such as Dynamic condition

p_R − p_L = 0 (or = f (κ_I)) Kinematic condition

(~u_R − ~u_L) · ~n = 0, (~u_R − ~u_L) · ~t = 0 (or 6= 0) Vacuum (or called cavitation) condition

(~u_R − ~u_L) · ~n <

Z _ρL

ρvL

c

ρ dρ +

Z _ρR

ρvR

c

ρ dρ

κI ρv c

curvature at interface density at vacuum sound speed

Aim of Talk

With the above compressible flow model & interface conditions in mind, the plan of this talk is to:

Discuss anomales of state-of-the-art methods for solving the following interface problems:

Slip line (shear flow) problem Material line problem

Cavitation problem

Discuss mapped grid approach for problems with complex geometries

Slip Line (Shear Flow) Problem

For simplicity, we consider a plane interface moving

vertically from the left to right in x¹-direction as shown in N = 2 below. Interface conditions for this problem are

Dynamic condition: p_R = p_L

Kinematic condition: u1,R = u1,L & u2,R − u2,L
6= 0

slip line

¯ u¹

¯

u² −u¯2

Slip Line Problem (Cont.)

Ignore gravitational-potential term, we may derive following transport equations for motion of ρ, ρu2, & ρe + ρu²₂/2,

respectively, as

∂ρ

∂t + ¯u1

∂ρ

∂x1

= 0 (from mass conservation law)

∂

∂t (ρu2) + ¯u1

∂

∂x1

(ρu2) = 0 (from momentum in x2-dir)

∂

∂t (ρe) + ¯u1

∂

∂x1

(ρe) +

∂

∂t

 1

2ρu²₂



+ ¯u1

∂

∂x1

 1

2ρu²₂



= 0 (from total energy)

Slip Line Problem (Cont.)

Assume a single phase ideal gas flow with EOS p(ρ, e) = (γ − 1)ρe; γ > 1 is ratio of specific heats

In this instance, to ensure pressure equilibrium, as it should be for this slip line problem, from

∂

∂t

 p γ − 1



+ ¯u1

∂

∂x1

 p γ − 1

 +

∂

∂t

 1

2ρu²₂



+ ¯u¹ ∂

∂x1

 1

2ρu²₂



= 0,

we should have

∂

∂t

 1

2ρu²₂



+ ¯u1

∂

∂x1

 1

2ρu²₂



= 0

Slip Line Problem (Cont.)

Note that when the problem is solved numerically by an interface-capturing method, it is common to compute pressure, p, based on conservative variables as

p = (γ − 1) E −

P²

i=1(ρu_i)² 2ρ

! ,

while generally (ρu2)²/2ρ 6= ρu²₂/2 when a slip line is smeared out

Immediate consequence of this is loss of pressure equilibrium, and so incorrect solution of other state variables (see below)

Slip Line Problem (Cont.)

Suppose a first order Godunov method of the form Qⁿ⁺¹_ij = Qⁿ_ij − ∆t

∆x1

F1 Qⁿ_i,j, Qⁿ_i+1,j
− F1 Qⁿ_i−1,j, Qⁿ_ij


is used for numerical discretization of this slip line problem Here Qⁿ_ij denotes approximate value of cell average of

solution q over cell C_ij at time t_n, ^i.e., Qⁿ_ij = 1

∆x1∆x2

Z Z

Cij

q(x1, x2, t_n)dx1dx2,

and F¹(Q_i±1,j, Q_ij) denotes numerical flux, say, defined by F1(Q_i±1,j, Q_ij) = f1

q_i±1/2,j^∗ 

(q_i±1/2,j^∗ Riemann solution)

Slip Line Problem (Cont.)

An example obtained by using a Godunov-type method Errors depend strongly on transverse velocity jump

0 0.2 0.4 0.6 0.8 1

1 2 3 4 5 6

0 0.2 0.4 0.6 0.8 1

1 1.05 1.1 1.15 1.2 1.25 1.3 1.35

0 0.2 0.4 0.6 0.8 1

9.95 10 10.05 10.1 10.15

0 0.2 0.4 0.6 0.8 1

−1

−0.5 0 0.5 1

x1

x1

u1 u 2

ρ p

Slip Line Problem (Cont.)

To devise a more accurate method for numerical resolution of slip lines, we may use

Diffuse interface approach

Include transverse kinetic energy equation in the model & use its solution for pressure update

p = (γ − 1)



E − (ρu1)²

2ρ + ρu²₂ 2



This transverse kinetic equation should be modified so that there is no difficulty to work with shock waves

Sharp interface approach

Front tracking or Lagrangian moving grid method

Slip Line Problem (Cont.)

Result obtained using a Lagrangian-like method Errors are on the order of machine epsilon

0 0.5 1 1.5

0 1 2 3 4 5 6

0 0.5 1 1.5

1 1 1 1 1

0 0.5 1 1.5

10 10 10 10 10 10 10 10

0 0.5 1 1.5

−1

−0.5 0 0.5 1

x1

x1

u1 u 2

ρ p

Slip Line Problem (Cont.)

Grid system for a Lagrangian run of slip line problem

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.05 0.1 0.15

Initial grid

0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2

0 0.05 0.1 0.15 0.2

Final grid

Euler Equations in Generalized Coord.

Lagrangian results shown above is based on solving Euler equations in generalized coordinates of the form

∂

∂τ









ρJ ρJu_i

JE







 +

N

X

j=1

∂

∂ξj

J









ρU_j

ρu_iU_j + p^∂ξ_∂x^j

i

EUj + pUj − p^∂ξ_∂t^j









=









0

−ρJ _∂x^∂φ

i

−ρJ~u · ∇φ









which is as a result of coordinate change (~x, t) 7→ (~ξ, τ ) via

dt = dτ,

dx¹ = u^g1 dτ + a¹ dξ¹ + a² dξ² + a³ dξ³, dx² = u^g₂ dτ + b¹ dξ¹ + b² dξ² + b³ dξ³, dx³ = u^g3 dτ + c¹ dξ¹ + c² dξ² + c³ dξ³.

Euler in Generalized Coord. (Cont.)

For convenience, we write

∂ ˜q

∂τ +

N

X

j=1

∂ ˜f_j(˜q)

∂ξ_j = ˜ψ(˜q), where

˜

q = [ρJ, ρu¹J, . . . , ρuNJ, EJ]^T ,

f˜_j = J [ρU_j, ρu¹U_j + p∂_x1ξ_j, . . . , ρu_NU_j + p∂_xNξ_j, EU_j + pU_j − p∂_tξ_j]^T , ψ = [0, −ρJ∂˜ x₁φ, . . . , −ρJ∂x_Nφ, −ρJ~u · ∇φ]^T

U_j = ∂_tξ_j + P_N

i=1 u_i∂_xiξ_j: contravariant velocity in ξ_j-direction

~u^g = h~u: speed of grid motion, h ≈ 1

Material Line Problem

Now let us move on to our second interface-only problem that concerns a material line, separating regions of two different fluid phases

Again for simplicity, we consider a plane interface moving vertically from the left to right in x1-direction.

Interface conditions for this problem are the same as in slip line case.

Dynamic condition: p_R = p_L (but there is jump in material quantities)

Kinematic condition: u1,R = u1,L & u2,R − u2,L
6= 0

Material Line Problem (Cont.)

Assume ideal gas law p(ρ, e) = (γ − 1)ρe for each fluid phase, where γ takes different quantities for each phase In this instance, to ensure pressure equilibrium, from

∂

∂t

 p γ − 1



+ ¯u1

∂

∂x¹

 p γ − 1

 +

∂

∂t

 1

2ρu²₂



+ ¯u1

∂

∂x¹

 1

2ρu²₂



= 0,

we should have

∂

∂t

 1 2ρu²2



+¯u¹ ∂

∂x¹

 1 2ρu²2



= 0 & ∂

∂t

 1 γ − 1



+ ¯u¹ ∂

∂x¹

 1 γ − 1



= 0

M a te ri a l L in e P ro b le m (C o n t.)

Comparisonwithdifferentγmodel(nou2 jump)

2.5 4.0 5.5

2.5 4.0 5.5

2.5 4.0 5.5

2.5 4.0 5.5

1.0 1.2

1.0 1.2

1.0 1.2

1.0 1.2

1.00 1.06

1.00 1.06

1.00 1.06

1.00 1.06

0.00.40.8

1.20 1.30 1.40

0.00.40.8

1.20 1.30 1.40

0.00.40.8

1.20 1.30 1.40

0.00.40.8

1.20 1.30 1.40

withργwithρ/(γ−1)withγwith1/(γ−1)

ρe ρe ρe ρe

u1

u1

u1

u1

p p p p

γ γ γ γ

DES2010,DepartmentofMathematics,NTU,8-9January,2010–p.20/44

Material Line Problem (Cont.)

Comments on problems with complex equation of state Even for single phase flow case, standard method will fail to attain pressure equilibrium, say, for example, van der Waals gas

p(ρ, e) = γ − 1

1 − bρ ρe + aρ²
− aρ²

For multiphase flow case, suitable transport equations for characterizing numerical fluid mixing can be derived (see in a later)

Cavitation Problem: 1D

Homogeneous flow with standard atmospheric condition

Due to opposite flow motion, pressure drop & formation of cavitation zone

u1 −u1

Cavitation test: 1D

Dynamic cavitation formation (gas-liquid mixture case)

−10 0 1

5

10x 10⁴

−1 0 1

−100

−50 0 50 100

u 1

−10 0 1

500 1000

ρ

−10 0 1

0.5 1

cavitation

p α

x x

Cavitation test: 1D

Cavitation formation after t = 0⁺ (elastic-plastic case)

−4 −2 0 2 4

7.8 7.85 7.9 7.95 8

ρ (Mg/m3 )

−4 −2 0 2 4

−1

−0.5 0 0.5 1

u (km/s)

−4 −2 0 2 4

0 2 4 6 8

x (m) σ x (GPa)

−4 −2 0 2 4

−0.6

−0.4

−0.2 0 0.2 0.4 0.6

x (m) s x (GPa)

Cavitating Problem (Cont.)

Isentropic relaxation model (Saurel ^{et al.} JCP ’09)

∂

∂t (α¹ρ¹) +

N

X

j=1

∂

∂xj

(α¹ρ¹uj) = 0,

∂

∂t (α²ρ²) +

N

X

j=1

∂

∂xj

(α²ρ²uj) = 0,

∂

∂t (ρu_i) +

N

X

j=1

∂

∂xj

(ρu_iu_j + pδ_ij) = 0, i = 1, . . . , N

∂α¹

∂t +

N

X

j=1

u_j ∂α¹

∂xj

= µ (p¹ (ρ¹) − p² (ρ²))

Each phasic pressure p_ι satisfies p_ι(ρ) = (p⁰ + B) (ρ/ρ⁰)^γ − B

Mixture pressure p satisfies p = α¹p¹ + α²p², µ : parameter

Cavitating Problem (Cont.)

There are some funadmental issues for numerical

resolution cavitating problems that I do not show here. But in general the following two things are important for the

problem

Devise of general cavitation-capturing method ?

Relaxation model to elastic-plastic & to other physical laws

Mapped Grid Method

To solve problem with complex geometries, here we consider body-fitted mapped grid method for its easy extension to three dimensions

Employ finite volume formulation of numerical solution

Qⁿ_ij ≈ 1

∆ξ¹∆ξ² Z

C_ij

q(ξ¹, ξ², τ_n) dξ¹ξ²

i − 1 i − 1

i

i j

j

j + 1 j + 1

C_ij Cˆ_ij

ξ¹ ξ²

mapping

∆ξ¹

∆ξ² logical domain physical domain

←−

x¹ = x¹(ξ¹, ξ²) x² = x²(ξ¹, ξ²) x¹

x²

Mapped Grid Method (Cont.)

In three dimensions N = 3, equations to be solved take

∂q

∂t +

N

X

j=1

∂ ˜fj

∂ξ_j = ψ(q)

A simple dimensional-splitting method based on wave propagation approach of LeVeque ^{et al.} is used, ^i.e.,

Solve one-dimensional Riemann problem normal 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 resolution

Shock waves over circular array

A Mach 1.42 shock wave in water over a circular array

t=0

shock

Shock waves over circular array

Grid system

−3 −2 −1 0 1 2 3 4 5

−4

−3

−2

−1 0 1 2 3 4

time=0

Shock waves over circular array

Contours for density

t=0.0024

Cavitation test: 2D steady state

Uniform gas-liquid mixture with speed 600m /s over a circular region

Cavitation test: 2D steady state

Pseudo colors of volume fraction

Formation of cavitation zone (Onset–shock induced, diffusion, . . . ?)

Cavitation test: 2D steady state

Pseudo colors of pressure

Smooth transition across liquid-gas phase boundary

Cavitation test: 2D steady state

Pseudo colors of volume fraction: 2 circular case

Convergence of solution as the mesh is refined ?

t=0.05

Moving cylindrical vessel

Initial setup

time=0

air helium

interface

Moving cylindrical vessel

Solution at time t = 0.25

Moving cylindrical vessel

Solution at time t = 0.5

Moving cylindrical vessel

Solution at time t = 0.75

Moving cylindrical vessel

Solution at time t = 1

Shock wave & Disperse phases

Shock wave & Disperse phases

Shock wave & Disperse phases

Shock wave & Disperse phases

Shock wave & Disperse phases

Shock wave & Disperse phases

Shock wave & Disperse phases

Shock wave & Disperse phases

Compressible Multiphase Flow

Homogeneous equilibrium pressure & velocity across material interfaces

Volume-fraction based model equations (Shyue JCP

’98, Allaire ^{et al.} JCP ’02)

∂

∂t (αiρi) + 1 J

N_d

X

j=1

∂

∂ξ_j (αiρiUj) = 0, i = 1, 2, . . . , mf

∂

∂t (ρu_i) + 1 J

N_d

X

j=1

∂

∂ξj

(ρu_iU_j + pJ_ji) = 0, i = 1, 2, . . . , N_d,

∂E

∂t + 1 J

N_d

X

j=1

∂

∂ξj

(EUj + pUj) = 0,

∂α_i

∂t + 1 J

N_d

XU_j ∂α_i

∂ξ = 0, i = 1, 2, . . . , m_f − 1;

Thank You