compressible multiphase flow

Full text

(1)

compressible multiphase flow

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

(2)

Objective

Discuss a class of volume of fluid (vs. level set, MAC,

particles) methods for interface problems with application to compressible multiphase flow

1. Adaptive moving grid approach (this lecture) Cartesian grid embedded volume tracking

Useful also for turbulent flow (Glimm & Abarzhi’s talk)

Moving mapped grid interface capturing

Simple method if mesh handling is robust

2. Eulerian interface sharpening approach (next lecture) Artificial interface compression method

Anti-diffusion method

(3)

Incompressible 2-phase flow

One model for free-surface (or 2-phase) flow can be described by incompressible Navier-Stokes 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 α-weighted average (α ∈ [0, 1]) as

ρ = αρ1 + (1 − α) ρ2, ǫ = αǫ1 + (1 − α) ǫ2, where

τ = ǫ ∇~u + ∇~uT

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

 ∇α

|∇α|



(4)

Incompressible 2-phase flow

Accurate advection scheme for volume fraction α is of fundamental importance. Typical methods employed are

Simple line interface calculation (SLIC) Donor acceptor scheme

Flux corrected transport method Youngs (PLIC) method

Level set method

CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes)

THINC (Tangent of Hyperbola for INterface Capturing) or its variant

(5)

Compressible flow in porous medium

One model for gas flow in variable porous medium can be described by compressible Euler equations:

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

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

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

Embeded discontinuous porosity α ∈ [0, 1]

Well-balanced & entropy solution due to non-conservative term

(6)

Baer-Nunziato 2-phase model

Known 2-phase mixture model for deflagration-to-detonation transition (DDT) in reactive granular materials

1ρ1)t + ∇ · (α1ρ1~u1) = 0

1ρ1~u1)t + ∇ · (α1ρ1~u1 ⊗ ~u1) + ∇(α1p1) = p0∇α1 + λ(~u2 − ~u1) (α1ρ1E1)t + ∇ · (α1ρ1E1~u1 + α1p1~u1) = p02)t + λ~u0 · (~u2 − ~u1) (α2ρ2)t + ∇ · (α2ρ2~u2) = 0

2ρ2~u2)t + ∇ · (α2ρ2~u2 ⊗ ~u2) + ∇(α2p2) = p0∇α2 − λ(~u2 − ~u1)

2ρ2E2) t + ∇ · (α2ρ2E2~u2 + α2p2~u2) = −p02)t − λ~u0 · (~u2 − ~u1) (α2)t + ~u0 · ∇α2 = µ (p2 − p1)

αk = Vk/V : volume fraction (α1 + α2 = 1), ρk: density,

~uk: velocity, pk: pressure, Ek = ek + ~u2k/2: specific total energy, ek: specific internal energy, k = 1, 2

(7)

Baer-Nunziato 2-phase model

Quantities p0 & ~u0 denote interfacial pressure & velocity Baer & Nunziato (1986)

p0 = p2, ~u0 = ~u1 Saurel & Abgrall (1999)

p0 = P2

k=1 αkpk, ~u0 = P2

k=1 αkρk~uk

P2

k=1 αkρk

Quantities λ & µ (> 0) are velocity & pressure relaxation parameters that give rates to equilibrium for velocitiy &

pressure, respectively

(8)

Reduced 2-phase model

For some problems, it is simple to use reduced 2-phase model of Kapila, Menikoff, Bdzil, Son, & Stewart (2001)

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 = ρ2c22 − ρ1c21 P2

k=1 ρkc2kk ∇ · ~u

Model is zero-order approximation of Baer-Nunziato model with stiff mechanical relaxation

Positivity & sharpness of α is important for accurate solution

(9)

Cartesian embedded boundary method

Towards volume tracking for sharp interface resolution, we consider conservation laws

tq + ∇ · ~f (q) = 0

with embedded boundary (i.e., interface with rigid BCs)

q ∈ lRm: vector of m conserved quantities

f = (f~ 1, f2, . . . , fd) ∈ lRm×d: flux matrix; fj ∈ lRm

Γ: embedded boundary

Γ →

(10)

