SLIDE 1
Simple Eulerian Methods for Compressible Fluids in Domains with - - PowerPoint PPT Presentation
Simple Eulerian Methods for Compressible Fluids in Domains with - - PowerPoint PPT Presentation
Simple Eulerian Methods for Compressible Fluids in Domains with Moving Boundaries Alina Chertock North Carolina State University chertock@math.ncsu.edu joint work with Armando Coco, Yuanzhen Cheng, Smadar Karni, Alexander Kurganov and Giovanni
SLIDE 2
SLIDE 3
Compressible Euler Equations
ρ ρu ρv E
t
+ ρu ρu2 + p ρuv u(E + p)
x
+ ρv ρuv ρv2 + p v(E + p)
y
= 0 E = p γ − 1 + ρ 2(u2 + v2)
- ρ is the fluid density;
- u, v are the velocities;
- E is the total energy;
- p is the pressure;
- c =
γp ρ is the speed of sound.
1
SLIDE 4
Two Types of Problems
- Multicomponent Flow: We assume that γ is a piecewise constant function
propagation with the fluid velocity according to γt + uγx + vγy = 0.
- Problems with Changing Geometries:
– The right boundary is a moving solid wall with a prescribed equation of motion x = xB(t)
B
Moving Boundary
x
– Flow around an (oscillating) solid circle.
I II
2
SLIDE 5
Conservative Methods
wt + F(w)x = 0 ⇔ ρ ρu E ργ
t
+ ρu ρu2 + p u(E + p) uργ
x
= 0 Godunov-type discretization: ¯ wn+1
j
= ¯ wn
j − ∆t
∆x
- Hj+1/2 − Hj−1/2
- ,
Hj+1/2 = H(. . . , ¯ wn
j , ¯
wn
j+1, . . .),
H(. . . , w, w, . . .) = F(w). The major difficulty: to ensure that the pressure equilibrium between fluid components is not disturbed by the numerical scheme.
3
SLIDE 6
What is going wrong? (S. Karni)
4
SLIDE 7
Multicomponent Flow
wt + f(w)x = 0 ⇐ ⇒ ρ ρu E
t
+ ρu ρu2 + p u(E + p)
x
= 0 EOS: p = (γ − 1)
- E − 1
2ρu2
- − γp∞
w = (ρ, ρu, E)T We assume that γ is a piecewise constant function propagation with the fluid velocity according to γt + uγx = 0.
5
SLIDE 8
Multi-Fluid Example
0.2 0.4 0.6 0.8 1 0.98 1 1.02 1.04 1.06 1.08 1.1 1.12 pressure
CONSERVATIVE METHOD EXACT SOLUTION
(ρ, u, p, γ)T =
- (1.0, 1.0, 1.0, 1.6)T,
if x < 0.25, (0.1, 1.0, 1.0, 1.4)T, if x > 0.25. Different EOS...
6
SLIDE 9
What is going wrong?
J-1 J J+1 tn tn+1 xJ1/2 xJ+1/2 material interface
- Fluxes are computed using the information in the mixed cell.
- No valid EOS in mixed cells.
7
SLIDE 10
Some Models/Schemes
- Conserve not only total mass and momentum but also total energy.
– Volume-of-Fluid Method (Noh & Woodward; Collela, Glaz & Ferguson). – Internal Energy Correction Algorithm (Jenny, Mueller & Thomann). – A Fluid-Mixture Type Algorithm (Shyue).
- Conserve only total mass and momentum but not total energy.
– Pressure Evolution Method (Karni). – Ghost-Fluid Method for the Poor (Abgrall & Karni).
- Nonconservative
– Ghost-Fluid Method (Fedkiw, Aslam, Merriman & Osher).
- – Conservative locally moving mesh method for multifluid flows (A.C.
& A. Kurganov) – A second-order finite-difference method for compressible fluids in domains with moving boundaries, (A. C., A. Coco, A. Kurganov and
- G. Russo)
– Interface tracking method for compressible multifluids, (A.C., S. Karni & A. Kurganov)
8
SLIDE 11
Semi-Discrete Central-Upwind Scheme – 1-D Case
wt + f(w)x = 0 ⇐ ⇒ ρ ρu E
t
+ ρu ρu2 + p u(E + p)
x
= 0 w := (ρ, ρu, E)T Computational cells: Cj := [xj−1
2, xj+1 2]
The cell averages: wj(t) := 1 ∆x
- Cj
w(x, t) dx Time evolution: d dtwj(t) = − Hj+1
2(t) − Hj−1 2(t)
∆x
- Hj+1
2
- : numerical fluxes
9
SLIDE 12
{wj(t)} → w(·, t) →
- w±
j+1
2(t)
- →
- Hj+1
2(t)
- → {wj(t + ∆t)}
(Discontinuous) piecewise-linear reconstruction:
- w(x, t) := wj(t) + (wx)j(x − xj),
x ∈ Cj It is conservative, second-order accurate, and non-oscillatory provided the slopes, {(wx)j}, are computed by a nonlinear limiter Example — Generalized Minmod Limiter (wx)j = minmod
- θwj − wj−1
∆x , wj+1 − wj−1 2∆x , θwj+1 − wj ∆x
- where
minmod(z1, z2, ...) := minj{zj}, if zj > 0 ∀j, maxj{zj}, if zj < 0 ∀j, 0,
- therwise,
and θ ∈ [1, 2] is a constant
10
SLIDE 13
{wj(t)} → w(·, t) →
- w±
j+1
2(t)
- →
- Hj+1
2(t)
- → {wj(t + ∆t)}
w±
j+1
2(t) are the point values of
- w(x, t) = wj + (wx)j(x − xj),
x ∈ Cj at xj+1
2:
w+
j+1
2 :=
w(xj+1
2 + 0, t) = wj+1 − ∆x
2 (wx)j+1 w−
j+1
2 :=
w(xj−1
2 − 0, t) = wj + ∆x
2 (wx)j
11
SLIDE 14
{wj(t)} → w(·, t) →
- w±
j (t)
- →
- Hj+1
2(t)
- → {wj(t + ∆t)}
Hj+1
2(t) = H(w−
j+1
2(t), w+
j+1
2(t))
Example — Central-Upwind Flux: [Kurganov, Lin, Noelle, Petrova, Tadmor, et al.; 2000–2007] Hj+1
2 =
a+
j+1
2f(w−
j+1
2) − a−
j+1
2f(w+
j+1
2)
a+
j+1
2 − a−
j+1
2
+ a+
j+1
2a−
j+1
2
a+
j+1
2 − a−
j+1
2
- w+
j+1
2 − w−
j+1
2
- a±
j+1
2(t) = a±
j+1
2
- w+
j+1
2, w−
j+1
2
- are one-sided local speeds, estimated using
the largest and the smallest eigenvalues of the Jacobian ∂f ∂w. d dtwj(t) = − Hj+1
2(t) − Hj−1 2(t)
∆x
12
SLIDE 15
Boundary Treatment in 1-D
- Assume that at some time level t the wall is located in cell J
- All computational cells are divided into 3 groups:
– Internal cells with reliable data – Boundary cells with unreliable data – External cells with no data
xJ5/2 xJ3/2 xJ1/2 xJ+1/2 xJ+3/2
Moving Boundary
xB
13
SLIDE 16
Boundary Treatment in 1-D
- The cell averages wJ(t) will not be evolved
- To evolve the cell averages wJ−1(t) we need:
– the numerical flux HJ−1
2(t)
– the point values w+
J−1
2(t)
xJ5/2 xJ3/2 xJ1/2 xJ+1/2 xJ+3/2
Moving Boundary
wJ1/2
J1
xB w
+
14
SLIDE 17
Solid Wall Extrapolation
The solid wall extrapolation of the solution values in cell (J − 1), ρJ−1, uJ−1 = (ρu)J−1 ρJ−1 , pJ−1 = (γ − 1)
- EJ−1 − ρJ−1u2
J−1
2
- To obtain the ”missing” ghost values:
ρ d2xB(t) dt2
- x=xB(t) = −px
- x=xB(t),
∂p ∂ρ
- x=xB(t) = γp
ρ
- x=xB(t)
xJ5/2 xJ3/2 xJ1/2 xJ+1/2 xJ+3/2 wJ1 wJ1/2 xB
+
wgh
15
SLIDE 18
Solid Wall Extrapolation
- Second-Order Approximation:
ρJ−1 + ρgh 2 · d2xB(t) dt2 = −pgh − pJ−1 ∆x pgh − pJ−1 ρgh − ρJ−1 = γ · pgh + pJ−1 ρgh + ρJ−1
- First-Order Approximation:
ρgh := ρJ−1, ugh := 2 ˙ xB(t) − uJ−1, pgh := pJ−1 ˙ xB(t): velocity of the wall at time t
16
SLIDE 19
Phase Space Interpolation
if u∗ > 0 then w+
J−1
2 =
w∗, if u∗ − c∗ < 0 wJ−1,
- therwise
else w+
J−1
2 =
w∗, if u∗ + c∗ > 0 wgh,
- therwise
gh
ρ u p
J−1 J−1 J−1
p
gh
ρ u p
* * *
ρ u
gh
Once w+
J−1
2
are available, we compute the slopes (wx)J−1 using w−
J−3
2, wJ−1 and w+
J−1
2
xJ5/2 xJ3/2 xJ1/2 xJ+1/2 xJ+3/2
Moving Boundary
- w
w wJ1/2
J3/2 J1 +
xB
17
SLIDE 20
Solution of the RP – A Very Simple Structure
- If ugh < uJ−1, the solution consists of three constant states, connected
by two shock waves. The intermediate state is equal to: p∗ = pJ−1 + (γ + 1)ρJ−1 ( ˙ xB(t) − uJ−1)2 4 1 +
- 1 +
- 4cJ−1
(γ + 1) ( ˙ xB(t) − uJ−1) 2 u∗ = ˙ xB(t), ρ∗ = ρJ−1 (γ − 1)pJ−1 + (γ + 1)p∗ (γ + 1)pJ−1 + (γ − 1)p∗
- If ugh > uJ−1, the solution consists of three constant states, connected
by two rarefaction waves.The intermediate state is given by: p∗ = pJ−1
- 1 − (γ − 1) ( ˙
xB(t) − uJ−1) 2cJ−1 2γ/(γ−1) u∗ = ˙ xB(t), ρ∗ = ρJ−1 p∗ pJ−1 1/γ
18
SLIDE 21
Time Evolution
d dtwj(t) = − Hj+1
2(t) − Hj−1 2(t)
∆x : wj(t) → wj(t + ∆t)
- If xB(t + ∆t) ∈ CJ, the evolution step is completed.
- If xB(t + ∆t) ∈ CJ−1, then
– Cell J − 1 becomes the new boundary cell – Cell J turns into an external one – The computed values wJ−1(t + ∆t) will not be used any more
J−1
xJ3/2 xJ1/2 xJ+1/2 t+ t WJ1 x t t
19
SLIDE 22
- If xB(t + ∆t) ∈ CJ+1, then
– Cell (J + 1) becomes a new boundary cell – Cell J turns into an internal one – We use the phase space interpolation between (ρJ−1, uJ−1, pJ−1) and (ρgh, ugh, pgh) to obtain (ρ∗, u∗, p∗) and then wJ(t + ∆t) := w∗
∈
xJ1/2 xJ+1/2 xJ+3/2 WJ t+ t x t t
20
SLIDE 23
Important Remarks
- The proposed boundary cell treatment procedures use the O(∆x)
accurate information on the location of the boundary instead of the exact one and thus the overall accuracy of the introduced method cannot be higher than the first order
- The main advantages of the proposed method are its simplicity and
robustness
- To achieve a higher resolution, one may use an adaptive mesh refinement
technique near the boundary
- The proposed boundary treatment procedure is not conservative.
However, the conservation error occurs only along the boundary surface
- f co-dimension one. Thus the total conservation error is expected to be
proportional to the mesh size, as confirmed by our numerical experiments
21
SLIDE 24
Example — Piston Problem
The tube is bounded by a stationary solid wall at x = 0 and a piston at x = xB(t), which is free to move without friction in the tube. The motion
- f the piston is governed by the Newton equation:
d2xB(t) dt2 = A m (p(xB(t), t) − pout) , where A is the area of the piston, m is its mass, A/m = 2, and pout = 2 is an external pressure. The velocity of the piston is obtained from d ˙ xB(t) dt = A m (p∗ − pout) , where p∗ is the solution of the corresponding RP. Initial Condition: u(x, 0) ≡ 0, ˙ xB(0) = 0, ρ(x, 0) ≡ p(x, 0) ≡ 1, xB(0) = 1.
22
SLIDE 25
What happens?
- Due to the difference between the interior and exterior pressure, the
piston starts moving into the tube, rising the pressure inside it.
- When the pressure in the tube becomes larger than the outside one, the
piston decelerates and turns back.
- The pattern repeats and the motion of the piston becomes oscillatory.
Remark: The piston problem is slightly different since the piston motion is not a-priori prescribed but depends on the difference between the pressures inside and outside the tube. However, when the boundary cell extrapolation procedure is used, the elastic properties of the piston are neglected and the wall (piston) is assumed to be solid.
23
SLIDE 26
0.5 1 1 2 3 4 5 6 7 8 9 10 x t 0.9 0.92 0.94 0.96 0.98 8 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 x t ∆ x=1/100 ∆ x=1/200 ∆ x=1/400 ∆ x=1/2000
The piston trajectory as a function of time (left) computed on four different uniform grids with ∆x = 1/100, 1/200, 1/400, and 1/2000. Zoom at [0.89, 0.98] × [7.95, 8.8] (right).
24
SLIDE 27
Example — Solid Oscillating Boundary
1-D :
wt + f(w)x = 0, w := (ρ, ρu, E)T, E = p γ − 1 + ρu2 2 The initial data correspond to a right moving shock: (ρ, u, p) =
- (4/3, 35/99, 1.5),
x < 0.5 (1, 0, 1), 0.5 < x < xB(0) γ = 1.4 xB(t) = 0.9 + 0.1 sin(10πt)
25
SLIDE 28
0.5 1 1 2 3 4 5
ρ (x,t=0.05)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200 0.5 1 1 2 3 4 5
ρ (x,t=0.15)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200 0.5 1 1 2 3 4 5
ρ (x,t=0.3)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200 0.5 1 1 2 3 4 5
ρ (x,t=0.4)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200
26
SLIDE 29
0.5 1 −3 −2 −1 1 2 3
u(x,t=0.05)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200 0.5 1 −3 −2 −1 1 2 3
u(x,t=0.15)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200 0.5 1 −3 −2 −1 1 2 3
u(x,t=0.3)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200 0.5 1 −3 −2 −1 1 2 3
u(x,t=0.4)
∆ x=1/100 ∆ x=1/200 ∆ x=1/51200
27
SLIDE 30
Mass conservation error t = 0.05 t = 0.4 ∆x Error Rate Error Rate 1/200 4.371 · 10−3 – 2.580 · 10−3 – 1/400 1.995 · 10−3 1.13 1.548 · 10−3 0.74 1/800 9.673 · 10−4 1.04 8.589 · 10−4 0.85 1/1600 4.565 · 10−4 1.08 4.279 · 10−4 1.00 1/3200 2.202 · 10−4 1.05 2.221 · 10−4 0.95
28
SLIDE 31
2-D compressible Euler equations
ρ ρu ρv E
t
+ ρu ρu2 + p ρuv u(E + p)
x
+ ρv ρuv ρv2 + p v(E + p)
y
= 0 EOS: E = p γ − 1 + ρ 2(u2 + v2), γ = const
- Boundary conditions:
– left: inflow – right: outflow – top, bottom, circle: solid wall
I II
29
SLIDE 32
Example — No Shock - Vertically Moving Ball
Initial Data: ρ(x, y, 0) ≡ p(x, y, 0) ≡ 1, u(x, y, 0) ≡ v(x, y, 0) ≡ 0 Vertical Movement of the Cylinder: xc(t) ≡ 0.5, yc(t) = 0.5 + 0.1 sin(10πt)
30
SLIDE 33
∆ x=∆ y=1/400, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/400, t=0.2
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.2
0.5 1 0.2 0.4 0.6 0.8 1
31
SLIDE 34
Example — No Shock - Diagonally Moving Ball
Initial Data: ρ(x, y, 0) ≡ p(x, y, 0) ≡ 1, u(x, y, 0) ≡ v(x, y, 0) ≡ 0 Diagonal Movement of the Cylinder: xc(t) = yc(t) = 0.5 + 0.1 √ 2 sin(10πt)
32
SLIDE 35
∆ x=∆ y=1/400, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/400, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
33
SLIDE 36
Example — Moving Shock - Moving Ball
Initial Data (Right-Moving Shock):
I:
ρ = 4
3,
u = 35
99,
v = 0, p = 3
2
II:
ρ = 1, u = 0, v = 0, p = 1 Vertical Movement of the Cylinder: xc(t) ≡ 0.5, yc(t) = 0.5 + 0.1 sin(10πt)
34
SLIDE 37
∆ x=∆ y=1/400, t=0.05
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/400, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.05
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.1
0.5 1 0.2 0.4 0.6 0.8 1
35
SLIDE 38
∆ x=∆ y=1/400, t=0.15
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/400, t=0.2
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.15
0.5 1 0.2 0.4 0.6 0.8 1
∆ x=∆ y=1/800, t=0.2
0.5 1 0.2 0.4 0.6 0.8 1
36
SLIDE 39
Mass conservation error
t = 0.05 t = 0.1 ∆x = ∆y Error Rate Error Rate 1/100 5.056 · 10−4 – 3.984 · 10−3 – 1/200 2.453 · 10−4 1.04 2.513 · 10−3 0.66 1/400 2.105 · 10−4 0.22 1.879 · 10−3 0.42 1/800 9.038 · 10−5 1.22 1.017 · 10−3 0.89 t = 0.15 t = 0.2 ∆x = ∆y Error Rate Error Rate 1/100 3.505 · 10−3 – 7.979 · 10−3 – 1/200 2.096 · 10−3 0.74 4.410 · 10−3 0.86 1/400 1.462 · 10−3 0.52 2.992 · 10−3 0.56 1/800 8.112 · 10−4 0.85 1.857 · 10−3 0.69
37
SLIDE 40
Second-Order Boundary Conditions
38
SLIDE 41
1-D Finite-Difference Scheme
dUj dt = − ˆ fj+1
2 − ˆ
fj−1
2
∆x , j = 1, . . . , J, Uj(t) :≈ U(xj, t) The numerical fluxes { ˆ fj+1
2} are computed as follows:
- Split the flux function f as
f(Uj) = f +
j + f − j ,
f ±
j := 1
2
- f(Uj) ± ajUj
- ,
aj := |uj| + cj
- Evaluate the values of f + and f − at the cell interfaces xj+1
2:
f E
j := f + j + ∆x
2 (fx)+
j ,
f W
j
:= f −
j − ∆x
2 (fx)−
j ,
where the slopes are computed using the generalized minmod limiter
- The numerical fluxes are given by
ˆ fj+1
2 = f E
j + f W j+1.
39
SLIDE 42
Second-Order Boundary Conditions – 1-D
xB x xJ xJ+1 x
J−1 J+2
More accurate solid wall boundary conditions: u(xB) = x′
B
∂p ∂x
- x=xB
= −ρ(xB)x′′
B
∂ρ ∂x
- x=xB
= 1 c2(xB) · ∂p ∂x
- x=xB
= − ρ(xB) c2(xB)x′′
B
, c := γp ρ Second-order discretization: αugh + (1 − α)uJ = x′
B
pgh − pJ ∆x = − [αρgh + (1 − α)ρJ] x′′
B
ρgh − ρJ ∆x = −αρgh + (1 − α)ρJ αc2
gh + (1 − α)c2 J
x′′
B
, α := xB − xJ ∆x
40
SLIDE 43
Second-Order Boundary Conditions – 2-D
- Internal points (xj, yk) ∈ Ω ((j, k) ∈ I).
- Ghost points (xj, yk) ((j, k) ∈ G).
- Inactive points.
These are grid points located outside of Ω, but not ghost points.
41
SLIDE 44
n
τ
Ω
n
τ
Ω
φ: level-set function, V (x, y, t): prescribed normal boundary velocity, ∂Ω : ∂φ ∂t + V |∇φ| = 0.
- (x, y) ∈ R2: φ(x, y) < 0
- : domain Ω
- (x, y) ∈ R2: φ(x, y) = 0
- : boundary ∂Ω
- (x, y) ∈ R2: φ(x, y) > 0
- : obstacle
Example – signed distance function: φ(x, y) =
- −dist
- (x, y), ∂Ω
- ,
if (x, y) ∈ Ω dist
- (x, y), ∂Ω
- ,
if (x, y) / ∈ Ω n and τ: normal and tangential unit vectors κ: signed curvature of ∂Ω
42
SLIDE 45
- Condition on the normal component of the velocity: un = V
- Condition on the pressure:
∂p ∂n = −ρu2
τκ − ρuτ
∂τ ∂t · n − ρDV Dt where DV Dt := ∂V ∂t + u ∂V ∂x + v ∂V ∂y and ∂τ ∂t · n = −∂φ ∂x · ∂2φ ∂y ∂t + ∂φ ∂y · ∂2φ ∂x ∂t
- Condition on the density:
∂ρ ∂n = ρ γp · ∂p ∂n
- Condition on the tangential component of the velocity:
∂uτ ∂n = ∂V ∂τ + uτκ
43
SLIDE 46
1-D Accuracy Test
- Computational domain: [0, 1]
- Solid wall boundary conditions imposed at both ends
- Left boundary is fixed, right boundary oscillates according to the a-priori
prescribed equation of motion xB(0) = 0.9, x′
B(t) = 0.25 sin3(2πt).
- The initial data: ρ(x, 0) ≡ 1,
u(x, 0) ≡ 0, p(x, 0) ≡ 1 ∆x ρ∆x − ρ2∆x1 r1 ρ∆x − ρ2∆x∞ r∞ 1/400 2.05 · 10−4 – 3.90 · 10−3 – 1/800 5.03 · 10−5 2.03 1.31 · 10−3 1.58 1/1600 1.29 · 10−5 1.96 3.82 · 10−4 1.77 1/3200 3.25 · 10−6 1.99 1.04 · 10−4 1.88
44
SLIDE 47
2-D Accuracy Test
- Computational domain: [0, 1] × [0, 1]
- Fluid domain is Ω = [0, 1] × [0, 1] \ BR(t), where BR(t) is the rigid ball
- f radius R = 0.1 and center (xc(t), yc(t))
- Inflow boundary conditions at the left boundary and outflow boundary
conditions at the right boundary, while the top and the bottom boundaries
- f the domain as well as the boundaries of the moving circle are assumed
to be solid walls. We compare:
- the second-order boundary conditions
- the first-order boundary conditions, which are obtained using
∂ρ ∂n = 0, un = V, ∂uτ ∂n = 0, ∂p ∂n = 0.
45
SLIDE 48
To better illustrate the advantages of the second-order boundary treatment, we will plot the 1-D slices of the computed solutions along the three directions:
46
SLIDE 49
2-D Accuracy Test
We consider a simple wave that propagates around a steady disk with xc(t) ≡ 0.6 and yc(t) ≡ 0.5 in a constant medium. The initial conditions are (ρ, u, v, p) =
- ρ(x, y),
u(x, y), 0, p(x, y)
- , if |x − 0.35| < 0.25,
- ρ0, 0, 0, p0
- ,
- therwise,
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 1 1.1 1.2 1.3 1.4 1.5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 1 1.1 1.2 1.3 1.4 1.5
47
SLIDE 50
∆x = ∆y ρ∆x − ρ2∆x1 rate u∆x − u2∆x1 rate 1/100 2.71 · 10−3 – 2.62 · 10−3 – 1/200 1.09 · 10−3 1.31 1.08 · 10−3 1.28 1/400 3.12 · 10−4 1.80 3.11 · 10−4 1.80
0.02 0.04 0.06 0.08 0.1 1.1 1.2 1.3 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.2 − 0.1 0.1
un
0.02 0.04 0.06 0.08 0.1 0.3 0.4 0.5 0.6
uτ
0.02 0.04 0.06 0.08 0.1 1.2 1.3 1.4 1.5 1.6
p
0.02 0.04 0.06 0.08 0.1 1.1 1.2 1.3 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.2 − 0.1 0.1
un
0.02 0.04 0.06 0.08 0.1 0.3 0.4 0.5 0.6
uτ
0.02 0.04 0.06 0.08 0.1 1.2 1.3 1.4 1.5 1.6
p
48
SLIDE 51
Example — Shock Hitting a Steady Ball
0.02 0.04 0.06 0.08 0.1 1.1 1.2 1.3 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.05 − 0.025 0.025 0.05
un
0.02 0.04 0.06 0.08 0.1 0.4 0.5 0.6 0.7 0.8 0.9
uτ
0.02 0.04 0.06 0.08 0.1 1.1 1.2 1.3 1.4 1.5
p
0.02 0.04 0.06 0.08 0.1 1.1 1.2 1.3 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.05 − 0.025 0.025 0.05
un
0.02 0.04 0.06 0.08 0.1 0.4 0.5 0.6 0.7 0.8 0.9
uτ
0.02 0.04 0.06 0.08 0.1 1.1 1.2 1.3 1.4 1.5
p
First-order BC (left) vs. second-order BC (right)
49
SLIDE 52
Example — Shock Hitting a Steady Ball
0.02 0.04 0.06 0.08 0.1 1.2 1.25 1.3 1.35 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.1 0.1 0.2
un
0.02 0.04 0.06 0.08 0.1 0.25 0.35 0.45 0.55
uτ
0.02 0.04 0.06 0.08 0.1 1.4
p
0.02 0.04 0.06 0.08 0.1 1.2 1.25 1.3 1.35 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.1 0.1 0.2
un
0.02 0.04 0.06 0.08 0.1 0.25 0.35 0.45 0.55
uτ
0.02 0.04 0.06 0.08 0.1 1.4
p
First-order BC (left) vs. second-order BC (right)
50
SLIDE 53
Example — Slow Moving Ball
0.02 0.04 0.06 0.08 0.1 0.8 1 1.2 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.5 − 0.25 0.25
un
0.02 0.04 0.06 0.08 0.1 − 0.05 − 0.025 0.025 0.05
uτ
0.02 0.04 0.06 0.08 0.1 0.7 0.9 1.1 1.3 1.5 1.7
p
0.02 0.04 0.06 0.08 0.1 0.8 1 1.2 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.5 − 0.25 0.25
un
0.02 0.04 0.06 0.08 0.1 − 0.05 − 0.025 0.025 0.05
uτ
0.02 0.04 0.06 0.08 0.1 0.7 0.9 1.1 1.3 1.5 1.7
p
First-order BC (left) vs. second-order BC (right)
51
SLIDE 54
Example — Slow Moving Ball
0.02 0.04 0.06 0.08 0.1 0.8 1 1.2 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.5 − 0.25 0.25
un
0.02 0.04 0.06 0.08 0.1 − 0.05 0.05 0.1
uτ
0.02 0.04 0.06 0.08 0.1 0.7 0.9 1.1 1.3 1.5
p
0.02 0.04 0.06 0.08 0.1 0.8 1 1.2 1.4
ρ
0.02 0.04 0.06 0.08 0.1 − 0.5 − 0.25 0.25
un
0.02 0.04 0.06 0.08 0.1 − 0.05 0.05 0.1
uτ
0.02 0.04 0.06 0.08 0.1 0.7 0.9 1.1 1.3 1.5
p
First-order BC (left) vs. second-order BC (right)
52
SLIDE 55
Example — Slow Moving Ball
0.02 0.04 0.06 0.08 0.1 0.975 1 1.025
ρ
0.02 0.04 0.06 0.08 0.1 − 0.025 0.025
un
0.02 0.04 0.06 0.08 0.1 − 0.05 − 0.025 0.025 0.05 0.075
uτ
0.02 0.04 0.06 0.08 0.1 0.95 0.975 1 1.025
p
0.02 0.04 0.06 0.08 0.1 0.975 1 1.025
ρ
0.02 0.04 0.06 0.08 0.1 − 0.025 0.025
un
0.02 0.04 0.06 0.08 0.1 − 0.05 − 0.025 0.025 0.05 0.075
uτ
0.02 0.04 0.06 0.08 0.1 0.95 0.975 1 1.025
p
First-order BC (left) vs. second-order BC (right)
53
SLIDE 56
Multifluid Example– Light Cylindrical Helium Bubble
−1.5 −1 −0.5 0.5 1 1.5 −0.4 −0.2 0.2 0.4
γ=7/5 γ=5/3
If x > 0.75 ρ = 4 3, u = − 707 2000, v = 0, p = 3 2. Inside the bubble ρ = 4 29, u = v = 0, p = 1. If x < 0.75 ρ = 1, u = 0, v = 0, p = 1.
54
SLIDE 57
−1.5 −1 −0.5 0.5 1 1.5 −0.4 −0.2 0.2 0.4
γ=7/5 γ=5/3
Upon hitting the bubble, the shock wave partly refracts, and partly reflects. The refracted shock is curved and propagates faster inside the helium bubble than in the air, as sound speed in helium is higher than sound speed in
- air. The bubble compresses under the impact of the shock and is set into
motion. Being lighter than air, the helium accelerates more under the impact of the shock, and shear is produced across the bubble interface. The refracted shock sweeps over the bubble, and continues to interacts with the bubble interface, giving rise to highly complicated refraction patterns. After the shock completes its sweep over the bubble, the bubble continues to evolve slowly into a familiar kidney shape as a result of the vorticity generated at the bubble interface, before losing its integrity and breaking
- up. Interface roll-up due to shear instability may be observed, as well as
wave reflection from the solid top/bottom boundaries.
55
SLIDE 58
56
SLIDE 59
57
SLIDE 60
Multifluid Example – Heavy Cylindrical Bubble
−1.5 −1 −0.5 0.5 1 1.5 −0.4 −0.2 0.2 0.4
γ=1.4 γ=1.249
If x > 0.75 ρ = 4 3, u = − 707 2000, v = 0, p = 3 2. Inside the bubble ρ = 3.1538, u = v = 0, p = 1. If x < 0.75 ρ = 1, u = 0, v = 0, p = 1.
58
SLIDE 61
−1.5 −1 −0.5 0.5 1 1.5 −0.4 −0.2 0.2 0.4
γ=1.4 γ=1.249
Upon being hit by the shock wave, the bubble compresses and undergoes a deformation, and the shock partly reflects and partly refracts. The speed
- f sound inside the R22 bubble is lower than outside, causing the refracted
shock wave to move more slowly than the shock outside the bubble. The bubble is heavier than the surrounding air, and accelerates less than the air under the impact of the shock. As a result, vorticity is generated at the bubble interface. The refracted shock can be observed to focus inside the bubble, causing a rise in pressure and resulting in a visible forward jet along the center of the bubble. The shear at the boundary interface causes the interface to roll-up, as the bubble continues to evolve under the induced vorticity field.
59
SLIDE 62
60
SLIDE 63
61
SLIDE 64
2-D Hughes Pedestrian Flow Model
ρt + f(ρ, φx)x + g(ρ, φy)y = 0 ||∇φ|| = 1 u
- ρ: pedestrian density
- u = umax
- 1 −
ρ ρmax
- : isotropic walking speed
- umax: free-flow speed
- ρmax: jam density
- φ: cost potential function satisfying the Eikonal equation
- (f, g)T: flux given by (f, g)T = −ρu ∇φ
||∇φ|| = −ρu2∇φ The pedestrian route choice strategy satisfies the reactive user equilibrium principle in which a pedestrian choose a route to minimize the instantaneous travel cost to the destination
62
SLIDE 65
Inside the obstacles and at solid wall boundaries: ρ ≡ 0 φ ≡ φmax, φmax: very large fixed positive number At exits: φ ≡ 0
63
SLIDE 66
Some References
- Introduced as a tool to design waliking facilities in
[Hughes; 2002] [Huang, Wong, Zhang, Shu, Lam; 2009]
- Used to study pedestrian flows in the domain with obstacles in
[Xia, Wong, Zhang, Shu, Lam; 2008] [Huang, Wong, Zhang, Shu, Lam; 2009]
- Numerical methods have been developed in
[Xia, Wong, Zhang, Shu, Lam; 2008] (DG on triangular meshes) [Huang, Wong, Zhang, Shu, Lam; 2009] (FD 5-order WENO)
64
SLIDE 67
Three Types of Cells: Interior, Mixed, Obstacle
65
SLIDE 68
Numerical Challenges
- We update the solution in time in the interior cells only, while the mixed
cell data required for numerical flux function evaluation are obtained using an interpolation procedure
- Interior cells: the solution is evolved using the second-order semi-discrete
central-upwind (CU) scheme
- The CU scheme, however, cannot be directly applied, since the flux
depends on the derivatives of the cost potential function φ ⇒ we take advantage of the fact that CU schemes are not based on (approximate) Riemann problem solvers and the only upwinding information required is the estimate of one-sided local speeds of propagation
- Additional source of difficulty: the flux is implicitly dependent on
the density ρ through the Eikonal equation ⇒ We apply the fast sweeping method to evolve the cost potential in the interior and mixed cells, while keeping its values large in obstacle cells
66
SLIDE 69
Numerical Examples
- Domain: [0, 100] × [0, 50]
- Upper and lower boundaries: solid walls, we add one layer of obstacle
cells, in which we set φ = φmax
- the initial density is ρ(x, y, 0) ≡ 0 and the pedestrians enter the domain
from the left, exit from the right
67
SLIDE 70
Numerical Examples
- 1. We take the solid obstacle to be the disk of radius 10 centered at (50, 20)
and the exit is on the right between y = 10 and y = 40.
- 2. The obstacle is shifted down by 8.5 so that the lower passage gets very
narrow, which is a problematic situation for the naive method, which extends the obstacle and thus almost blocks the lower passage.
- 3. We make two substantial modifications to the data from Example 1. First,
we take a larger circular obstacle of radius 15 centered at (81.5, 25). Second, we take a more narrow exit on the right being now between y = 15 and y = 35. This makes the pedestrians to get substantially blocked when the obstacle is extended by the naive method.
68
SLIDE 71
Example 1
t = 30 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 60 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 120 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 180 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
69
SLIDE 72
- Setting for the “naive” approach
70
SLIDE 73
Example 2
t = 30 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 30 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 60 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 60 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
We compare the proposed method (left) with the “naive” one (right)
71
SLIDE 74
t = 90 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 90 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 120 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 120 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
72
SLIDE 75
t = 150 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 150 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 180 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 180 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
73
SLIDE 76
Example 3
t = 60 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 60 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 120 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 120 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
We compare the proposed method (left) with the “naive” one (right)
74
SLIDE 77
t = 180 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 180 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 240 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 240 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
75
SLIDE 78
t = 300 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 300 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 360 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10 t = 360 20 40 60 80 100 10 20 30 40 50 2 4 6 8 10
76
SLIDE 79
THANK YOU!
77