• 沒有找到結果。

interface-sharpening methods for

N/A
N/A
Protected

Share "interface-sharpening methods for"

Copied!
46
0
0

(1)

Application to compressible multiphase flow

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

(2)

Outline

Computing monotone sharp resolution of interfaces is of fundamental importance in many practical problems of interest

Discuss a simple Eulerian interface sharpening approach (vs. Lagrangian, interface tracking, or adaptive moving mesh) for hyperbolic problems Review two PDE-based interface sharpening techniques for solving volume-fraction linear

transport equation that arises, for example, from viscous incompressibe 2-phase flow

Extend method for computing material lines or free surfaces arising from compressible multiphase flow

(3)

Incompressible 2-phase flow: Review

Consider unsteady, incompressible, viscous, immiscible 2-phase flow with governing equations

∇ · ~u = 0 (Continuity)

t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = ∇ · τ + ρ~g + f~σ (Momentum)

tα + ~u · ∇α = 0 (Volume-fraction transport ) Material quantities in 2-phase coexistent region are often

computed by α-based weighted average as

ρ = αρ1 + (1 − α) ρ2, ǫ = αǫ1 + (1 − α) ǫ2, α ∈ [0, 1], where source terms are volume-fraction dependent

τ = ǫ ∇~u + ∇~uT

, f~σ = −σκ∇α with κ = ∇ ·

 ∇α

|∇α|



(4)

Interface sharpening techniques

Typical interface sharpening methods for computing volume fraction in incompressible 2-phase flow include:

Algebraic based approach

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: Xiao et al.

PDE based approah

Artificial compression: Harten CPAM 1977, Olsson &

Kreiss JCP 2005

Anti-diffusion: So, Hu & Adams JCP 2011

(5)

Artificial interface compression

Our first interface-sharpening model concerns nonlinear artificial compression of form proposed by Olsson & Kreiss

tα + ~u · ∇α = 1

µ∇ · [(D (∇α · ~n) − α (1 − α))~n]

where ~n = ∇α/|∇α|, D > 0, µ ≫ 1

Numerical method based on fractional step may apply 1. Advection step over a time step ∆t to solve

tα + ~u · ∇α = 0 or ∂tα + ∇ · (α~u) = 0 since by assumption ∇ · ~u = 0

2. Interface compression step towards τ-steady state

∂ α = ∇ · [(D (∇α · ~n) − α (1 − α))~n] , τ = t/µ

(6)

Square-wave pluse moving with u = 1 after 4 periodic cycle

0 0.2 0.4 0.6 0.8 1

No compress With compress Exact

α

(7)

Zalesak’s rotating disc

Contours α = 0.5 at 4 different times within 1 period in that

~u = (1/2 − y, x − 1/2)

No compression With compression

(8)

Vortex in cell

Contours α = (0.05, 0.5, 0.95) at 6 different times in 1 period

~u = − sin2 (πx) sin (2πy), sin (2πx) sin2 (πy)

cos (πt/8) No compression

With compression

(9)

Interface compression: Remarks

Consider 1D model problem with u > 0 of form



tα + u∂xα = 1

µ∂x [D (∂xα · ~n) − α (1 − α)]

α(x, 0) = α0(x) = 1/ (1 + exp (−x/D)) , x ∈ lR, t > 0 Exact solution for this problem is simply α(x, t) = α0(x − ut) When α(x, 0) is perturbed to α0(x) + δ(x), δ(x) ≪ 1, we have

τα + ∂˜ ξ α˜2/2

= ∂ξ D∂ξα˜

, α(ξ, 0) = ˜˜ α0(ξ) with ξ = x − ut, τ = t/µ, & α(ξ, τ ) = 1 − α(ξ, τ )˜ , yielding steady state solution α(ξ, τ ) = ˜α0 (ξ + ξ0) as τ → ∞ for some suitably chosen shift ξ0, see Sattinger (1976)

(10)

Interface compression: Remarks

If perturbation is zero mass, i.e., R

−∞ δ(ξ, 0)dξ = 0 we have true solution with ξ0 = 0, see Goodman (1986)

When model is solved by a conservative method, truncation errors will be of zero mass, yielding convergence of

numerical solution to exact one in time we want

In multi-D case, let Kα = D∇α · ~n − α (1 − α). We solve

τα = ∇ · [(D (∇α · ~n) − α (1 − α))~n] = Kα∇ · ~n + ~n · ∇Kα yielding τ-steady state solution as µ → ∞, when Kα = 0 &

1D profile in coordinate normal to interface

When µ finite, Kα∇ · ~n + ~n · ∇Kα 6= 0, strength & accuracy of curvature ∇ · ~n plays important role in interface resolution