Finite-volume discretization

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

QnS ≈ 1 M(S)

Z

S

q(x, tn) dV = 1

κS (∆x)d Z

S

q(x, tn) dV M(S): measure (area in 2D

or volume in 3D) of cell S κS = M(S)/ (∆x)d: volume fraction

hkS = M(∂Sk)/∆x(d−1): area fraction

Assume uniform discretiza- tion ∆x = ∆x , j = 1, . . . , d

A B

C D

(11)

Finite-volume discretization

Divergence theorem over each control volume leads to finite volume approximation for

∇ · ~f (q)

S ≈ 1

M(S) Z

S

∇ · ~f (q) dV = 1 M(S)

Z

∂S

f (q) · ~n dA~

= 1

κS∆x

X

k∈∂S

hkSf (q)~ k · ~nk = 1 κS∆x

X

k∈∂S

hkSk

≡ DS · ~f (q)

Away from boundaries,

method reduces to standard conservative finite difference discretization

τ∆x = 

∇ · ~f

S − DS · ~f = O(∆xp−1/κ ) if F˘k = O(∆p)

A B

C D

(12)

Finite-volume discretization

For conservation laws, On each control volume S, a fully discretized finite volume method can be written as

Qn+1S = QnS − ∆tDS · ~f (q) = QnS − ∆t κS∆x

X

k∈∂S

hkSk

which satisfies discrete conservation identity X

S

κSQn+1 = X

S

κSQn − ∆t

∆x

X

S

X

k∈∂S

hkSk

Method is subject to stringent CFL stability constraint

∆t = O

(κmin)1/d∆x

|λ|max



≪ 1 if κmin ≪ 1 (Small-cell problem)

(13)

Small-cell problem

Methods proposed for small-cell problem include Non-conservation approach: Glimm & coworker

Flux-redistribution approach: Chern, Colella, & Modiano Cell-merging approach: Quirk, Powell, Causon

H-box approach: Berger, Helzel, LeVeque (see below)

R( ¯Q) Q¯R

L

(14)

Flux calculation: Wave propagation

On 2D uniform grid, high-resolution in fluctutation form Qn+1ij = Qnij − ∆t

∆x1

A+1 ∆Qi−1

2,j + A1 ∆Qi+1

2,j



− ∆t

∆x2

A+2 ∆Qi,j−1

2 + A2 ∆Qi,j+1

2



− ∆t

∆x1

F˜i+1 1

2,j − ˜Fi−1 1 2,j

 − ∆t

∆x2

F˜i,j+2 1 2

− ˜Fi,j−2 1 2



A+1 ∆Q

A1 ∆Q

A2 A+1 ∆Q

A+2 A1 ∆Q

A+2 A+1 ∆Q

A2 A1 ∆Q

(15)

Flux calculation: Wave propagation

Fluctuation terms in wave propagation form are obtained from solving Riemann problems at each cell edge, yielding speeds λj,m & waves Wj,k, j = 1, 2, k = 1, . . . , mw,

A±1 ∆Q =

mw

X

k=1

λ1,k±

W1,k, A±2 ∆Q =

mw

X

k=1

λ2,k±

W2,k

A±2 A±1 ∆Q =

mw

X

p=1

λ1,p±

A±2 W1,p =

mw

X

p=1

λ1,p± Xmw

q=1

λ2,q±

W2,q

1 := F˜1 + 1 2

mw

X

k=1

λ1,k



1 − ∆t

∆x1

λ1,k



Wf1,k

To aviod oscillations near discontinuities, a wave limiter is applied leading to limited waves Wf

(16)

Wave propagation with boundary

First order upwind with normal-wave update

Solve Riemann problem in direction normal at cell edge,

& use resulting waves to update nearby cell averages Qn+1S := Qn+1S − M (Zp ∩ S)

M(S) Wp

Z↓p

A B

C

D Z↓p

A B

C D

(17)

Wave propagation with boundary

First order upwind with transverse-wave update

Use transverse portion of equation, solve Riemann

problem in transverse direction, & use resulting waves to update neighboring cell averages as usual

Stability of method is typically improved, while conservation of method is maintained

