Shock in air & R22 bubble interaction

56  Download (0)

Full text

(1)

Eulerian

interface-sharpening methods for

compressible flow

Keh-Ming Shyue

Department of Mathematics National Taiwan University

(2)

Objective

Derive consistent interface-sharpening model & method for compressible (single & multiphase) flow problems

PDE-based approach

tq + ∇ · f (q) + B∇q = 1 µDq

Algebraic-based approach (with F. Xiao, Tokyo Tech.) 1. Sharp cell edge solution reconstruction

2. Solution update: MUSCL or semi-discretize scheme

(3)

Shock in air & R22 bubble interaction

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

(4)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(5)

Shock in air-R22 bubble (cont.)

With anti-diffusion WENO 5

(6)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(7)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(8)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(9)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(10)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(11)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(12)

Shock air-R22 bubble (cont.)

With anti-diffusion WENO 5

(13)

Shock air-R22 bubble (cont.)

With anti-diffusion With THINC

(14)

Shock air-R22 bubble (cont.)

Volume tracking (Shyue 2006) 2nd order

(15)

PDE-based interface sharpening

Incompressible 2-phase flow: PDE-based interface sharpening for volume-fraction transport

tα + u · ∇α = 1

µDα, µ ∈ lR ≫ 1

Artificial compression: Harten CPAM 1977, Olsson &

Kreiss JCP 2005

Dα := ∇ · [(D(∆x)∇α · ~n − α (1 − α)) ~n]

Anti-diffusion: So, Hu & Adams JCP 2011 Dα := −∇ · (D(u)∇α)

(16)

Model interface-only problem

Interface-only problem for compressible 1-phase Euler equations with constant pressure p, velocity u, & jump in density ρ across interfaces

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

t (ρu) + ∇ · (ρu ⊗ u + pIN) = 0 (Momentum)

tE + ∇ · (Eu + pu) = 0 (Energy) yielding basic transport equations for interfaces as

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

2 (∂tρ + u · ∇ρ) + ∂t (ρe) + u · ∇ (ρe) = 0

(17)

Interface-only anti-diffusion model

Shyue (2011) proposed anti-diffusion model for density

tρ + u · ∇ρ = 1

µDρ, Dρ := −∇ · (D∇ρ)

To ensure velocity equilibrium for momentum

u(∂tρ + u · ∇ρ) = 1

µDρu, Dρu := u Dρ

& velocity-pressure equilibrium, for total energy

u · u

2 (∂tρ + u · ∇ρ) + ∂t (ρe) + u · ∇ (ρe) = 1

µDE, DE := hu · u

2 + ∂ρ(ρe)i Dρ

Mie-Grüneisen EOS p(ρ, e) = p(ρ) + Γ(ρ)ρ [e − e(ρ)]

(18)

1 -phase anti-diffusion model

To deal with shock waves, anti-diffusion model for compressible 1-phase flow

tρ + ∇ · (ρu) = 1 µDρ

t (ρu) + ∇ · (ρu ⊗ u + pIN) = 1

µDρu

tE + ∇ · (E~u + pu) = 1

µDE Dρ := −HI∇ · (D∇ρ) , Dρu := u Dρ, DE := hu · u

2 + ∂ρ(ρe)i Dρ HI : Interface indicator =

