Time Advancement of the Navier-Stokes Equations: p-Adaptive Exponential Methods

Shujie Li

Beijing Computational Science Research Center, Beijing, China.

DOI: 10.4236/jfcmv.2020.82004                

Abstract

An adaptive exponential time advancement framework is developed for solving the multidimensional Navier-Stokes equations with a variable-order discontinuous Galerkin (DG) discretization on hybrid unstructured curved grids. The adaptive framework is realized with cell-wise, variable-order DG refinements and a dynamic assembly of elemental Jacobian matrices. The accuracy and performance gain are investigated for several benchmark cases up to a realistic, three-dimensional rotor flow. Numerical results are shown to be more efficient than the use of uniform-order exponential DG for simulating viscous flows.

Keywords

Exponential Time Discreitzation, Navier-Stokes Equation, Discontinuous Galerkin, Curved Grids

1. Introduction

Computational efficiency is one of the most concerned problems towards industrial applications of high-order methods. While various high-order methods have been gaining popularity, their computational efficiency can still be enhanced by employing a fast time marching method. To this end, a class of exponential time integration methods has been developed for fast time stepping of high-order discretized compressible Navier-Stokes equations, demonstrated in an arbitrarily high-order discontinuous Galerkin framework HA3D [1] - [6] on hybrid curved elements, although it is generally applicable to any spatial discretization. The developed predictor-corrector exponential (PCEXP) time scheme [3] is shown to be able to eliminate the Courant-Friedrichs-Lewy (CFL) restriction with low absolutely temporal errors, showing a special capability of achieving fast time stepping for steady and unsteady, inviscid and viscous flows.

In this work, we exploit the possibility of reducing the computational cost associated with the use of uniform-order exponential DG methods in a variable-order solution adaptation frame, and to our best knowledge, it has not been reported yet. In the current adaptation framework, cell order of accuracy is dynamically determined with a local cell error indicator in a modal discontinuous Galerkin method, resulting in a cell-wise adaptive distribution of cell orders and variable-order elemental Jacobian matrices.

The remaining parts of this abstract are organized as follows. Section 2 introduces the spatial discretization method used. Section 3 gives the overview of the exponential time-marching methods including the accuracy (p) adaptation methods. Finally, the last section presents some numerical tests and the Caradonna-Tung rotor in hover [7] is simulated with the adaptive method for demonstrating its applicability to real-world problems.

2. Spatial Discretization

2.1. Governing Equations

Governing equations consider the three-dimensional compressible Navier-Stokes equations with a source term in the rotating reference frame

$U_t+\nabla \cdot F^c\left(U\right)+\nabla \cdot F^v\left(U\mathrm{,}\nabla U\right)=S\left(U\right)\mathrm{,}$ (1)

where $U,F^c$ are the conventional conservative vector and viscous flux. For rotational reference frame, the source term $S$ and inviscid flux $F^c$ are defined as

$F^c=\left(\begin{array}{c}\rho {\left(v-v_r\right)}^{\text{T}}\\ \rho \left(v-v_r\right)v^{\text{T}}+pI\\ \rho H{\left(v-v_r\right)}^{\text{T}}\end{array}\right)\mathrm{,}S=\left(\begin{array}{c}0\\ -\rho \omega \times v\\ 0\end{array}\right)\mathrm{,}$ (2)

where $v={\left(u\mathrm{,}v\mathrm{,}w\right)}^{\text{T}}$ is the absolute velocity, $\omega ={\left({\omega }_x\mathrm{,}{\omega }_y\mathrm{,}{\omega }_z\right)}^{\text{T}}$ is the angular velocity of the rotating frame of reference, $v_r=\omega \times x$ the relative velocity; $\tau$, $q$ the viscous stress tensor and the heat flux vector. $\rho$, p, and e denote the

flow density, pressure, and the specific internal energy; $E=e+\frac{1}{2}{\Vert v\Vert }^2$ and