Z↓pq

A B

C

D Z↓pq

A B

C D

(18)

Wave propagation with boundary

High resolution with piecewise linear wave update

wave before propagation after propagation

a) b)

c) d)

αprp/2

αprp/2

λp∆ t

λp∆ t

(19)

Moving boundary tracking

With a prescribed moving-object velocity ~ub at current time, tracked boundaries are advanced in a Lagrangian manner over time step ∆t

Final moving boundary locations are set by where new segments intersect grid lines; interpolation is required with care to initialize cut-cell values

DC BA new object↓

old object

(20)

Embedded boundary condition

Rigid embedded boundary (stationary or moving), reflection principle is used to assign states for fictitious subcells, i.e.,

z := zE (z = ρ, p, α)

~u := ~uE − 2(~uE · ~n)~n + 2(~ub · ~n), where ~ub is velocity of object

If motion of object is governed by Newton’s second law, velocity of object ~un+1b at time step tn+1 is computed by

~un+1b = ~unb + ∆t

MbF~bn,

where Mb is mass of object & F~bn is force on object

(21)

Moving piston problem

Numerical results agree well with exact solution Piston moves with speed ~ub = (0, 300)m/s

0 0.5 1

0.5 1 1.5 2

x (m) ρ (kg/m3 )

0 0.5 1

0 0.05 0.1 0.15 0.2 0.25 0.3

x (m)

u (km/s)

0 0.5 1

0.5 1 1.5 2 2.5

x (m)

p (Bar)

piston

(22)

Cylinder lift off

Cylinder motion 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.05 0.1 0.15 0.2

t = 0

t = 0.1641s

t = 0.30085s

(23)

Cylinder lift off

A convergence study of center of cylinder & relative mass loss 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

(24)

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

Density Pressure 3 Volume fraction air

water t = 1ms

−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

air

water t = 3ms

(25)

VOF interface tracking

Given volume fractions, VOF tracking consists of:

1. Volume fraction update

Take a time step on current grid using wave propagation method to update volume fractions α governed by

tα + ~u · ∇α = 0

at next time step; ~u underlying velocity field

2. Interface reconstruction (by PLIC or MOF, see Z. Jia’s talk at 17:40-18:00, 9/15)

Standard PLIC method does:

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

(26)

Interface reconstruction: Example

Cell-averaged volume fractions (left) & reconstructed interface (right); a set of disconnected line segments

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↓

(27)

Interface reconstruction: Example

With α in previous slide & ~u = (1, 1) updated volume fraction over ∆t = 0.06 shown here ( κ1 = 5.7 × 103 &

κ2 = 1.3 × 103)

New reconstructed interface location

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 κ1 κ2

old interface

new interface

(28)

Solution update

Solution is double-valued on cells subdivided by VOF tracked interface; Again use finite volume formulation of solution QnS for approximate value of solution q

QnS ≈ 1 M(S)

Z

S

q(x, tn) dV

C E

D F

G H

(29)

Wave propagation method

Normal- & transverse-edge wave propagation

Zp Zp

↓ ↓

Zpq Zpq

(30)

Stability issues

Choose time step ∆t based on uniform mesh size ∆x as ν = |λ|max∆t

∆x ≤ 1

Use wave propagation method described above which expands range of influence of solution to maintain

stability of method even in presence of small Cartesian cut cells

Use smoothing operator (such as, h-box approach of Berger et al. ) locally for cell averages in irregular cells for accuracy

(31)

(32)

(33)

(34)

Shock in air & R22 bubble

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)

(35)

(36)

Flying Aluminum plate

Track vacuum-Al boundary Impact velocity 400m/s

Pressure contours are shown

Al target Al flyer

vacuum

time = 0.5µs time = 1µs time = 2µs time = 4µs

(37)

Moving mapped grid method

Towards sharp interface resolution, one may also consider moving mapped grid in that hyperbolic conservation laws

∂q

∂t +

Xd j=1

∂fj(q)

∂xj = 0

in transformed coordinates (x, t) 7→ (ξ, τ ) becomes conservation laws with patially varying fluxes:

