# Predictor-Corrector LU-SGS Discontinuous Galerkin Finite Element Method for Conservation Laws

Authors: Xinrong Ma; Sanyang Liu; Gongnan Xie
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Efficient implicit predictor-corrector LU-SGS discontinuous Galerkin (DG) approach for compressible Euler equations on unstructured grids is investigated by adding the error compensation of high-order  term. The original LU-SGS and GMRES schemes for DG method are discussed. Van Albada limiter is employed to make the scheme monotone. The numerical experiments performed for the transonic inviscid flows around NACA0012 airfoil, RAE2822 airfoil, and ONERA M6 wing indicate that the present algorithm has the advantages of low storage requirements and high convergence acceleration. The computational efficiency is close to that of GMRES scheme,  nearly 2.1 times greater than that of LU-SGS scheme on unstructured grids for 2D cases, and almost 5.5 times greater than that of RK4 on unstructured grids for 3D cases.


---

## Body

## 1. Introduction

High-order discontinuous Galerkin (DG) finite element methods were developed based on weighted residual theory; they maintain advantages of both the traditional high resolution finite difference method and the finite volume method while overcoming their shortcomings. Indeed, the DG method can be considered as a mixture of classic finite element method (FEM) and finite volume method (FVM), which is a better solution strategy for solving problems in the presence of strong shocks and discontinuities because the solution across each element can be discontinuous. DG methods can easily deal with complex boundary-value problem and flexibly handle discontinuity, which have a low requirement of the regularity of grids. And high accuracy can be achieved by selecting appropriate basis functions by improving the order of the piecewise interpolation polynomials functions. In addition, the methods are highly parallelizable as each element is independent and the interelement communications are minimal. And they have several useful mathematical properties.While DG method was introduced by Reed and Hill [1] for solving the neutron transport equation back in 1973, nowadays, DG methods have been widely used in the computational fluid dynamics, computational aeroacoustics, and computational electromagnetics. See [2–25].In recent years, significant progress has been made in developing numerical algorithms for solving the compressible flow problems. Many numerical methods are based on the semidiscrete approach: DG methods are used for the spatial discretization, rendering the original partial differential equations (PDE) into a system of ordinary differential equations (ODE) in time. Usually, for time-dependent problems, DG methods have been used in conjunction with explicit high-order accurate time-integration methods, such as nonlinear stable Runge-Kutta DG methods in the literatures [7–13]. In general, explicit schemes are easy to implement and parallelize and require only limited memory storage. Such methods are well suited for problems with similar spatial and temporal scales, while being notoriously time-consuming and inefficient for problems with disparate temporal and spatial scales, such as low reduced frequency phenomena and steady-state problems. As a consequence, implicit time-integration strategies should be developed exclusively in order to avoid the stability restrictions of explicit methods, which are unconditionally stable; for details, refer to [5, 6, 26]. Implicit solvers, which do allow large time steps, are widely used in the computational fluid dynamics community for the steady solution of nonlinear conservation laws in [27]. The Newton-Krylov-Schwarz method has recently emerged as a promising technique for the parallel implicit solution of large-scale aerodynamics problems in [28], which is specially well suited for the discontinuous spectral Galerkin method, since each subdomain can be treated separately.The lower-upper symmetric Gauss-Seidel (LU-SGS) time-marching scheme, which was originally given for structured grids, has been established in [29] and has been applied to tetrahedral/prism unstructured grids. Another attractive implicit scheme is the generalized minimum residual scheme (GMRES), which was introduced by Saad and Schultz [2] firstly. Then, Bassi and Rebay showed the efficiency of GMRES method in [14] and used a simple block Jacobi preconditioner for the implicit solution of the compressible Navier-Stokes equations. Xia and Luo [3] presented a matrix-free GMRES algorithm with an LU-SGS preconditioner reconstructed discontinuous Galerkin method on tetrahedron grids for compressible flow problems. Then, they proposed implicit reconstructed discontinuous Galerkin (IRDG) method based on the automatic differentiation technique [4].In this work, we focus exclusively on a predictor-corrector LU-SGS (PCLU-SGS) strategy for discontinuous Galerkin method in conjunction with Van Albada limiter [30] to solve the Euler equations on unstructured grids. The governing equations are listed in Section 2. DG method is presented and the limiter and numerical flux are described in detail in Section 3. The implicit time-marching procedures are given including TVD Runge-Kutta, LU-SGS, and GMRES schemes, and the PCLU-SGS scheme is established in Section 4. The numerical experiments are discussed in Section 5. Conclusion is given in Section 6.


## 2. Governing Equations