$H=E+p/\rho$ denote the total energy and total enthalpy, respectively; $I$ denotes the $3\times 3$ unit matrix; and the pressure p is given by the equation of state for a perfect gas

$p=\rho \left(\gamma -1\right)e\mathrm{,}$ (3)

where $\gamma =7/5$ is the ratio of specific heats for perfect gas.

2.2. Modal Discontinuous Galerkin Method

The discontinuous Galerkin discretization of Equation (1) is defined on a computational domain $\Omega$ divided into curved elements of arbitrary shape. The present adaptive discontinuous Galerkin method seeks an variable-order approximation $U$ in each element $E\in \Omega$ with a p-order polynomial $\left\{P_p\left(E\right)\in L_2\left(\Omega \right)\mathrm{,}\forall E\in \Omega \right\}$, namely

$U\left(x\mathrm{,}t\right)={\displaystyle \underset{j=1}{\overset{N\left(p\right)}{\sum }}}{u}_j\left(t\right){\psi }_j\left(x\right)\mathrm{.}$ (4)

For 3-D problems $N\left(p\right)=\left(p+1\right)\left(p+2\right)\left(p+3\right)/6$. The cell order p is a local variable for each cell, thus an adaptive approximation is simply obtained by adding or removing terms of (18), accordingly. An orthonormal basis set $\psi \left(x\right)$ expressed in the Cartesian coordinates is used to facilitate the implementation of the adaptation framework [1] [3]. By multiplying Equation (1) with the adaptive-order basis functions, the weak form is obtained

${\displaystyle {\int }_E}{\psi }_i{\psi }_j\text{d}x\frac{\text{d}{u}_j}{\text{d}t}=-{\displaystyle {\int }_{\partial E}}{\psi }_i\tilde{F}\cdot \hat{n}\text{d}\sigma +{\displaystyle {\int }_E}\left(F\cdot \nabla {\psi }_i+{\psi }_{i}S\right)\text{d}x\mathrm{<mo>\; :\; </mo>}=R_i\mathrm{,}$ (5)

where the Einstein summation convention is used. Here $\hat{n}$ denotes the outward unit normal of the surface element $\partial E$ of the element E. The flux terms are defined as

$\tilde{F}={\tilde{F}}^c\left(u_h^{\pm }\right)+{\tilde{F}}^v\left(u_h^{\pm }\mathrm{,}{\left({\nabla }_h{u}_h+{\delta }_f\right)}^{\pm }\right)\mathrm{<mo>\; ;\; </mo>}F=F^c\left(u_h\right)+F^v\left(u_h\mathrm{,}\left({\nabla }_h{u}_h+\delta \right)\right)\mathrm{,}$ (6)

where ${\tilde{F}}^c$ is computed by a Riemann solver, and the viscous flux ${\tilde{F}}^v$ is computed with the second approach of Bassi and Rebay (BR2) which introduces the local and global lifting operators ${\delta }_f$ and $\delta$ [8]. The BR2 viscous flux discretization is compact where the global lifting operator is used for volume integration while the local lifting operator is used face-wised

$\delta ={\displaystyle \underset{\partial E}{\sum }}{\delta }_f\mathrm{,}$ (7)

and the local lifting operator is solved by Galerkin projection on each surface $\sigma$,

where we introduce the average operator $\left\{a\right\}=\frac{1}{2}\left(a^{+}+a^{-}\right)$

${\displaystyle {\int }_{\partial E}}\psi \cdot {\delta }_f={\displaystyle {\int }_{\partial E}}\psi \cdot \left(\left\{u\right\}-u_h^{+}\right)\text{d}\sigma$ (8)

The implementation of BR2 formulation to exponential schemes uses analytical global Jacobian [5] [6], which is composed of variable order cell Jacobians in the current framework.

3. Exponential Time Integration

The PCEXP scheme that originally developed for fast time stepping of semi-discretized equations is independent of the choice of spatial discretization. And this feature makes it applicable as a general time integration tool. The previous results show it permits the use of very large time steps [1] - [6], and are more efficient than the third-order explicit TVD Runge-Kutta method and the second-order implicit backward difference formula [3]. In this section, the PCEXP scheme is briefly overviewed and its variable-order implementation is discussed.

