A generalized Lagrangian scheme for hyperbolic balance laws
Application to compressible multifluid flow
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Taiwan
Outline
Model scientific problem
Previous work & current work motivation Mathematical formulation
Hyperbolic balance law in generalized coordinate Grid movement conditions
Examples to shallow water/compressible flow model Numerical discretization
Generalized Riemann problem
Godunov-type flux-based wave propagation method
Sample examples
Model Scientific Problem
Asteroid impact problem (a geophysical example)
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 solver
Discretization based on structured, unstructured, or
Previous Work
Fluid-mixture interface-capturing method (JCP 1998, 1999, 2001, 2004, Shock Waves 2006)
Shock over MORB (Mid-Ocean Ridge Basalt) liquid Falling liquid drop problem
Volume-of-fluid interface tracking method (JCP 2006) Shock-bubble interaction
Surface tracking for moving boundaries (Hyp 2006) Falling rigid object in water tank
Moving cylindrical vessel
Unified coordinate method (2006, with Hui & Hu)
Supersonic NACA0012 over heavier gas
Current Work
Motivated by well-known facts that
Lagrangian method can resolve material or slip lines sharply if there is no grid tangling
Generalized curvilinear grid is often superior to
Cartesian when employed in numerical methods for
complex fixed or moving geometries
Current Work
Motivated by well-known facts that
Lagrangian method can resolve material or slip lines sharply if there is no grid tangling
Generalized curvilinear grid is often superior to
Cartesian when employed in numerical methods for complex fixed or moving geometries
Aim is to devise Lagrange-like moving grid approach for nonlinear hyperbolic system of balance law
∂
∂t q (~x, t) +
N d
X
j=1
∂
∂x j f j (q, ~x) = ψ (q, ~x)
in general N d ≥ 1 space dimension that is more robust
than aforementioned Eulerian-based method
Mathematical Formulation
Begin by considering canonical hyperbolic balance law
∂
∂t q (~x, t) +
N d
X
j=1
∂
∂x j f j (q, ~x) = ψ (q, ~x)
in Cartesian coordinate system
Hyperbolic if P N j=1 d α j (∂f j /∂q) is diagonalizable with real eigenvalues, α j ∈ lR
q ∈ lR m : vector of m state quantities
f j ∈ lR m : flux vector, j = 1, 2, · · · , N d , ψ ∈ lR m : sources
~x = (x 1 , x 2 , · · · , x N d ) : spatial vector, t ≥ 0 : time
Mathematical Formulation
Now consider a general non-rectangular domain Ω in
N d = 2 & introduce coordinate change (~x, t) 7→ (~ ξ, τ ) via
ξ = (ξ ~ 1 , ξ 2 , · · · , ξ N d ) , ξ j = ξ j (~x, t), τ = t,
that maps Ω to a logical domain Ω ˆ & also transforms balance law to a new form
−1 −0.5 0 0.5 1
−1.5
−1
−0.5 0
−1 −0.5 0 0.5
−2
−1.8
−1.6
−1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
x
y η Ω ˆ Ω
logical domain
physical domain
Mathematical Formulation
That is, using chain rule of partial differentiation, derivatives in physical space become
∂
∂t = ∂
∂τ +
N d
X
i=1
∂ξ i
∂t
∂
∂ξ i , ∂
∂x j =
N d
X
i=1
∂ξ i
∂x j
∂
∂ξ i for j = 1, 2, · · · , N d ,
yielding strong form of balance law in transformed space
∂
∂τ q ˜ ~ξ, τ + X N d
j=1
∂
∂ξ j f ˜ j
˜
q, ∇~ ξ
= ˜ ψ
˜
q, ∇~ ξ
with
˜ ∂ξ j X N d ∂ξ j !
˜ −1
Mathematical Formulation
Assume existence of inverse transformation
t = τ, x j = x j (~ ξ, t) for j = 1, 2, · · · , N d ,
To find basic geometric-metric relations between different coordinates, employ elementary differential rule
∂(τ, ~ ξ)
∂(t, ~x) = ∂(t, ~x)
∂(τ, ~ ξ)
−1
,
yielding in N d = 3 case, for example, as
1 0 0 0
∂ t ξ 1 ∂ x 1 ξ 1 ∂ x 2 ξ 1 ∂ x 3 ξ 1
∂ t ξ 2 ∂ x 1 ξ 2 ∂ x 2 ξ 2 ∂ x 3 ξ 2
∂ t ξ ∂ x ξ ∂ x ξ ∂ x ξ
= 1 J
J 0 0 0
J 01 J 11 J 21 J 31
J 02 J 12 J 22 J 32
J J J J
Mathematical Formulation
Here
J =
∂(x 1 , x 2 , x 3 )
∂(ξ 1 , ξ 2 , ξ 3 )
= det ∂(x 1 , x 2 , x 3 )
∂(ξ 1 , ξ 2 , ξ 3 )
, J 11 =
∂(x 2 , x 3 )
∂(ξ 2 , ξ 3 )
, J 21 =
∂(x 1 , x 3 )
∂(ξ 3 , ξ 2 )
, J 31 =
∂(x 1 , x 2 )
∂(ξ 2 , ξ 3 )
, J 12 =
∂(x 2 , x 3 )
∂(ξ 3 , ξ 1 )
, J 22 =
∂(x 1 , x 3 )
∂(ξ 1 , ξ 3 )
, J 32 =
∂(x 1 , x 2 )
∂(ξ 3 , ξ 1 )
, J 13 =
∂(x 2 , x 3 )
∂(ξ 1 , ξ 2 )
, J 23 =
∂(x 1 , x 3 )
∂(ξ 2 , ξ 1 )
, J 33 =
∂(x 1 , x 2 )
∂(ξ 1 , ξ 2 )
,
J 0 j = −
N d
X
i=1
J ij ∂ τ x i , j = 1, 2, 3,
and so numerically computatable ∇ξ j , j = 1, 2, 3 , as
Mathematical Formulation
Note that to complete the model When ∂ τ ~x = 0 (stationary case)
∂ t ξ = 0 ~ & ∇ ~ x ξ ~ time-independent; no more condition
Mathematical Formulation
Note that to complete the model When ∂ τ ~x = 0 (stationary case)
∂ t ξ = 0 ~ & ∇ ~ x ξ ~ time-independent; no more condition When ∂ τ ~x = ~u 0 , ~u 0 is constant (quasi-stationary case)
Both ∂ t ~ ξ & ∇ ~ x ξ ~ time-independent; no more condition
Mathematical Formulation
Note that to complete the model When ∂ τ ~x = 0 (stationary case)
∂ t ξ = 0 ~ & ∇ ~ x ξ ~ time-independent; no more condition When ∂ τ ~x = ~u 0 , ~u 0 is constant (quasi-stationary case)
Both ∂ t ~ ξ & ∇ ~ x ξ ~ time-independent; no more condition Otherwise ∂ τ ~x 6= 0 (genuine moving case)
Both ∂ t ~ ξ & ∇ ~ x ξ ~ time-dependent; require N d degree of freedom for ∂ τ ~x & numerical approach to compute
∇ ξ ~ ~x over time
Grid Movement Conditions
Here we are interested in
Lagrange-like condition ∂ τ ~x = h 0 ~u , ~u velocity, h 0 ∈ [0, 1]
Compatibility conditions for ∂ τ ∂ ξ j x i & ∂ ξ j ∂ τ x i , i.e. ,
∂
∂τ
∂x i
∂ξ j
+ ∂
∂ξ j
− ∂x i
∂τ
= 0 for i, j = 1, 2, · · · , N d
Grid Movement Conditions
Here we are interested in
Lagrange-like condition ∂ τ ~x = h 0 ~u , ~u velocity, h 0 ∈ [0, 1]
Compatibility conditions for ∂ τ ∂ ξ j x i & ∂ ξ j ∂ τ x i , i.e. ,
∂
∂τ
∂x i
∂ξ j
+ ∂
∂ξ j
− ∂x i
∂τ
= 0 for i, j = 1, 2, · · · , N d
General 1 -parameter ∂ τ ~x = h~u with h ∈ [0, 1] chosen by J or grid-angle preserving, yielding
A 0 h +
N d
X
j=1
A j ∂ ξ j h = 0 (1st order PDE constraint )
General k -parameter, k > 1 , · · ·
Shallow Water Equations
With bottom topography included, for example Cartesian coordinate case
∂
∂t
h hu i
+
N d
X
j=1
∂
∂x j
hu j
hu i u j + 1 2 gh 2
=
0
−gh ∂x ∂B
i
, i = 1, · · · , N d
Generalized coordinate case
∂
∂τ
hJ hJu i
+
N d
X
j=1
∂
∂ξ j J
hU j
hu i U j + 1 2 gh 2 ∂ξ ∂x j
i
=
0
−ghJ ∂x ∂B
i
with U j = ∂ t ξ j + P N d
i=1 u i ∂ x i ξ j , j = 1, 2, · · · , N d
h : water height, u : velocity in x -direction
Shallow Water Equations
With ∂ τ ~x = h 0 ~u & grid-metric conditions, complete model system in N d = 2 transformed space, for example, takes
∂
∂τ
hJ hJu 1
hJu 2
A B C D
+ ∂
∂ξ 1
hJU 1
hJu 1 U 1 + 1 2 gJh 2 ∂x ∂ξ 1 1 hJu 2 U 1 + 1 2 gJh 2 ∂x ∂ξ 1
2
−h 0 u 1
−h 0 u 2
0 0
+ ∂
∂ξ 2
hJU 2
hJu 1 U 2 + 1 2 gJh 2 ∂x ∂ξ 2 1 hJu 2 U 2 + 1 2 gJh 2 ∂x ∂ξ 2
2
0 0
−h 0 u 1
−h 0 u 2
= ˜ ψ
Here A = ∂ ξ 1 x 1 , B = ∂ ξ 1 x 2 , C = ∂ ξ 2 x 1 , D = ∂ ξ 2 x 2
Remarks
Hyperbolicity
In Cartesian coordinates, model is hyperbolic In generalized coord., model is hyperbolic when
h 0 6= 1 , & is weakly hyperbolic when h 0 = 1 & N d > 1 Canonical form
In Cartesian coordinates
∂
∂t q (~x, t) +
N d
X
j=1
∂
∂x j f j (q) = ψ (q)
In generalized coordinates: spatially varying fluxes
∂ N d ∂
Compressible Euler Equations
With gravity effect included, for example Cartesian coordinate case
∂
∂t
ρ ρu i
E
+
N d
X
j=1
∂
∂x j
ρu j
ρu i u j + pδ ij Eu j + pu j
=
0
−ρ ∂x ∂φ
i
−ρ~u · ∇φ
, i = 1, · · · , N d
Generalized coordinate case
∂
∂τ
ρJ ρJu i
JE
+
N d
X
j=1
∂
∂ξ j J
ρU j
ρu i U j + p ∂ξ ∂x j
i
EU j + pU j − p ∂ξ ∂t j
=
0
−ρJ ∂x ∂φ
i
−ρJ~u · ∇φ
ρ : density, p = p(ρ, e) : pressure , e : internal energy E = ρe + ρ P N d
j=1 u 2 j /2 : total energy, φ : gravitational potential
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
Extension to Multifluid
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
Model barotropic 2 -phase flow problem with
fluid component 1 & 2 characterized by Tait EOS
p(ρ) = (p 0 ι + B ι )
ρ ρ 0 ι
γ ι
− B ι , ι = 1, 2
Barotropic 2-Phase Flow
Define mixture pressure law (Shyue, JCP 2004)
p(ρ, ρe) =
(p 0 ι + B ι )
ρ ρ 0 ι
γ ι
− B ι if α = 0 or 1 (γ − 1)
ρe + ρB ρ 0
− γB if α ∈ (0, 1)
Derived from de = −pd(1/ρ) using
p(ρ, S) = A(S) (p 0 + B) ρ ρ 0
γ
− B
Here A(S) = e [(S−S 0 )/C V ] , S , C V : specific entropy & heat
at constant vol. α : volume fraction of one fluid component
Barotropic 2-Phase Flow
Transport equations for material quantities γ , B , & ρ 0 γ -based equations
∂
∂τ
1 γ − 1
+
N d
X
j=1
U j ∂
∂ξ j
1 γ − 1
= 0
∂
∂τ
γB γ − 1
+
N d
X
j=1
U j ∂
∂ξ j
γB γ − 1
= 0
∂
∂τ
J B
ρ 0
ρ
+
N d
X
j=1
∂
∂ξ j
J B
ρ 0
ρU j
= 0
Above equations are derived from energy equation &
make use of homogeneous equilibrium flow
assumption together with mass conservation law
Barotropic 2-Phase Flow
α -based equations
∂ α
∂τ +
N d
X
j=1
U j ∂α
∂ξ = 0, with z =
2
X
ι=1
α ι z ι , z = 1
γ − 1 & γB γ − 1
∂
∂τ
J B
ρ 0
ρ
+
N d
X
j=1
∂
∂ξ j
J B
ρ 0
ρU j
= 0
α -based equations (Allaire et al. , JCP 2002)
∂ α
∂τ +
N d
X
j=1
U j ∂α
∂ξ j = 0 with z =
2
X
ι=1
α ι z ι , z = 1
γ − 1 & γB γ − 1
∂ N d ∂
Multifluid Model
With (x τ , y τ ) = h 0 (u, v) & sample EOS described above, our α -based model for multifluid flow is
∂
∂τ
Jρ Jρu Jρv JE
x ξ y ξ x η y η Jρ 1 α
+ ∂
∂ξ
JρU
JρuU + y η p JρvU − x η p
JEU + (y η u − x η v)p
−h 0 u
−h 0 v 0 0 Jρ 1 αU
+ ∂
∂η
JρV
JρuV − y ξ p JρvV + x ξ p
JEV + (x ξ v − y ξ u)p 0
0
−h 0 u
−h 0 v Jρ 1 αV
= ˜ ψ
∂ α
∂τ + U ∂α
∂ξ + V ∂α
∂η = 0, plus α-averaged material quantities
Multifluid Model
For convenience, our multifluid model is written into
∂q
∂τ + f 1
∂
∂ξ , q, ∇~ ξ
+ f 2
∂
∂η , q, ∇~ ξ
= ˜ ψ with
q = [Jρ, Jρu, Jρv, JE, x ξ , y ξ , x η , y η , Jρ 1 α, α] T f 1 = ∂
∂ξ (JρU ), ∂
∂ξ (JρuU + y η p), ∂
∂ξ (JρvU − x η p), ∂
∂ξ (JEU + (y η u − x η v)p),
∂
∂ξ (−h 0 u), ∂
∂ξ (−h 0 v), 0, 0, ∂
∂ξ (Jρ 1 αU ), U ∂α
∂ξ
T
f 2 = ∂
∂η (JρV ), ∂
∂η (JρuV − y ξ p), ∂
∂η (JρvV + x ξ p), ∂
∂η (JEV + (x ξ v − y ξ u)p),
Multifluid model: Remarks
As before, under thermodyn. stability condition, our
multifluid model in generalized coordinates is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1
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 general non-barotropic
multifluid flow can be made in an analogous manner
Multifluid model: Remarks
As before, under thermodyn. stability condition, our
multifluid model in generalized coordinates is hyperbolic when h 0 6= 1 , & is weakly hyperbolic when h 0 = 1
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 general non-barotropic
multifluid flow can be made in an analogous manner
Numerical discretization?
Numerical Discretization
In 2 D, equations to be solved takes the form
∂q
∂τ + f 1
∂
∂ξ , q, ∇~ ξ
+ f 2
∂
∂η , q, ∇~ ξ
= ˜ ψ
A simple dimensional-splitting approach based on f -wave formulation of LeVeque et al. is used
Solve one-dimensional generalized Riemann problem (defined below) at each cell interfaces Use resulting jumps of fluxes (decomposed into
each wave family) of Riemann solution to update cell averages
Introduce limited jumps of fluxes to achieve high
resolution
Numerical Discretization
Employ finite volume formulation of numerical solution
Q n ij ≈ 1
∆ξ∆η Z
C ij
q(ξ, η, τ n ) dA
that gives approximate value of cell average of solution q over cell C ij = [ξ i , ξ i+1 ] × [η j , η j+1 ] at time τ n
−1.5
−1
−0.5 0
−1.6
−1.4
−1.2
−1
−0.8
−0.6
−0.4
−0.2
y η
computational grid
physical grid
Generalized Riemann Problem
Generalized Riemann problem of our multifluid model at cell interface ξ i−1/2 consists of the equation
∂q
∂τ + F i− 1
2 ,j
∂ ξ , q, ∇~ ξ
= 0 together with flux function
F i− 1
2 ,j =
f i−1,j
∂ ξ , q, ∇~ ξ
for ξ < ξ i−1/2 f ij
∂ ξ , q, ∇~ ξ
for ξ > ξ i−1/2 and piecewise constant initial data
q(ξ, 0) =
Q n i−1,j for ξ < ξ i−1/2
Q n ij for ξ > ξ i−1/2
Generalized Riemann Problem
Generalized Riemann problem at time τ = 0
ξ τ
Q n i−1,j Q n ij
Generalized Riemann Problem
Exact generalized Riemann solution: basic structure
q τ + f i−1,j
∂ ξ , q, ∇~ ξ
= 0 q τ + f i,j
∂ ξ , q, ∇~ ξ
= 0 ξ τ
Q n i−1,j Q n ij
Generalized Riemann Problem
Shock-only approximate Riemann solution: basic structure
ξ
τ
Q n i−1,j Q n ij
λ 1 λ 2
λ 3 q mL − q mL +
q mR
Z 1 = f L (q mL − ) − f L (Q n i−1,j ) Z 2 = f R (q mR ) − f R (q mL + )
Z 3 = f R (Q n ij ) − f R (q mR )
Numerical Discretization
Basic steps of a dimensional-splitting scheme ξ -sweeps: solve
∂q
∂τ + f 1
∂
∂ξ , q, ∇~ ξ
= 0
updating Q n ij to Q ∗ i,j η -sweeps: solve
∂q
∂τ + f 2
∂
∂η , q, ∇~ ξ
= 0
updating Q ∗ ij to Q n+1 i,j
Numerical Discretization
That is to say,
ξ -sweeps: we use
Q ∗ ij = Q n ij − ∆τ
∆ξ
F i+ − 1
2 ,j − F i− + 1
2 ,j
− ∆τ
∆ξ ˜ Z i+ 1
2 ,j − ˜ Z i− 1
2 ,j
with Z ˜ i− 1
2 ,j = 1 2
m w
X
p=1
sign
λ p i− 1
2 ,j
1 − ∆τ
∆ξ λ
p
i− 1 2 ,j
Z ˜ p i− 1
2 ,j
η -sweeps: we use
Q n+1 ij = Q ∗ ij − ∆τ
∆η
G i,j+ − 1
2
− G +
i,j− 1 2
− ∆τ
∆η ˜ Z i,j+ 1
2 − ˜ Z i,j− 1
2
with Z ˜ i,j− 1
2 = 1 2
m w
X sign λ p
i,j− 1
1 − ∆τ
∆η λ
p i,j− 1
Z ˜ p
i,j− 1
Numerical Discretization
Flux-based wave decomposition
f i,j − f i−1,j =
m w
X
p=1
Z i−1/2 p =
m w
X
p=1
λ p i−1/2 W i−1/2 p
Some care should be taken on the limited jump of fluxes W ˜ p , for p = 2 (contact wave), in particular to ensure
correct pressure equilibrium across material interfaces MUSCL-type (slope limited) high resolution extension is not simple as one might think of for multifluid problems Splitting of discontinuous fluxes at cell interfaces:
significance ?
First order or high resolution method for geometric
conservation laws: significance to grid uniformity ?
Lax’s Riemann problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 0.5 1 1.5 2
x
ρ
Exact h 0 =0 h 0 =0.99
0 0.5 1
0 0.5 1 1.5 2
x
u
0 0.5 1
0 1 2 3 4
x
p
Lax’s Riemann problem
Physical grid coordinates at selected times
Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system
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
x
time
Woodward-Colella’s problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
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
Woodward-Colella’s problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 5 10 15 20
ρ
Fine grid h 0 =0 h 0 =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
Woodward-Colella’s problem
h 0 = 0 Eulerian result
h 0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 2 4 6 8
ρ
0 0.5 1
−5 0 5 10 15
0 0.5 1
0 200 400 600
u p
x x
x
t = 0.038 t = 0.038
t = 0.038
Woodward-Colella’s problem
Physical grid coordinates at selected times
Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.01 0.02 0.03 0.04
x
time
2 D Riemann problem
With initial 4 -shock wave pattern
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ρ u v p
=
1.5
0 0 1.5
0.532 1.206
0 0.3
0.138 1.206 1.206 0.029
0.532 0 1.206
0.3
2 D Riemann problem
With initial 4 -shock wave pattern Lagrange-like result
Occurrence of simple Mach reflection
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
Density Pressure Physical grid
2 D Riemann problem
With initial 4 -shock wave pattern Eulerian result
Poor resolution around simple Mach reflection
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
Density Pressure Physical grid
More Examples
Two-dimensional case
Radially symmetric problem Underwater explosion
Shock-bubble interaction Helium bubble case
Refrigerant bubble case Three-dimensional case
Underwater explosion Shock-bubble interaction
Helium bubble case
Refrigerant bubble case
Conclusion
Have described fluid-mixture type algorithm in generalized moving-curvilinear grid
Have shown results in 1 , 2 & 3 D to demonstrate
feasibility of method for practical problems
Conclusion
Have described fluid-mixture type algorithm in generalized moving-curvilinear grid
Have shown results in 1 , 2 & 3 D to demonstrate feasibility of method for practical problems
Future direction
Efficient & accurate grid movement strategy Static & Moving 3 D geometry problems
Weakly compressible flow Viscous flow extension
· · ·
Conclusion
Have described fluid-mixture type algorithm in generalized moving-curvilinear grid
Have shown results in 1 , 2 & 3 D to demonstrate feasibility of method for practical problems
Future direction
Efficient & accurate grid movement strategy Static & Moving 3 D geometry problems
Weakly compressible flow Viscous flow extension
· · ·
Thank You
Euler in General. Coord.
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
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 2 (α 1 c − up E )/2c 2 (α 2 c − vp E )/2c 2 p E /2c 2
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 2 (β 1 c − up E )/2c 2 (β 2 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
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 η )
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
A = uy η − vx η , B = vx ξ − uy ξ , C = u ξ y η + v η x ξ − u η y ξ − v ξ x η
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 θ)
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 ?
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)
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
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
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] ?
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:
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 1 = uy η , B 1 = uy ξ , C 1 = u ξ y η − u η y ξ
A 2 = −vx η , B 2 = vx ξ , C 2 = v η x ξ − v ξ x η
Shock in molybdenum over MORB
Density Pressure
Molybdenum
MORB time = 50µ s
time = 100µ s
Falling Liquid Drop Problem
Interface capturing with gravity
100 200 300 400 500 600 700 800 900 1000
air
t = 0
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900 1000
air
t = 0.2s
Falling Liquid Drop Problem
Interface diffused badily
100 200 300 400 500 600 700 800 900
air
t = 0.4s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.6s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.8s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 1s
Shock-Bubble Interaction
Volume tracking for material interface
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Shock-Bubble Interaction
Small moving irregular cells: stability & accuracy
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
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
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
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
Moving Cylindrical Vessel
initial condition interface
air helium
Moving Cylindrical Vessel
Two sample grid systems used in computation
−1
−0.5 0 0.5 1
−0.6
−0.4
−0.2 0 0.2 0.4 0.6
Cartesian body-fitted
Moving Cylindrical Vessel
Cartesian grid results with embedded moving boundary
time=2
Moving Cylindrical Vessel
Cartesian grid results with embedded stationary boundary
time=2
