**Wave-propagation based methods** **for**

**compressible homogeneous** **two-phase flow**

### Keh-Ming Shyue

### Department of Mathematics

### National Taiwan University

**Outline**

### Eulerian formulation

### Mathematical models

### Wave-propagation based volume tracking method Sample examples

### Generalized Lagrangian formulation Mathematical models

### Flux-based wave decomposition method Sample examples

### Future work

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 2/72

**Two-Phase Flow Model (I)**

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

**Two-Phase Flow Model (II)**

### Saurel & Gallouet (1998) (α 1 ρ 1 ) _{t} + ∇ · (α 1 ρ 1 ~u 1 ) = m ˙

### (α 1 ρ 1 ~u 1 ) _{t} + ∇ · (α 1 ρ 1 ~u 1 ⊗ ~u 1 ) + ∇(α 1 p 1 ) = p 0 ∇α 1 + m~u ˙ 0 + F _{d}

### (α 1 ρ 1 E 1 ) _{t} + ∇ · (α 1 ρ 1 E 1 ~u 1 + α 1 p 1 ~u 1 ) = p 0 (α 2 ) _{t} + mE ˙ 0 + F _{d} ~u 0 + Q 0

### (α 2 ρ 2 ) _{t} + ∇ · (α 2 ρ 2 ~u 2 ) = − m ˙

### (α 2 ρ 2 ~u 2 ) _{t} + ∇ · (α 2 ρ 2 ~u 2 ⊗ ~u 2 ) + ∇(α 2 p 2 ) = p 0 ∇α 2 − m~u ˙ 0 − F _{d}

### (α ^{2} ρ ^{2} E ^{2} ) _{t} + ∇ · (α ^{2} ρ ^{2} E ^{2} ~u ^{2} + α ^{2} p ^{2} ~u ^{2} ) = −p ^{0} (α ^{2} ) _{t} − mE ˙ ^{0} − F _{d} ~u ^{0} − Q ^{0} (α 2 ) _{t} + ~u 0 · ∇α 2 = µ (p 2 − p 1 )

### ˙

### m : mass transfer, F _{d} : drag force Q 0 : convective heat exchange

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 4/72

**Two-Phase Flow 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

**Two-Phase Flow Model: Derivation**

### Standard way to derive these equations is based on averaging theory of Drew (Theory of Multicomponent Fluids, D.A. Drew & S. L. Passman, Springer, 1999) Namely, introduce indicator function χ _{k} as

### χ _{k} (M, t) =

### ( 1 if M belongs to phase k

### 0 otherwise

### Denote < ψ > as volume averaged for flow variable ψ , hψi = 1

### V Z

### V

### ψ dV Gauss & Leibnitz rules

### hχ _{k} ∇ψi = h∇(χ _{k} ψ)i−hψ∇χ _{k} i & hχ _{k} ψ _{t} i = h(χ _{k} ψ) _{t} i−hψ(χ _{k} ) _{t} i

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 6/72

**Two-Phase Flow Model (cont.)**

### Take product of each conservation law with χ _{k} & perform averaging process. In case of mass conservation equation, for example, we have

### hχ _{k} ρ _{k} i _{t} + ∇· < χ _{k} ρ _{k} ~u _{k} >= hρ _{k} (χ _{k} ) _{t} + ρ _{k} ~u _{k} · ∇χ _{k} i Since χ _{k} is governed by

### (χ _{k} ) _{t} + ~u 0 · ∇χ _{k} = 0 (~u 0 : interface velocity ), this leads to mass averaged equation for phase k

### hχ ρ i + ∇· < χ ρ ~u >= hρ (~u − ~u ) · ∇χ i

**Two-Phase Flow Model (cont.)**

### In summary, averaged model system, we have, are

### hχ _{k} ρ _{k} i _{t} + ∇· < χ _{k} ρ _{k} ~u _{k} >= hρ _{k} (~u _{k} − ~u 0 ) · ∇χ _{k} i hχ _{k} ρ _{k} ~u _{k} i _{t} + ∇· < χ _{k} ρ _{k} ~u _{k} ⊗ ~u _{k} > +∇ hχ _{k} p _{k} i = hp _{k} ∇χ _{k} i +

### hρ _{k} ~u _{k} (~u _{k} − ~u 0 ) · ∇χ _{k} i hχ _{k} ρ _{k} E _{k} i _{t} + ∇· < χ _{k} ρ _{k} E _{k} ~u _{k} + χ _{k} p _{k} ~u _{k} >= hp _{k} ~u _{k} · ∇χ _{k} i +

### hρ _{k} E (~u _{k} − ~u 0 ) · ∇χ _{k} i hχ _{k} i _{t} + h~u _{k} · ∇χ _{k} i = h(~u _{k} − ~u 0 ) · ∇χ _{k} i

