• 沒有找到結果。

Keh-Ming Shyue

N/A
N/A
Protected

Academic year: 2022

Share "Keh-Ming Shyue"

Copied!
170
0
0

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

全文

(1)

Towards

a unified coordinate method for

compressible multifluid problems

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

(2)

Talk Outline

Model scientific problem Mathematical formulation

Basic compressible flow model Grid movement conditions

Multifluid extension Numerical discretization

Generalized Riemann problem Godunov-type f -wave method Sample examples

Future direction

(3)

Model Scientific Problem

Asteroid impact problem

(4)

Fundamental Challenges

Mathematical model aspect

Incompressible or compressible flow modelling ? Equations of motion & constitutive laws

· Gas phase: air

· Liquid phase: ocean water

· Solid phase: asteroid, basalt crust, mantle Interface conditions ?

Mass transfer, cavitation, fracture, · · · Numerical method aspect

Multiphase, multiscale, Eulerian or Lagrangian type

method ?

(5)

Preliminary

Basic facts

Lagrangian method resolve material or slip lines sharply if no grid tangling

Generalized curvilinear grid is often superior to

Cartesian when employed in numerical methods for

complex fixed or moving geometries

(6)

Preliminary

Basic facts

Lagrangian method resolve material or slip lines sharply if no grid tangling

Generalized curvilinear grid is often superior to

Cartesian when employed in numerical methods for complex fixed or moving geometries

Previous work done by Eulerian compressible solver Falling liquid drop problem

Shock-bubble interaction

Flying Aluminum-plate problem

Falling rigid object in water tank

(7)

Preliminary

Basic facts

Lagrangian method resolve material or slip lines sharply if no grid tangling

Generalized curvilinear grid is often superior to

Cartesian when employed in numerical methods for complex fixed or moving geometries

Previous work done by Eulerian compressible solver Falling liquid drop problem

Shock-bubble interaction

Flying Aluminum-plate problem

Falling rigid object in water tank

(8)

Falling Liquid Drop Problem

Interface capturing with gravity

100 200 300 400 500 600 700 800 900 1000

air

t = 0

(9)

Falling Liquid Drop Problem

100 200 300 400 500 600 700 800 900 1000

air

t = 0.2s

(10)

Falling Liquid Drop Problem

Interface diffused badily

100 200 300 400 500 600 700 800 900

air

t = 0.4s

(11)

Falling Liquid Drop Problem

100 200 300 400 500 600 700 800 900

air

t = 0.6s

(12)

Falling Liquid Drop Problem

100 200 300 400 500 600 700 800 900

air

t = 0.8s

(13)

Falling Liquid Drop Problem

100 200 300 400 500 600 700 800 900

air

t = 1s

(14)

Shock-Bubble Interaction

Volume tracking for material interface

(15)

Shock-Bubble Interaction

(16)

Shock-Bubble Interaction

(17)

Shock-Bubble Interaction

(18)

Shock-Bubble Interaction

(19)

Shock-Bubble Interaction

(20)

Shock-Bubble Interaction

(21)

Shock-Bubble Interaction

(22)

Shock-Bubble Interaction

(23)

Shock-Bubble Interaction

Small moving irregular cells: stability & accuracy

(24)

Flying Aluminum-plate problem

Vacuum-Al interface tracking

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

(25)

Flying Aluminum-plate problem

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

(26)

Flying Aluminum-plate problem

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

(27)

Flying Aluminum-plate problem

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

(28)

Flying Aluminum-plate problem

Small moving irregular cells: stability & accuracy

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

(29)

Falling Rigid Object in Water Tank

Moving boundary tracking & interface capturing

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

Density Pressure Volume fraction

air

water

t = 0ms

(30)

Falling Rigid Object in Water Tank

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

Density Pressure Volume fraction

air

water

t = 1ms

(31)

Falling Rigid Object in Water Tank

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

Density Pressure Volume fraction

air

water

t = 2ms

(32)

Falling Rigid Object in Water Tank

Small moving irregular cells: stability & accuracy

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