(11)

Interface compression runs

Methods used here are very elementary, i.e., 1. Use Clawpack for advection in Step 1

2. Use simple forward Euler in time, second order in space for interface compression in Step 2

Diffusion coefficient D = ε min∀i ∆xi Time step ∆τ

∆τ ≤ 1 2D

Xd

i=1

∆x2i

Stopping criterion: simple 1-norm error measure

(12)

Extension to compressible flow

Shukla, Pantano & Freund (JCP 2010) proposed extension of interface-compression method for incompressible flow to compressible flow governed by reduced 2-phase model as

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

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

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

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

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

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

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

µH(α1)~n · (∇ (D∇ρ · ~n) − (1 − 2α1) ∇ρ) Mixture pressure is computed based on isobaric closure

(13)

Compressible flow: Density correction

To see how density compression term comes from, we assume ∇ρ · ~n ∼ ∇α1 · ~n & consider case when

Kα1 = D∇α1 · ~n − α1 (1 − α1) ≈ 0 =⇒ D∇α1 · ~n ≈ α1(1 − α1) yielding density diffusion normal to interface at τ-steady as

∇ (D∇ρ · ~n) · ~n ≈ ∇ (α1(1 − α1)) · ~n = (1 − 2α1) ∇α1 · ~n

∼ (1 − 2α1) ∇ρ · ~n

Define Kρ = ∇(D∇ρ · ~n) − (1 − 2α1) ∇ρ & form

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

µH(α1)~nKρ

H(α1) = tanh (α1(1 − α1)/D)2 is localized-interface function

(14)

Shuklaet al.interface compression

In each time step, Shukla’s interface-compression algorithm for compressible 2-phase flow consists of following steps:

1. Solve model equation without interface-compression terms by WENO method

2. Compute primitive variable w = (ρ1, ρ2, ρ, ~u, p, α1) from conservative variables q = (α1ρ1, α2ρ2, ρ~u, ρE, α1)

3. Iterate interface & density compression equations to τ-steady state until convergence

4. Update conserved variables at end of time step from

primitive variables in step 2 & new values of ρ , α1 from step 3

(15)

Underwater explosion (UNDEX)

Solution adpated from Shukla’s paper (JCP 2010)

(16)

Underwater explosion

Solution adpated from Shyue’s paper (JCP 2006)

(17)

Shuklaet al.algorithm: Remarks

In Shukla’s results there are noises in pressure contours for UNDEX means poor calculation of pressure near interface To understand method better, consider simple interface only problem where p & ~u are constants in domain, while ρ &

material quantities in EOS have jumps across interfaces Assume consistent approximation in step 1 for model

equation without interface-compression, yielding

smeared (α1ρ1, α2ρ2, α1) & retain (~u, p) = (~u, p) In step 3, ρ = (α1ρ1) + (α2ρ2) & α1 are compressed to ρ˜ &

˜

α1, which in step 4, total mass & momentum are set

(ρ, ρu)n+1 = (˜ρ, ˜ρ~u) =⇒ ~un+1 = ˜ρ~u/˜ρ = ~u as expected

(18)

Shuklaet al.algorithm: Remarks

In addition, for total energy, we set

(ρE)n+1 =

1

2ρ|~u|2 + ρe

n+1

= 1

2ρ|~˜ u|2 + ρe(?)e

Consider stiffened gas EOS for phasic pressure pk = (γk − 1) (ρe)k − γkBk, k = 1, 2. We then have

e ρe =

X2 k=1

αkρkek =

X2 k=1

˜

αk p + γkBk γk − 1

= p

X2 k=1

˜ αk

γk − 1 +

X2 k=1

˜

αk γkBk γk − 1

yielding equilibrium pressure pn+1 = p if

 1 n+1

=

X2 α˜k

&

 γB n+1

=

X2

˜

αk γkBk

(19)

Shuklaet al.algorithm: Remarks

Next example concerns linearized Mie-Grüneisen EOS for phasic pressure pk = (γk − 1) (ρe)k + (ρk − ρ0k) Bk

e ρe =

X2 k=1

αkρkek =

X2 k=1

˜

αkp

γk − 1 − (α˜kρk − α˜kρ0k) Bk γk − 1

= p

X2 k=1

˜ αk

γk − 1 −

X2 k=1

(α˜kρk − α˜kρ0k) Bk γk − 1 yielding equilibrium pressure pn+1 = p if

 1 γ − 1

n+1

=

X2 k=1

˜ αk

γk − 1 &

(ρ − ρ0)B γ − 1

n+1

=

X2 k=1

(α˜kρk α˜kρ0k) Bk γk − 1

(20)

Shuklaet al.algorithm: Remarks