The unsteady, compressible inviscid 2D Euler equations can be expressed in the following conservative form:(1)
∂
U
x
,
t
∂
t
+
∂
F
U
x
,
t
∂
x
+
∂
G
U
x
,
t
∂
y
=
0
i
n
Ω
×
0
,
T
,where T
>
0 is the length of time interval and Ω is a two-dimensional bounded domain. The conservative state vector U and the inviscid flux component vectors F
(
U
) and G
(
U
) are defined by(2)
U
=
ρ
ρ
u
ρ
v
ρ
e
,
F
U
=
ρ
u
ρ
u
2
+
p
ρ
u
v
ρ
e
+
p
u
,
G
U
=
ρ
v
ρ
u
v
ρ
v
2
ρ
e
+
p
v
,where the notations ρ, p, and e denote the density, pressure, and specific total energy per unit mass of the fluid, respectively. u and v are the velocity components of the flow in the x and y coordinate directions. This set of equations is completed by the perfect gas equation of state given as follows:(3)
p
=
γ
-
1
ρ
e
-
1
2
u
2
+
v
2
,where γ is defined as the ratio of specific heat of the fluid (γ
=
1.4 for air). Then the equations are applied with the initial and boundary conditions denoted by (4), where ∂
Ω represents the boundary of domain Ω. Consider the following:(4)
U
x
,
0
=
U
0
x
,
B
U
=
0
o
n
∂
Ω
×
0
,
T
.


## 3. Discontinuous Galerkin Finite Element Method

### 3.1. DG Spatial Discretization

The computational domainΩ is partitioned into an ensemble of nonoverlapping elements, triangles in 2D; that is, Ω
=
e
1
∪
e
2
∪
⋯
∪
e
n, where n denotes the number of elements in the domain. We consider possible choices of the piecewise basis functions and then obtain the following weak formulation of (5) by multiplying a test function Φ and integrating by parts over the e
m:(5)
d
d
t
∫
e
m
Φ
h
U
h
d
V
-
∫
e
m
∂
Φ
h
∂
x
F
U
h
+
∂
Φ
h
∂
y
G
U
h
d
V
+
∫
∂
e
m
Φ
h
F
U
h
n
x
+
G
U
h
n
y
d
S
=
0
∀
Φ
,where U
h and Φ
h represent the finite element approximations of the analytical solution U and the test function Φ, respectively. n
=
(
n
x
,
n
y
) is the unit normal vector of outward to the boundary. Let the approximate solution and test function U
h and Φ
h be expressed as(6)
U
h
=
∑
j
=
1
N
U
j
t
φ
j
p
x
,
Φ
h
=
∑
j
=
1
N
Φ
j
φ
j
p
x
,where φ
j
p
(
x
) is the shape function of the polynomials of degree p. Equation (5) must be satisfied for any test function Φ
h, so by substituting (6) to (5), we obtain the following system of N equations:(7)
d
U
j
t
d
t
∫
e
m
φ
i
p
φ
j
p
d
V
-
∫
e
m
∂
φ
i
p
∂
x
F
∑
j
=
1
N
U
j
φ
j
p
+
∂
φ
i
p
∂
y
G
∑
j
=
1
N
U
j
φ
j
p
d
V
+
∫
∂
e
m
φ
i
p
F
∑
j
=
1
N
U
j
φ
j
p
n
x
+
G
∑
j
=
1
N
U
j
φ
j
p
n
y
d
S
=
0
i
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
i
i
i
f
f
f
f
i
i
=
1,2
,
…
,
N
.The interface flux function F
(
U
h
)
n
x
+
G
(
U
h
)
n
y can be treated as a numerical Riemann flux function H
(
U
h
L
,
U
h
R
,
n
), where U
h
L and U
h
R represent the internal element interface solution and neighboring element interface solution, respectively. In the present work, the Roe, LLF, and HLLC approximate Riemann solvers are employed. The domain and boundary integrals in (7) are calculated by use of 2
p and 2
p
+
1 order accurate Gauss quadrature formulas [31, 32] with a number of quadrature points corresponding to the degree of interpolating polynomials.By grouping together all the elemental time-dependent and spatial contributions, (7) can be written as a system of ordinary differential equations:(8)
M
d
U
d
t
=
R
U
,where the mass matrix M has identical diagonal blocks M
i
j, U is the global vector of the degrees of freedom, and R
(
U
) represents the steady state residual vector. As a result, the inverse of the mass matrix M can be easily computed, especially, using the orthogonal basis functions, and stored in advance due to the fact that it remains unchanged during the process.In the present paper, we explore the orthogonal basis functions through Gram-Schmidt orthogonalization method, and high accuracy can be achieved by improving the order of the piecewise interpolation polynomials functions. For 2D problems,φ
0
=
1, φ
1
=
ξ, φ
2
=
η, φ
3
=
ξ
2
,
…, and, for 3D problems, φ
0
=
1, φ
1
=
ξ, φ
2
=
η, φ
3
=
ζ
,
….