−2 0 2

−3

−2

−1 0 1 2 3

Density Pressure Volume fraction

air

water

t = 3ms

(33)

Euler Eqs. in Generalized Coord.

With gravity effect included, for example, 2 D compressible Euler eqs. in Cartesian coordinates take

∂q

∂t + ∂f (q)

∂x + ∂g(q)

∂y = ψ(q)

where

q =

 ρ ρu ρv E

, f (q) =

ρu ρu 2 + p

ρuv Eu + pu

, g(q) =

ρv ρuv ρv 2 + p Ev + pv

, ψ =

 0 0 ρg ρgv

ρ : density , (u, v) : vector of particle velocity

2 2

(34)

Euler in General. Coord. (Cont.)

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 ξ

J = x ξ y η − x η y ξ : grid Jacobian

(35)

Euler in General. Coord. (Cont.)

With these notations, Euler eqs. in generalized coord. are

∂ ˜ q

∂τ + ∂ ˜ f

∂ξ + ∂˜ g

∂η = ˜ ψ

where

˜

q = J

 ρ ρu ρv E

, ˜ f = J

ρU

ρuU + ξ x p ρvU + ξ y p EU + pU − ξ t p

, ˜ g = J

ρV

ρuV + η x p ρvV + η y p E V + pV − η t p

, ˜ ψ = J

 0 0 ρg ρgv

with contravariant velocities U & V defined by

(36)

Euler in General. Coord. (Cont.)

Model system in quasi-linear form

∂ ˜ q

∂τ + A ∂ ˜ q

∂ξ + B ∂ ˜ q

∂η = ˜ ψ

A = ∂ ˜ f

∂ ˜ q =

ξ t ξ x ξ y 0

ξ x p ρ − uU ξ x u(1 − p E ) + U ξ y u − ξ x vp E ξ x p E ξ y p ρ − vU ξ x v − ξ y up E ξ y v(1 − p E ) + U ξ y p E (p ρ − H)U ξ x H − uU p E ξ y H − vU p E U + p E U

B = ∂˜ g

∂ ˜ q =

η t η x η y 0

η x p ρ − uV η x u(1 − p E ) + V η y u − η x vp E η x p E η y p ρ − vV η x v − η y up E η y v(1 − p E ) + V η y p E (p ρ − H)V η x H − uV p E η y H − vV p E V + p E V

with H = (E + p)/ρ, U = U − ξ t = ξ x u + ξ y v, V = V − η t = η x u + η y v

(37)

Euler in General. Coord. (Cont.)

Eigen-structure of matrix A is

Λ A = diag 

U − c q

ξ x 2 + ξ y 2 , U, U, U + c q

ξ x 2 + ξ y 2 

R A =

1 1 0 1

u − α 1 c u α 2 u + α 1 c v − α 2 c v −α 1 v + α 2 c H − U 1 c H − c 2 /p E −U 2 H + U 1 c

L A =

(p ρ + cU 1 )/2c 2 −(α 1 c + up E )/2c 2 −(α 2 c + vp E )/2c 2 p E /2c 2 1 − p ρ /c 2 up E /c 2 vp E /c 2 −p E /c 2

U 2 α 2 −α 1 0

(p ρ − cU 1 )/2c 21 c − up E )/2c 22 c − vp E )/2c 2 p E /2c 2

(38)

Euler in General. Coord. (Cont.)

Eigen-structure of matrix B is

Λ B = diag 

V − c q

η x 2 + η y 2 , V, V, V + c q

η x 2 + η y 2 

R B =

1 1 0 1

u − β 1 c u β 2 u + β 1 c v − β 2 c v −β 1 v + β 2 c H − V 1 c H − c 2 /p E −V 2 H + V 1 c

L B =

(p ρ + cV 1 )/2c 2 −(β 1 c + up E )/2c 2 −(β 2 c + vp E )/2c 2 p E /2c 2 1 − p ρ /c 2 up E /c 2 vp E /c 2 −p E /c 2