∂τ (Jq) +

Xd j=1

∂ξj Jq ∂ξj

∂t +

Xd k=1

Jfk(q) ∂ξj

∂xk

!

= 0

To close system, conditions on grid motion ∂tξ & so grid metrics ∂ξj/∂xk with J = det (∂ξ/∂x)1 should be set

(38)

Computing geometric quantities

To set ∂tξ, ∂ξj/∂xk, & J, one may use Hui’s unified coordinate approach

∂xi

∂τ = hui & ∂

∂τ

∂xi

∂ξj



+ ∂

∂ξj



−∂xi

∂τ



= 0

for i, j = 1, 2, . . . , d & compute J = det (∂ξ/∂x)1 directly.

Quantity h may be kept as a constant or varied with time by geometric measure

Variant ALE approach, see next talk by B. Chen or B.

Tian’s talk at 11:50-12:10, 9/15

∂J

∂τ +

Xd

∂ξj



J ∂ξj

∂t



= 0

(39)

Reduced 2-phase model (revisit)

Physical balance laws

∂τ





α1ρ1J α2ρ2J ρuiJ

ρEJ



 +

Xd j=1

∂ξj J





α1ρ1Uj α2ρ2Uj ρuiUj + p∂x∂ξj

i

ρEUj + pUj − p∂ξ∂tj



 = 0

∂α1

∂τ +

Xd j=1

Uj ∂α1

∂ξj = ρ2c22 − ρ1c21 P2

k=1 ρkc2kk

Xd j=1

∂Uj

∂ξj Uj = ∂tξj + Pd

i=1 uixiξj: ξj’s contravariant velocity

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

(40)

Numerical approximation

To ensure free-stream (constant states) preservation, which is one of basic properties that a moving mesh method

should have, a consistent discretization of

∂J

∂τ +

Xd j=1

∂ξj



J ∂ξj

∂t



= 0

is required which may be difficult for 2D quadrilateral or 3D hexahedral mesh

In this regards, non-conservative form is often employed

∂q

∂τ + 1 J

Xd j=1

∂ξj

Xd

k=1

Jfk(q) ∂ξj

∂xk + Jq∂ξj

∂t

!

= 0

(41)

We remark that Menciger & ˘Zun (JCP 2011) developed a PLIC-VOF adaptive moving grid method in that monitor function D employed in elliptic grid generator, say in form proposed by Ceniceros & Hou (JCP 2001),

ξ · (D∇ξxj) = 0, j = 1, . . . , d

for x(ξ), make use of volume fraction α explicitly, i.e.,

∂D

∂τ = ε∇2D with D(x, 0) = 1 + βM (∇α) & D|∂Ω = 0 where M (z) = |z|R

dV / R

|z|dV

This can be used readily in moving-mesh method developed by T. Tang & many others

(42)

Woodward-Colella’s problem

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

0 0.5 1

0 5 10 15 20

ρ

Fine grid h0=0 h0=0.99

0 0.5 1

−5 0 5 10 15

0 0.5 1

0 200 400 600

u p

t = 0.032 t = 0.032

t = 0.032

0 2 4 6 8

ρ

−5 0 5 10 15

0 200 400 600

u p

t = 0.038 t = 0.038

t = 0.038

(43)

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 →

(44)

Oscillating water column

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

time

oscilw grid

0 1 2 3 4 5 6 7 8 9 10

0.8 1 1.2 1.4 1.6 1.8

pL

0.8 1 1.2 1.4 1.6 1.8

pR

mmesh mesh0

mmesh mesh0

(45)

Water-vapor cavitation

Initially, in closed shock tube, flow is homogeneous (contains α = 106 gas in bulk liquid) at standard

atmospheric condition & exists interface separating flow with opposite motion (u = 100 m/s)

Result in pressure drop & formation of cavitation zone in middle; shocks form also from both ends

Eventually, shock-cavitation collision occurs

−u u

← Inteface

(46)

Water-vapor cavitation

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1x 10−3

time

cavi grid

(47)

Underwater explosions

t=0.2ms

t=0.4ms

t=0.8ms

t=1.2ms

Moving mesh Eulerian mesh air