### 3.2. Flux Functions

The numerical flux function can be evaluated using any upwind flux functions. This is exactly similar to FVM because discontinuities can be allowed across the interface. Therefore, approximate Riemann solvers can be used to compute the flux function. In the present work, three flux functions have been employed including HLLC Riemann numerical solver by Toro [33] which has easier and lower computational cost in comparison with many other available Riemann solvers, such as local Lax-Friedrich (LLF) scheme. HLLC flux function not only maintains the advantages of the HLL solver but also resolves isolated contact discontinuities exactly, which has been extended in conjunction with time-derivative preconditioning to compute flow problems at all speeds. HLL flux can be expressed as(9)
F
H
L
L
=
F
L
i
f
S
L
>
0
F
L
R
∗
i
f
S
L
≤
0
≤
S
R
F
R
i
f
S
R
<
0
,where S
L and S
R represent the fastest wave speed for the left and right states, respectively. F
L
R
∗ is written as(10)
F
L
R
∗
=
S
R
F
L
-
S
L
F
R
+
S
L
S
R
U
R
-
U
L
S
R
-
S
L
.The HLLC flux is a modification of HLL flux, which can be written as(11)
F
H
L
L
C
=
F
L
0
<
S
L
F
L
∗
=
F
L
+
S
L
U
L
∗
-
U
L
S
L
≤
0
,
S
M
>
0
F
R
∗
=
F
R
+
S
R
U
R
∗
-
U
R
S
M
≤
0
,
S
R
≥
0
F
R
0
>
S
R
,where S
M is constant between the two acoustic waves:(12)
S
M
=
p
R
-
p
L
+
ρ
L
q
L
S
L
-
q
L
-
ρ
R
q
R
S
R
-
q
R
ρ
L
S
L
-
q
L
-
ρ
R
S
R
-
q
R
.For details, refer to [33, 34].The second flux function implemented in this paper is LLF solver [35]. It is more dissipate than both the HLLC flux function and the Roe flux function, but it is more robust. The LLF flux can be written as(13)
F
L
L
F
=
1
2
F
L
+
F
R
-
λ
′
+
c
′
U
R
-
U
L
,where λ
′ is the velocity normal to the interface and c
′ is the speed of sound at the interface. λ
′
+
c
′ is the largest wave speed in the direction normal to the interface.The third flux function is Roe numerical flux by Roe [36]. Roe format is a typical flux differential splitting scheme, which contains more feature information, and therefore has a strong ability to capture shock. The Roe flux can be expressed as(14)
F
R
o
e
=
1
2
F
L
+
F
R
-
A
~
U
R
-
U
L
,where(15)
A
~
U
R
-
U
L
=
α
4
u
α
4
+
n
x
α
5
+
α
6
v
α
4
+
n
y
α
5
+
α
7
h
α
4
+
λ
′
α
5
+
u
α
6
+
v
α
7
-
c
′
2
γ
-
1
α
1
,
α
1
=
λ
′
Δ
ρ
-
Δ
p
c
′
2
,
α
2
=
λ
′
+
c
′
Δ
p
+
ρ
c
′
Δ
λ
′
2
c
′
2
,
α
3
=
λ
′
-
c
′
Δ
p
-
ρ
c
′
Δ
λ
′
2
c
′
2
,
α
4
=
α
1
+
α
2
+
α
3
,
α
5
=
c
′
α
2
-
α
3
,
α
6
=
ρ
λ
′
Δ
u
-
n
x
Δ
λ
′
,
α
7
=
ρ
λ
′
Δ
v
-
n
y
Δ
λ
′
,
λ
′ is the velocity normal to the interface, and c
′ is the speed of sound at the interface.In this paper, Van Albada limiter [30] is employed to make the scheme monotone for 2D problems and Barth-Jespersen limiter [37] is used for 3D problems.


## 3.1. DG Spatial Discretization