V 2 β 2 −β 1 0

(p ρ − cV 1 )/2c 21 c − up E )/2c 22 c − vp E )/2c 2 p E /2c 2

with (β 1 , β 2 ) = (η x , η y )/ q

η x 2 + η y 2 , V 1 = β 1 u + β 2 v, V 2 = −β 2 u + β 1 v

(39)

Grid Movement Conditions

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

Lagrangian case: (x τ , y τ ) = (u, v)

(40)

Grid Movement (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 ξ )

(41)

Grid Movement (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

(42)

Grid Movement (Cont.)

To ensure h ∈ [0, 1] , transformed variable ˜h = κ(h) is used,

e.g. , Hui et al. employed κ = ln (εh|~u|) , ε normalized constant, yielding

A˜ ˜ h ξ + ˜ B˜ h η + ˜ C = 0

Grid-angle preserving case

A = ˜ q

x 2 η + y η 2 (x ξ sin θ − y ξ cos θ) , B = ˜ q

x 2 ξ + y ξ 2 (y η cos θ − x η sin θ) C = ˜ q

x 2 ξ + y ξ 2 [y η (cos θ) η − x η (sin θ) η ] − q

x 2 η + y η 2 [y ξ (cos θ) ξ − x ξ (sin θ) ξ ]

Mesh-area preserving case

A = y ˜ η cos θ − x η sin θ, B = x ˜ ξ sin θ − y ξ cos θ

C = y ˜ η (cos θ) ξ − x η (sin θ) ξ + x ξ (sin θ) η − y ξ (cos θ) η

where ~u = (u, v) = |~u|(cos θ, sin θ)

(43)

Grid Movement: Remarks

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 ?

(44)

Grid Movement: 2 Free Degrees

2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions

1. Grid-angle preserving

2. Specialized grid-material line matching (see next)

(45)

Grid Movement: 2 Free Degrees

2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions

1. Grid-angle preserving

2. Specialized grid-material line matching (see next)

Good results are shown for steady-state problems

(46)

Grid Movement: 2 Free Degrees

2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions

1. Grid-angle preserving

2. Specialized grid-material line matching (see next)

Good results are shown for steady-state problems

Little results for time-dependent problems with rapid

transient solution structures

(47)

Grid Movement: 2 Free Degrees

2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions

1. Grid-angle preserving

2. Specialized grid-material line matching (see next) Good results are shown for steady-state problems Little results for time-dependent problems with rapid transient solution structures

Other 2 -parameter case: (x τ , y τ ) = (hu, k v)

Novel imposed conditions for h ∈ [0, 1] & k ∈ [0, 1] ?

(48)

Grid Movement: 2 Free Degrees

2 -parameter case of Hui et al. (2005): (x τ , y τ ) = (U g , V g ) Imposed conditions

1. Grid-angle preserving

2. Specialized grid-material line matching (see next) Good results are shown for steady-state problems Little results for time-dependent problems with rapid transient solution structures

Other 2 -parameter case: (x τ , y τ ) = (hu, k v)

Novel imposed conditions for h ∈ [0, 1] & k ∈ [0, 1] ? Roadmap of current work:

(x τ , y τ ) = h 0 (u, v) → (x τ , y τ ) = h(u, v) → · · ·

(49)

Novel Conditions for h & k

Mesh-area preserving case

∂J

∂τ = ∂

∂τ (x ξ y η − x η y ξ )

= x ξτ y η + x ξ y ητ − x ητ y ξ − x η y ξτ

= · · ·

= (A 1 h ξ + B 1 h η + C 1 h) + (A 2 k ξ + B 2 k η + C 2 k) = 0, yielding, for example,

A 1 h ξ + B 1 h η + C 1 h = 0 A 2 k ξ + B 2 k η + C 2 k = 0 with

A = uy , B = uy , C = u y − u y

(50)

Single-Fluid Model

With (x τ , y τ ) = h 0 (u, v) , our model system for single-phase flow reads

∂τ

 Jρ Jρu Jρv JE

x ξ y ξ x η y η

+ ∂

∂ξ

ρU

ρuU + y η p ρvU − x η p

EU + (y η u − x η v)p

−h 0 u

−h 0 v 0 0

+ ∂

∂η

ρV

ρuV − y ξ p ρvV + x ξ p

EV + (x ξ v − y ξ u)p 0

0

−h 0 u

−h 0 v

= ˜ ψ

where U = (1 − h 0 )(y η u − x η v) & V = (1 − h 0 )(x ξ v − y ξ u)

(51)

Single-Fluid Model: Remarks

Hyperbolicity (under thermodyn. stability cond.)

In Cartesian coordinates, model is hyperbolic

(52)

Single-Fluid Model: Remarks

Hyperbolicity (under thermodyn. stability cond.) In Cartesian coordinates, model is hyperbolic

In generalized-moving coord., model is hyperbolic

when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1

(53)

Single-Fluid Model: Remarks

Hyperbolicity (under thermodyn. stability cond.) In Cartesian coordinates, model is hyperbolic

In generalized-moving coord., model is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1 Canonical form

In Cartesian coordinates

∂q

∂t + ∂f (q)

∂x + ∂g(q)

∂y = ψ(q) In generalized coordinates

∂q ∂f (q, Ξ) ∂g(q, Ξ)

(54)

Single-Fluid Model: Remarks

Hyperbolicity (under thermodyn. stability cond.) In Cartesian coordinates, model is hyperbolic

In generalized-moving coord., model is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1 Canonical form

In Cartesian coordinates

∂q

∂t + ∂f (q)

∂x + ∂g(q)

∂y = ψ(q)

In generalized coordinates : spatially varying fluxes

∂q

∂τ + ∂ f (q, Ξ)

∂ξ + ∂ g(q, Ξ)

∂η = ψ(q), Ξ : grid metrics

(55)

Three Space Dimensions

Euler equations for inviscid compressible flow

∂t

 ρ ρu ρv ρw ρE

+ ∂

∂x

ρu ρu 2 + p

ρuv ρuw ρEu + pu

+ ∂

∂y

ρv ρuv ρv 2 + p

ρvw ρEv + pv

+ ∂

∂z

ρw ρuw ρvw ρw 2 + p ρEw + pw

= ψ

E = e + (u 2 + v 2 + w 2 )/2 , e(ρ, p) : internal energy

ψ : source terms (geometrical, gravitational, & so on)

(56)

Three Space Dimensions (Cont.)

Introduce transformation (t, x, y, z) → (τ, ξ, η, ζ) via

 dt dx dy dz

=

1 0 0 0

U A 1 B 1 C 1 V A 2 B 2 C 2 W A 3 B 3 C 3

 dτ dξ dη dζ

where

Q = (U, V, W ) ~ : grid velocity Q = 0 ~ Eulerian case

Q = (u, v, w) ~ Lagrangian case

A i , B i , C i : geometric variables, i = 1, 2, 3

(57)

Three Space Dimensions (Cont.)

Inverse transformation (τ, ξ, η, ζ) → (t, x, y, z) reads

 dτ dξ dη dζ

= 1 J

J 0 0 0

J 01 J 11 J 21 J 31 J 02 J 12 J 22 J 32 J 03 J 13 J 23 J 33

 dt dx dy dz

, J =

A 1 B 1 C 1 A 2 B 2 C 2 A 3 B 3 C 3

where

J 11 = B 2 C 3 − B 3 C 2 , J 21 = C 1 B 3 − B 1 C 3 , J 31 = B 1 C 2 − C 1 B 2

J 12 = C 2 A 3 − A 2 C 3 , J 22 = A 1 C 3 − C 1 A 3 , J 32 = C 1 A 2 − A 1 C 2

J 13 = A 2 B 3 − B 2 A 3 , J 23 = B 1 A 3 − A 1 B 3 , J 33 = A 1 B 2 − B 1 A 2

(58)

Three Space Dimensions (Cont.)

Euler equations in generalized curvilinear coordinates

∂τ

 ρJ ρJu ρJv ρJw ρJE

+ ∂

∂ξ

ρU

ρuU + pJ 11 ρvU + pJ 21 ρwU + pJ 31

ρEU + pX

+ ∂

∂η

ρV

ρuV + pJ 12 ρvV + pJ 22 ρwV + pJ 32

ρEV + pY

+ ∂

∂ζ

ρW

ρuW + pJ 13 ρvW + pJ 23 ρwW + pJ 33

ρEW + pZ

= ψ

where

U = (u − U )J 11 + (v − V )J 21 + (w − W )J 31 , X = uJ 11 + vJ 21 + wJ 31

V = (u − U )J 12 + (v − V )J 22 + (w − W )J 32 , Y = uJ 12 + vJ 22 + wJ 32

W = (u − U )J 13 + (v − V )J 23 + (w − W )J 33 , Z = uJ 13 + vJ 23 + wJ 33

(59)

Three Space Dimensions (Cont.)

Geometrical conservation laws

∂τ

 A 1 A 2 A 3 B 1 B 2 B 3 C 1 C 2 C 3

+ ∂

∂ξ

−U

−V

−W 0 0 0 0 0 0

+ ∂

∂η

 0 0 0

−U

−V

−W 0 0 0

+ ∂

∂ζ

 0 0 0 0 0 0

−U

−V

−W

= 0

(60)

Grid Movement Conditions

Pseudo-Lagrangian like

(U, V, W ) = h 0 (u, v, w), h 0 ∈ (0, 1) Mesh-volume preserving: ∂J/∂t = 0

Grid-angle preserving

Other novel approach

(61)

Three Space Dimensions (Cont.)

In summary, Euler equations in generalized coord. takes

∂q

∂t + ∂f (q, Ξ)

∂ξ + ∂g(q, Ξ)

∂η + ∂h(q, Ξ)

∂ζ = ψ where

q = (ρJ, ρJu, ρJv, ρJw, ρJE, A i , B i , C i )

f (q, Ξ) = (ρU, ρuU + pJ 11 , ρvU + pJ 21 , ρwU + pJ 31 , ρEU + pX , · · · ) g(q, Ξ) = (ρV, ρuV + pJ 12 , ρvV + pJ 22 , ρwV + pJ 32 , ρEV + pY, · · · )

h(q, Ξ) = (ρW, ρuW + pJ 13 , ρvW + pJ 23 , ρwW + pJ 33 , ρEW + pZ, · · · )

with Ξ : grid metrics & equation of state p = p(ρ, e)

(62)

Hyperbolic Balance Laws

Canonical form

∂t q (~x, t) +

N d

X

j=1

∂x j f j (q, ~x) = ψ (q, ~x)

q ∈ lR m : vector of m state quantities f j ∈ lR m : flux vector, j = 1, 2, · · · , N d ψ ∈ lR m : source term

~x = (x 1 , x 2 , · · · , x N d ) : spatial vector, t : time Hyperbolicity

Assume real e-values & complete e-vectors for

Jacobian matrix P N j=1 d α j (∂f j /∂q), α j ∈ lR

(63)

Hyperbolic Balance Laws (Cont.)

Introduce coordinate transformation (~x, t) 7→ (~ ξ, τ ) via

ξ = (ξ ~ 1 , ξ 2 , · · · , ξ N d ) , ξ j = ξ j (~x, t), τ = t,

Generalized coordinate form

∂τ q ˜ ~ξ, τ + X N d

j=1

∂ξ j f ˜ j 

˜ q, ~ ξ 

= ˜ ψ 

˜ q, ~ ξ 

,

where

˜

q = Jq, f ˜ j = J q ∂ξ j

∂t +

N d

X

k=1

f k ∂ξ j

∂x k

!

, ψ = Jψ ˜

(64)

Extension to Multifluid

Assume homogeneous ( 1 -pressure & 1 -velocity) flow;

i.e. , across interfaces p ι = p & ~u ι = ~u , ∀ fluid phase ι

gas

gas gas

gas

gas gas

gas

gas gas

gas

liquid

(65)

Extension to Multifluid

Assume homogeneous ( 1 -pressure & 1 -velocity) flow;

i.e. , across interfaces p ι = p & ~u ι = ~u , ∀ fluid phase ι Mathematical model: Fluid-mixture type

Use basic conservation (or balance) laws for single

& multicomponent fluid mixtures

Introduce additional transport equations for

problem-dependent material quantities near

numerically diffused interfaces, yielding direct

computation of pressure from EOS

(66)

Extension to Multifluid

Assume homogeneous ( 1 -pressure & 1 -velocity) flow;

i.e. , across interfaces p ι = p & ~u ι = ~u , ∀ fluid phase ι Mathematical model: Fluid-mixture type

Use basic conservation (or balance) laws for single

& multicomponent fluid mixtures

Introduce additional transport equations for problem-dependent material quantities near numerically diffused interfaces, yielding direct computation of pressure from EOS

Sample examples

Barotropic 2 -phase flow

Hybrid barotropic & non-barotropic 2 -phase flow

(67)

Barotropic 2-Phase Flow

Equations of state

Fluid component 1 & 2 : Tait EOS

p(ρ) = (p + B ι )

 ρ ρ

 γ ι

− B ι , ι = 1, 2

(68)

Barotropic 2-Phase Flow

Equations of state

Fluid component 1 & 2 : Tait EOS

p(ρ) = (p + B ι )

 ρ ρ

 γ ι

− B ι , ι = 1, 2

Mixture pressure law (Shyue, JCP 2004)

p(ρ, ρe) =

 

 

 

 

(p + B ι )

 ρ ρ

 γ ι

− B ι if α = 0 or 1 (γ − 1)



ρe + ρB ρ 0



− γB if α ∈ (0, 1)

Here α denotes volume fraction of one chosen fluid component

(69)

Barotropic 2-Phase Flow

Equations of state

Fluid component 1 & 2 : Tait EOS

p(ρ) = (p + B ι )

 ρ ρ

 γ ι

− B ι , ι = 1, 2

Mixture pressure law (Shyue, JCP 2004)

p(ρ, ρe) =

 

 

 

 

(p + B ι )

 ρ ρ

 γ ι

− B ι if α = 0 or 1 (γ − 1)



ρe + ρB ρ 0



− γB if α ∈ (0, 1)

variant form of p(ρ, S) = A(S) (p 0 + B)  ρ ρ

 γ

− B

(70)

Barotropic 2-Phase Flow (Cont.)

Transport equations for material quantities γ , B , & ρ 0

γ -based equations

∂τ

 1 γ − 1



+ U ∂

∂ξ

 1 γ − 1



+ V ∂

∂η

 1 γ − 1



= 0

∂τ

 γB γ − 1



+ U ∂

∂ξ

 γB γ − 1



+ V ∂

∂η

 γB γ − 1



= 0

∂τ

 J B

ρ 0 ρ



+ ∂

∂ξ

 J B

ρ 0 ρU



+ ∂

∂η

 J B

ρ 0 ρV



= 0

(71)

Barotropic 2-Phase Flow (Cont.)

Transport equations for material quantities γ , B , & ρ 0

γ -based equations

∂τ

 1 γ − 1



+ U ∂

∂ξ

 1 γ − 1



+ V ∂

∂η

 1 γ − 1



= 0

∂τ

 γB γ − 1



+ U ∂

∂ξ

 γB γ − 1



+ V ∂

∂η

 γB γ − 1



= 0

∂τ

 J B

ρ 0 ρ



+ ∂

∂ξ

 J B

ρ 0 ρU



+ ∂

∂η

 J B

ρ 0 ρV



= 0

Above equations are derived from energy equation

& make use of homogeneous equilibrium flow

assumption together with mass conservation law

(72)

Barotropic 2-Phase Flow (Cont.)

Transport equations for material quantities γ , B , & ρ 0

γ -based equations

∂τ

 1 γ − 1



+ U ∂

∂ξ

 1 γ − 1



+ V ∂

∂η

 1 γ − 1



= 0

∂τ

 γB γ − 1



+ U ∂

∂ξ

 γB γ − 1



+ V ∂

∂η

 γB γ − 1



= 0

∂τ

 J B

ρ 0 ρ



+ ∂

∂ξ

 J B

ρ 0 ρU



+ ∂

∂η

 J B

ρ 0 ρV



= 0

If we ignore JBρ/ρ 0 term, they are essentially

equations proposed by Saurel & Abgrall (SISC

1999), but are written in generalized coord.

(73)

Barotropic 2-Phase Flow (Cont.)

Transport equations for material quantities γ , B , & ρ 0

γ -based equations

∂τ

 1 γ − 1



+ U ∂

∂ξ

 1 γ − 1



+ V ∂

∂η

 1 γ − 1



= 0

∂τ

 γB γ − 1



+ U ∂

∂ξ

 γB γ − 1



+ V ∂

∂η

 γB γ − 1



= 0

∂τ

 J B

ρ 0 ρ



+ ∂

∂ξ

 J B

ρ 0 ρU



+ ∂

∂η

 J B

ρ 0 ρV



= 0

α -based equations

∂ α

∂τ + U ∂α

∂ξ + V ∂α

∂η = 0, with z =

2

X

ι=1

α ι z ι , z = 1

γ − 1 & γB

γ − 1

(74)

Barotropic 2-Phase Flow (Cont.)

Transport equations for material quantities γ , B , & ρ 0

γ -based equations

∂τ

 1 γ − 1



+ U ∂

∂ξ

 1 γ − 1



+ V ∂

∂η

 1 γ − 1



= 0

∂τ

 γB γ − 1



+ U ∂

∂ξ

 γB γ − 1



+ V ∂

∂η

 γB γ − 1



= 0

∂τ

 J B

ρ 0 ρ



+ ∂

∂ξ

 J B

ρ 0 ρU



+ ∂

∂η

 J B

ρ 0 ρV



= 0

α -based equations (Allaire et al. , JCP 2002)

∂ α

∂τ + U ∂α

∂ξ + V ∂α

∂η = 0 with z =

2

X

ι=1

α ι z ι , z = 1

γ − 1 & γB γ − 1

∂τ (Jρ 1 α) + ∂

∂ξ (Jρ 1 αU ) + ∂

∂η (Jρ 1 αV ) = 0 with z = B

ρ ρ

(75)

Barotropic 2-Phase Flow (Cont.)

Transport equations for material quantities γ , B , & ρ 0

γ -based equations

∂τ

 1 γ − 1



+ U ∂

∂ξ

 1 γ − 1



+ V ∂

∂η

 1 γ − 1



= 0

∂τ

 γB γ − 1



+ U ∂

∂ξ

 γB γ − 1



+ V ∂

∂η

 γB γ − 1



= 0

∂τ

 J B

ρ 0 ρ



+ ∂

∂ξ

 J B

ρ 0 ρU



+ ∂

∂η

 J B

ρ 0 ρV



= 0

α -based equations (Kapila et al. , Phys. Fluid 2001)

∂ α

∂τ + U ∂α

∂ξ + V ∂α

∂η = α 1 α 2 ρ 1 c 2 1 − ρ 2 c 2 2 P 2

k=1 α k ρ k c 2 k

!

∇ · ~u

(76)

Barotropic & Non-Barotropic Flow

Equations of state

Fluid component 1 : Tait EOS

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

 γ

− B

Fluid component 2 : Noble-Abel EOS

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



ρe (0 ≤ b ≤ 1/ρ)

(77)

Barotropic & Non-Barotropic Flow

Equations of state

Fluid component 1 : Tait EOS

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

 γ

− B

Fluid component 2 : Noble-Abel EOS

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



ρe (0 ≤ b ≤ 1/ρ)

Mixture pressure law (Shyue, Shock Waves 2006)

 

  (p 0 + B)  ρ ρ 0

 γ

− B if α = 1 (fluid 1)

(78)

Baro. & Non-Baro. Flow (Cont.)

Equations of state

Fluid component 1 : Tait EOS

p(V ) = A(S 0 ) (p 0 + B)  V 0 V

 γ

− B, V = 1/ρ

Fluid component 2 : Noble-Abel EOS

p(V, S) = A(S)p 0  V 0 − b V − b

 γ

Mixture pressure law

p(V, S) = A(S) (p 0 + B)  V 0 − b V − b

 γ

− B

(79)

Baro. & Non-Baro. Flow (Cont.)

Equations of state

Fluid component 1 : Tait EOS

p(V ) = A(S 0 ) (p 0 + B)  V 0 V

 γ

− B, V = 1/ρ

Fluid component 2 : Noble-Abel EOS

p(V, S) = A(S)p 0  V 0 − b V − b

 γ

Mixture pressure law

p(V, S) = A(S) (p 0 + B)  V 0 − b V − b

 γ

− B

(80)

Baro. & Non-Baro. Flow (Cont.)

Transport equations for material quantities γ , b , & B α -based equations

∂ α

∂τ + U ∂α

∂ξ + V ∂α

∂η = 0

∂τ (Jρ 1 α) + ∂

∂ξ (Jρ 1 αU ) + ∂

∂η (Jρ 1 αV ) = 0 with z = P 2

ι=1 α ι z ι , z = γ−1 1 , γ−1 , & γ−bρ γ−1 B

(81)

Baro. & Non-Baro. Flow (Cont.)

Transport equations for material quantities γ , b , & B α -based equations

∂ α

∂τ + U ∂α

∂ξ + V ∂α

∂η = 0

∂τ (Jρ 1 α) + ∂

∂ξ (Jρ 1 αU ) + ∂

∂η (Jρ 1 αV ) = 0 with z = P 2

ι=1 α ι z ι , z = γ−1 1 , γ−1 , & γ−bρ γ−1 B

Note: 1 − bρ

γ − 1 p + γ − bρ

γ − 1 B = ρe =

2

X

ι=1

α ι ρ ι e ι

2

X  1 − b ι ρ ι γ ι − b ι ρ ι 

(82)

Baro. & Non-Baro. Flow (Cont.)

Transport equations for material quantities γ , b , & B α -based equations

∂ α

∂τ + U ∂α

∂ξ + V ∂α

∂η = 0

∂τ (Jρ 1 α) + ∂

∂ξ (Jρ 1 αU ) + ∂

∂η (Jρ 1 αV ) = 0 with z = P 2

ι=1 α ι z ι , z = γ−1 1 , γ−1 , & γ−bρ γ−1 B

Note: 1 − bρ

γ − 1 p + γ − bρ

γ − 1 B = ρe =

2

X

ι=1

α ι ρ ι e ι

=

2

X

ι=1

α ι  1 − b ι ρ ι

γ ι − 1 p ι + γ ι − b ι ρ ι γ ι − 1 B ι



參考文獻

相關文件

Have shown results in 1 , 2 & 3 D to demonstrate feasibility of method for inviscid compressible flow problems. Department of Applied Mathematics, Ta-Tung University, April 23,

Our model system is written in quasi-conservative form with spatially varying fluxes in generalized coordinates Our grid system is a time-varying grid. Extension of the model to

Lagrangian method resolve material or slip lines sharply if no grid tangling.. Generalized curvilinear grid is often

Lagrangian method can resolve material or slip lines sharply if there is not too much grid tangling.. Generalized curvilinear grid is often

GMRES: Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems..

A diamagnetic material placed in an external magnetic field B ext develops a magnetic dipole moment directed opposite B ext. If the field is nonuniform, the diamagnetic material

A diamagnetic material placed in an external magnetic field B ext develops a magnetic dipole moment directed opposite B ext.. If the field is nonuniform, the diamagnetic material

2 Center for Theoretical Sciences and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan..