• 沒有找到結果。

Illustrative Example 1

N/A
N/A
Protected

Academic year: 2022

Share "Illustrative Example 1"

Copied!
61
0
0

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

全文

(1)

Numerical Methods for

Compressible Multicomponent Flow with

Moving Interfaces and Boundaries

Keh-Ming Shyue

Department of Mathematics National Taiwan University

(2)

Illustrative Example 1

Flying Aluminum-plate problem Vacuum-Al boundary

0 1 2 3

−2

−1 0 1 2

Density

Al target Al flyer

vacuum

0 1 2 3

−2

−1 0 1 2

Pressure

0 1 2 3

−2

−1 0 1 2

Volume fraction

t = 0µs

(3)

Illustrative Example 1

Flying Aluminum-plate problem Vacuum-Al boundary

0 1 2 3

−2

−1 0 1 2

Density

Al target Al flyer

vacuum

0 1 2 3

−2

−1 0 1 2

Pressure

0 1 2 3

−2

−1 0 1 2

Volume fraction

t = 0.5µs

(4)

Illustrative Example 1

Flying Aluminum-plate problem Vacuum-Al boundary

0 1 2 3

−2

−1 0 1 2

Density

Al target Al flyer

vacuum

0 1 2 3

−2

−1 0 1 2

Pressure

0 1 2 3

−2

−1 0 1 2

Volume fraction

t = 1µs

(5)

Illustrative Example 1

Flying Aluminum-plate problem Vacuum-Al boundary

0 1 2 3

−2

−1 0 1 2

Density

Al target Al flyer

vacuum

0 1 2 3

−2

−1 0 1 2

Pressure

0 1 2 3

−2

−1 0 1 2

Volume fraction

t = 2µs

(6)

Illustrative Example 1

Flying Aluminum-plate problem Vacuum-Al boundary

0 1 2 3

−2

−1 0 1 2

Density

Al target Al flyer

vacuum

0 1 2 3

−2

−1 0 1 2

Pressure

0 1 2 3

−2

−1 0 1 2

Volume fraction

t = 4µs

(7)

Illustrative Example 2

Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary

−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

(8)

Illustrative Example 2

Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary

−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

(9)

Illustrative Example 2

Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary

−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

(10)

Illustrative Example 2

Flying rigid-cylinder in partially air-filled water tank Air-water interface & moving rigid boundary

−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

(11)

Overview

Mathematical model

Fluid-mixture type equations of motion for homogeneous two-phase flow

Mie-Grüneisen EOS for real materials Numerical techniques

Finite volume method based on wave propagation Surface tracking for moving boundaries

Volume tracking for moving interfaces Numerical results

Future work

(12)

Two Phase Flow Problem

Ignore physical effects such as viscosity, surface tension, mass diffusion, and so on

Each fluid component k, k = 1, 2, satisties Eulerian conservation laws

ρt + ∇ · (ρ~u) = 0

(ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (ρE)t + ∇ · (ρE~u + p~u) = 0 Mie-Grüneisen equation of state

p(ρ, e) = pref(ρ) + ρ Γ(ρ)

he − eref(ρ) i

(13)

Mie-Grüneisen Equations of State

Typical examples are:

(pref, eref) lies along an isentrope

(14)

Mie-Grüneisen Equations of State

Typical examples are:

(pref, eref) lies along an isentrope

1. Jones-Wilkins-Lee EOS for gaseous explosives Γ(V ) = γ − 1

eref(V ) = e0 +

A V0

R1

exp  −R1V V0



+ B V0

R2

exp  −R2V V0



pref(V ) = p0 + A exp

 −R1V V0



+ B exp −R2V V0



Here V = 1/ρ

(15)

Mie-Grüneisen Equations of State

Typical examples are:

(pref, eref) lies along an isentrope

2. Cochran-Chan EOS for solid explosives Γ(V ) = Γ0 = γ − 1

eref(V ) = e0 +

−A V0

1 − E1

"

 V V0

1−E1

− 1

#

+ B V0

1 − E2

"

 V V0

1−E2

− 1

#

pref(V ) = p0 + A

 V V0

−E1

− B  V V0

−E2

(16)

Mie-Grüneisen Equations of State

(pref, eref) lies along a Hugoniot locus

Assume linear shock speed us & particle velocity up us = c0 + s up

We may derive the relations Γ(V ) = Γ0

 V V0

α

pref(V ) = p0 +

c02

(V0 − V ) [V0 − s(V0 − V )]2 eref(V ) = e0 +

1 2

hpref(V ) + p0

i(V0 − V )

(17)

Material Quantities for Model EOS

JWL EOS ρ0(kg/m3) A(GPa) B(GPa) R1 R2 Γ

TNT1 1630 371.2 3.23 4.15 0.95 0.30

TNT2 1630 548.4 9.375 4.94 1.21 1.28

Water 1004 1582 −4.67 8.94 1.45 1.17

CC EOS ρ0(kg/m3) A(GPa) B(GPa) E1 E2 Γ

TNT 1840 12.87 13.42 4.1 3.1 0.93

Copper 8900 145.67 147.75 2.99 1.99 2

Shock EOS ρ0(kg/m3) c0(m/s) s Γ0 α

Aluminum 2785 5328 1.338 2.0 1

Copper 8924 3910 1.51 1.96 1

(18)

Two-Phase Flow Model

Model derivation 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ψ∇χki & hχkψti = h(χkψ)ti−hψ(χk)ti

(19)

Two-Phase Flow Model (Cont.)

Take product of each conservation (or balance) law with χk

& perform averaging process. In case of mass conservation equation, for example, we have

kρkit + ∇· < χkρk~uk >= hρkk)t + ρk~uk · ∇χki Since χk is governed by