### Note: existence of various interfacial source terms

### Mathematical as well as numerical modelling of these terms are important (but difficult) for general multiphase flow

### problems

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 8/72

**Reduced Two-Phase Flow Model**

### Murrone & Guillard (JCP 2005)

### Assume λ = λ ^{′} /ε & µ = µ ^{′} /ε , λ ^{′} = O(1) & µ ^{′} = O(1) Apply formal asymptotic analysis to Baer &

### Nunziato’s model, as ε → 0 , gives leading order approximation

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

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

### (ρE) _{t} + ∇ · (ρE~u + p~u) = 0 ( mixture total energy )

**Reduced Two-Phase Model (cont.)**

### Remarks:

### 1 . In this case, p 1 → p 2 & ~u 1 → ~u 2 , as ε → 0 , which means the flow is homogeneous ( 1 -pressure & 1 -velocity) with p _{ι} = p & ~u _{ι} = ~u , ι = 0, 1, 2 , across interfaces

### 2 . Mixture equation of state: p = p(α 2 , α 1 ρ 1 , α 2 ρ 2 , ρe) 3 . Isobaric closure: p 1 = p 2 = p

### For some EOS, explicit formula for p is available (examples are given next)

### For some other EOS, p is found by solving coupled equations

### p 1 (ρ 1 , ρ 1 e 1 ) = p 2 (ρ 2 , ρ 2 e 2 ) & α 1 ρ 1 e 1 + α 2 ρ 2 e 2 = ρe

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 10/72

**Reduced Two-Phase Model (cont.)**

### Polytropic ideal gas: p _{k} = (γ _{k} − 1)ρ _{k} e _{k}

### ρe =

### 2

### X

### k=1

### α _{k} ρ _{k} e _{k} =

### 2

### X

### k=1

### α _{k} p

### γ _{k} − 1 =⇒

### p = ρe

### ^{2} X

### k=1

### α _{k}

### γ _{k} − 1

**Reduced Two-Phase Model (cont.)**

### Polytropic ideal gas: p _{k} = (γ _{k} − 1)ρ _{k} e _{k}

### ρe =

### 2

### X

### k=1

### α _{k} ρ _{k} e _{k} =

### 2

### X

### k=1

### α _{k} p

### γ _{k} − 1 =⇒

### p = ρe

### ^{2} X

### k=1

### α _{k} γ _{k} − 1

### Van der Waals gas: p _{k} = ( _{1−b} ^{γ} ^{k} ^{−1}

### k ρ k )(ρ _{k} e _{k} + a _{k} ρ ^{2} _{k} ) − a _{k} ρ ^{2} _{k} ρe =

### 2

### X

### k=1

### α _{k} ρ _{k} e _{k} =

### 2

### X

### k=1

### α _{k} 1 − b _{k} ρ _{k} γ _{k} − 1

### (p + a _{k} ρ ^{2} _{k} ) − a _{k} ρ ^{2} _{k}

### =⇒

### p =

### "

### ρe −

### 2

### X

### k=1

### α _{k} 1 − b _{k} ρ _{k}

### γ _{k} − 1 − 1

### a _{k} ρ ^{2} _{k}

### # ^{2} X

### k=1

### α _{k} 1 − b _{k} ρ _{k} γ _{k} − 1

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 11/72

**Reduced Two-Phase Model (cont.)**

### Two-molecular vibrating gas: p _{k} = ρ _{k} R _{k} T (e _{k} ) , T satisfies

### e = RT

### γ − 1 + RT vib

### exp T vib /T − 1

### As before, we now have

### ρe =

### 2

### X

### k=1

### α _{k} ρ _{k} e _{k} =

### 2

### X

### k=1

### α _{k}

###

###

### ρ _{k} R _{k} T _{k} γ _{k} − 1

### + ρ _{k} R _{k} T vib ,k

### exp

### T vib ,k /T _{k}

### − 1

###

###

### =

### 2

### X

### k=1

### α _{k}

###

###

### p γ _{k} − 1

### + p vib ,k

### exp

### p vib ,k /p

### − 1

###

### (Nonlinear eq.)

**Reduced Model: Remarks**

### 4 . It can be shown entropy of each phase S _{k} now satisfies DS _{k}

### Dt = ∂S _{k}

### ∂t + ~u · ∇S _{k} = 0, for k = 1, 2 5 . Model system is hyperbolic under suitable

### thermodynamic stability condition

### 6 . When α _{k} = 0 , ρ _{k} can not be recovered from α _{k} & α _{k} ρ _{k} , and so take α _{k} ∈ [ε, 1 − ε] , ε ≪ 1

### 7 . Other model systems exist in the literature that are more robust for homogeneous flow (examples)

### 8 . When individual pressure law differs in form (see

### below), new mixture pressure law should be devised first & construct model equations based on that

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 13/72