water

(48)

Underwater explosions

Moving mesh (coarsen by factor 5)

−1 0 1

−1

−0.5 0 0.5

−1 0 1

−1

−0.5 0 0.5

time = 0.2ms time = 0.4ms

−1

−0.5 0 0.5

−1

−0.5 0 0.5

time = 0.8ms time = 1.2ms

(49)

Accuracy test in 2D

Consider 2D compressible Euler equations with ideal gas law as governing equations

Take smooth vortex flow with initial condition

ρ =



1 − 25(γ − 1)

8γπ2 exp (1 − r2)

1/(γ−1)

, p = ργ,

u1 = 1 − 5

exp ((1 − r2)/2) (x2 − 5) , u2 = 1 + 5

exp ((1 − r2)/2) (x1 − 5) ,

& periodic boundary conditions as an example, where r = p

(x1 − 5)2 + (x2 − 5)2

(50)

Accuracy test in 2D

Grids used for this smooth vortex flow test

0 5 10

0 2 4 6 8 10

0 5 10

0 2 4 6 8 10

0 5 10

0 2 4 6 8

Grid 1 Grid 2 10 Grid 3

kEzk1,∞ = kzcomput − zexactk1,∞ denotes discrete 1- or maximum-norm error for state variable z

Results shown below are at time t = 10 on N × N mesh

(51)

Accuracy results in 2D: Grid 1

N E1(ρ) Order E1(u1) Order E1(u2) Order E1(p) Order

40 0.6673 2.3443 1.7121 0.8143

80 0.1792 1.90 0.6194 1.92 0.4378 1.97 0.2128 1.94 160 0.0451 1.99 0.1537 2.01 0.1104 1.99 0.0536 1.99 320 0.0113 2.00 0.0384 2.00 0.0276 2.00 0.0134 2.00

N E(ρ) Order E(u1) Order E(u2) Order E(p) Order

40 0.1373 0.3929 0.1810 0.1742

80 0.0377 1.87 0.1014 1.95 0.0502 1.85 0.0482 1.85 160 0.0093 2.02 0.0248 2.03 0.0123 2.03 0.0119 2.02 320 0.0022 2.07 0.0062 2.00 0.0030 2.04 0.0029 2.04

(52)

Accuracy results in 2D: Grid 2

N E1(ρ) Order E1(u1) Order E1(u2) Order E1(p) Order

40 0.9298 2.6248 2.1119 1.2104

80 0.2643 1.81 0.7258 1.85 0.5296 2.00 0.3277 1.89 160 0.0674 1.97 0.1833 1.99 0.1309 2.02 0.0845 1.96 320 0.0169 2.00 0.0458 2.00 0.0327 2.00 0.0212 1.99

N E(ρ) Order E(u1) Order E(u2) Order E(p) Order

40 0.1676 0.4112 0.2259 0.2111

80 0.0471 1.83 0.1242 1.73 0.0645 1.79 0.0586 1.85 160 0.0126 1.91 0.0333 1.90 0.0162 2.02 0.0149 1.97 320 0.0033 1.93 0.0085 1.97 0.0040 2.00 0.0038 1.98

(53)

Accuracy results in 2D: Grid 3

N E1(ρ) Order E1(u1) Order E1(u2) Order E1(p) Order

40 4.8272 4.7734 5.3367 5.4717

80 1.5740 1.62 1.5633 1.61 1.5660 1.77 1.5634 1.81 160 0.4536 1.79 0.4559 1.78 0.4537 1.79 0.4560 1.78 320 0.1215 1.90 0.1221 1.90 0.1222 1.89 0.1221 1.90

N E(ρ) Order E(u1) Order E(u2) Order E(p) Order

40 0.4481 0.4475 0.4765 0.4817

80 0.1170 1.94 0.1181 1.92 0.1196 1.99 0.1191 2.02 160 0.0434 1.43 0.0431 1.45 0.0442 1.43 0.0440 1.44 320 0.0117 1.89 0.0119 1.86 0.0119 1.89 0.0118 1.89

(54)

Future perspectives