( 1 if near interface, 0 otherwise

(19)

Anti-diffusion method

Anti-diffusion model in compact form

tq + ∇ · f = 1 µDq with q, f, & Dq defined (not shown)

Fractional step method:

1. Solve homogeneous equation without source terms

tq + ∇ · f = 0

2. Iterate model equation with source terms

τq = Dq

to sharp layer; τ = µt (pseudo time)

(20)

Numerical interface-only problem

Consistency of numerical solution for interface-only problem Step 1: assume consistent approximation model equation without anti-diffusion,

smeared ρ & equilibrium u = un, p = pn

(may be difficult with highly nonlinear MG EOS)

Step 2: assume anti-diffusion for ρ is done consistently &

stably, yielding update ρ to ρn+1 & for momentum,

(ρu)n+1 := (ρu) + Dρu = (ρu) + uDρ

= (ρu) + u ρn+1 ρ = ρn+1u =⇒ un+1 = ρn+1u

ρn+1 = u

(21)

Numerical interface-only (cont.)

For total energy, En+1 := E + DE, i.e.,

(ρK + ρe)n+1 := (ρK + ρe) + (K + ∂ρ(ρe)) Dρ, K = u · u 2

Splitting kinetic & internal energy with MG EOS

(ρK)n+1 := (ρK) + K ρn+1 ρ = ρn+1K =⇒ Kn+1 = K

=⇒ un+1 = u (not true if u has transverse jump)

 p − p

Γ + ρe

n+1

:=  p − p

Γ + ρe



+ ∂ρ  p − p

Γ + ρe



ρn+1 ρ

=⇒ pn+1 = p (for linearized MG EOS only)

(22)

Numerical shock wave problems

Weak solution for problems with shock & rarefaction waves Interface indicator HI takes value zero away from interfacs, yielding standard compressible Euler equations in

conservation form

Step 1, use state-of-the-art shock capturing method for entropy-satisfying weak solutions

(23)

Implementation issues

Methods used here are very elementary, i.e.,

Step 1: Clawpack (or Wenoclaw/Sharpclaw) for

homogeneous equation over CFL-constrained ∆t Step 2: explicit 1-step method for anti-diffusion

Numerical regularization to ρ & Dρ

Diffusion coefficient D = D(u) (local in space & time) Time step (forward Euler) ∆τ ≤ min 

∆t, min

N

i=1 x2i 2N Dmax



Stopping criterion: Run 1 − 2 iteration currently HI is chosen based on checking jumps in ρ & p

(24)

Sod Riemann problem

High-resolution result with anti-diffusion

Ideal gas: p(ρ, e) = (γ − 1)ρe

0 0.5 1

0 0.2 0.4 0.6 0.8 1

0 0.5 1

0 0.2 0.4 0.6 0.8 1

0 0.5 1

0 0.2 0.4 0.6 0.8

Density Velocity 1 Pressure

x x

x

(25)

Aluminum impact problem

rest Al moving Al IC:

(ρ, u, p)L = 

4 × 103 kg/m3, 0 m/s, 7.93 × 109 Pa (ρ, u, p)R = 

2785 kg/m3, −2 × 103 m/s, 0 Pa

Mie-Grüneisen EOS with

Γ(ρ) = Γ0(1 − η)α, p(ρ) = ρ0c20η

(1 − sη)2 , e(ρ) = η 0

(p0 + p(ρ)) ρ0(kg/m3) c0 (m/s) s Γ0 α p0 e0

2785 5328 1.338 2.0 1 0 0

(26)

Aluminum impact (cont.)

High-resolution result with anti-diffusion at time t = 50µs

0 0.5 1

2.5 3 3.5 4 4.5

sharpen no sharpen exact

0 0.5 1

−2

−1.5

−1

−0.5 0

0 0.5 1

−5 0 5 10 15 20 25

Density (Mg/m3) Velocity (km/s) 30Pressure (GPa)

x x

x

(27)

1 -fluid multiphase anti-diffusion

Compressible multiphase: One-fluid anti-diffusion model

tρ + ∇ · (ρu) = 1 µDρ

t (ρu) + ∇ · (ρu ⊗ u + pIN) = 1

µDρu

tE + ∇ · (E~u + pu) = 1 µDE

tφj + u · ∇φj = ρdρφj∇ · u + 1

µDφj

Mixture pressure modeled by Mie-Grüneisen EOS:

φ1 = 1

Γ(ρ), φ2 = p(ρ)

Γ(ρ) , φ3 = ρe(ρ)

µ → ∞ reduces to fluid-mixture model (Shyue JCP 2001)

(28)

5 -equation 2-phase flow model

Unsteady, inviscid, compressible homogeneous 2-phase flow governed by 5-equation model

t1ρ1) + ∇ · (α1ρ1u) = 0 (Phasic 1 continuity)

t2ρ2) + ∇ · (α2ρ2u) = 0 (Phasic 2 continuity)

t (ρu) + ∇ · (ρu ⊗ u + pIN) = 0 (Momentum)

tE + ∇ · (Eu + pu) = 0 (Energy)

tα1 + u · ∇α1 = 0 (Volume fraction) Phasic pressure pk follows Mie-Grüneisen EOS

pkk, ek) = p,kk) + Γkkk 

ek − e,kk)

k = 1, 2

(29)

5 -equation model (cont.)

Isobaric closure leads to mixture pressure , i.e., Substitute p = p1 = p2 in ρe =

2

X

k=1

ρkek, yielding

p = ρe −

2

X

k=1

αkρke,kk) +

2

X

k=1

αk p,kk) Γkk)

