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.
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:where is the length of time interval and is a two-dimensional bounded domain. The conservative state vector and the inviscid flux component vectors and are defined bywhere the notations , , and denote the density, pressure, and specific total energy per unit mass of the fluid, respectively. and are the velocity components of the flow in the and coordinate directions. This set of equations is completed by the perfect gas equation of state given as follows:where is defined as the ratio of specific heat of the fluid ( 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:
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, , where 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 :where and represent the finite element approximations of the analytical solution and the test function , respectively. is the unit normal vector of outward to the boundary. Let the approximate solution and test function and be expressed aswhere is the shape function of the polynomials of degree . Equation (5) must be satisfied for any test function , so by substituting (6) to (5), we obtain the following system of equations:The interface flux function can be treated as a numerical Riemann flux function , where and 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 and 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:where the mass matrix has identical diagonal blocks , is the global vector of the degrees of freedom, and represents the steady state residual vector. As a result, the inverse of the mass matrix 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, , , , , and, for 3D problems, , , , .
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 aswhere and represent the fastest wave speed for the left and right states, respectively. is written asThe HLLC flux is a modification of HLL flux, which can be written aswhere is constant between the two acoustic waves: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 aswhere is the velocity normal to the interface and is the speed of sound at the interface. 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 aswhere is the velocity normal to the interface, and 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:where represents residual term.
Then the coefficient matrix using the decomposition method can be written as ; we obtainwhere represents the diagonal matrix and and represent the lower and upper matrices. Ignoring the infinitesimal quantity , (17) is then solved using one sweep of symmetric Gauss-Seidel iteration as shown in the following: forward sweep: backward sweep: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.Use the original LU-SGS scheme to solve , forward sweep and backward sweep:Compute the high-order infinitesimal quantity using the value of to correct the residual .Use the original LU-SGS scheme to compute :
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)Denote , where is the projection operator on .(2)For , compute and denote , and, for , compute the intermediate function, :(3)Denote .The scheme is linearly stable for a Courant number less than or equal to . 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 state and 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.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Furthermore, the Sod shock tube problem with the initial conditions is given as follows:where , and the mesh consists of 100 elements in Figure 7.
(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 state and angle of attack . 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.
(a)
(b)
5.3. Transonic Flow over ONERA M6 Wing
This case is about a transonic flow at Mach number around the ONERA M6 wing with 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].
(a)
(b)
(a)
(b)
(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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was supported by the National Natural Science Funds of China under Grants nos. 11002117 and 61373174 and the Xianyang Normal University Research Funds under Grants nos. 09XSYK209 and 09XSYK204.