k)t + ~u0 · ∇χk = 0, ~u0: interface velocity yielding averaged equation for mass

kρkit + ∇· < χkρk~uk >= hρk (~uk − ~u0) · ∇χki Analogously, we may derive averaged equation for momentum, energy, & entropy (not shown here)

(20)

Two-Phase Flow Model (Cont.)

In summary, averaged model system, we have, are

kρkit + ∇· < χkρk~uk >= hρk (~uk − ~u0) · ∇χki

kρk~ukit + ∇· < χkρk~uk ⊗ ~uk > +∇ hχkpki = hpk∇χki + hρk~uk (~uk − ~u0) · ∇χki

kρkEkit + ∇· < χkρkEk~uk + χkpk~uk >= hpk~uk · ∇χki + hρkE (~uk − ~u0) · ∇χki

kit + h~uk · ∇χki = h(~uk − ~u0) · ∇χki Existence of interfacial source terms

Mathematical as well as numerical modelling these terms are essential & difficult for multiphase flow problems

(21)

Homogeneous Two-Phase Flow Model

Assume homogeneous flow (i.e., across interfaces:

pk = p & ~uk = ~u, k = 0, 1, 2)

Introduce volume fraction αk = Vk/V (α1 + α2 = 1)

By dropping all interfacial terms, we may obtain a simplified model as

kρk)t + ∇ · (αkρk~u) = 0

kρk~u)t + ∇ · (αkρk~u ⊗ ~u) + ∇ (αkp) = p∇αkkρkEk)t + ∇ · (αkρkEk~u + αkp~u) = p~u · ∇αkk)t + ~u · ∇αk = 0

for k = 1, 2

(22)

Homogeneous Two-Phase Flow Model

Assume homogeneous flow (i.e., across interfaces:

pk = p & ~uk = ~u, k = 0, 1, 2)

Introduce volume fraction αk = Vk/V (α1 + α2 = 1) Alternatively, a simplified model as

kρk)t + ∇ · (αkρk~u) = 0

(ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (ρE)t + ∇ · (ρE~u + p~u) = 0 (αk)t + ~u · ∇αk = 0

Here ρ = P2

k=1 αkρk, ρE = P2

k=1 αkρkEk

(23)

Homogeneous Flow Model (Cont.)

Mixture equation of state: p = p(α2, α1ρ1, α2ρ2, ρe) Isobaric closure: p1 = p2 = p

For a class of EOS, explicit formula for p is available (examples are given next)

For some complex EOS, from (α2, ρ1, ρ2, ρe) in model equations we recover p by solving

p11, ρ1e1) = p22, ρ2e2) &

2

X

k=1

αkρkek = ρe This homogeneous two-phase model was called a five-equation model by Allaire, Clerc, & Kokh (JCP

2002) or a volume-fraction model by Shyue (JCP 1998)

(24)

Homogeneous Flow Model (Cont.)

Polytropic ideal gas: pk = (γk − 1)ρkek ρe =

2