**Barotropic & Non-Barotropic Flow**

### Fluid component 1 : Tait EOS

### p(ρ) = (p 0 + B) ρ ρ 0

### γ

### − B

### Fluid component 2 : Noble-Abel EOS

### p(ρ, e) = γ − 1 1 − bρ

### ρe

### Mixture pressure law (Shyue, Shock Waves 2006)

###

### (p 0 + B) ρ γ

### − B if α = 1

**Barotropic Two-Phase Flow**

### Fluid component ι : Tait EOS

### p(ρ) = (p 0 ι + B _{ι} )

### ρ ρ 0 ι

### γ ι

### − B _{ι} , ι = 1, 2

### Mixture pressure law (Shyue, JCP 2004)

### p =

###

###

###

###

###

### (p 0 ι + B _{ι} )

### ρ ρ 0 ι

### γ ι

### − B _{ι} if α = α _{ι} (0 or 1) (γ − 1) ρ

### e + B ρ 0

### − γB if α ∈ (0, 1)

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 15/72

**Homogeneous Two-Phase Model**

### In summary, mathematical model for compressible homogeneous two-phase flow:

### Equations of motion

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

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

### (α ) + ~u · ∇α = α α ρ 1 c ^{2} _{1} − ρ 2 c ^{2} _{2} !

### ∇ · ~u

**Wave Propagation Method**

### Finite volume formulation of wave propagation method, Q ^{n} _{S} gives approximate value of cell average of solution q over cell S at time t _{n}

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

### Z

### S

### q(X, t _{n} ) dV

### M(S) : measure (area in 2 D or volume in 3 D) of cell S

### C E

### D F

### G H

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 17/72

**Wave Propagation Method (cont.)**

### First order version: Piecewise constant wave update Godunov-type method: Solve Riemann problem at each cell interface in normal direction & use resulting waves to update cell averages

### Q ^{n+1} _{S} := Q ^{n+1} _{S} − M (W _{p} ∩ S)

### M(S) R _{p} , R _{p} being jump from RP

### W ↓ _{p}

**Wave Propagation Method (cont.)**

### First order version: Transverse-wave included

### Use transverse portion of equation, solve Riemann problem in transverse direction, & use resulting

### waves to update cell averages as usual

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

### ↓ ↓

### W _{pq} W _{pq}

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 19/72

**Wave Propagation Method (cont.)**

### High resolution version: Piecewise linear wave update wave before propagation after propagation

a) b)

c) d)

α_{p}r

p/2

λ_{p}∆ t

**Volume Tracking Algorithm**

### 1. Volume moving procedure *(a) Volume fraction update*

### Take a time step on current grid to update cell averages of volume fractions at next time step *(b) Interface reconstruction*

### Find new interface location based on volume fractions obtained in (a) using an interface reconstruction scheme. Some cells will be subdivided & values in each subcell must be initialized.

### 2. Physical solution update

### Take same time interval as in (a), but use a method to update cell averages of multicomponent model on new grid created in (b)

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 21/72

**Interface Reconstruction Scheme**

### Given volume fractions on current grid, piecewise linear interface reconstruction (PLIC) method does:

### 1. Compute interface normal

### Gradient method of Parker & Youngs Least squares method of Puckett

### 2. Determine interface location by iterative bisection

### 0 0

### 0 0

### 0 0

### 0 0

### 0

### 0 .09 0 .51 0 .29

### ↓

### ↓ interface

### interface

### Data set Parker & Youngs Puckett

**Volume Moving Procedure**

### (a) Volume fractions given in previous slide are updated with uniform (u, v) = (1, 1) over ∆t = 0.06

### (b) New interface location is reconstructed

### 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 5(−3) 1(−3)

### ↑

### ↓

### old interface

### new interface

### (a) (b)

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 23/72

**Surface Moving Procedure**

### Solve Riemann problem at tracked interfaces & use

### resulting wave speed of the tracked wave family over ∆t to find new location of interface at the next time step

### o

### o

### o o

### o

### ↑

### ↑ old front

### old front

### new front

**Boundary Conditions**

### For tracked segments representing rigid (solid wall)

### boundary (stationary or moving), reflection principle is used to assign states for fictitious subcells in each time step:

### z _{C} := z _{E} (z = ρ, p, α)

### ~u _{C} := ~u _{E} − 2(~u _{E} · ~n)~n + 2(~u 0 · ~n)

### ~u 0 : moving boundary velocity

### C

### E F

### G H

### ~n

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 25/72

**Interface Conditions**

### For tracked segments representing material interfaces, pressure equilibrium as well as velocity continuity

### conditions across interfaces are fulfilled by 1. Devise of the wave-propagation method

### 2. Choice of Riemann solver used in the method

**Stability Issues**

### Choose time step ∆t based on uniform grid mesh size

### ∆x , ∆y as

