Mapped grid methods

Wave propagation methods for hyperbolic problems on mapped

grids

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

Outline

Review a high-resolution wave propagation method for solving hyperbolic problems on mapped grids (which is basic integration scheme implemented in CLAWPACK) Describe an interpolation-type adaptive moving mesh approach that is an easy generalization of the above method for improving numerical resolution

Show sample results

Mapped grid methods

To begin with, we consider system of conservation laws

∂_tq + ∇ · f (q) = 0 ⁽¹⁾ with suitable initial condition for q in N_d ≥ 1 spatial domain Here q ∈ lR^m denotes vector of m conserved quantities, f = (f¹, f², . . . , f_Nd) ∈ lR^m×Nd denotes flux matrix; f_j ∈ lR^m It is known integral form of (1) over any control volume C is

d dt

Z

C

q dx = − Z

∂C

f(q) · n ds,

where n is outward-pointing normal vector at boundary ∂C

Mapped grid methods

Using above integral conservation law, a finite volume method on a control volume C can be written as

Qⁿ⁺¹ = Qⁿ − ∆t M(C)

N^s

X

j=1

h_jF˘_j,

where M(C) is measure (area in 2D or volume in 3D) of C, N_s is number of sides,

h_j is length of j-th side (in 2D) or area of cell edge (in 3D) measured in physical space,

F˘_j is numerical approximation to normal flux in average across j-th side of grid cell

Mapped grid methods

In the following we assume that our mapped grids are logically rectangular, & will restrict our consideration to two-dimensional case N_d = 2 as illustrated below

i − 1 i − 1

i

i j

j

j + 1 j + 1

Cˆ_ij C_ij

ξ₁ ξ₂

mapping

∆ξ₁

∆ξ₂

computational grid physical grid

←−

x₁ = x₁(ξ₁, ξ₂) x₂ = x₂(ξ₁, ξ₂) x₁

x₂

Then on a curvilinear grid, a finite volume method takes

∆t   ∆t  

Mapped grid methods

On a curvilinear grid, a finite volume method takes Qⁿ⁺¹_ij = Qⁿ_ij− ∆t

κ_ij∆ξ1

F_i+¹ 1

2,j − F_i−¹ 1 2,j

− ∆t κ_ij∆ξ2

F_i,j+² 1 2

− F_i,j−² 1 2



where ∆ξ1, ∆ξ2 denote equidistant discretization of computational domain,

κ_ij = M(C_ij)/∆ξ¹∆ξ² is area ratio between area of grid cell in physical space & area of a computational grid,

F_i−¹ ₁

2,j = γ_i−¹

2,jF˘_i−¹

2,j, F_i,j−² ₁

2

= γ_i,j−¹

2

F˘_i,j−¹

2 are fluxes per unit length in computational space with γ_i−¹

2,j = h_i−¹

2,j/∆ξ1 &

γ_i,j−¹

2 = h_i,j−¹

2 /∆ξ2 representing length ratios

Mapped grid methods

In this setup, first order wave propagation method is a Godunov-type finite volume method that takes form

Qⁿ⁺¹_ij = Qⁿ_ij− ∆t κ_ij∆ξ₁

A⁺₁ ∆Q_i−¹

2,j + A⁻₁ ∆Q_i+¹

2,j

−

∆t κ_ij∆ξ₂

A⁺₂ ∆Q_i,j−¹

2 + A⁻₂ ∆Q_i,j+¹

2



with right-, left-, up-, and down-moving fluctuations

A⁺₁ ∆Q_i−¹

2,j, A⁻₁ ∆Q_i+¹

2,j, A⁺₂ ∆Q_i,j−¹

2, & A⁻₂ ∆Q_i,j+¹

2 that are entering into grid cell

To determine these fluctuations, we need to solve

one-dimensional Riemann problems normal to cell edges

Computing fluctuations

Let ⁿ_i−¹

2,j = (n¹, n²) & ^t_i−¹

