Well-balanced DG scheme for Euler equations with gravity Praveen. C - - PowerPoint PPT Presentation
Well-balanced DG scheme for Euler equations with gravity Praveen. C praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065 Oberseminar Mathematische Str omungsmechanik Institut f
praveen@tifrbng.res.in
Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065
Oberseminar Mathematische Str¨
Institut f¨ ur Mathematik-der Julius-Maximilians-Universit¨ at W¨ urzburg 16 April 2015 Supported by Airbus Foundation Chair at TIFR-CAM/ICTS
1 / 57
Flow properties ρ = density, u = velocity p = pressure, E = total energy Gravitational potential Φ; force per unit volume of fluid −ρ∇Φ System of conservation laws ∂ρ ∂t + ∂ ∂x(ρu) = ∂ ∂t(ρu) + ∂ ∂x(p + ρu2) = −ρ∂Φ ∂x ∂E ∂t + ∂ ∂x(E + p)u = −ρu∂Φ ∂x
2 / 57
Perfect gas assumption p = (γ − 1)
2ρu2
γ = cp cv > 1 In compact notation ∂q ∂t + ∂f ∂x = − ρ ρu ∂Φ ∂x where q = ρ ρu E , f = ρu p + ρu2 (E + p)u
3 / 57
ue = 0
dpe dx = −ρe dΦ dx (1)
p = ρRT, R = gas constant
4 / 57
pe(x) exp Φ(x) RTe
(2) Density ρe(x) = pe(x) RTe
peρ−ν
e
= const, peT
−
ν ν−1
e
= const, ρeT
−
1 ν−1
e
= const (3) where ν > 1 is some constant. From (1) and (3), we obtain νRTe(x) ν − 1 + Φ(x) = const (4) E.g., pressure is pe(x) = C1 [C2 − Φ(x)]
ν−1 ν 5 / 57
◮ Isothermal case: Xing and Shu [3], well-balanced WENO scheme ◮ If ν = γ we are in isentropic case
h(x) + Φ(x) = const has been considered by Kappeli and Mishra [1].
◮ Desveaux et al: Relaxation schemes, general hydrostatic states ◮ PC and Klingenberg: isothermal and polytropic, ideal gas
◮ Yulong Xing: ideal gas, isothermal well-balanced 6 / 57
∂q ∂t = R(q)
R(qe) = 0
q(x, 0) = qe(x) + ε˜ q(x, 0), ε ≪ 1
∂ ˜ q ∂t = R′(qe)˜ q
7 / 57
∂qh ∂t = Rh(qh)
Rh(qh,e) = 0 = ⇒ ∂qh ∂t = 0
qh(x, t) = qh,e(x) + ε˜ qh(x, t)
8 / 57
∂ ∂t(qh,e + ε˜ qh) = Rh(qh,e + ε˜ qh) = Rh(qh,e) + εR′
h(qh,e)˜
qh
∂ ˜ qh ∂t = 1 εRh(qh,e) + R′
h(qh,e)˜
qh
∂ ˜ qh ∂t = 1 εChrqh,e + R′
h(qh,e)˜
qh
∂ ˜ qh ∂t = R′
h(qh,e)˜
qh
9 / 57
◮ Gauss-Lobatto-Legendre nodes ◮ Arbitrary quadrilateral cells in 2-D
10 / 57
Consider conservation law with source term qt + f(q)x = s(q) which has a stationary solution qe f(qe)x = s(qe) Weak formulation: Find q(t) ∈ V such that d dt (q, φ) + a(q, φ) = (s(q), φ), ∀φ ∈ V Finite element scheme: Find qh(t) ∈ Vh such that d dt (qh, φh) + a(qh, φh) = (s(qh), φh), ∀φh ∈ Vh
11 / 57
This is not true in general since we need to use quadratures. FEM with quadrature: Find qh(t) ∈ Vh such that d dt (qh, φh)h + ah(qh, φh) = (sh(qh), φh)h, ∀φh ∈ Vh Let qe,h = Πh(qe), Πh : V → Vh, interpolation or projection For the above scheme to be well-balanced, we require that ah(qe,h, φh) = (sh(qe,h), φh)h, ∀φh ∈ Vh because if qh(0) = qe,h = ⇒ qh(t) = qe,h ∀t
12 / 57
Ci = (xi− 1
2 , xi+ 1 2 ),
∆xi = xi+ 1
2 − xi− 1 2
x = ξ∆xi + xi− 1
2
(5)
Gauss-Lobatto-Legendre nodes, roots of (1 − ξ2)P ′
N(ξ)
the interpolation property ℓj(ξk) = δjk, 0 ≤ j, k ≤ N
13 / 57
φj(x) = ℓj(ξ), 0 ≤ j ≤ N
differentiation d dxφj(x) = ℓ′
j(ξ) dξ
dx = 1 ∆xi ℓ′
j(ξ)
xj = ξj∆xi + xi− 1
2 ,
0 ≤ j ≤ N
14 / 57
Consider the single conservation law with source term ∂q ∂t + ∂f ∂x = s Solution inside cell Ci is polynomial of degree N qh(x, t) =
N
qj(t)φj(x), qj(t) = qh(xj, t) Also approximate flux fh(x, t) =
N
f(qh(xj, t))φj(x) =
N
fj(t)φj(x)
15 / 57
Gauss-Lobatto-Legendre quadrature
φ(x)ψ(x)dx ≈ (φ, ψ)h = (φ, ψ)N,Ci = ∆xi
N
ωqφ(xq)ψ(xq) GLL weights ωq correspond to the reference interval [0, 1].
Semi-discrete DG: For 0 ≤ j ≤ N
d dt (qh, φj)h + (∂xfh, φj)h +[ ˆ fi+ 1
2 − fh(x−
i+ 1
2 )]φj(x−
i+ 1
2 )
−[ ˆ fi− 1
2 − fh(x+
i− 1
2 )]φj(x+
i− 1
2 ) = (sh, φj)h
(6) where ˆ fi+ 1
2 = ˆ
f(q−
i+ 1
2 , q+
i+ 1
2 ) is a numerical flux function. 16 / 57
Numerical flux is consistent: ˆ f(q, q) = f(q)
Def: Contact property
The numerical flux ˆ f is said to satisfy contact property if for any two states qL = [ρL, 0, p/(γ − 1)] and qR = [ρR, 0, p/(γ − 1)] we have ˆ f(qL, qR) = [0, p, 0]⊤
contact discontinuity.
⇒ numerical flux exactly support a stationary contact discontinuity.
17 / 57
Let ¯ Ti = temperature corresponding to the cell average value in cell Ci Rewrite the source term in the momentum equation as (Xing & Shu) s = −ρ∂Φ ∂x = ρR ¯ Ti exp Φ R ¯ Ti ∂ ∂xexp
R ¯ Ti
sh(x) = ρh(x)R ¯ Ti exp Φ(x) R ¯ Ti ∂ ∂x
N
exp
R ¯ Ti
(7) Source term in the energy equation 1 ρh (ρu)hsh
18 / 57
Well-balanced property
Let the initial condition to the DG scheme (6), (7) be obtained by interpolating the isothermal hydrostatic solution corresponding to a continuous gravitational potential Φ. Then the scheme (6), (7) preserves the initial condition under any time integration scheme. Proof: For continuous hydrostatic solution, by flux consistency ˆ fi+ 1
2 − fh(x−
i+ 1
2 ) = 0,
ˆ fi− 1
2 − fh(x+
i− 1
2 ) = 0
Above is true even if density is discontinuous, provided flux satisfies contact property. = ⇒ density and energy equations are well-balanced
19 / 57
The flux fh has the form fh(x, t) =
N
pj(t)φj(x), pj = pressure at the GLL point xj Isothermal initial condition, ¯ Ti = Te = const. The source term evaluated
20 / 57
at any GLL node xk is given by sh(xk) = ρh(xk)RTe
exp Φ(xk) RTe N
exp
RTe ∂ ∂xφj(xk) = pk exp Φ(xk) RTe N
exp
RTe ∂ ∂xφj(xk) =
N
pk exp Φ(xk) RTe
RTe ∂ ∂xφj(xk) =
N
pj exp Φ(xj) RTe
RTe ∂ ∂xφj(xk) =
N
pj ∂ ∂xφj(xk) = ∂ ∂xfh(xk)
21 / 57
Since ∂xfh(xk) = sh(xk) at all the GLL nodes xk we can conclude that (∂xfh, φj)h = (sh, φj)h , 0 ≤ j ≤ N and hence the scheme is well-balanced for the momentum equation also.
22 / 57
∂q ∂t + ∂f ∂x + ∂g ∂y = s where q is the vector of conserved variables, (f, g) is the flux vector and s is the source term, given by q = ρ ρu ρv E , f = ρu p + ρu2 ρuv (E + p)u , g = ρv ρuv p + ρv2 (E + p)v s = −ρ∂Φ
∂x
−ρ ∂Φ
∂y
−
∂x + v ∂Φ ∂y
23 / 57
Momentum equation ∇pe = −ρe∇Φ Assuming the ideal gas equation of state p = ρRT and a constant temperature T = Te = const, we get dpe = −ρedΦ = − pe RTe dΦ Integrating this equations gives the condition pe(x, y) exp Φ(x, y) RTe
(8) We will exploit the above property of the hydrostatic state to construct the well-balanced scheme.
24 / 57
A quadrilateral cell K and reference cell ˆ K = [0, 1] × [0, 1]
1 2 3 4 K x y ξ η 1 2 3 4 ˆ K TK
1-D GLL points ξr ∈ [0, 1], 0 ≤ r ≤ N
25 / 57
Tensor product of GLL points N = 1 N = 2 N = 3 Basis functions in cell K: i ↔ (r, s) ϕK
i (x, y) = ℓr(ξ)ℓs(η),
(x, y) = TK(ξ, η) Computation of ∂xfh requires derivatives of the basis functions ∂ ∂xφK
i (x, y) = ℓ′ r(ξ)ℓs(η) ∂ξ
∂x + ℓr(ξ)ℓ′
s(η)∂η
∂x
26 / 57
Consider a single conservation law with source term ∂q ∂t + ∂f ∂x + ∂g ∂x = s Solution inside cell K qh(x, y, t) =
M
qK
i (t)ϕK i (x, y),
M = (N + 1)2 Approximate fluxes inside the cell by interpolation at the same GLL nodes fh(x, y) =
M
f(qh(xi, yi))φK
i (x, y) = M
fiφK
i (x, y)
gh(x, y) =
M
g(qh(xi, yi))φK
i (x, y) = M
giφK
i (x, y)
27 / 57
Quadrature on element K using GLL points (φ, ψ)K =
N
N
φ(ξr, ξs)ψ(ξr, ξs)ωrωs|JK(ξr, ξs)| Quadrature on the faces of the cell K (φ, ψ)∂K =
(φ, ψ)e where (·, ·)e are one dimensional quadrature rules using the subset of GLL nodes located on the boundary of the cell.
28 / 57
Semi-discrete DG scheme: For 1 ≤ i ≤ M
d dt
i
i
i
+
Fh − F −
h , φK i
i
Numerical flux function ˆ Fh = ˆ F(q−
h , q+ h , n),
n = (nx, ny), unit outward normal to ∂K Trace of flux on ∂K from cell K F −
h = f− h nx + g− h ny
29 / 57
Let ¯ TK = temperature corresponding to the cell average value in cell K Rewrite source term in the x momentum equation s = −ρ∂Φ ∂x = ρR ¯ TK exp
R ¯ TK ∂ ∂x exp
Φ R ¯ TK
sh(x, y) = ρh(x, y)R ¯ TK exp Φ(x, y) R ¯ TK ∂ ∂x
M
exp
R ¯ TK
j (x, y)
(10)
30 / 57
Well-balanced property
Let the initial condition to the DG scheme (9), (10) be obtained by interpolating the isothermal hydrostatic solution corresponding to a continuous gravitational potential Φ. Then the scheme preserves the initial condition under any time integration scheme.
31 / 57
◮ If residual in cell K is zero, then dont apply limiter in that cell. 32 / 57
dt = R(t, q)
Second order accurate SSP Runge-Kutta scheme q(1) = qn + ∆tR(tn, qn) q(2) = 1 2qn + 1 2[q(1) + ∆tR(tn + ∆t, q(1))] qn+1 = q(2) Third order accurate SSP Runge-Kutta scheme q(1) = qn + ∆tR(tn, qn) q(2) = 3 4qn + 1 4[q(1) + ∆tR(tn + ∆t, q(1))] q(3) = 1 3qn + 2 3[q(2) + ∆tR(tn + ∆t/2, q(2))] qn+1 = q(3)
33 / 57
34 / 57
Potential Φ = x, Φ = sin(2πx) Initial state u = 0, ρ = p = exp(−Φ(x)) Mesh ρu ρv ρ E 25x25 1.03822e-13 6.68114e-15 2.72604e-14 9.53913e-14 50x50 1.04783e-13 5.92391e-15 2.67559e-14 9.36725e-14 100x100 1.05019e-13 5.6383e-15 2.66323e-14 9.34503e-14 200x200 1.05088e-13 5.54862e-15 2.66601e-14 9.33861e-14
Table: Well-balanced test on Cartesian mesh using Q1 and potential Φ = y
35 / 57
Mesh ρu ρv ρ E 25x25 1.04518e-13 7.29936e-15 2.7548e-14 9.64205e-14 50x50 1.04983e-13 6.03994e-15 2.69317e-14 9.43158e-14 100x100 1.05069e-13 5.68612e-15 2.69998e-14 9.39126e-14 200x200 1.05089e-13 5.69125e-15 2.68828e-14 9.462e-14
Table: Well-balanced test on Cartesian mesh using Q2 and potential Φ = x
Mesh ρu ρv ρ E 25x25 9.23424e-13 1.16432e-13 2.31405e-13 8.16645e-13 50x50 9.36459e-13 1.04921e-13 2.28315e-13 8.04602e-13 100x100 9.39613e-13 1.00384e-13 2.28001e-13 8.03005e-13 200x200 9.40422e-13 9.89098e-14 2.2792e-13 8.02653e-13
Table: Well-balanced test on Cartesian mesh using Q1 and Φ = sin(2πx)
36 / 57
Mesh ρu ρv ρ E 25x25 9.34536e-13 1.32134e-13 2.35173e-13 8.30316e-13 50x50 9.39556e-13 1.08172e-13 2.29908e-13 8.10055e-13 100x100 9.40442e-13 1.00923e-13 2.28538e-13 8.04638e-13 200x200 9.40668e-13 9.90613e-14 2.28051e-13 8.0357e-13
Table: Well-balanced test on Cartesian mesh using Q2 and Φ = sin(2πx)
37 / 57
Potential Φ(x) = x Initial condition u = 0, ρ = exp(−x), p = exp(−x) + η exp(−100(x − 1/2)2)
38 / 57
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.002 0.000 0.002 0.004 0.006 0.008 0.010
Pressure perturbation Initial 100 cells 1000 cells
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.002 0.000 0.002 0.004 0.006 0.008 0.010
Pressure perturbation Initial 50 cells 500 cells
(a) (b)
Figure: Evolution of perturbations for η = 10−2: (a) Q1 (b) Q2
39 / 57
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.00002 0.00000 0.00002 0.00004 0.00006 0.00008 0.00010
Pressure perturbation Initial 100 cells 1000 cells
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.00002 0.00000 0.00002 0.00004 0.00006 0.00008 0.00010
Pressure perturbation Initial 50 cells 500 cells
(a) (b)
Figure: Evolution of perturbations for η = 10−4: (a) Q1 (b) Q2
40 / 57
Potential Φ(x) = x and initial condition (ρ, u, p) =
x < 1
2
(0.125, 0, 0.1) x > 1
2
0.0 0.2 0.4 0.6 0.8 1.0
x
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Density Q1, 100 cells Q2, 100 cells FVM, 2000 cells
0.0 0.2 0.4 0.6 0.8 1.0
x
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Density Q1, 200 cells Q2, 200 cells FVM, 2000 cells
41 / 57
Potential Φ = x + y Hydrostatic solution is given by ρ = ρ0 exp
p0 (x + y)
p = p0 exp
p0 (x + y)
u = v = 0 ρu ρv ρ E Q1, 25 × 25 9.85926e-14 9.85855e-14 5.32357e-14 1.55361e-13 Q1, 50 × 50 9.94493e-14 9.94451e-14 5.37084e-14 1.56669e-13 Q1, 100 × 100 9.96481e-14 9.96474e-14 5.38404e-14 1.57062e-13 Q2, 25 × 25 9.9256e-14 9.92682e-14 5.39863e-14 1.57435e-13 Q2, 50 × 50 9.961e-14 9.96538e-14 5.41091e-14 1.57521e-13 Q2, 100 × 100 9.95889e-14 9.97907e-14 5.43145e-14 1.57728e-13
42 / 57
p = p0 exp
p0 (x + y)
p0 [(x − 0.3)2 + (y − 0.3)2]
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
(a) (b)
Figure: Pressure perturbation on 50 × 50 mesh at time t = 0.15 (a) Q1 (b) Q2. Showing 20 contours lines between −0.0002 to +0.0002
44 / 57
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
(a) (b)
Figure: Pressure perturbation on 200 × 200 mesh at time t = 0.15 (a) Q1 (b) Q2. Showing 20 contours lines between −0.0002 to +0.0002
45 / 57
Potential Φ(x, y) = y Temperature is discontinuous T =
y < 0 Tu y > 0 Hydrostatic pressure and density p =
y ≤ 0 p0e−y/RTu, y > 0 , ρ =
y < 0 p/RTu, y > 0 Density is discontinuous at y = 0.
46 / 57
Case 1, Tl = 1, Tu = 2 Q1, 25x100 5.13108e-13 1.10971e-13 2.56836e-13 8.91784e-13 Q1, 50x200 5.15744e-13 1.12234e-13 2.84309e-13 8.97389e-13 Q2, 25x100 5.15726e-13 1.1341e-13 3.22713e-13 9.01183e-13 Q2, 50x200 5.16397e-13 1.13707e-13 3.93623e-13 9.02503e-13
Table: Well-balanced test for Rayleigh-Taylor problem
Case 2, Tl = 2, Tu = 1 Q1, 25x100 3.487e-13 1.0989e-13 1.98949e-13 7.42265e-13 Q1, 50x200 3.50153e-13 1.1121e-13 2.34991e-13 7.4747e-13 Q2, 25x100 3.50149e-13 1.12159e-13 2.80645e-13 7.5104e-13 Q2, 50x200 3.50548e-13 1.12533e-13 3.63806e-13 7.52151e-13
Table: Well-balanced test for Rayleigh-Taylor problem
47 / 57
Exact solution of Euler equations with gravity given by ρ = 1 + 0.2 sin(π(x + y − t(u0 + v0))), u = u0, v = v0 p = p0 + t(u0 + v0) − x − y + 0.2 cos(π(x + y − t(u0 + v0)))/π Compute solution error in L2 norm at time t = 0.1
1/h ρu ρv ρ E Error Rate Error Rate Error Rate Error Rate 50 0.00134154 – 0.00134154 – 0.0012837 – 0.00161287 – 100 0.000335446 1.99 0.000335446 1.99 0.00032044 2.00 0.000411141 1.97 200 8.35627e-05 2.00 8.35627e-05 2.00 7.97842e-05 2.00 0.00010335 1.99 400 2.08348e-05 2.00 2.08348e-05 2.00 1.98754e-05 2.00 2.58109e-05 2.00
Table: Convergence of error for degree N = 1
48 / 57
1/h ρu ρv ρ E Error Rate Error Rate Error Rate Error Rate 25 7.7019e-05 – 7.7019e-05 – 7.80868e-05 – 9.32865e-05 – 50 9.68863e-06 2.99 9.68863e-06 2.99 9.76471e-06 2.99 1.16849e-05 2.99 100 1.21506e-06 2.99 1.21506e-06 2.99 1.22031e-06 3.00 1.46256e-06 2.99 200 1.52134e-07 2.99 1.52134e-07 2.99 1.52503e-07 3.00 1.8247e-07 3.00
Table: Convergence of error for degree N = 2
49 / 57
30 × 30 cells 956 cells
50 / 57
Radial potential Φ(x, y) = r =
Hydrostatic solution ρ = p = exp(−r) ρu ρv ρ E Q1, 30x30 6.08113e-13 6.06081e-13 2.11611e-14 4.13127e-14 Q1, 50x50 5.45634e-13 5.44973e-13 1.29319e-14 2.66376e-14 Q1, 100x100 5.4918e-13 5.49012e-13 9.63433e-15 2.11463e-14 Q2, 30x30 6.29776e-13 6.27278e-13 1.9777e-14 5.08689e-14 Q2, 50x50 5.5376e-13 5.52903e-13 1.5635e-14 3.85134e-14 Q2, 100x100 5.51645e-13 5.5142e-13 2.173e-14 5.25578e-14
Table: Well balanced test for radial Rayleigh-Taylor problem on Cartesian mesh
51 / 57
ρu ρv ρ E Q1, 956 3.03033e-16 3.23738e-16 6.66245e-16 2.11449e-16 Q1, 2037 5.20653e-16 5.10865e-16 9.82565e-16 4.06402e-16 Q1, 10710 2.00037e-15 1.57078e-15 2.21284e-15 4.66515e-15 Q2, 956 2.69429e-14 3.31591e-14 4.50499e-14 1.61455e-13 Q2, 2037 7.68549e-14 1.1834e-13 1.01632e-13 3.84886e-13 Q2, 10710 2.92633e-13 2.44323e-13 2.9449e-13 4.10069e-13
Table: Well balanced test for radial Rayleigh-Taylor problem on unstructured mesh
52 / 57
Initial pressure and density p =
r ≤ r0 e− r
α +r0 (1−α) α
r > r0 , ρ =
r ≤ ri
1 αe− r
α +r0 (1−α) α
r > ri where ri = r0(1 + η cos(kθ)), α = exp(−r0)/(exp(−r0) + ∆ρ) The density jumps by an amount ∆ρ > 0 at the interface defined by r = ri whereas the pressure is continuous. Following [2], we take ∆ρ = 0.1, η = 0.02, k = 20
53 / 57
Cartesian mesh of 240 × 240 and Q1 basis
54 / 57
◮ General equation of state ◮ Hydrostatic solution may not be known explicitly 55 / 57
56 / 57
[1] R. K¨ appeli and S. Mishra, Well-balanced schemes for the Euler equations with gravitation, J. Comput. Phys., 259 (2014),
[2] Randall. J. LeVeque and Derek. S. Bale, Wave propagation methods for conservation laws with source terms, in Hyperbolic Problems: Theory, Numerics, Applications, Rolf Jeltsch and Michael Fey, eds., vol. 130 of International Series of Numerical Mathematics, Birkh¨ auser Basel, 1999, pp. 609–618. [3] Yulong Xing and Chi-Wang Shu, High order well-balanced WENO scheme for the gas dynamics equations under gravitational fields, J. Sci. Comput., 54 (2013), pp. 645–662.
57 / 57