### ∆t max _{p,q} (λ _{p} , µ _{q} )

### min(∆x, ∆y) ≤ 1,

### λ _{p} , µ _{q} : speed of p -wave, q -wave from Riemann

### problem solution in normal-, transverse-directions Use large time step method of LeVeque ( ^{i.e.} , wave

^{i.e.}

### interactions are assumed to behave in linear manner) to maintain stability of method even in the presence of

### small Cartesian cut cells

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

^{et al.}

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 27/72

**Shock-Bubble Interaction**

**Shock-Bubble Interaction**

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72

**Shock-Bubble Interaction**

**Shock-Bubble Interaction**

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72

**Shock-Bubble Interaction**

**Shock-Bubble Interaction**

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 28/72

**Shock-Bubble Interaction**

**Shock-Bubble Interaction**

**Shock-Bubble Interaction**

**Shock-Bubble Interaction**

**Shock-Bubble Interaction (cont.)**

### Approximate locations of interfaces

### time=55µs

### air R22

### time=115µs time=135µs

### time=187µs time=247µs time=200µs

### time=342µs time=417µs time=1020µs

**Shock-Bubble Interaction (cont.)**

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

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 30/72

**Underwater Explosions**

### Numerical schlieren images for density

**Underwater Explosions (cont.)**

### Approximate locations of interfaces

### time=0.2

### air

### water

### time=0.4

### air

### water

### time=0.8

### air

### water

### time=1.2

### air

### water

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 32/72

**Falling Rigid Object in Water Tank**

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Density air

### water

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Pressure

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Volume fraction

### t = 0ms

**Falling Rigid Object in Water Tank**

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Density air

### water

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Pressure

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Volume fraction

### t = 1ms

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 33/72

**Falling Rigid Object in Water Tank**

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Density air

### water

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Pressure

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Volume fraction

### t = 2ms

**Falling Rigid Object in Water Tank**

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Density air

### water

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Pressure

### −2 0 2

### −3

### −2

### −1 0 1 2 3

### Volume fraction

### t = 3ms

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 33/72

**Generalized Lagrangian Model**

### Introduce transformation (t, x, y) ↔ (τ, ξ, η) via

###

###

###

### dt dx dy

###

###

###

###

### =

###

###

###

###

### 1 0 0

### x _{τ} x _{ξ} x _{η} y _{τ} y _{ξ} y _{η}

###

###

###

###

###

###

###

### dτ dξ dη

###

###

###

###

### or

###

###

###

### dτ dξ dη

###

###

###

###

### =

###

###

###

###

### 1 0 0 ξ _{t} ξ _{x} ξ _{y} η _{t} η _{x} η _{y}

###

###

###

###

###

###

###

### dt dx dy

###

###

###

###

### Basic grid-metric relations:

###

###

###

###

### 1 0 0 ξ _{t} ξ _{x} ξ _{y} η _{t} η _{x} η _{y}

###

###

###

###

### =

###

###

###

###

### 1 0 0

### x _{τ} x _{ξ} x _{η} y _{τ} y _{ξ} y _{η}

###

###

###

###

### −1

### = 1 J

###

###

###

###

### x _{ξ} y _{η} − x _{η} y _{ξ} 0 0

### −x _{τ} y _{η} + y _{τ} x _{η} y _{η} −x _{η} x _{τ} y _{ξ} − y _{τ} x _{ξ} −y _{ξ} x _{ξ}

###

###

###

###

**G. Lagrangian Model (cont.)**

### Homogeneous two-phase model in N _{d} generalized coord.:

### ∂

### ∂τ (α 1 ρ 1 J) +

### N d

### X

### j=1

### ∂

### ∂ξ _{j} (α 1 ρ 1 JU _{j} ) = 0,

### ∂

### ∂τ (α 2 ρ 2 J) +

### N d

### X

### j=1

### ∂

### ∂ξ _{j} (α 2 ρ 2 JU _{j} ) = 0,

### ∂

### ∂τ (ρJu _{i} ) +

### N d

### X

### j=1

### ∂

### ∂ξ _{j} J

### ρu _{i} U _{j} + p ∂ξ _{j}

### ∂x _{i}

### = 0 for i = 1, 2, . . . , N _{d} ,

### ∂

### ∂τ (JE) +

### N d

### X

### j=1

### ∂

### ∂ξ _{j} J

### EU _{j} + pU _{j} − p ∂ξ _{j}

### ∂t

### = 0,

### ∂α 2

### ∂τ +

### N d

### X

### j=1

### U _{j} ∂α 2

### ∂ξ _{j} = 0, U _{j} = ∂ _{t} ξ _{j} +

### N d

### X

### i=1

### u _{i} ∂ _{x} _{i} ξ _{j}

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 35/72

**G. Lagrangian Model (cont.)**