We start with the following semi-discrete system of ordinary differential equations which may be obtained from a spatial discretization with problem associated boundary conditions:

$\frac{\text{d}u}{\text{d}t}=R\left(u\right)\mathrm{,}$ (9)

where $u=u\left(t\right)\in {ℝ}^K$ denotes the vector of the solution variables and $R\left(u\right)\in {ℝ}^K$ the right-hand-side term which may be the spatially discretized residual terms of the discontinuous Galerkin method used in this work. The dimension K is the degrees of freedom which can be very large for 3-D problems. Without loss of generality, we consider $u\left(t\right)$ in the interval of one time step, i.e., $t\in \left[t_n\mathrm{,}{t}_{n+1}\right]$. Splitting the right hand side leads to a different exact expression

$\frac{\text{d}u}{\text{d}t}=J_{n}u+N\left(u\right)\mathrm{,}$ (10)

where the subscript n indicates the value evaluated at $t=t_n$, $J_n$ denotes the Jacobian matrix $J_n={\partial R\left(u\right)/\partial u|}_{t=t_n}=\partial R\left(u_n\right)/\partial u$ and $N\left(u\right)=R\left(u\right)-J_{n}u$ denotes the remainder, which in general is nonlinear. Equation (10) admits the following formal solution:

$u_{n+1}=exp\left(\Delta t{J}_n\right)u_n+{\displaystyle {\int }_0^{\Delta t}}exp\left(\left(\Delta t-\tau \right)J_n\right)N\left(u\left(t_n+\tau \right)\right)\text{d}\tau \mathrm{,}$ (11)

Here, the matrix exponential is defined similarly to its scalar version

$exp\left(-t{J}_n\right)={\displaystyle \underset{m=0}{\overset{\infty }{\sum }}}\frac{{\left(-t{J}_n\right)}^m}{m\mathrm{<mo>\; !\; </mo>}}$ (12)

is the integrating factor. The formal solution (11) is the starting point to derive the proposed exponential scheme in which the stiff part is computed analytically whereas the nonlinear term is approximated numerically. The PCEXP approximation of (11) is read as

$u^{\mathrm{<mo>\; *\; </mo>}}=u_n+\Delta t{\Phi }_1\left(\Delta t{J}_n\right)R\left(u_n\right)\mathrm{.}$ (13)

$u_{n+1}=u^{\mathrm{<mo>\; *\; </mo>}}+\frac{1}{2}\Delta t{\Phi }_1\left(\Delta t{J}_n\right)\left[N\left(u^{\mathrm{<mo>\; *\; </mo>}}\right)-N\left(u_n\right)\right]\mathrm{.}$ (14)

where a new matrix function is defined as

${\Phi }_1\left(\Delta t{J}_n\right)\mathrm{<mo>\; :\; </mo>}=\frac{J^{-1}}{\Delta t}\left[exp\left(\Delta tJ\right)-I\right]\mathrm{,}$ (15)

and $I$ denotes the $K\times K$ identity matrix.

The physical nature of such type of exponential schemes relies on the global coupling feature via the global Jacobian matrix $J$, so that flow transportation information can be broadcasted to the whole computational domain without a CFL restriction. That is why the exponential schemes behave like a fully implicit method but only depends on the current solution, i.e., in an explicit way of formulation (14). While the second-order PCEXP scheme is more appropriate for computing unsteady problems, its first-order version, EXP1 only takes the first stage of (14), thus is more efficient for steady flows as shown in reference [2] [3] [5]. In this paper, we used PCEXP for all the test cases for both steady and unsteady flows.

3.1. p-Adaptive Time Steps

The time step of the adaptive exponential DG is dynamically determined via a residual monitoring strategy [3]. The local time step $\Delta t$ is chosen by the minimum of the convection time step $\Delta t_{\text{c}}$ and the diffusion time step $\Delta t_{\text{d}}$, which are given below