2,j = (t¹, t²) be normalized normal and tangential vectors to cell edge (i − ¹₂, j) between cells (i − 1, j) & (i, j). Then rotation matrix which rotates velocity components (^e.g., for shallow water or acoustic equations) has form

R_i−¹

2,j =







1 0 0

0 n¹ n² 0 t¹ t²







Fluctuations A^±₁ ∆Q_i−¹

2,j can be calculated by steps:

1. Determine Q^˘_L & Q^˘_R by rotating velocity components,

Q˘_L = R_i−¹

2,jQ_i−1,j, Q˘_R = R_i−¹

2,jQ_i,j

Computing fluctuations

2. Solve Riemann problem for 1D conservation law

∂_tq + ∇ · f₁ = 0

with initial data

q(ξ₁, 0) =





Q˘_L if x < x_i−¹

2

Q˘_R if x > x_i−¹

2

This results in waves W^˘_i−^1,k₁

2,j that are moving with speeds

λ˘^1,k_i−1

2,j, k = 1, . . . , N_w, N_w is the number of waves. Here jumps are decomposed into waves as Q^˘_R − ˘Q_L = PNw

k=1 W˘_i−^1,k₁

2,j

& solution at cell edge can be expressed

X

Computing fluctuations

3. Define scaled speeds

λ^1,k_i−1

2,j = γ_i−¹

2,jλ˘^1,k_i−1

2,j

and rotate waves back to Cartesian coordinates by

W_i−^1,k₁

2,j = R^T_i−1

2,jW˘_i−^1,k₁

2,j, k = 1, . . . , N_w

4. Determine left- and right-moving fluctuations in the form

A^±₁ ∆Q_i−¹

2,j =

Nw

X

k=1

λ^1,k_i−1

2,j

^±

W_i−^1,k₁

2,j

5. In an analogous way up and down-moving fluctuations at cell edge (i, j − ¹₂) can be calculated

High resolution corrections

The speeds and limited versions of waves are used to

calculate second order correction terms. These terms are added to the method in flux difference form as

Qⁿ⁺¹_ij := Qⁿ⁺¹_ij − 1 κ_ij

∆t

∆ξ₁

Fe_i+¹ 1

2,j − eF_i−¹ 1 2,j

 − 1 κ_ij

∆t

∆ξ₂

Fe_i,j+² 1

2 − eF_i,j−² 1 2



For example, at cell edge (i − ¹₂, j) correction flux takes

Fe_i−¹ 1

2,j = 1 2

Nw

X

k=1

λ^1,k_i−¹

2,j

 1 − ∆t κ_i−¹

2,j∆ξ₁ 

λ^1,k_i−¹

2,j

!

Wf_i−^1,k₁

2,j

where κ_i−¹

2,j = (κ_i−1,j + κ_ij)/2. To aviod oscillations near discontinuities, a wave limiter is applied leading to limited

f

High resolution corrections

To ensure second order accuracy & also improve stability, a transverse wave propagation is included in the algorithm

that left- & right-going fluctuations A^±₁ ∆Q_i−¹

2,j are each split into two transverse fluctuations: up- & down-going

A^±₂ A⁺₁ ∆Q_i−¹

2,j & A^±₂ A⁻₁ ∆Q_i−¹

2,j

This wave propagation method can be shown to be conservative & stable under a variant of CFL

(Courant-Friedrichs-Lewy) condition of form

ν = ∆t max

i,j,k



 

λ^1,k_i−¹

2,j

  J_ip,j∆ξ₁ ,

λ^2,k_i,j−¹

2

  J_i,jp∆ξ₂



 ≤ 1,

where i_p = i if λ¹_i−^,k₁

2,j > 0 & i − 1 if λ¹_i−^,k₁

2,j < 0

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γπ² exp (1 − r²)

1/(γ−1)

, p = ρ^γ,

u₁ = 1 − 5

2π exp ((1 − r²)/2) (x₂ − 5) , u₂ = 1 + 5

2π exp ((1 − r²)/2) (x₁ − 5) ,

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

(x1 − 5)² + (x2 − 5)²

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

kE_zk1,∞ = kzcomput ^{− z}exact^k1,∞ denotes discrete 1- or maximum-norm error for state variable z

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