### Continuity on mixed derivatives of grid coordinates gives geometrical conservation laws

### ∂

### ∂τ

###

###

###

###

###

###

### x _{ξ}

### y _{ξ} x _{η} y _{η}

###

###

###

###

###

###

###

### + ∂

### ∂ξ

###

###

###

###

###

###

###

### −x _{τ}

### −y _{τ} 0 0

###

###

###

###

###

###

###

### + ∂

### ∂η

###

###

###

###

###

###

### 0 0

### −x _{τ}

### −y _{τ}

###

###

###

###

###

###

###

### = 0

### with (x _{τ} , y _{τ} ) to be specified as, for example,

### Eulerian case: (x _{τ} , y _{τ} ) = ~0

**G. Lagrangian Model (cont.)**

### General 1 -parameter case: (x _{τ} , y _{τ} ) = h(u, v) , h ∈ [0, 1]

### At given time instance, h can be chosen based on

### Grid-angle preserving condition (Hui ^{et al.} JCP 1999)

^{et al.}

### ∂

### ∂τ cos ^{−1} ∇ξ

### |∇ξ| · ∇η

### |∇η|

### = ∂

### ∂τ cos ^{−1}

###

###

### −y _{η} x _{η} − y _{ξ} x _{ξ} q y _{ξ} ^{2} + y _{η} ^{2} q

### x ^{2} _{ξ} + x ^{2} _{η}

###

###

### = · · ·

### = Ah _{ξ} + Bh _{η} + Ch = 0 (1st order PDE )

### with

### A = q

### x ^{2} _{η} + y _{η} ^{2} (vx _{ξ} − uy _{ξ} ) , B = q

### x ^{2} _{ξ} + y _{ξ} ^{2} (uy _{η} − vx _{η} ) C = q

### x ^{2} _{ξ} + y _{ξ} ^{2} (u _{η} y _{η} − v _{η} x _{η} ) − q

### x ^{2} _{η} + y _{η} ^{2} (u _{ξ} y _{ξ} − v _{ξ} x _{ξ} )

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 37/72

**G. Lagrangian Model (cont.)**

### General 1 -parameter case: (x _{τ} , y _{τ} ) = h(u, v) , h ∈ [0, 1]

### Or alternatively, based on

### Mesh-area preserving condition

### ∂J

### ∂τ = ∂

### ∂τ (x _{ξ} y _{η} − x _{η} y _{ξ} )

### = x _{ξτ} y _{η} + x _{ξ} y _{ητ} − x _{ητ} y _{ξ} − x _{η} y _{ξτ}

### = · · ·

### = Ah _{ξ} + Bh _{η} + Ch = 0 (1st order PDE )

### with

**G. Lagrangian Model (cont.)**

### Numerics: h - or ˜h -equation constraint geometrical laws

### ∂

### ∂τ

###

###

###

###

###

###

### x _{ξ} y _{ξ} x _{η} y _{η}

###

###

###

###

###

###

###

### − ∂

### ∂ξ

###

###

###

###

###

###

### hu hv 0 0

###

###

###

###

###

###

###

### − ∂

### ∂η

###

###

###

###

###

###

### 0 0 hu hv

###

###

###

###

###

###

###

### = 0

### Usability: Mesh-area evolution equation

### ∂J

### ∂τ − ∂

### ∂ξ [h (uy _{η} − vx _{η} )] − ∂

### ∂η [h (vx _{ξ} − uy _{ξ} )] = 0

### Initial & boundary conditions for h - or ˜h -equation ?

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 38/72

**G. Lagrangian Model (cont.)**

### In summary, with (x _{τ} , y _{τ} ) = h 0 (u, v) & EOS, model system for homogeneous two-phase flow reads

### ∂

### ∂τ

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

### Jα 1 ρ 1

### Jα 2 ρ 2

### Jρu Jρv JE

### x _{ξ} y _{ξ} x _{η}

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

### + ∂

### ∂ξ

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

### Jα 1 ρ 1 U Jα 2 ρ 2 U JρuU + y _{η} p JρvU − x _{η} p

### JEU + (y _{η} u − x _{η} v)p

### −h 0 u

### −h 0 v 0

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

### + ∂

### ∂η

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

### Jα 1 ρ 1 V Jα 2 ρ 2 V JρuV − y _{ξ} p JρvV + x _{ξ} p

### JEV + (x _{ξ} v − y _{ξ} u)p 0

### 0

### −h u

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

###

### = 0

**G. Lagrangian Model (cont.)**

### Under thermodyn. stability condition, our multifluid model in generalized coordinates is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1

### Model system is written in quasi-conservative form with spatially varying fluxes in generalized coordinates

### Grid system is a time-varying grid

### Extension of the model to general non-barotropic

### multifluid flow can be made in an analogous manner

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 40/72

**Flux-based Wave Decomposition**