The computational domainΩ is partitioned into an ensemble of nonoverlapping elements, triangles in 2D; that is, Ω
=
e
1
∪
e
2
∪
⋯
∪
e
n, where n denotes the number of elements in the domain. We consider possible choices of the piecewise basis functions and then obtain the following weak formulation of (5) by multiplying a test function Φ and integrating by parts over the e
m:(5)
d
d
t
∫
e
m
Φ
h
U
h
d
V
-
∫
e
m
∂
Φ
h
∂
x
F
U
h
+
∂
Φ
h
∂
y
G
U
h
d
V
+
∫
∂
e
m
Φ
h
F
U
h
n
x
+
G
U
h
n
y
d
S
=
0
∀
Φ
,where U
h and Φ
h represent the finite element approximations of the analytical solution U and the test function Φ, respectively. n
=
(
n
x
,
n
y
) is the unit normal vector of outward to the boundary. Let the approximate solution and test function U
h and Φ
h be expressed as(6)
U
h
=
∑
j
=
1
N
U
j
t
φ
j
p
x
,
Φ
h
=
∑
j
=
1
N
Φ
j
φ
j
p
x
,where φ
j
p
(
x
) is the shape function of the polynomials of degree p. Equation (5) must be satisfied for any test function Φ
h, so by substituting (6) to (5), we obtain the following system of N equations:(7)
d
U
j
t
d
t
∫
e
m
φ
i
p
φ
j
p
d
V
-
∫
e
m
∂
φ
i
p
∂
x
F
∑
j
=
1
N
U
j
φ
j
p
+
∂
φ
i
p
∂
y
G
∑
j
=
1
N
U
j
φ
j
p
d
V
+
∫
∂
e
m
φ
i
p
F
∑
j
=
1
N
U
j
φ
j
p
n
x
+
G
∑
j
=
1
N
U
j
φ
j
p
n
y
d
S
=
0
i
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
f
i
i
i
i
f
f
f
f
i
i
=
1,2
,
…
,
N
.The interface flux function F
(
U
h
)
n
x
+
G
(
U
h
)
n
y can be treated as a numerical Riemann flux function H
(
U
h
L
,
U
h
R
,
n
), where U
h
L and U
h
R represent the internal element interface solution and neighboring element interface solution, respectively. In the present work, the Roe, LLF, and HLLC approximate Riemann solvers are employed. The domain and boundary integrals in (7) are calculated by use of 2
p and 2
p
+
1 order accurate Gauss quadrature formulas [31, 32] with a number of quadrature points corresponding to the degree of interpolating polynomials.By grouping together all the elemental time-dependent and spatial contributions, (7) can be written as a system of ordinary differential equations:(8)
M
d
U
d
t
=
R
U
,where the mass matrix M has identical diagonal blocks M
i
j, U is the global vector of the degrees of freedom, and R
(
U
) represents the steady state residual vector. As a result, the inverse of the mass matrix M can be easily computed, especially, using the orthogonal basis functions, and stored in advance due to the fact that it remains unchanged during the process.In the present paper, we explore the orthogonal basis functions through Gram-Schmidt orthogonalization method, and high accuracy can be achieved by improving the order of the piecewise interpolation polynomials functions. For 2D problems,φ
0
=
1, φ
1
=
ξ, φ
2
=
η, φ
3
=
ξ
2
,
…, and, for 3D problems, φ
0
=
1, φ
1
=
ξ, φ
2
=
η, φ
3
=
ζ
,
….


## 3.2. Flux Functions