X

k=1

αkρkek =

2

X

k=1

αk p

γk − 1 =⇒ p = ρe

 2 X

k=1

αk γk − 1 Van der Waals gas: pk = (1−γkb−1kρk )(ρkek + akρ2k) − akρ2k

ρe =

2

X

k=1

αkρkek =

2

X

k=1

αk  1 − bkρk γk − 1



(p + akρ2k) − akρ2k



=⇒

p =

"

ρe −

2

X

k=1

αk  1 − bkρk

γk − 1 − 1



akρ2k

#  2 X

k=1

αk  1 − bkρk γk − 1



(25)

Homogeneous Flow Model (Cont.)

Two-molecular vibrating gas: pk = ρkRkT (ek), T satisfies e = RT

γ − 1 + RTvib

exp Tvib/T − 1 As before, we now have

ρe =

2

X

k=1

αkρkek =

2

X

k=1

αk

 ρkRkTk γk − 1



+ ρkRkTvib,k exp 

Tvib,k/Tk

 − 1

=

2

X

k=1

αk

 p

γk − 1



+ pvib,k exp 

pvib,k/p

 − 1

 (Nonlin. eq.)

(26)

Homogeneous Flow Model (Cont.)

It can be shown entropies, Sk, k = 1, 2, satisfy

 ∂p1

∂S1



ρ1

DS1

Dt −  ∂p2

∂S2



ρ2

DS2

Dt = ρ1c21 − ρ2c22 ∇ · ~u Murrone & Guillard (JCP 2005) propsed a reduced two-phase flow model in which

2)t + ~u · ∇α2 = α1α2

ρ1c21 − ρ2c22 P2

k=1 αkρkc2k

!

and now phase entropies satisfy DSk

Dt = ∂Sk

∂t + ~u · ∇Sk = 0, for k = 1, 2

(27)

Homogeneous Flow Model (Cont.)

Model system is hyperbolic under suitable

thermodynamic stability condition (see below)

In the model, when α2 = 0 (or = 1), ρ2 (or ρ1) can not be recovered from α2 & α2ρ2 (or α1 & α1ρ1).

It is not absolutely clear in the model how to compute nonlinear term ρι, ι > 1 from αk & αkρk

This formulation of model equation would not work when one fluid component is adiabatic, but the other fluid component is not

Surely, there are other set of model systems proposed in the literature that are robust for homogeneous flow

(28)

Thermodynamic Stability

Fundamental derivative of gas dynamics

G = −V 2

(∂2p/∂V 2)S

(∂p/∂V )S , S : specific entropy

(29)

Thermodynamic Stability

Fundamental derivative of gas dynamics

G = −V 2

(∂2p/∂V 2)S

(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,

(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0

(30)

Thermodynamic Stability

Fundamental derivative of gas dynamics

G = −V 2

(∂2p/∂V 2)S

(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,

(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0 (∂2p/∂V 2)S > 0 means convex EOS

(31)

Thermodynamic Stability

Fundamental derivative of gas dynamics

G = −V 2

(∂2p/∂V 2)S

(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,

(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0 (∂2p/∂V 2)S > 0 means convex EOS

(∂p/∂V )S < 0 means real speed of sound, for c2 =  ∂p

∂ρ



S

= −V 2  ∂p

∂V



S

> 0

(32)

Finite Volume Wave Propagation Method

Finite volume method, QnS gives approximate value of cell average of solution q over cell S at time tn

QnS ≈ 1 M(S)

Z

S

q(X, tn) dV

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

C E

D F

G H

(33)

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

Qn+1S := Qn+1S −M (Wp ∩ S)

M(S) Rp, Rp being jump from RP

Wp Wp

(34)

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

↓ ↓

Wpq Wpq

(35)

Wave Propagation Method (cont.)

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

a) b)

c) d)

αpr

p/2

αpr

p/2

λp∆ t

λp∆ t

(36)

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)

(37)

Interface Reconstruction Scheme

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

1. Compute interface normal

Gradient method of Parker and Youngs Least squares method of Puckett

2. Determine interface location by iterative bisection

0 0

0 0

0 0

0 0

0 0

0

0 0

0 0

0 0

0 0

0

1 0.29 0.68 0.09

0.51 0.51

0.09 0.68 0.29

↓ interface

interface

Data set Parker & Youngs Puckett

(38)

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)

(39)