Extension of Cartesian embedded boundary method in 3D will be based on either FronTier or Chombo

framework developed by Stony Brook, or ANAG, LBL Free-stream perserving PLIC-VOF moving mesh

method is under development

Low Mach solver for cavitation problems is ongoing

(55)

Thank you

(56)

Basic grid-metric relations

Assume existence of inverse transformation

t = τ, xj = xj(~ξ, t) for j = 1, 2, . . . , d,

To find basic grid-metric relations between different coordinates, employ elementary differential rule

∂(τ, ~ξ)

∂(t, ~x) = ∂(t, ~x)

∂(τ, ~ξ)

1

yielding in d = 3 case, for example, as

1 0 0 0

tξ1 x1ξ1 x2ξ1 x3ξ1

tξ2 x1ξ2 x2ξ2 x3ξ2

= 1 J

J 0 0 0

J01 J11 J21 J31 J02 J12 J22 J32

(57)

Basic grid-metric relations

Here

J =

∂(x1, x2, x3)

∂(ξ1, ξ2, ξ3)

= det

∂(x1, x2, x3)

∂(ξ1, ξ2, ξ3)

 , J11 =

∂(x2, x3)

∂(ξ2, ξ3)

, J21 =

∂(x1, x3)

∂(ξ3, ξ2)

, J31 =

∂(x1, x2)

∂(ξ2, ξ3) ,

J12 =

∂(x2, x3)

∂(ξ3, ξ1)

, J22 =

∂(x1, x3)

∂(ξ1, ξ3)

, J32 =

∂(x1, x2)

∂(ξ3, ξ1) ,

J13 =

∂(x2, x3)

∂(ξ1, ξ2)

, J23 =

∂(x1, x3)

∂(ξ2, ξ1)

, J33 =

∂(x1, x2)

∂(ξ1, ξ2) ,

J0j = −

Nd

X

i=1

Jijτxi, j = 1, 2, 3,

and so grid-metric relations between different coordinates

∇ξj = (∂tξj, ∇~xξj) = (∂tξj, ∂x1ξj, ∂x2ξj, ∂x3ξj) = 1

J (J0j, J1j, J2j, J3j)

(58)

Basic grid-metric relations

Note in two dimensions d = 2, we have

∂ξ1

∂t , ∂ξ1

∂x1 , ∂ξ1

∂x2



= 1 J



∂x1

∂τ

∂x2

∂ξ2 + ∂x2

∂τ

∂x1

∂ξ2 , ∂x2

∂ξ2 , −∂x1

∂ξ2



∂ξ2

∂t , ∂ξ2

∂x1 , ∂ξ2

∂x2



= 1 J

∂x1

∂τ

∂x2

∂ξ1 ∂x2

∂τ

∂x1

∂ξ1 , −∂x2

∂ξ1 , ∂x1

∂ξ1



J = ∂x1

∂ξ1

∂x2

∂ξ2 ∂x1

∂ξ2

∂x2

∂ξ1

Thus to have G = 0 fulfilled, grid-metrics should obey

∂J

∂τ +

∂ξ1



J ∂ξ1

∂t



+

∂ξ2



J ∂ξ2

∂t



= 0

∂ξ1



J ∂ξ1

∂x1



+

∂ξ2



J ∂ξ2

∂x1



=

∂ξ1

∂x2

∂ξ2



+

∂ξ2



∂x2

∂ξ1



= 0



J ∂ξ1 

+ 

J ∂ξ2 

= 

−∂x1 

+ 

∂x1 

= 0

(59)

Basic grid-metric relations

Here we have

G = q

∂J

∂τ +

Xd j=1

∂ξj



J ∂ξj

∂t



 +

Xd j=1

fj

" d X

k=1

∂ξk



J ∂ξk

∂xj

#

Note with the use of basic grid-metric relations (see next), it is known that

∂J

∂τ +

Xd j=1

∂ξj



J ∂ξj

∂t



= 0 (geometrical conservation law) Xd

k=1

∂ξk



J ∂ξk

∂xj



= 0 ∀ j = 1, 2, . . . , d (compatibility condition)

and hence G = 0

Updating...

References

Related subjects :