In Shukla et al. algorithm, there is a consistent problem as X2

k=1

kρk)n+1 =

X2 k=1

˜

αkρk 6= ρ˜ = ρn+1

One way to remove this inconsistency is to include

compression terms in partial density αkρk directly, k = 1, 2,

tkρk) + ∇ · (αkρk~u) = 1

µH(αk)~n · (∇ (D∇ (αkρk) · ~n) − (1 − 2αk) ∇ (αkρk)) We then set ρn+1 = P2

k=1kρk)n+1 = P2

k=1 α˜kρ˜k Validation of this approach is required

(21)

Positivity & accuracy

In compressible multiphase flow, positivity of volume fraction, i.e., αk ≥ 0, ∀k, is important due to provision of 1. information on interface location

2. information on thermodynamic states such as ρe & p in numerical “mixture” region & so ρk from αkρk

It is known that devise of oscillation-free higher-order

method (WENO, DG, or variant) for multiphase flow is still an open problem

In this regards, interface-sharpening of some kind should be a useful tool as opposed to Eulerian higher-order

methods or other adaptive mesh methods

(22)

Anti-diffusion interface sharpening

Our second interface-sharpening model concerns

anti-diffusion proposed by So, Hu & Adams (JCP 2011)

tα + ~u · ∇α = − 1

µ∇ · (D∇α), D > 0, µ ≫ 1 Standard fractional step method may still apply

1. Advection step over a time step

tα + ~u · ∇α = 0 2. Anti-diffusion step towards sharp layer

τα = −∇·(D∇α) or ∂τα = −∇·(D∇α · ~n) ~n, τ = t/µ Numerical regularization is required such as employ MINMOD

(23)

Square-wave pluse moving with u = 1 after 4 periodic cycle

0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8 1

No sharpen Compress Anti−diffus Exact

α

time=4

(24)

Vortex in cell (revisit)

Contours α = (0.05, 0.5, 0.95) at 6 different times in 1 period No interface sharpening (second order)

With interface compression

With anti-diffusion

(25)

Deformation flow in 3D

In this test, consider velocity field

~u = 2 sin2 (πx) sin (2πy) sin (2πz), − sin (2πx) sin2 (πy) sin (2πz),

− sin (2πx) sin (2πy) sin2 (πz)

cos (πt/3)

(26)

Deformation flow in 3D

No anti-diffusion With anti-diffusion

(27)

Deformation flow in 3D

No anti-diffusion With anti-diffusion

(28)

Anti-diffusion runs

Methods used here are essentially the same as artificial interface compression runs, i.e.,

1. Use Clawpack for advection in Step 1

2. Use first order explicit method for anti-diffusion in Step 2 Diffusion coefficient D = max |~u|

Time step ∆τ

∆τ ≤ 1 2D

Xd i=1

∆x2i

Stopping criterion: some measure of interface sharpness

(29)

Anti-diffusion to compressible flow

Reduced 2-phase model with anti-diffusion (Shyue 2011)

tα1 + ~u · ∇α1 = − 1

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

µKα1

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

µHI∇ · (D∇α1ρ1) = −1

µHIKα1ρ1

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

µHI∇ · (D∇α2ρ2) = −1

µHIKα2ρ2

t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = − 1

µHI~u ∇ · (D∇ρ) = − 1

µHIKρ~u

t (ρE) + ∇ · (ρE~u + p~u) = − 1

µHI Kρ|~u|2/2 + Kρe

Isobaric closure for mixture pressure is used as usual

HI denotes interface indicator, K denotes “diffusion” term

(30)

Anti-diffusion to compressible flow

To find Kρ|~u|2/2 assuming |~u|2 is constant, we observe

1

2ρ|~u|2



= 1

2|~u|2∇ρ yielding Kρ|~u|2/2 = 1

2|~u2|∇ · (D∇ρ) To find Kρe, we need to know equation of state. Now in

stiffened gas case with pk = (γk − 1) (ρe)k − γkBk,

∇(ρe) = ∇

X2 k=1

αkρkek

!

= ∇

X2 k=1

αk p + γkBk γk − 1

!

=

X2 k=1

p + γkBk γk − 1



∇αk =

p + γ1B1

γ1 − 1 − p + γ2B2 γ2 − 1



∇α1

(31)

Anti-diffusion to compressible flow

We next consider case with linearized Mie-Grüneisen EOS pk = (γk − 1) (ρe)k + (ρk − ρ0k) Bk k = 1, 2, & proceed same procedure as before

∇(ρe) = ∇

X2 k=1

αkρkek

!

= ∇

X2 k=1

αk p − (ρk − ρ0k)Bk γk − 1

!

=