The numerical flux function can be evaluated using any upwind flux functions. This is exactly similar to FVM because discontinuities can be allowed across the interface. Therefore, approximate Riemann solvers can be used to compute the flux function. In the present work, three flux functions have been employed including HLLC Riemann numerical solver by Toro [33] which has easier and lower computational cost in comparison with many other available Riemann solvers, such as local Lax-Friedrich (LLF) scheme. HLLC flux function not only maintains the advantages of the HLL solver but also resolves isolated contact discontinuities exactly, which has been extended in conjunction with time-derivative preconditioning to compute flow problems at all speeds. HLL flux can be expressed as(9)
F
H
L
L
=
F
L
i
f
S
L
>
0
F
L
R
∗
i
f
S
L
≤
0
≤
S
R
F
R
i
f
S
R
<
0
,where S
L and S
R represent the fastest wave speed for the left and right states, respectively. F
L
R
∗ is written as(10)
F
L
R
∗
=
S
R
F
L
-
S
L
F
R
+
S
L
S
R
U
R
-
U
L
S
R
-
S
L
.The HLLC flux is a modification of HLL flux, which can be written as(11)
F
H
L
L
C
=
F
L
0
<
S
L
F
L
∗
=
F
L
+
S
L
U
L
∗
-
U
L
S
L
≤
0
,
S
M
>
0
F
R
∗
=
F
R
+
S
R
U
R
∗
-
U
R
S
M
≤
0
,
S
R
≥
0
F
R
0
>
S
R
,where S
M is constant between the two acoustic waves:(12)
S
M
=
p
R
-
p
L
+
ρ
L
q
L
S
L
-
q
L
-
ρ
R
q
R
S
R
-
q
R
ρ
L
S
L
-
q
L
-
ρ
R
S
R
-
q
R
.For details, refer to [33, 34].The second flux function implemented in this paper is LLF solver [35]. It is more dissipate than both the HLLC flux function and the Roe flux function, but it is more robust. The LLF flux can be written as(13)
F
L
L
F
=
1
2
F
L
+
F
R
-
λ
′
+
c
′
U
R
-
U
L
,where λ
′ is the velocity normal to the interface and c
′ is the speed of sound at the interface. λ
′
+
c
′ is the largest wave speed in the direction normal to the interface.The third flux function is Roe numerical flux by Roe [36]. Roe format is a typical flux differential splitting scheme, which contains more feature information, and therefore has a strong ability to capture shock. The Roe flux can be expressed as(14)
F
R
o
e
=
1
2
F
L
+
F
R
-
A
~
U
R
-
U
L
,where(15)
A
~
U
R
-
U
L
=
α
4
u
α
4
+
n
x
α
5
+
α
6
v
α
4
+
n
y
α
5
+
α
7
h
α
4
+
λ
′
α
5
+
u
α
6
+
v
α
7
-
c
′
2
γ
-
1
α
1
,
α
1
=
λ
′
Δ
ρ
-
Δ
p
c
′
2
,
α
2
=
λ
′
+
c
′
Δ
p
+
ρ
c
′
Δ
λ
′
2
c
′
2
,
α
3
=
λ
′
-
c
′
Δ
p
-
ρ
c
′
Δ
λ
′
2
c
′
2
,
α
4
=
α
1
+
α
2
+
α
3
,
α
5
=
c
′
α
2
-
α
3
,
α
6
=
ρ
λ
′
Δ
u
-
n
x
Δ
λ
′
,
α
7
=
ρ
λ
′
Δ
v
-
n
y
Δ
λ
′
,
λ
′ is the velocity normal to the interface, and c
′ is the speed of sound at the interface.In this paper, Van Albada limiter [30] is employed to make the scheme monotone for 2D problems and Barth-Jespersen limiter [37] is used for 3D problems.


## 4. Time-Marching Schemes

In order to resolve the time-dependent problem, the semidiscrete system can be integrated in time in this paper. The implicit time-integration schemes have been widespread for DG discretization. Iterative algorithms such as GMRES and CGS are often used to approximately solve the sparse linear equations due to the enormous computational cost and the large memory requirement of direct methods. Another implicit scheme, the LU-SGS scheme originally developed for structured grids, has been extended to unstructured and hybrid grids, which does not require any extra storage compared to explicit methods. The LU-SGS procedures are described as follows.


### 4.1. PCLU-SGS Scheme

In the original LU-SGS approach [29], (8) can be translated into the following system:(16)
A
Δ
Q
=
R
e
s
n
,where R
e
s
n represents residual term.Then the coefficient matrix using the decomposition method can be written asA
=
D
+
L
+
U; we obtain(17)
D
+
L
D
-
1
D
+
U
Δ
Q
n
=
R
e
s
n
+
L
D
-
1
U
Δ
Q
n
,where D represents the diagonal matrix and L and U represent the lower and upper matrices. Ignoring the infinitesimal quantity (
L
D
-
1
U
)
Δ
Q
n, (17) is then solved using one sweep of symmetric Gauss-Seidel iteration as shown in the following:forward sweep:(18)
D
+
L
Δ
Q
~
=
R
e
s
n
,
backward sweep:(19)
D
-
1
D
+
U
Δ
Q
-
=
Δ
Q
~
.While ignoring the higher-order infinitesimal quantity does not affect the accuracy of the method, increasing the truncation error will affect the rate of convergence. Therefore, we obtain the PCLU-SGS algorithm through the compensation of high-order term for the original LU-SGS scheme. The computational procedure is shown in detail as follows.(
1
)
Use the original LU-SGS scheme to solveΔ
Q
-,
forward sweep and backward sweep:(20)
D
+
L
Δ
Q
~
=
R
e
s
n
,
D
-
1
D
+
U
Δ
Q
-
=
Δ
Q
~
.
(
2
)
Compute the high-order infinitesimal quantity(
L
D
-
1
U
)
Δ
Q
n using the value of Δ
Q
- to correct the residual R
e
s
n.
(
3
)
Use the original LU-SGS scheme to computeΔ
Q:(21)
D
+
L
D
-
1
D
+
U
Δ
Q
n
=
R
e
s
¯
n
=
R
e
s
n
+
L
D
-
1
U
Δ
Q
-
.


### 4.2. TVD Runge-Kutta Scheme