Surface Moving Procedure

Solve Riemann problem at tracked interfaces & use resulting waves to find new location of interface at the next time step

o

o

o o

o

o

↑ old front

old front

new front

(40)

Interface Conditions

For tracked segments representing rigid (solid wall)

boundary (stationary or moving), appropriate boundary states are assigned for fictitious subcells in each time step

For tracked segments representing material interfaces, jump conditions across interfaces are satisfied only in an approximate manner, & would not be imposed

explicitly in each time step

(41)

Stability Issues

Choose time step ∆t based on uniform grid mesh size

∆x, ∆y as

∆t maxp,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 interpolation operator (such as, h-box approach of Berger et al. ) locally for cell averages in irregular cells

(42)

References

(JCP 1998) An efficient shock-capturing algorithm for compressible multicomponent problems

(JCP 1999, 2001) A fluid-mixture type algorithm for

compressible multicomponent flow with van der Waals (Mie-Grüneisen) equation of state

(JCP 2004) A fluid-mixture type algorithm for barotropic two-fluid flow Problems

(JCP 2006) A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions

(Shock Waves 2006) A volume-fraction based algorithm for hybrid barotropic & non-barotropic two-fluid flow

problems

(43)

Flying Projectile Problem

0 1 2 3

−1

−0.5 0 0.5 1

Stationary Density

0 1 2 3

−1

−0.5 0 0.5 1

Pressure

0 1 2 3

−1

−0.5 0 0.5 1

Entropy

0 1 2 3

−1

−0.5 0 0.5 1

Moving Density

0 1 2 3

−1

−0.5 0 0.5 1

Pressure

0 1 2 3

−1

−0.5 0 0.5 1

Entropy

(44)

Cylinder lift-off Problem

Moving speed of cylinder is governed by Newton’s law Pressure contours are shown with a 1000 × 200 grid

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15 0.2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15 0.2

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.05 0.1 0.15 0.2

t = 0

t = 0.1641s

t = 0.30085s

(45)

Cylinder lift-off Problem

A convergence study of center of cylinder & relative mass loss for at final stopping time t = 0.30085s

Mesh size Center of cylinder Relative mass loss 250 × 50 (0.618181, 0.134456) −0.257528

500 × 100 (0.620266, 0.136807) −0.131474 1000 × 200 (0.623075, 0.138929) −0.066984

Results are comparable with numerical appeared in literature

(46)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(47)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(48)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(49)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(50)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(51)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(52)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(53)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(54)

Shock in Air & R22 Bubble Interaction

Leftward-going Mach 1.22 shock wave in air over heavier R22 bubble

Numerical schlieren images for density

(55)

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

(56)

Shock-Bubble Interaction (cont.)

Quantitative assessment of prominent flow velocities:

Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf 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

Vs (VR, VT) Incident (refracted, transmitted) shock speed t ∈ [0, 250]µs (t ∈ [0, 202]µs, t ∈ [202, 250]µs ) Vui (Vuf) Initial (final) upstream bubble wall speed t ∈ [0, 400]µs (t ∈ [400, 1000]µs)

Vdi (Vdf) Initial (final) downstream bubble wall speed t ∈ [200, 400]µs (t ∈ [400, 1000]µs)

(57)

Shock wave in molybdenum over MORB

Numerical schlieren images for density

(58)

Shock-MORB Interaction (cont.)

Numerical schlieren images for pressure

(59)

Shock-MORB Interaction (cont.)

Approximate locations of interfaces

molybdenum

MORB

molybdenum

MORB

time = 50µs time = 100µs

(60)

Future Work

Extension to low Mach number flow

Remove sound-speed stiffness by preconditioning techniques or pressure-based method

Extension to include more physics towards real applications

Such as capillary, diffusion, or elastic-plastic effect Extension to 3D volume tracking method

Surface reconstruction

Finite volume method with moving interfaces Stability in presence of small cut cells

Extension to unified coordinates of W.-H. Hui

(61)

Thank You

參考文獻

相關文件

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

He proposed a fixed point algorithm and a gradient projection method with constant step size based on the dual formulation of total variation.. These two algorithms soon became

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

Idea: condition the neural network on all previous words and tie the weights at each time step. Assumption: temporal

Model checking: residuals can be obtained by 1-step ahead forecast errors at each time point, (observations minus the fitted value) Then if the model is adequate, residuals should