### In 2 D, equations to be solved takes the form

### ∂q

### ∂τ + f 1

### ∂

### ∂ξ , q, ∇~ ξ

### + f 2

### ∂

### ∂η , q, ∇~ ξ

### = ˜ ψ

### A simple dimensional-splitting approach based on f -wave formulation of LeVeque ^{et al.} is used

^{et al.}

### Solve one-dimensional generalized Riemann problem (defined below) at each cell interfaces Use resulting jumps of fluxes (decomposed into

### each wave family) of Riemann solution to update cell

### averages

**F-Waves Decomposition (cont.)**

### Employ finite volume formulation of numerical solution

### Q ^{n} _{ij} ≈ 1

### ∆ξ∆η Z

### C ij

### q(ξ, η, τ _{n} ) dA

### that gives approximate value of cell average of solution q over cell C _{ij} = [ξ _{i} , ξ _{i+1} ] × [η _{j} , η _{j+1} ] at time τ _{n}

−1 −0.5 0 0.5 1

−1.5

−1

−0.5 0

−1 −0.5 0 0.5

−2

−1.8

−1.6

−1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2

### x

### y

### ξ

### η

### computational grid physical grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 42/72

**Generalized Riemann Problem**

### Generalized Riemann problem of our multifluid model at cell interface ξ _{i−1/2} consists of the equation

### ∂q

### ∂τ + F _{i−} ^{1}

### 2 ,j

### ∂ _{ξ} , q, ∇~ ξ

### = 0 together with flux function

### F _{i−} ^{1}

### 2 ,j =

###

###

###

### f _{i−1,j}

### ∂ _{ξ} , q, ∇~ ξ

### for ξ < ξ _{i−1/2} f _{ij}

### ∂ _{ξ} , q, ∇~ ξ

### for ξ > ξ _{i−1/2}

### and piecewise constant initial data

**Generalized Riemann Problem**

### Generalized Riemann problem at time τ = 0

### q _{τ} + f _{i−1,j}

### ∂ _{ξ} , q, ∇~ ξ

### = 0 q τ + f _{i,j}

### ∂ _{ξ} , q, ∇~ ξ

### = 0 ξ τ

### Q ^{n} _{i−1,j} Q ^{n} _{ij}

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 44/72

**Generalized Riemann Problem**

### Exact generalized Riemann solution: basic structure

### τ

### n Q ^{n}

**Generalized Riemann Problem**

### Shock-only approximate Riemann solution: basic structure

### q _{τ} + f _{i−1,j}

### ∂ _{ξ} , q, ∇~ ξ

### = 0 q τ + f _{i,j}

### ∂ _{ξ} , q, ∇~ ξ

### = 0 ξ

### τ

### Q ^{n} _{i−1,j} Q ^{n} _{ij}

### λ ^{1} λ ^{2}

### λ ^{3} q _{mL} ^{−} q _{mL} ^{+}

### q _{mR}

### Z ^{1} = f _{L} (q _{mL} ^{−} ) − f _{L} (Q ^{n} _{i−1,j} ) Z ^{2} = f _{R} (q _{mR} ) − f _{R} (q _{mL} ^{+} )

### Z ^{3} = f _{R} (Q ^{n} _{ij} ) − f _{R} (q _{mR} )

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 44/72

**F-Waves Decomposition (cont.)**

### Basic steps of a dimensional-splitting scheme ξ -sweeps: solve

### ∂q

### ∂τ + f 1

### ∂

### ∂ξ , q, ∇~ ξ

### = 0

### updating Q ^{n} _{ij} to Q ^{∗} _{i,j} η -sweeps: solve

### ∂q

### ∂τ + f 2

### ∂

### ∂η , q, ∇~ ξ

### = 0

**F-Waves Decomposition (cont.)**

### That is to say,

### ξ -sweeps: we use

### Q ^{∗} _{ij} = Q ^{n} _{ij} − ∆τ

### ∆ξ

### F _{i+} ^{−} _{1}

### 2 ,j − F _{i−} ^{+} _{1}

### 2 ,j

### − ∆τ

### ∆ξ ˜ Z _{i+} ^{1}

### 2 ,j − ˜ Z _{i−} ^{1}

### 2 ,j

### with Z ˜ _{i−} ^{1}

### 2 ,j = 1 2

### m w

### X

### p=1

### sign

### λ ^{p} _{i−} _{1}

### 2 ,j

### 1 − ∆τ

### ∆ξ λ

### p

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

### Z ˜ ^{p} _{i−} _{1}

### 2 ,j

### η -sweeps: we use

### Q ^{n+1} _{ij} = Q ^{∗} _{ij} − ∆τ

### ∆η

### G _{i,j+} ^{−} _{1}

### 2

### − G ^{+}

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

### − ∆τ

### ∆η ˜ Z _{i,j+} ^{1}