The explicit time-integration schemes have been widespread for DG discretization. The TVD Runge-Kutta scheme of Cockburn [10] for Euler equations can be expressed as follows.(1)
Denoteu
h
0
=
P
v
h
(
u
0
), where P
v
h is the projection operator on L
2.
(2)
Forn
=
0,1
,
…
,
N
-
1, compute u
h
n
+
1 and denote u
h
0
=
u
h
n, and,
forj
=
1
,
…
,
p
+
1, compute the intermediate function, u
h
(
j
):(22)
u
h
j
=
∑
m
=
0
j
-
1
α
j
m
u
h
m
+
β
j
m
Δ
t
n
M
-
1
L
h
u
h
m
,
γ
h
t
n
+
d
m
Δ
t
n
.
(3)
Denoteu
h
n
+
1
=
u
h
p
+
1.The scheme is linearly stable for a Courant number less than or equal to 1
/
(
2
p
+
1
). In this paper, the RK4 scheme is employed to compute the Euler equations.


## 4.1. PCLU-SGS Scheme

In the original LU-SGS approach [29], (8) can be translated into the following system:(16)
A
Δ
Q
=
R
e
s
n
,where R
e
s
n represents residual term.Then the coefficient matrix using the decomposition method can be written asA
=
D
+
L
+
U; we obtain(17)
D
+
L
D
-
1
D
+
U
Δ
Q
n
=
R
e
s
n
+
L
D
-
1
U
Δ
Q
n
,where D represents the diagonal matrix and L and U represent the lower and upper matrices. Ignoring the infinitesimal quantity (
L
D
-
1
U
)
Δ
Q
n, (17) is then solved using one sweep of symmetric Gauss-Seidel iteration as shown in the following:forward sweep:(18)
D
+
L
Δ
Q
~
=
R
e
s
n
,
backward sweep:(19)
D
-
1
D
+
U
Δ
Q
-
=
Δ
Q
~
.While ignoring the higher-order infinitesimal quantity does not affect the accuracy of the method, increasing the truncation error will affect the rate of convergence. Therefore, we obtain the PCLU-SGS algorithm through the compensation of high-order term for the original LU-SGS scheme. The computational procedure is shown in detail as follows.(
1
)
Use the original LU-SGS scheme to solveΔ
Q
-,
forward sweep and backward sweep:(20)
D
+
L
Δ
Q
~
=
R
e
s
n
,
D
-
1
D
+
U
Δ
Q
-
=
Δ
Q
~
.
(
2
)
Compute the high-order infinitesimal quantity(
L
D
-
1
U
)
Δ
Q
n using the value of Δ
Q
- to correct the residual R
e
s
n.
(
3
)
Use the original LU-SGS scheme to computeΔ
Q:(21)
D
+
L
D
-
1
D
+
U
Δ
Q
n
=
R
e
s
¯
n
=
R
e
s
n
+
L
D
-
1
U
Δ
Q
-
.


## 4.2. TVD Runge-Kutta Scheme

The explicit time-integration schemes have been widespread for DG discretization. The TVD Runge-Kutta scheme of Cockburn [10] for Euler equations can be expressed as follows.(1)
Denoteu
h
0
=
P
v
h
(
u
0
), where P
v
h is the projection operator on L
2.
(2)
Forn
=
0,1
,
…
,
N
-
1, compute u
h
n
+
1 and denote u
h
0
=
u
h
n, and,
forj
=
1
,
…
,
p
+
1, compute the intermediate function, u
h
(
j
):(22)
u
h
j
=
∑
m
=
0
j
-
1
α
j
m
u
h
m
+
β
j
m
Δ
t
n
M
-
1
L
h
u
h
m
,
γ
h
t
n
+
d
m
Δ
t
n
.
(3)
Denoteu
h
n
+
1
=
u
h
p
+
1.The scheme is linearly stable for a Courant number less than or equal to 1
/
(
2
p
+
1
). In this paper, the RK4 scheme is employed to compute the Euler equations.


## 5. Numerical Experiments

### 5.1. Transonic Flow around NACA0012 Airfoil