! 2 X

k=1

αk Γkk)

Model is hyperbolic with mixture acoustic impedance Allaire et al. (JCP 2002)

ρc2 =

2

X

k=1

αkρkc2k, ck : phasic sound speed

(30)

5 -equation model (cont.)

In cavitated regions, p < psat, cutoff model Non-conservative energy correction

E := Esat =

2

X

k=1

αkkek)sat + ρK

cutoff phasic internal energy is

kek)sat = psat p,ksat)

Γksat) + ρsate,ksat)

α-based energy-preserving correction (Shyue 2012)

α1 := (α1)sat = ρe − 2e2)sat 1e1)sat 2e2)sat

(31)

5 -equation anti-diffusion model

Proposed 5-equation anti-diffusion model

t1ρ1) + ∇ · (α1ρ1u) = 1

µDα

1ρ1

t2ρ2) + ∇ · (α2ρ2u) = 1

µDα2ρ2

t (ρu) + ∇ · (ρu ⊗ u + pIN) = 1

µDρu

tE + ∇ · (Eu + pu) = 1

µDE

tα1 + u · ∇α1 = 1

µDα1

Exist two ways to set Dz, z = α1ρ1, . . . , α1, in literature

(32)

5 -equation anti-diffusion (cont.)

α-ρ based (Shyue 2011)

Dα1 := −∇ · (D∇α1), Dαkρk := HI∇ · (D∇αkρk), k = 1, 2 Dρ :=

2

X

k=1

Dαkρk, Dρu := uDρ,

DE := KDρ +

2

X

k=1

αkρkkek)Dαkρk +

2

X

k=1

ρkekDαk

α-based only (So, Hu, & Adams JCP 2012)

Dα1 := −∇ · (D∇α1), Dαkρk := ρkDαk, k = 1, 2, Dα2 := −Dα1 Dρ :=

2

X

k=1

Dαkρk, Dρu := uDρ, DE := KDρ +

2

X

k=1

ρkekDαk

Use D := ∇ · [(D∇α ·~n − α (1 − α )) ~n] (Shyue 2012)

(33)

5 -equation interface-compression

Shukla, Pantano & Freund (JCP 2010)

t 1ρ1) + ∇ · (α1ρ1u) = 0

t 2ρ2) + ∇ · (α2ρ2u) = 0

t (ρu) + ∇ · (ρu ⊗ u + pIN) = 0

tE + ∇ · (Eu + pu) = 0

tα1 + u · ∇α1 = 1

µn · ∇(D∇α1 · n α1 (1 − α1))

tρ + ∇ · (ρu) = 1

µHI1)n · (∇ (D∇ρ · n) − (1 − 2α1) ∇ρ)

Nonliner compression & linear diffusion for interface sharpening & stability

Method proposed there is unstable numerically

(34)

PDE-based interface sharpening

PDE-based interface-sharpening model in compact form

tq + ∇ · f + B∇q = 1 µDq with q, f, B, & Dq defined (not shown)

Fractional step method is used

1. Solve homogenous equation without source terms

tq + ∇ · f + B∇q = 0 2. Iterate model equation with source terms

τq = Dq

to sharp layer; τ = µt (pseudo time)

(35)

Numerical interface-only problem

Consistency of numerical solution for interface-only problem Step 1: assume consistent approximation model equation without anti-diffusion,

smeared (α1ρ1), (α2ρ2), α1 & equilibrium u = un, p = pn (can be done even with highly nonlinear MG EOS)

Step 2: update (α1ρ1), (α2ρ2), α1 to (α1ρ1)n+1, (α2ρ2)n+1, α1n+1 via anti-diffusion, & for momentum,

(ρu)n+1 := (ρu) + Dρu = (ρu) + uDρ =⇒ un+1 = u

(36)

Numerical interface-only (cont.)

For total energy, En+1 := E + DE, i.e., α-based

(ρK + ρe)n+1 := (ρK + ρe) + KDρ +

2

X

k=1

kek) Dαkρk

Splitting kinetic & internal energy with MG EOS

(ρK)n+1 := (ρK) + K ρn+1 ρ = ρn+1K =⇒ Kn+1 = K

2

X

k=1

kρkek)n+1 :=

2

X

k=1

kρkek) +

2

X

k=1

kek) αn+1k αk

2

X

k=1



αk  p − p,k

Γk + ρke,k

n+1

:=

2

X

k=1

 p − p,k

Γk + ρke,k



αn+1k