Accuracy results in 2D: Grid 1

N E₁(ρ) Order E₁(u₁) Order E₁(u₂) Order E₁(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^∞(u₁) Order E^∞(u₂) 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

Accuracy results in 2D: Grid 2

N E₁(ρ) Order E₁(u₁) Order E₁(u₂) Order E₁(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^∞(u₁) Order E^∞(u₂) 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

Accuracy results in 2D: Grid 3

N E₁(ρ) Order E₁(u₁) Order E₁(u₂) Order E₁(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^∞(u₁) Order E^∞(u₂) 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

Accuracy test in 3D

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

Take smooth radially-symmetric flow with the flow condition that is at rest initially with density

ρ(r) = 1 + exp(−30(r − 1)²)/10 & pressure p(r) = ρ^γ Grids used for smooth radially-symmetric flow test

1 1 2

2 0 1 2

1 1 2

2 0 1 2

1 1 2

2 0 1 2

Grid 1 Grid 2 Grid 3

Accuracy results in 3D: Grid 1

N E₁(ρ) Order E₁(|~u|) Order E₁(p) Order 20 7.227 · 10⁻³ 8.920 · 10⁻³ 1.019 · 10⁻²

40 2.418 · 10⁻³ 1.58 2.558 · 10⁻³ 1.80 3.415 · 10⁻³ 1.58 80 6.356 · 10⁻⁴ 1.93 6.754 · 10⁻⁴ 1.92 8.980 · 10⁻⁴ 1.93 160 1.616 · 10⁻⁴ 1.98 1.718 · 10⁻⁴ 1.97 2.282 · 10⁻⁴ 1.98

N E^∞(ρ) Order E^∞(|~u|) Order E^∞(p) Order 20 1.096 · 10⁻² 1.200 · 10⁻² 1.569 · 10⁻²

40 4.085 · 10⁻³ 1.42 4.381 · 10⁻³ 1.45 5.848 · 10⁻³ 1.42 80 1.235 · 10⁻³ 1.73 1.263 · 10⁻³ 1.79 1.765 · 10⁻³ 1.73 160 3.517 · 10⁻⁴ 1.81 3.349 · 10⁻⁴ 1.91 5.030 · 10⁻⁴ 1.81

Accuracy results in 3D: Grid 2

N E₁(ρ) Order E₁(|~u|) Order E₁(p) Order 20 7.227 · 10⁻³ 8.920 · 10⁻³ 1.019 · 10⁻²

40 2.418 · 10⁻³ 1.58 2.558 · 10⁻³ 1.80 3.415 · 10⁻³ 1.58 80 6.356 · 10⁻⁴ 1.93 6.754 · 10⁻⁴ 1.92 8.980 · 10⁻⁴ 1.93 160 1.616 · 10⁻⁴ 1.98 1.718 · 10⁻⁴ 1.97 2.282 · 10⁻⁴ 1.98

N E^∞(ρ) Order E^∞(|~u|) Order E^∞(p) Order 20 7.227 · 10⁻³ 8.920 · 10⁻³ 1.019 · 10⁻²

40 2.418 · 10⁻³ 1.58 2.558 · 10⁻³ 1.80 3.415 · 10⁻³ 1.58 80 6.356 · 10⁻⁴ 1.93 6.754 · 10⁻⁴ 1.92 8.980 · 10⁻⁴ 1.93 160 1.616 · 10⁻⁴ 1.98 1.718 · 10⁻⁴ 1.97 2.282 · 10⁻⁴ 1.98

Accuracy results in 3D: Grid 3

N E₁(ρ) Order E₁(|~u|) Order E₁(p) Order 20 1.290 · 10⁻² 1.641 · 10⁻² 1.816 · 10⁻²

40 4.694 · 10⁻³ 1.46 4.999 · 10⁻³ 1.71 6.623 · 10⁻³ 1.46 80 1.257 · 10⁻³ 1.90 1.379 · 10⁻³ 1.86 1.774 · 10⁻³ 1.90 160 3.209 · 10⁻⁴ 1.97 3.546 · 10⁻⁴ 1.96 4.527 · 10⁻⁴ 1.97

N E^∞(ρ) Order E^∞(|~u|) Order E^∞(p) Order 20 1.632 · 10⁻² 1.984 · 10⁻² 2.316 · 10⁻²

40 5.819 · 10⁻³ 1.49 6.745 · 10⁻³ 1.56 8.307 · 10⁻³ 1.48 80 1.823 · 10⁻³ 1.67 4.290 · 10⁻³ 0.65 2.710 · 10⁻³ 1.67 160 5.053 · 10⁻⁴ 1.85 3.271 · 10⁻³ 0.39 7.237 · 10⁻⁴ 1.85

Extension to moving mesh

One simple way to extend mapped grid method described above to solution adaptive moving grid method is to take approach proposed by

H. Tang & T. Tang, Adaptive mesh methods for one- and two-dimensional hyperbolic conservation laws, SIAM J.

Numer. Anal., 2003

In each time step, this moving mesh method consists of three basic steps:

1. Mesh redistribution

2. Conservative interpolation of solution state 3. Solution update on a fixed mapped grid

Mesh redistribution scheme

Winslow’s approach (1981)

Solve ∇ · (D∇ξ_j) = 0, j = 1, . . . , N_d

for ξ(x). Coefficient D is a positive definite matrix which may depend on solution gradient

Variational approach (Tang & many others)

Solve ∇_ξ · (D∇_ξx_j) = 0, j = 1, . . . , N_d for x(ξ) that minimizes the “energy” functional

E(x(ξ)) = 1 2

Z

Ω Nd

X

j=1

∇^T_ξ D∇x_jdξ

Lagrangian (ALE)-type approach ( , CAVEAT code)

Mesh redistribution

Dashed lines represent initial mesh & solid lines represent new mesh after a redistribution step

0 0.5 1 1.5 2 2.5 3 3.5 4

0 0.5 1 1.5 2 2.5 3 3.5 4

Conservative interpolation

Numerical solutions need to be updated conservatively, ^i.e.

XM C^k+1

Q^k+1 = X

M C^k
Q^k

after each mesh redistribution iterate k. This can be done Finite-volume approach (Tang & Tang, SIAM 03)

M C^k+1

Q^k+1 = M C^k

Q^k −

Ns

X

j=1

h_jG˘_j, G = ( ˙x · n)Q˘

Geometric approach (Shyue 2010 & others)

"

X

S

M C_p^k+1 ∩ S_p^k
#

Q^k+1_C = X

S

M C_p^k+1 ∩ S_p^k
Q^k_S

Conservative interpolation

Geometric way to perform conservative interpolation

0 0.5 1 1.5 2 2.5 3 3.5 4

0 0.5 1 1.5 2 2.5 3 3.5 4

Interpolation-free moving mesh

One way to dervise an interpolation-free moving mesh

method is to consider coordinate change of equations via (x, t) 7→ (ξ, t). In this case, the transformed conservation law reads

∂_tq˜ + ∇_ξ · ˜f = G ⁽²⁾ with

˜

q = Jq, J = det (∂ξ/∂x)⁻¹ f˜_j = J (q ∂_tξ_j + ∇ξ_j · f )

G = q [∂_tJ + ∇_ξ · (J∂_tξ_j)] +

XN j=1

f_j∇_ξ · J∂_xjξ_k

= 0 (if GCL & SCL are satisfied)

We may then devise numerical method to solve (2) (Shyue:

2 D Riemann Problem

4-shock wave pattern initially

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1 



 ρ u¹ u2

p





 =





 1.1

0 0 1.1













0.507 0.893

0 0.35













1.1 0.893 0.893

1.1













0.507 0 0.893

0.35







2 D Riemann Problem

Mesh redistribution scheme: variational approach Grid system

time=0.25

2 D Riemann Problem

Mesh redistribution scheme: variational approach Density contours

t=0.25

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Sedov problem

Mesh redistribution scheme: Lagrangian approach 30 × 30 mesh grid

0.2 0.4 0.6 0.8 1 1.2

y

Grid system

1 2 3 4 5 6

Density

Lagrangian Eulerian

Sedov problem

Mesh redistribution scheme: variational approach

“Monitor” function D based on density gradient

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2

y

Grid system

0 0.2 0.4 0.6 0.8 1

0 1 2 3 4 5 6

Density

Variational Eulerian

Shock waves over circular array

A Mach 1.42 shock wave in water over a circular array

t=0

shock

Shock waves over circular array

Grid system

−3 −2 −1 0 1 2 3 4 5

−4

−3

−2

−1 0 1 2 3 4

time=0

Shock waves over circular array

Contours for density

t=0.0024

Moving cylindrical vessel

Impose uniform flow velocity (u1, u2) = (−1, 0) (^i.e., in the frame of vessel moving with speed one in x1-direction)

air helium interface

Moving cylindrical vessel

Solution at time t = 0.25

Moving cylindrical vessel

Solution at time t = 0.5

Moving cylindrical vessel

Solution at time t = 0.75

Moving cylindrical vessel

Solution at time t = 1

Cavitation test: 2D steady state

Uniform gas-liquid mixture with speed 600m /s over a circular region

Cavitation test: 2D steady state

Pseudo colors of volume fraction

Formation of cavitation zone (Onset–shock induced, diffusion, . . . ?)

Cavitation test: 2D steady state

Pseudo colors of pressure

Smooth transition across liquid-gas phase boundary

Cavitation test: 2D steady state

Pseudo colors of volume fraction: 2 circular case

Convergence of solution as the mesh is refined ?

t=0.05

Thank You

Compressible Multiphase Flow

Homogeneous equilibrium pressure & velocity across material interfaces

Volume-fraction based model equations (Shyue JCP

’98, Allaire ^{et al.} JCP ’02)

∂

∂t (α_iρ_i) + 1 J

Nd

X

j=1

∂

∂ξ_j (α_iρ_iU_j) = 0, i = 1, 2, . . . , m_f

∂

∂t (ρu_i) + 1 J

Nd

X

j=1

∂

∂ξ_j (ρu_iU_j + pJ_ji) = 0, i = 1, 2, . . . , N_d,

∂E

∂t + 1 J

Nd

X

j=1

∂

∂ξ_j (EU_j + pU_j) = 0,

∂α_i

∂t + 1 J

Nd

XU_j ∂α_i

∂ξ_j = 0, i = 1, 2, . . . , m_f − 1

Barotropic two-fluid flow

Two-phase flow model

∂_t (α1ρ1) + ∇ · (α1ρ1u) = 0,

∂_t (α2ρ2) + ∇ · (α2ρ2u) = 0,

∂_t (ρu) + ∇ · (ρu ⊗ u + pδ) = 0 α1ρ1

ρ1(p) + α2ρ2

ρ2(p) = 1 (Iterative solve for p)

Barotropic two-fluid flow

A relaxation model of Saurel ^{et al.} (JCP ’090

∂

∂t (α₁ρ₁) +

XN j=1

∂

∂x_j (α₁ρ₁u_j) = 0,

∂

∂t (α₂ρ₂) +

XN j=1

∂

∂x_j (α₂ρ₂u_j) = 0,

∂

∂t (ρu_i) +

XN j=1

∂

∂x_j (ρu_iu_j + pδ_ij) = 0, i = 1, . . . , N

∂α₁

∂t +

XN j=1

u_j ∂α₁

∂x_j = 1

µ (p₁ (ρ₁) − p₂ (ρ₂))

Each phasic pressure p_ι is a function of density only Mixture pressure p satisfies p = α₁p₁ + α₂p₂, µparameter

Mixture speed of sound

Wood formula (stiffness in c¯ vs. α)

1

ρ¯c² = α ρwc²w

+ 1 − α ρgc²g

ρw = 10³kg/m³, cw = 1449.4m/s, ρg = 1.0kg/m³, cg = 374.2m/s

500 1000 1500

¯c