### 2 − ˜ Z _{i,j−} ^{1}

### 2

### with Z ˜ _{i,j−} ^{1}

### 2 = 1 2

### m w

### X

### p=1

### sign λ ^{p}

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

### 1 − ∆τ

### ∆η λ

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

### Z ˜ ^{p}

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

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 46/72

**F-Waves Decomposition (cont.)**

### Flux-based wave decomposition

### f _{i,j} − f _{i−1,j} =

### m w

### X

### p=1

### Z _{i−1/2} ^{p} =

### m w

### X

### p=1

### λ ^{p} _{i−1/2} W _{i−1/2} ^{p}

### Some care should be taken on the limited jump of fluxes W ˜ ^{p} , for p = 2 (contact wave), in particular to ensure

### correct pressure equilibrium across material interfaces

### MUSCL-type (slope limited) high resolution extension is

### not simple as one might think of for multifluid problems

### Splitting of discontinuous fluxes at cell interfaces:

**Lax’s Riemann Problem**

### h 0 = 0 Eulerian result

### h 0 = 0.99 Lagrangian-like result

### sharper resolution for contact discontinuity

### 0 0.5 1

### 0 0.5 1 1.5 2

### x

### ρ

### Exact h 0 =0 h 0 =0.99

### 0 0.5 1

### 0 0.5 1 1.5 2

### x

### u

### 0 0.5 1

### 0 1 2 3 4

### x

### p

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 48/72

**Lax’s Riemann Problem**

### Physical grid coordinates at selected times

### Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system

### 0.05 0.1 0.15

### time

**Woodward-Colella’s Problem**

### h 0 = 0 Eulerian result

### h 0 = 0.99 Lagrangian-like result

### sharper resolution for contact discontinuity

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

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 50/72

**Woodward-Colella’s Problem**

### h 0 = 0 Eulerian result

### h 0 = 0.99 Lagrangian-like result

### sharper resolution for contact discontinuity

### 0 5 10 15 20

### ρ

### Fine grid h 0 =0 h 0 =0.99

### −5 0 5 10 15

### 0 200 400 600

### u p

### t = 0.032 t = 0.032

### t = 0.032

**Woodward-Colella’s Problem**

### h 0 = 0 Eulerian result

### h 0 = 0.99 Lagrangian-like result

### sharper resolution for contact discontinuity

### 0 0.5 1

### 0 2 4 6 8

### ρ

### 0 0.5 1

### −5 0 5 10 15

### 0 0.5 1

### 0 200 400 600

### u p

### x x

### x

### t = 0.038 t = 0.038

### t = 0.038

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 52/72

**Woodward-Colella’s Problem**

### Physical grid coordinates at selected times

### Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system

### 0.01 0.02 0.03 0.04

### time

### 2 **D Riemann Problem**

### With initial 4 -shock wave pattern

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

###

###

###

###

### ρ u v p

###

###

###

###

###

### =

###

###

###

###

### 1.5

### 0 0 1.5

###

###

###

###

###

###

###

###

###

###

### 0.532 1.206

### 0 0.3

###

###

###

###

###

###

###

###

###

###

### 0.138 1.206 1.206 0.029

###

###

###

###

###

###

###

###

###

###

### 0.532 0 1.206

### 0.3

###

###

###

###

###

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 54/72

### 2 **D Riemann Problem**

### With initial 4 -shock wave pattern Lagrange-like result

### Occurrence of simple Mach reflection

### 0.6 0.8

### 0.6 0.8

### 0.6 0.8

### Density Pressure Physical grid

### 2 **D Riemann Problem**

### With initial 4 -shock wave pattern Eulerian result

### Poor resolution around simple Mach reflection

### 0.2 0.4 0.6 0.8 0.2

### 0.4 0.6 0.8

### 0.2 0.4 0.6 0.8 0.2

### 0.4 0.6 0.8

### 0.2 0.4 0.6 0.8 0.2

### 0.4 0.6 0.8

### Density Pressure Physical grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 56/72

**More Examples**

### Two-dimensional case

### Radially symmetric problem Underwater explosion

### Shock-bubble interaction Helium bubble case

### Refrigerant bubble case Three-dimensional case

### Underwater explosion Shock-bubble interaction

### Helium bubble case

**Radially Symmetric Problem**

### 0 0.2 0.4

### 0 0.1 0.2 0.3 0.4 0.5

### 0 0.2 0.4

### 0 0.1 0.2 0.3 0.4 0.5

### 0 0.2 0.4

### 0 0.1 0.2 0.3 0.4

### Density Pressure 0.5 Physical grid

### gas

### liquid a) h 0 = 0.99

### 0 0.2 0.4

### 0 0.1 0.2 0.3 0.4 0.5

### 0 0.2 0.4

### 0 0.1 0.2 0.3 0.4 0.5

