• 沒有找到結果。

Department of Mathematics

N/A
N/A
Protected

Academic year: 2022

Share "Department of Mathematics"

Copied!
138
0
0

加載中.... (立即查看全文)

全文

(1)

Wave-propagation based methods for

compressible homogeneous two-phase flow

Keh-Ming Shyue

Department of Mathematics

National Taiwan University

(2)

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

(3)

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 02 ) t − λ~u 0 · (~u 2 − ~u 1 )

(α 2 ) t + ~u 0 · ∇α 2 = µ (p 2 − p 1 )

(4)

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

(5)

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

(6)

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

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

(7)

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

k ρ k i t + ∇· < χ k ρ k ~u k >= hρ kk ) 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

(8)

Two-Phase Flow Model (cont.)

In summary, averaged model system, we have, are

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 +

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 +

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

(9)

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 )

(10)

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

(11)

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

(12)

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

(13)

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

(14)

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

(15)

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

(16)

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

(17)

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

(18)

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

(19)

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

(20)

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

(21)

Wave Propagation Method (cont.)

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

a) b)

c) d)

αpr

p/2

λp∆ t

(22)

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

(23)

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

(24)

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

(25)

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

(26)

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

(27)

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

(28)

Stability Issues

Choose time step ∆t based on uniform grid mesh size

∆x , ∆y as

∆t max p,qp , µ 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

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

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

(29)

Shock-Bubble Interaction

(30)

Shock-Bubble Interaction

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

(31)

Shock-Bubble Interaction

(32)

Shock-Bubble Interaction

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

(33)

Shock-Bubble Interaction

(34)

Shock-Bubble Interaction

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

(35)

Shock-Bubble Interaction

(36)

Shock-Bubble Interaction

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

(37)

Shock-Bubble Interaction

(38)

Shock-Bubble Interaction

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

(39)

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

(40)

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

(41)

Underwater Explosions

Numerical schlieren images for density

(42)

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

(43)

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

(44)

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

(45)

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

(46)

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

(47)

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 ξ

(48)

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 ix i ξ j

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

(49)

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

(50)

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)

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

(51)

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

(52)

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

(53)

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

(54)

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

(55)

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

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

(56)

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

(57)

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

(58)

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

(59)

Generalized Riemann Problem

Exact generalized Riemann solution: basic structure

τ

n Q n

(60)

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

(61)

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

(62)

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

(63)

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:

(64)

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

(65)

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

(66)

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

(67)

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

(68)

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

(69)

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

(70)

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

(71)

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

(72)

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

(73)

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

(74)

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

(75)

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

(76)

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

(77)

Underwater Explosions

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

(78)

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

(79)

Underwater Explosions

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

(80)

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

(81)

Underwater Explosions (Cont.)

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

−1

−0.5 0 0.5

time=0ms

(82)

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

(83)

Underwater Explosions (Cont.)

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

−1

−0.5 0 0.5

time=0.4ms

(84)

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

(85)

Underwater Explosions (Cont.)

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

−1

−0.5 0 0.5

time=1.2ms

(86)

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

(87)

Shock-Bubble (Helium)

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

(88)

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

(89)

Shock-Bubble (Helium)

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

(90)

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

(91)

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

(92)

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

(93)

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

(94)

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

(95)

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

(96)

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

(97)

Shock-Bubble (Refrigerant)

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

(98)

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

(99)

Shock-Bubble (Refrigerant)

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

(100)

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

(101)

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

(102)

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

(103)

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

(104)

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

(105)

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

(106)

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

(107)

Underwater Explosions

Numerical schlieren images h 0 = 0.6 , 100 3 grid

(108)

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

(109)

Underwater Explosions

Numerical schlieren images h 0 = 0.6 , 100 3 grid

(110)

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

(111)

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

(112)

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

(113)

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

(114)

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

(115)

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

(116)

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

(117)

3 D Shock-Bubble (Helium)

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

(118)

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

(119)

3 D Shock-Bubble (Helium)

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

(120)

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

(121)

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

(122)

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

(123)

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

(124)

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

參考文獻

相關文件

Now given the volume fraction for the interface cell C i , we seek a reconstruction that mimics the sub-grid structure of the jump between 0 and 1 in the volume fraction

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

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

Compute interface normal, using method such as Gradient or least squares of Youngs or Puckett Determine interface location by iterative bisection..

Ireland, Kenneth F./Rosen, Michael I.: A Classical Introduction to Modern Number Theory, Volume 84 of Graduate Texts in Mathematics, Springer-Verlag, New York, Second Edition,

Jing Yu, NTU, Taiwan Values at Algebraic Points.. Thiery 1995) Suppose the Drinfeld module ρ is of rank 1. Drinfeld modules corresponding to algebraic points can be defined over ¯

Based on [BL], by checking the strong pseudoconvexity and the transmission conditions in a neighborhood of a fixed point at the interface, we can derive a Car- leman estimate for

Reinforcement learning is based on reward hypothesis A reward r t is a scalar feedback signal. ◦ Indicates how well agent is doing at