Consider the calculation stateM
a
∞
=
0.8 and 1.2
5
° angle of attack. The pressure coefficient distribution using different time-marching schemes compared with the experiment results for transonic flow around NACA0012 airfoil is given in Figure 2. Good agreement can be seen in terms of the location and strength of shocks (see Figure 1). The flood contours are shown in Figure 3. From Figure 4 it can be seen that it takes only 1500 iterations to obtain the result by using the present algorithm, which are far fewer than the 4500 iterations required to obtain the same results by using the original LU-SGS method. The computational efficiency is close to that of GMRES algorithm and nearly 2.1 times greater than that of the LU-SGS one. Figure 5 shows that LLF solver similar to the HLL solver is more dissipate than both the HLLC flux function and the Roe flux function. Figure 6 shows the effects of convergence performance with different CFL number; it can be clearly seen that the results are almost the same when CFL number is greater than 100. The results of the test cases verify the effectiveness and the ability to capture discontinuous of PCLU-SGS DG method.Figure 1
Mesh around NACA0012 airfoil: (a) whole and (b) close-up mesh near the airfoil.
(a)
(b)Figure 2
Pressure coefficient distribution with different time-marching format for transonic flow around NACA0012 airfoil (p
=
2).Figure 3
Flood contours.Figure 4
Convergence history for transonic flow around NACA0012 airfoil: (a) CPU-Time and (b) iteration number.
(a)
(b)Figure 5
Pressure coefficient distribution with different numerical flux for transonic flow around NACA0012 airfoil: (a)p
=
2 and (b) p
=
3.
(a)
(b)Figure 6
Convergence history with different CFL number for transonic flow around NACA0012 airfoil: (a) CPU-Time and (b) iteration number.
(a)
(b)Furthermore, the Sod shock tube problem with the initial conditions is given as follows:(23)
U
=
ρ
L
,
u
L
,
p
L
T
=
1,0
,
1
T
,
i
f
-
5
≤
x
≤
0
ρ
R
,
u
R
,
p
R
T
=
0.125,0
,
0.1
T
,
i
f
0
≤
x
≤
5
,where t
=
2.0, and the mesh consists of 100 elements in Figure 7.Figure 7
Density with HLLC flux for Sod shock tube problem: (a)p
=
1 and (b) p
=
2.
(a)
(b)Obviously, with the improvement of the accuracy, the method with HLLC flux resolves better contact discontinuity.In the present paper, though the LLF flux function is more dissipate than both the HLLC flux function and the Roe flux function, it is more robust and more economical. Then we use it in the following examples.


### 5.2. Transonic Flow around RAE2822 Airfoil

Consider the calculation stateM
a
=
0.725 and angle of attack α
=
2.5
4
°. The pressure coefficient distribution using different time-marching schemes for transonic flow around RAE2822 airfoil is given in Figure 8. The numerical solutions demonstrate that it takes only 1200 iterations to obtain the result by using the present algorithm, which are far fewer than 4000 iterations required to obtain the same results by using the original LU-SGS method. And the results show that the convergence acceleration is nearly 2.3 times that of the original LU-SGS one from Figure 9.Figure 8
Pressure coefficient distribution with different time-marching format flow around RAE2822 airfoil(
p
=
2
).Figure 9
Convergence history with different time-marching format for transonic flow around RAE2822 airfoil: (a) CPU-Time and (b) iteration number.
(a)
(b)


### 5.3. Transonic Flow over ONERA M6 Wing

This case is about a transonic flow at Mach numberM
a
=
0.84 around the ONERA M6 wing with α
=
3.0
6
° angle of attack. The unstructured mesh consists of 582752 triangles in Figure 10. The surface pressure coefficient distribution with different spanwise location is given in Figure 12. From Figure 11 it can be seen that it takes only 5050 iterations to obtain the result by using the present algorithm, which are far fewer than 44021 iterations required to obtain the same results by using the RK4 method. The convergence acceleration is nearly 5.5 times that of the RK4 one and nearly half that of GMRES one. The convergence performance is similar to that in [38].Figure 10
Surface mesh around ONERA M6 wing and flood contours: (a) surface mesh and (b) flood contours.
(a)
(b)Figure 11
Convergence history with different time-marching format for transonic flow around ONERA M6 wing: (a) CPU-Time and (b) iteration number.
(a)
(b)Figure 12
Surface pressure coefficient distribution with different spanwise location for transonic flow over ONERA M6 wing: (a) 20% and (b) 80%(
p
=
1
).
(a)
(b)The results are not as good as those in the literature [39], because the shock detector is used, which directly affect the accuracy of solutions on smooth region. How to accurately judge problem units will be our future efforts.


## 5.1. Transonic Flow around NACA0012 Airfoil