X2 k=1

p + ρ0kBk

γk − 1 ∇αk +

X2 k=1

Bk

γk − 1∇ (αkρk)

= β0∇α1 +

X2 k=1

βk∇(αkρk)

We choose Kρe = β0 ∇ · (D∇α1) + P2

1 βk ∇ · (D∇αkρk)

(32)

Anti-diffusion to compressible flow

Write anti-diffusion model in compact form

tq + ∇ · ~f + B∇q = −1

µψ(q) with q, f~, B, & ψ defined (not shown)

In each time step, proposed anti-diffusion algorithm for compressible 2-phase flow consists of following steps:

1. Solve model equation without anti-diffusion terms

tq + ∇ · ~f + B∇q = 0

2. Iterate model equation with anti-diffusion terms

τq = −ψ(q)

(33)

Circular water column in uniform flow

Density surface plot (moving speed ~u = (1, 1/10) )

(34)

Single-phase Riemann problem

0 0.2 0.4 0.6 0.8 1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

No sharpen Anti−diffu Exact

ρ

x

Sod problem

(35)

Single-phase Riemann problem

0 0.2 0.4 0.6 0.8 1

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

No sharpen Anti−diffu Exact

ρ

x

Lax problem

(36)

Oscillating water column

Initially, in closed shock tube, water column moves at u = 1 from left to right, yielding air compression at right

& air expansion at left

Subsequently, pressure difference built up across water column resulting deceleration of column of water to

right, makes a stop, & then acceleration to left; a reverse pressure difference built up across water column redirecting flow from left to right again

Eventually, water column starts to oscillate

air

air water

u →

(37)

Oscillating water column

0.6 0.8 1 1.2 1.4 1.6 1.8 2

pL

With anti−diffusion No anti−diffusion

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2

pR

(38)

2 D Riemann Problem

With initial 4-slip lines wave pattern

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1



 ρ u v p



 =



 1 0.75

−0.5 1







 2 0.75

0.5 1









1

−0.75 0.5

1









3

−0.75

−0.5 1





(39)

2 D Riemann Problem

Density on top & pressure on bottom

(40)

Shock in air & water cylinder

Solution adpated from Shukla’s paper (JCP 2010)

(41)

Shock in air & water cylinder

Schlieren images for anti-diffusion results at times t = 0.15, 0.4, 0.65 (volume fraction on right)

(42)

Anti-diffusion sharpening: Remarks

1. Local interface identification Algebraic

H(α) = tanh (α(1 − α)/D)2 (Shukla et al. ) Physical-jump (density, volume fraction, . . . )

2. Diffusion coefficient definition

Global D = max |~u| or Dj = max |uj|

Local D = maxM (C) |~u| or DjC = maxM (C) |uj| 3. Stopping criterion

Run 1 − 2 anti-diffusion iteration currently Interface-shapreness measure (?)

(43)

Effect of diffusion coefficient

Anti-diffusion results for Zalesak’s rotating disc

Global Dj Local Dj

(44)

Effect of stopping criterion

Anti-diffusion results for Zalesak’s rotating disc

Local Dj: 1 step Local Dj: 2 step

(45)

Future perspectives

Extend method to mapped grid & model with phase

transition proposed by Saurel, Petipas, Abgrall (JFM 2008)

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

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

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

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

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

KS1 ∇ · ~u + K˜Γ

S Q1 + K˜cS ρ ˙Y K¯S =

 α1

KS1 + α2 KS2

−1

, K˜S =

KS1

α1 + KS2 α2



, K˜Γ =

1

α1 + Γ2 α2



c =

 c21

+ c22 

, Kι = ριc2ι , Q1 = H(T2 − T1), ˙Y (mass trans.)

(46)

Thank you

Weak solution for problems with shock &amp; rarefaction waves Interface indicator H I takes value zero away from interfacs, yielding standard compressible Euler equations

A constant state u − is formed on the left side of the initial wave train followed by a right facing (with respect to the velocity u − ) dispersive shock having smaller

 If I buy a call option from you, I am paying you a certain amount of money in return for the right to force you to sell me a share of the stock, if I want it, at the strike price,

Interface positions at different instants: experimental (left) and numerical results computed without (Simulation 1, middle) and with (Simulation 2, right)

Eulerian interface sharpening approach (this lecture) Artificial interface compression method..

Objectives  To introduce the Learning Progression Framework LPF for English Language as a reference tool to identify students’ strengths and weaknesses, and give constructive

Learning elements of the knowledge contexts at junior secondary level in the TEKLA Curriculum Guide was enriched to give students a broad and balanced. foundation on

• Figure 26.26 at the right shows why it is safer to use a three-prong plug for..