=⇒ pn+1 = p (for general MG EOS)

(37)

Piston-induced liquid drops depression

Liquid & vapor governed by stiffened gas EOS Piston velocity u = −100m/s

vapor

liquid liquid

← piston

(38)

Liquid drops depression

Vapor mass fraction Mixture pressure

t=2ms

(39)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=4ms

(40)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=6ms

(41)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=8ms

(42)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=10ms

(43)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=12ms

(44)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=14ms

(45)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=16ms

(46)

Liquid drops depression (cont.)

Vapor mass fraction Mixture pressure

t=18ms

(47)

Underwater explosions: cylindrical wall

High pressure gaseous explosive in water (cylindrical case)

water

explosive

(48)

UNDEX: cylindrical wall (cont.)

t = 30µs

shock

reflected rarefaction reflected shock

t = 45µs

cavitation

interface

(49)

UNDEX: cylindrical wall (cont.)

t = 60µs

cavitation

t = 120µs

cavitation

(50)

UNDEX: cylindrical wall (cont.)

0 20 40 60 80 100 120

0 2000 4000 6000 8000

With anti−difusion Without−anti−diffusion

0 20 40 60 80 100 120

−10 0 10 20 30 40 50 60

With anti−difusion Without−anti−diffusion

Pressure(atm) Velocity(m/s)

Time (µs)

cavi cutoff cavi collapse & reload

(51)

Algebraic-based interface sharpening

Incompressible flow: Algebraic-based interface sharpening CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes): Ubbink & Issa JCP 1999

THINC (Tangent of Hyperbola for INterface Capturing):

Xiao, Honma & Kono Int. J. Numer. Meth. Fluids 2005 Improved THINC

Compressible flow: THINC-type interface sharpening 1. Cell edge solution reconstruction: THINC

2. Solution update: MUSCL or semi-discretize scheme

(52)

THINC interface sharpening

THINC reconstruction assumes

αi(x) = αmin + αmax αmin

2



1 + γ tanh



βx − xi−1/2

∆x x¯i



αmin = min (αi−1, αi+1) , αmax = max (αi−1, αi+1) , γ = signi+1 αi−1)

β measures sharpeness (given constant) & x¯i chosen by

αi = 1

∆x

Z xi+1/2 xi−1/2

αi(x) dx

Cell edges determined by

αi+1/2,L = αmin + αmax αmin

2



1 + γ tanh β + C 1 + C tanh β



αi−1/2,R = αmin + αmax αmin

(1 + γC) (C not shown)

(53)

Homogeneous-equilibrium THINC

With αi+1/2,L & αi−1/2,R obtained by THINC reconstruction Homogeneous-equilibrium reconstruction (analgously to α-based 5-equation anti-diffusion model)

1ρ1)i+1/2,L = ρ1,i α1,i+1/2,L − α1,i (α2ρ2)i+1/2,L = ρ2,i α2,i+1/2,L − α2,i

ρi+1/2,L = (α1ρ1)i+1/2,L + (α2ρ2)i+1/2,L (ρu)i+1/2,L = uiρi+1/2,L

Ei+1/2,L = Kiρi+1/2,L +

2

X

k=1

kek) i αk,i+1/2,L − αk,i

zi−1/2,R can be made in a similar manner

(54)

6 -equation anti-diffusion model

Anti-diffusion to 6-equation model (Pelanti & Shyue 2012)

t 1ρ1) + ∇ · (α1ρ1~u) = 1

µDα1ρ1

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

µDα2ρ2

t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1 + α2p2) = 1

µDρu

t 1E1) + ∇ · (α1E1~u + α1p1~u) + w (q, ∇q) = −νpI (p1 p2) + 1

µDα1E1

t 2E2) + ∇ · (α2E2~u + α2p2~u) − w (q, ∇q) = νpI (p1 p2) + 1

µDα2E2

tα1 + ~u · ∇α1 = ν (p1 p2) + 1

µDα1

w = −~u ((Y2p1 + Y1p2) ∇α1 + α1Y2p1 α2Y1p2)

(55)

Future perspective

Shape-preserving interface shapening on mapped grid Interface shapening to flow with physical sources

(usefulness ?)

tq + ∇ · f (q) + B∇q = ψ(q) + 1 µDq

Extension to low Mach weakly compressible flow Hybrid interface-sharpening & WENO towards high order for compressible turbulent flow

(56)

Thank you

Figure

Updating...

References

Related subjects :