$\begin{array}{l}\Delta t_{\text{c}}=\frac{\text{CFL}{h}_{\text{3D}}}{\left(2p+1\right)\left(\Vert v\Vert +c\right)}\mathrm{<mo>\; ;\; </mo>}\\ \Delta t_{\text{d}}=\frac{\text{CFL}{h}_{\text{3D}}^2}{{\left(p+1\right)}^2\left(\frac{2\mu M_{\infty }}{\rho Re}\max \left(\frac{4}{3}\mathrm{,}\frac{\gamma }{Pr}\right)\right)}\mathrm{.}\end{array}$ (16)

Here, p denotes the cell order, $v$ and c the velocity vector and sound speed evaluated at the cell center, d the spatial dimension, and $h_{\text{3D}}$ represents a characteristic size of a 3-D cell defined by the ratio of volume and surface area. 2-D computations are also supported by the HA3D solver developed by the author which uses a quasi-3D mesh obtained by extruding the 2-D mesh by one layer of cells. In the 2-D computations, $h_{\text{2d}}$ is used instead for eliminating the effect of z dimension. Given the cell size $\Delta z$ in the z direction, $h_{\text{2D}}$ is determined by

$\frac{2}{h_{\text{2D}}}=\frac{3}{h_{\text{3D}}}-\frac{1}{\Delta z}\mathrm{.}$ (17)

A proven efficient CFL formula is employed in (16). For the current adaptive computations, we take the global minimal time step for time advancement. Using local time steps is possible, but might suffer from stability issues with very large time steps.

3.2. Adaptation Strategy for the p-Refinement

An adaptive, p-refinement strategy defines DG polynomial order p in a cell-wise way so that higher-order approximations can be placed locally in key flow regions such as those of near body and wake flows. To identify these regions, a so-called spectral delay error (SDE) indicator [9] is developed for locating shock waves originally. The SDE indicator, relating to the numerical resolution of the approximation, can also be used for p-adaptive refinement. Recalling the modal DG expansion of (18), we have

$U\left(x\mathrm{,}t\right)={\displaystyle \underset{j=1}{\overset{N\left(p\right)}{\sum }}}{u}_j\left(t\right){\psi }_j\left(x\right)\mathrm{.}$ (18)

The truncated expansion $U_{tr}\left(x\mathrm{,}t\right)$ for a lower degree ( $p-1$ ) is computed by cutting the summation

$U_{tr}\left(x\mathrm{,}t\right)={\displaystyle \underset{j=1}{\overset{N\left(p-1\right)}{\sum }}}{u}_j\left(t\right){\psi }_j\left(x\right)$ (19)

Applying the truncated expansion to the total energy variable

$\rho E=\rho e+\frac{1}{2}\rho v^2$ leads to

$\beta =\frac{{\displaystyle \int }{\left(U-U_{tr}\right)}^2\text{d}x}{{\displaystyle \int }{U}^2\text{d}x}=\frac{{\displaystyle \int }{\left(\rho E-{\left(\rho E\right)}_{tr}\right)}^2\text{d}x}{{\displaystyle \int }{\left(\rho E\right)}^2\text{d}x}$ (20)

The variable $\beta$ is finally used as the cell error indicator of the p-adaptation. The adaptation process starts with $p=0$ globally, namely first-order approximation. And then the cell error indicator $\beta$ is computed in each cell after a full convergence is obtained in each level of adaptation. The cell order p is updated in each cell with the following criterion