Consider the calculation stateM
a
∞
=
0.8 and 1.2
5
° angle of attack. The pressure coefficient distribution using different time-marching schemes compared with the experiment results for transonic flow around NACA0012 airfoil is given in Figure 2. Good agreement can be seen in terms of the location and strength of shocks (see Figure 1). The flood contours are shown in Figure 3. From Figure 4 it can be seen that it takes only 1500 iterations to obtain the result by using the present algorithm, which are far fewer than the 4500 iterations required to obtain the same results by using the original LU-SGS method. The computational efficiency is close to that of GMRES algorithm and nearly 2.1 times greater than that of the LU-SGS one. Figure 5 shows that LLF solver similar to the HLL solver is more dissipate than both the HLLC flux function and the Roe flux function. Figure 6 shows the effects of convergence performance with different CFL number; it can be clearly seen that the results are almost the same when CFL number is greater than 100. The results of the test cases verify the effectiveness and the ability to capture discontinuous of PCLU-SGS DG method.Figure 1
Mesh around NACA0012 airfoil: (a) whole and (b) close-up mesh near the airfoil.
(a)
(b)Figure 2
Pressure coefficient distribution with different time-marching format for transonic flow around NACA0012 airfoil (p
=
2).Figure 3
Flood contours.Figure 4
Convergence history for transonic flow around NACA0012 airfoil: (a) CPU-Time and (b) iteration number.
(a)
(b)Figure 5
Pressure coefficient distribution with different numerical flux for transonic flow around NACA0012 airfoil: (a)p
=
2 and (b) p
=
3.
(a)
(b)Figure 6
Convergence history with different CFL number for transonic flow around NACA0012 airfoil: (a) CPU-Time and (b) iteration number.
(a)
(b)Furthermore, the Sod shock tube problem with the initial conditions is given as follows:(23)
U
=
ρ
L
,
u
L
,
p
L
T
=
1,0
,
1
T
,
i
f
-
5
≤
x
≤
0
ρ
R
,
u
R
,
p
R
T
=
0.125,0
,
0.1
T
,
i
f
0
≤
x
≤
5
,where t
=
2.0, and the mesh consists of 100 elements in Figure 7.Figure 7
Density with HLLC flux for Sod shock tube problem: (a)p
=
1 and (b) p
=
2.
(a)
(b)Obviously, with the improvement of the accuracy, the method with HLLC flux resolves better contact discontinuity.In the present paper, though the LLF flux function is more dissipate than both the HLLC flux function and the Roe flux function, it is more robust and more economical. Then we use it in the following examples.


## 5.2. Transonic Flow around RAE2822 Airfoil

Consider the calculation stateM
a
=
0.725 and angle of attack α
=
2.5
4
°. The pressure coefficient distribution using different time-marching schemes for transonic flow around RAE2822 airfoil is given in Figure 8. The numerical solutions demonstrate that it takes only 1200 iterations to obtain the result by using the present algorithm, which are far fewer than 4000 iterations required to obtain the same results by using the original LU-SGS method. And the results show that the convergence acceleration is nearly 2.3 times that of the original LU-SGS one from Figure 9.Figure 8
Pressure coefficient distribution with different time-marching format flow around RAE2822 airfoil(
p
=
2
).Figure 9
Convergence history with different time-marching format for transonic flow around RAE2822 airfoil: (a) CPU-Time and (b) iteration number.
(a)
(b)


## 5.3. Transonic Flow over ONERA M6 Wing

This case is about a transonic flow at Mach numberM
a
=
0.84 around the ONERA M6 wing with α
=
3.0
6
° angle of attack. The unstructured mesh consists of 582752 triangles in Figure 10. The surface pressure coefficient distribution with different spanwise location is given in Figure 12. From Figure 11 it can be seen that it takes only 5050 iterations to obtain the result by using the present algorithm, which are far fewer than 44021 iterations required to obtain the same results by using the RK4 method. The convergence acceleration is nearly 5.5 times that of the RK4 one and nearly half that of GMRES one. The convergence performance is similar to that in [38].Figure 10
Surface mesh around ONERA M6 wing and flood contours: (a) surface mesh and (b) flood contours.
(a)
(b)Figure 11
Convergence history with different time-marching format for transonic flow around ONERA M6 wing: (a) CPU-Time and (b) iteration number.
(a)
(b)Figure 12
Surface pressure coefficient distribution with different spanwise location for transonic flow over ONERA M6 wing: (a) 20% and (b) 80%(
p
=
1
).
(a)
(b)The results are not as good as those in the literature [39], because the shock detector is used, which directly affect the accuracy of solutions on smooth region. How to accurately judge problem units will be our future efforts.


## 6. Conclusion

An improved implicit time-marching scheme based on the original LU-SGS scheme is developed and applied for the discontinuous Galerkin method on unstructured grids. The developed new algorithm has been used to compute the transonic flows around NACA0012 airfoil, RAE2822 airfoil, and ONERA M6 wing. The implicit PCLU-SGS scheme for the DG method on unstructured grids is significantly more efficient and robust than the original LU-SGS scheme. The convergence performance of the present scheme can compete with the GMRES scheme.



---

*Source: 940257-2015-06-23.xml*