### 0 0.2 0.4

### 0 0.1 0.2 0.3 0.4

### Density Pressure 0.5 Physical grid

### gas

### liquid b) h 0 = 0

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 58/72

**Radially Symmetric Prob. (Cont.)**

### 0 0.1 0.2 0.3 0.4 0.5

### 0.2 0.4 0.6 0.8 1 1.2

### one−d h 0 =0 h 0 =0.99

### 0 0.1 0.2 0.3 0.4 0.5

### 0 0.1 0.2 0.3 0.4

### 0 0.1 0.2 0.3 0.4 0.5

### 0 0.5 1 1.5

### 0 0.1 0.2 0.3 0.4 0.5

### 0.5 1 1.5 2

### ρ (M g /m 3 ) ¯u (k m /s )

### p (G P a ) J

### gas liquid

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.9 , 800 × 500 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 60/72

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.9 , 800 × 500 grid

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.9 , 800 × 500 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 60/72

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.9 , 800 × 500 grid

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.9 , 800 × 500 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 60/72

**Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.9

### −1

### −0.5 0 0.5

### time=0ms

**Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.9

### −2 −1.5 −1 −0.5 0 0.5 1 1.5 2

### −1

### −0.5 0 0.5

### time=0.2ms

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 61/72

**Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.9

### −1

### −0.5 0 0.5

### time=0.4ms

**Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.9

### −2 −1.5 −1 −0.5 0 0.5 1 1.5 2

### −1

### −0.5 0 0.5

### time=0.8ms

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 61/72

**Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.9

### −1

### −0.5 0 0.5

### time=1.2ms

**Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 62/72

**Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

**Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 62/72

**Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

**Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 62/72

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

### −0.2

### −0.15

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.02

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 63/72

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.08

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

### −0.2

### −0.15

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.16

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 63/72

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.35

**Shock-Bubble (Refrigerant)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 64/72

**Shock-Bubble (Refrigerant)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

**Shock-Bubble (Refrigerant)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 64/72

**Shock-Bubble (Refrigerant)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

**Shock-Bubble (Refrigerant)**

### Numerical schlieren images: h 0 = 0.5 , 600 × 400 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 64/72

**Shock-Bubble (R22) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0

**Shock-Bubble (R22) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

### −0.2

### −0.15

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.02

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 65/72

**Shock-Bubble (R22) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.08

**Shock-Bubble (R22) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

### −0.2

### −0.15

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.16

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 65/72

**Shock-Bubble (R22) (Cont.)**

### Grid system (coarsen by factor 5 ) with h 0 = 0.5

### −0.1

### −0.05 0 0.05 0.1 0.15 0.2

### time=0.35

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.6 , 100 ^{3} grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 66/72

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.6 , 100 ^{3} grid

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.6 , 100 ^{3} grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 66/72

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.6 , 100 ^{3} grid

**Underwater Explosions**

### Numerical schlieren images h 0 = 0.6 , 100 ^{3} grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 66/72

### 3 **D Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

### −0.4

### −0.2 0 0.2 0.4 0.6 0.8 1

### time = 0

### 3 **D Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

−1

−0.5 0

0.5 1

−1

−0.5 0

0.5 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

### time = 0.25 ms

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 67/72

### 3 **D Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

### time = 0.5 ms

### 3 **D Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

−1

−0.5 0

0.5 1

−1

−0.5 0

0.5 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

### time = 1.0 ms

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 67/72

### 3 **D Underwater Explosions (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

### time = 1.5 ms

### 3 **D Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 68/72

### 3 **D Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid

### 3 **D Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 68/72

### 3 **D Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid

### 3 **D Shock-Bubble (Helium)**

### Numerical schlieren images: h 0 = 0.6 , 150 × 50 × 50 grid

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 68/72

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

### 0.3

### 0.4

### 0.5

### 0.6

### 0.7

### 0 0.05 0.1 0.15 0.2 0.25

### z

### time = 0

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

### 0

### 0.1

### 0.2

### 0.3

### 0.4

### 0.5

### 0.6

### 0.7

### 0.1 0 0.2

### 0 0.05 0.1 0.15 0.2 0.25

### x y

### z

### time = 0.02

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 69/72

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

### 0.3

### 0.4

### 0.5

### 0.6

### 0.7

### 0 0.05 0.1 0.15 0.2 0.25

### z

### time = 0.08

**Shock-Bubble (Helium) (Cont.)**

### Grid system (coarsen by factor 2 ) with h 0 = 0.6

### 0

### 0.1

### 0.2

### 0.3

### 0.4

### 0.5

### 0.6

### 0.7

### 0.1 0 0.2

### 0 0.05 0.1 0.15 0.2 0.25

### x y

### z

### time = 0.16

Department of Aeronautics and Astronautics, National Cheng Kung University, March 21, 2008 – p. 69/72