$p=\{\begin{array}{l}\mathrm{0,}{10}^{-1}\le \beta \\ p+\mathrm{1,}{10}^{-9}\le \beta \le {10}^{-1}\\ p\mathrm{,}\beta \le {10}^{-9}\end{array}$ (21)

The above criterion is applied only when the variation of global maximal Mach number is within 10⁻⁵ to prevent occurring premature adaptive solutions that are still in the process of dynamic evaluations.

3.3. Variable Order Residual Jacobians

The global residual Jacobian $J=\partial R/\partial u$ is required by the PCEXP scheme, which is composed of cell-wise, variable order Jacobians in the current adaptation framework. Unlike traditional nodal based high-order methods that require special treatments of interface order coupling [10] when different order approximations exist in two adjacent cells, the current modal based method is naturally compatible with variable accuracy without any treatment. The dimension of global residual Jacobian dynamically depends on the distribution of cell orders so that the dimension of elemental residual Jacobian depends on $N\left(p\right)$.

4. Numerical Results

4.1. Rotating Flow between Two Concentric Cylinders

The implementation of high-order discontinuous Galerkin discretization for the Navier-Stokes equations is verified on the rotating flow between two concentric cylinders or Taylor Couette flow [6]. The fluid is driven by two concentric cylinders which are rotating with constant angular velocities of ${\omega }_0$ and ${\omega }_1$. A low Reynolds number $Re=10$ is used for maintaining a laminar state, and the Reynolds number is defined by the tangential velocity and the radius of the inner cylinder. The analytical solution of this problem is given below

$u_{\theta }=r_0{\omega }_0\frac{r_1/r-r/r_1}{r_1/r_0-r_0/r_1}+r_1{\omega }_1\frac{r/r_0-r_0/r}{r_1/r_0-r_0/r_1}$ (22)

where $u_{\theta }$ is the tangential velocity, $r_0=1$ and $r_1=2$ are the inner radius and the outer radius. The isothermal boundary condition is set on the inner cylinder with angular velocity ${\omega }_0$ of Mach number $M_0=0.2$. The outer cylinder is stationary with ${\omega }_1=0$ and uses the adiabatic wall boundary condition. The front and the back faces of z-direction use the symmetric boundary condition for conducting quasi-2D computations. The order of accuracy is computed on a sequence of three meshes, using first- to fifth-order DG schemes ( ${\text{P}}^{0\to 4}$ ). Quadratic curved elements are used on the boundary surfaces. For this case, we focus on the error convergences. The $L_2$ norms of velocity errors are detailed in

Table 1. They are computed as $\sqrt{{\displaystyle \int }{\left(u-u_{\theta }\right)}^2\text{d}\Omega /{\displaystyle \int }\text{d}\Omega }$ integrating in the entire

computational domain $\Omega$. The expected order of convergence is observed for all the p levels, thereby verifying the high-order implementation of DG viscous discretization and also the use of curved elements.

4.2. Lid-Driven Cavity Flows

The developed adaptive framework is evaluated on the lid-driven cavity flows. The accuracy-adaptive solutions with $Re={10}^4$ are computed and compared with the baseline results of Ghia [11]. The top boundary is moving at a velocity at $M_{\infty }=0.2$, and the left, right, and bottom walls are set as the nonslip boundary condition. The front and the back faces of z-direction use the symmetric boundary condition for conducting quasi-2D computations. For the flow condition of $Re={10}^4$, secondary vortices show up in the corners of the cavity and high-resolution simulations are usually required. In this work, instead of using fine mesh, the adaptive frame is used so that a very coarse mesh ( $20\times 20$ ) can be used. In this case, we start by studying the impacts of discrete accuracy on the flow solution. Figure 1 presents the streamline plots using the ${\text{P}}_{0\to 1}$, ${\text{P}}_{0\to 4}$ solutions. It is observed that the ${\text{P}}_{0\to 1}$ solution is over-damped, where corner vortices are nearly disappeared due to the presence of high numerical viscosity of low-order approximations. The situation is improved when using higher-order adaptations where the corner vortices exhibit increased resolution of flow structures. It is also verified that the ${\text{P}}_{0\to 4}$ solution is essentially identical to the uniform ${\text{P}}_4$ solution of Figure 2(a). This agreement confirms the effectiveness of

Table 1. Rotating flow between two concentric cylinders: uniform and adaptive order results of the $L_2$ error expressed in the $lo{g}_{10}$ scale. $\mathrm{<mo>\; \#\; </mo>}{\text{P}}_4\left(\mathrm{<mi>\; \%\; </mi>}\right)$ denotes the percentage of the number of P₄ cells of all the cells obtained with the adaptive method using cellwise polynomial refinement from P₀ to P₄ order.

Figure 1. Lid-driven cavity flow at $Re={10}^4$ with polynomial adaptation: ${\text{P}}_{0\to 1}$, ${\text{P}}_{0\to 2}$, ${\text{P}}_{0\to 3}$, ${\text{P}}_{0\to 4}$ ( from left to right, from top to bottom).

(a) (b)

Figure 2. Lid-driven cavity flow at: (a) The baseline uniform reference solution; (b) order distribution of the adaptive solution.

the procedure of variable accuracy refinement. Interestingly, the polynomial refinement occurs primarily in the boundary layer and strong shear regions where energy cascade occurs that transfers energy from large-scale vortices to small-scale ones, as shown in Figure 2(b). A quantitative comparison is shown in Figure 3, where horizontal and vertical velocity profiles along the central lines of x and y axes are given. We notice that the uniformly and adaptively refined solutions are in good agreement with the reference solutions of Ghia [11] while the low-order solutions fail to match the results. In Table 2, we compare the

(a)(b)

Figure 3. Lid-driven cavity flow at: Comparison of the velocity profiles. (a) x-velocity profile along the line; (b) y-velocity profile along the line.

Table 2. Computational cost of the lid-driven cavity flow at. The CPU time is normalized by the one of. denotes the number of the k-order cells in the percentage of total cells. The adaptive solution achieves a nearly three-fold speedup compared to the P₄ one. The memory usage of the global Jacobian is expressed in Megabytes (M).

computational cost of all the runs, all the adaptation procedures offer cost reduction in terms of the total degree of freedoms (DOFs) and work unit. Especially, the solution that starts up at and ends up at offers a twofold reduction in DOFs and three times speedup in work units compared with the uniform solution. This case shows that the adaptive strategy is effective and accurate for viscous flows.

4.3. Caradonna-Tung Rotor in Hover

A real-world problem is considered in this case corresponding to the experimental model hover test conditions of Caradonna and Tung [7]. The experimental model consists of a two-bladed rigid rotor with rectangular planform blades with no twist or taper. The blades are made of NACA0012 airfoil sections with an aspect ratio of 6 as shown in Figure 4. The computational condition uses the case of tip Mach number, collective pitch, and the Reynolds number based on the blade tip speed and chord. More computational parameters are listed in Table 3. For this case, adaptive solutions are computed and compared with the experimental data [7]. The adaptive order distribution is shown in Figure 5(a) and the cell order is locally refined just around the trajectory of wake vortex shown in Figure 5(b). Comparisons of pressure coefficients at different span-wise locations defined by r/R are given in Figure 6, where r is the span-wise radius from the rotation axis and R is the tip radius defined in Table 3. The results demonstrate the feasibility of using p-adaptation framework of exponential time marching for 3-D problems.

Figure 4. The experimental configuration of Caradonna-Tung rotor in hover [7].

Table 3. Computational parameters of the Caradonna-Tung rotor in hover.

(a)(b)

Figure 5. Caradonna-Tung rotor in hover: (a), and Re = 1.92 × 10⁶; adaptive order distribution; (b) vorticity contour of the rotor wake flow.

Figure 6. Caradonna-Tung rotor in hover:, and; surface pressure coefficient (top); pressure coefficient comparisons at different span-wise, sectional locations r/R (computational: 96,222 cells; experimental [7]:); trajectory of vortex core comparison (right bottom).

5. Conclusions

An adaptive high-order DG framework has been developed with the PCEXP exponential time marching scheme. By using the adaptive control strategy, the performance of PCEXP scheme is shown to be enhanced in terms of computational cost and memory usage. The correctness of the variable order implementation of PCEXP is firstly validated in the case of rotating flow between two concentric cylinders. While uniform order solver can deliver formal convergence rates, the adaptive solution can even give a comparable error to the uniform solution. In the cavity flow case, it is found that the use of spectral decay error indicator is effective in capturing key flow structures such that the secondary corner vortices can be captured. Performance statistic shows that a three-fold reduction in CPU time and a two-fold reduction in DOFs are gained. Finally, the adaptive solver is applied to the Caradonna-Tung rotor in hover, where three-dimensional wake flows are adaptively captured with the coarse mesh. The results are in good agreement to the experimental data [7] thus demonstrates its applicability to practical 3D complex flows.

Acknowledgements

This work is funded by the National Natural Science Foundation of China (NSFC) under the Grant U1930402. The computational resources are provided by Beijing Computational Science Research Center (CSRC).

Nomenclatures

: specific heat ratio

: Prandtl number

: k-order DG polynomial

: cell number of the k-order cells

: quasi-2D characteristic cell size

: 3D characteristic cell size

: mesh dimension of coordinate z

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Li, S.-J. (2013) A Parallel Discontinuous Galerkin Method with Physical Orthogonal Basis on Curved Elements. Procedia Engineering, 61, 144-151. https://doi.org/10.1016/j.proeng.2013.07.107
[2]	Li, S.-J., Wang, Z.J., Ju, L. and Luo, L.-S. (2017) Explicit Large Time Stepping with a Second-Order Exponential Time Integrator Scheme for Unsteady and Steady Flows. 55th AIAA Aerospace Sciences Meeting, Grapevine, 9-13 January 2017, AIAA Paper 2017-0753. https://doi.org/10.2514/6.2017-0753
[3]	Li, S.-J., Luo, L.-S., Wang, Z.J. and Ju, L. (2018) An Exponential Time-Integrator Scheme for Steady and Unsteady Inviscid Flows. Journal of Computational Physics, 365, 206-225. https://doi.org/10.1016/j.jcp.2018.03.020
[4]	Li, S.-J. (2018) Efficient p-Multigrid Method Based on an Exponential Time Discretization for Compressible Steady Flows. arXiv:1807.0115.
[5]	Li, S.-J., Wang, Z.J., Ju, L. and Luo, L.-S. (2018) Fast Time Integration of Navier-Stokes Equations with an Exponential-Integrator Scheme. AIAA Aerospace Sciences Meeting, Kissimmee, 8-12 January 2018, AIAA Paper 2018-0369. https://doi.org/10.2514/6.2018-0369
[6]	Li, S.-J. and Ju, L. (2019) Exponential Time-Marching Method for the Unsteady Navier-Stokes Equations. AIAA Scitech 2019 Forum, San Diego, 7-11 January 2019, AIAA Paper 2019-0907. https://doi.org/10.2514/6.2019-0907
[7]	Caradonna, F.X. and Tung, C. (1981) Experimental and Analytical Studies of a Model Helicopter Rotor in Hover. NASA Technical Memorandum 1981-81232, NASA, Ames Research Center, Moffett Field.
[8]	Bassi, F. and Rebay, S. (1997) A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations. Journal of Computational Physics, 131, 267-279. https://doi.org/10.1006/jcph.1996.5572
[9]	Persson, P.-O. and Peraire, J. (2006) Sub-Cell Shock Capturing for Discontinuous Galerkin Methods. 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, 9-12 January 2006, AIAA Paper 2006-112. https://doi.org/10.2514/6.2006-112
[10]	Cagnone, J.S., et al. (2013) A p-Adaptive LCP Formulation for the Compressible Navier Stokes Equations. Journal of Computational Physics, 233, 324-338. https://doi.org/10.1016/j.jcp.2012.08.053
[11]	Ghia, U., Ghia, K.N. and Shin, C.T. (1982) High-Resolutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method. Journal of Computational Physics, 48, 387-411. https://doi.org/10.1016/0021-9991(82)90058-4