# An Implicit High-Order Preconditioned Flux Reconstruction Method for Low-Mach-Number Flow Simulation with Dynamic Meshes

## Abstract

A fully implicit high-order preconditioned flux reconstruction/correction procedure via reconstruction (FR/CPR) method is developed to solve the compressible Navier–Stokes equations at low Mach numbers. A dual-time stepping approach with the second-order backward differentiation formula (BDF2) is employed to ensure temporal accuracy for unsteady flow simulation. When dynamic meshes are used to handle moving/deforming domains, the geometric conservation law (GCL) is implicitly enforced to eliminate errors due to the resolution discrepancy between BDF2 and the spatial FR/CPR discretization. The large linear system resulting from the spatial and temporal discretizations is tackled with the restarted Generalized Minimal Residual (GMRES) solver in the PETSc (Portable, Extensible Toolkit for Scientific Computation) library. Through several benchmark steady and unsteady numerical tests, the preconditioned FR/CPR methods have demonstrated good convergence and accuracy for simulating flows at low Mach numbers. The new flow solver is then used to study the effects of Mach number on unsteady force generation over a plunging airfoil when operating in low-Mach-number flows. It is observed that weak compressibility has a significant impact on thrust generation but a negligible effect on lift generation of an oscillating airfoil.

## Key Words

Flux Reconstruction/Correction Procedure via Reconstruction; Low-Mach-Number Preconditioning; High-Order Methods; Dual-Time Stepping; Unsteady Flows; Dynamic Meshes; Weak Compressibility

## 1 Introduction

Analysis of unsteady fluid flows at low Mach numbers is an indispensable component in many mechanical engineering applications, such as design and control of wind turbines, micro air vehicles, autonomous underwater vehicles, and medical devices, just to name a few. There are typically two approaches to handle low-Mach-number flows. In the first one, fluids are considered incompressible, and their motion is governed by the incompressible Navier–Stokes equations. The pressure-correction method (e.g., SIMPLE (Semi-Implicit Method for Pressure Linked Equations) families [1], and the projection method [2]), and artificial compressibility method [3, 4, 5, 6, 7, 8] are among the most popular approaches for solving incompressible flows. An alternative is to directly solve the compressible Navier–Stokes equations at low Mach numbers. In this study, we pursue the second approach to solve low-Mach-number flows. We note that the compressible flow equations become stiff to solve numerically at low Mach numbers due to the large disparity between the flow (or entropy characteristic) velocity and the speed of sound. This disparity not only degrades the convergence of numerical solvers for the compressible Navier–Stokes equations, but also results in loss of numerical accuracy [9].

Various local preconditioning methods have been proposed to decrease the stiffness of the compressible Navier–Stokes equations at low Mach numbers. In general, after local preconditioning, the physical acoustic wave speeds are modified to be of the same order of magnitude as the flow speeds. Thus, the convergence rate of the flow solver can be greatly improved. Three major groups of preconditioning methods have been reported in the literatures [10]. Inspired by the artificial compressibility method developed by Chorin [11], Turkel adopted , where is the pressure, is the velocity vector, and is the change of entropy, as working variables in a series of work [9, 12, 13, 14]. Choi and Merkle [15], and Weiss and Smith [16] employed as working variables. Herein, is the temperature. Van Leer et al. [17, 18] developed a symmetric Van Leer–Lee–Roe preconditioner using as working variables, where is the fluid density, and is the speed of sound. We note that different sets of working variables have been assessed by Turkel et al. [13]. Hauke and Hughes in Ref. [19] pointed out that the conservative incompressible formulation is well defined only for the entropy variables and the primitive variables. As reported in Refs. [15, 16, 20, 21, 22], direct approximations of the gradients of primitive variables can be more accurate than those computed from gradients of conservative variables. Therefore, in this study, the primitive variables are selected as working variables, and the preconditioning approach by Weiss and Smith [16] is adopted.

Contributions. Preconditioned high-order methods have been employed to solve the compressible Navier–Stokes equations at low Mach numbers in past decades. For example, the discontinuous Galerkin (DG) method has been developed to solve low-Mach-number flows in Refs. [23, 24, 22, 25, 26]. In this paper, an implicit preconditioned FR/CPR method is developed for the first time to solve the unsteady compressible Navier–Stokes equations at low Mach numbers involving dynamic meshes. The dual-time stepping method with BDF2 is employed to achieve second-order accuracy in time for unsteady flow simulations. Within the framework of dual-time stepping, the GCL is enforced analytically following the approach proposed in Refs. [27, 28]. In a previous work [29], we have compared the numerical performance of the preconditioned FR/CPR formulation with the incompressible FR/CPR method with artificial compressibility [8]. In this study, we conduct thorough numerical experiments to test accuracy, convergence and efficiency of the preconditioned FR/CPR formulation for solving steady and unsteady flow problems.

Article overview. The remainder of the paper is organized as follows. In 2, the non-dimensionalized compressible Navier–Stokes equations are presented first ( 2.1). The preconditioning methodology for steady flow problems is then explained in 2.2. After that, preconditioning for unsteady flow problems on moving grids with implicit GCL enforcement is elaborated in 2.3. Numerical methods are introduced in 3. A brief review of the FR/CPR method is given in 3.1, and 3.2 formulates the concrete discrete forms of time marching methods. Boundary conditions for preconditioned systems are then discussed in 3.3. In 4, numerical results from both steady and unsteady simulations using the preconditioned FR/CPR methods are presented. The present study is summarized in 5.

## 2 Low-Mach-Number Preconditioning

### 2.1 Governing Equations

The non-dimensionalized compressible Navier–Stokes equations can be written as

(1) |

with the following non-dimensionalizations [30],

(2) |

Herein, is the variable in the physical domain, is the reference variable, are the non-dimensionalized conservative variables, is the non-dimensionalized density, is the non-dimensionalized velocity vector, is the non-dimensionalized pressure, is the non-dimensionalized temperature, is the non-dimensionalized total energy, the heat capacity ratio , and is the Mach number. The fluxes and consist of inviscid and viscous parts, i.e., and . The components in , , , can be found below,

(3) |

(4) |

Herein, the non-dimensionalized thermal conductivity , the Prandtl number and is the non-dimensionalized dynamic viscosity, which is treated as a constant in this work. is the Reynolds number of the reference flow at far fields. Non-dimensionalized viscous stresses are given by

(5) |

At low-Mach-number regions, the eigenvalues of the convection part of Eq. (1) are significantly different from each other. This disparity makes the hyperbolic system stiff to solve. Preconditioning methods replace the characteristic wave speeds with modified ones to decrease the stiffness of the hyperbolic system.

### 2.2 Preconditioning for Steady Flow Problems

A pseudo-transient continuation approach is used to handle steady flow problems. In this approach, on applying the chain-rule, Eq. (1) can be written as

(6) |

where is the pseudo-time, are the non-dimensionalized primitive variables, is the Jacobian matrix of the two sets of variables, i.e., the conservative and primitive variables, which can be expressed as

(7) |

Herein, is the non-dimensionalized specific total enthalpy defined as .

In the preconditioning method adopted here, is replaced with a matrix called preconditioning matrix, which can cluster the eigenvalues of the hyperbolic system. The preconditioning matrix is given by

(9) |

Note that in is substituted with in the preconditioning matrix . The definition of follows that in Ref. [16]

(10) |

where , is the local speed of sound, and is a free parameter related to the local Mach number , and the global cut-off Mach number which is usually chosen as the free stream Mach number .

After introducing the preconditioning matrix into the compressible Navier–Stokes equations, the governing equations can be re-expressed as

(11) |

The preconditioned Jacobian matrix related to inviscid fluxes in the normal direction of any surface becomes , where , and is the unit normal vector of the surface. The eigenvalues of this system are

(12) |

where

(13) |

For an ideal gas, . At low speed, when approaches zero, will approach . All the eigenvalues will then have the same magnitude as . Thus, the stiffness of the compressible Navier–Stokes equations is significantly decreased.

In the present study, when no moving grids are involved, the parameter in is defined as

(14) |

where is a free parameter. Typically, . Turkel suggested that can be up to in Ref. [14]. The global cut-off parameter [20, 14] is employed to prevent robustness deterioration instabilities near stagnation points.

### 2.3 Preconditioning for Unsteady Flow Problems with Dynamic Meshes

The preconditioning matrix in Eq. (6) destroys the temporal accuracy of Eq. (1). The dual-time stepping method [31] formulated as

(15) |

can be an intuitive solution for preconditioned unsteady flow problems. In this work, the BDF2 method is employed to discretize the physical time derivative term, thus achieving a second-order accuracy temporally. When dynamic meshes are involved, one can first transfer Eq. (15) from the physical domain into the computational domain as

(16) |

where

(17) |

The GCL of the transformation can be formulated as

(18) |

In this study, , thus,

(19) |

and

(20) |

Herein, is the grid velocity vector. Since grid velocities are only related to the physical time or , when marching along the pseudo-time , GCL can be enforced analytically to eliminate the error discrepancy between BDF2 and the spatial FR/CPR schemes following the approach in Ref. [27, 28]. Specifically, Eq. (16) can be written as

(21) |

As a result, the GCL error only depends on the accuracy of the spatial discretization. Finally, the governing equation in the physical domain reads

(22) |

When dynamic meshes are involved, the preconditioned Jacobian matrix in the normal direction of any surface for the convection terms becomes , where

(23) |

The eigenvalues of this system preserve the same form as those in Eq. (12) or Eq. (13) except a minor change, i.e., in Eq. (13) is modified to

(24) |

due to the presence of the grid velocity. The preconditioning parameter for dynamic meshes is still defined as [32]

(25) |

where is the Mach number calculated from .

## 3 Numerical Methods

### 3.1 Spatial Discretization Method

For completeness, a brief review of the FR/CPR method [33, 34, 35, 36] is presented in this section. In FR/CPR, fluxes can be divided into two parts, namely, local fluxes constructed directly from local solutions, and correction fluxes accounting for the differences between the local fluxes and common fluxes constructed from those on element interfaces. In the current study, the spatial domain is discretized into non-overlapping quadrilateral elements. When the FR/CPR method is employed for the spatial discretization, Eq. (22) can be reformulated as

(26) | ||||

For quadrilateral elements, and are given by

(27) |

where and are the common (or numerical) fluxes on the interfaces of standard elements. Herein, subscripts ‘L’ and ‘R’ stand for left and right edges of an element in a dimension by dimension sense. are correction functions which map the differences between the common and local fluxes on the element interfaces into the entire standard element. In this study, and are the right and left Radau polynomials, which can recover the nodal DG methods. Note that common fluxes are always calculated in the physical domain. Thus and can be expresses as

(28) |

where are the common fluxes in the normal direction of a physical surface. On element interfaces, the inviscid fluxes along can be written as

(29) |

Therefore, can be formulated as, taking the Rusanov approximate Riemann solver as an example,

(30) |

where is the local preconditioning matrix calculated by averaging primitive variables on elements interfaces, and superscripts ‘’ and ‘’ denote the left and right side of the interface, respectively. stands for the eigenvalue of the Jacobian matrix or .

The preconditioning does not affect the viscous terms and . To calculate the common viscous fluxes, we need to define common and common on element interfaces. Simply taking average of the primitive variables yields

(31) |

The common gradient is calculated following the the second approach of Bassi and Rebay (BR2) [37] as

(32) |

where and are the corrections to the gradients on the interfaces. See Ref. [38] for more information.

### 3.2 Time Marching Method

For steady flow problems, define and substitute it into Eq. (11) to obtain

(33) |

A backward Euler method is used to discretize the pseudo-time derivative, and the residual is linearized around at the pseudo-time step m. As a result, the linear system to be solved becomes

(34) |

where . The initial pseudo-time step size , and it will grow as the residual drops following

(35) |

where is the relative residual of pressure at the pseudo-time step , and . This approach is a simple modification of the switched evolution relaxation (SER) developed by Mulder and van Leer [39]. We have observed a good convergence with this aggressive pseudo-time step size adaptation strategy through numerical experiments.

For unsteady flow problems with the dual-time stepping method, define . On plugging it into Eq. (22), we obtain

(36) |

If denotes the current physical-time step, and denotes the pseudo-time step, then we have

(37) |

where and when , . To achieve a second-order temporal accuracy, BDF2 is used to discretize the derivative with respect to . Thus, Eq. (37) is discretized as

(38) |

Note that . After linearizing with respect to , Eq. (38) can be further written as

(39) | ||||

where . For unsteady flow problems tested here, is always fixed as for pseudo-time iterations.

### 3.3 Boundary Conditions

The characteristics of the inviscid parts of the compressible Navier–Stokes equation are changed due to the low-Mach-number preconditioning. No-reflection far-field boundary conditions or simplified far-field boundary conditions have been proposed in Refs. [13, 22]. However, the far field in this study is always located more than 100 characteristic lengths (e.g., the diameter of a cylinder, the chord length of an airfoil) away from the solid wall. Therefore, we simply use on boundaries at the far field, and enforce boundary conditions using the approximate Riemann solver. For inviscid flows, the slip wall boundary condition is employed. The adiabatic wall boundary condition is enforced for no-slip viscous walls. To ensure high-order accuracy on wall boundaries, mesh elements are employed to handle the curvatures of wall boundaries.

## 4 Numerical Results

Several benchmark tests have been conducted to verify the accuracy and convergence of the preconditioned FR/CPR method in this section. The error is used for the accuracy study. Specifically, for any solution variables on domain , the error is defined as

(40) |

Exclusively, the entropy error is defined as

(41) |

To monitor the convergence of nonlinear systems, such as Eqs. (33) and (36), two convergence criteria are set up for the pseudo-time iteration: one is for , denoted as ; the other is for the restarted GMRES linear solver, denoted as . For steady flow simulation using a pseudo-transient continuation (see Eq. (33)), the goal is to decrease the steady relative residual as much as possible with a reasonable criterion for ; for unsteady flow simulation (see Eq. (36)), the goal is to decrease the unsteady relative residual to a satisfactory level with a reasonable criterion for .

The drag coefficient , lift coefficient and thrust coefficient are defined as

(42) |

where is the area that the force acts on, and is the free stream velocity. The pressure is normalized as

(43) |

for better visualization, where and are the minimum and maximum pressure in the entire flow field.

A parallel code is developed to accelerate the simulations. The block Jacobi preconditioner is employed for parallelization, and ILU(0) is used as the local preconditioner for the restarted GMRES solver in PETSc. Each process possesses only one mesh block [41]. The dimension of the Krylov subspace is 150 for steady flow problems and 30 for unsteady ones. The Jacobian matrix is updated every pseudo-time iteration for steady flow problems and every two pseudo-time iterations for unsteady flow problems. All the steady flow problems are simulated with one process, and all the unsteady cases are solved with 16 processes.

### 4.1 Inviscid Flow over a Circular Cylinder

The solver developed in the present study is tested for the inviscid flow over a circular cylinder (smooth curvature) in this section. The physical domain of the problem is within a square of width . The diameter of the circular cylinder is one. The baseline mesh used in this section is illustrated in Figure 1, which has 20 elements in the normal direction of the wall and 24 elements in the circumferential direction. On splitting the baseline mesh once and twice, one can obtain the and meshes, respectively.

(a) | (b) |

Grid refinements for the third- () and fourth-order () FR/CPR schemes have been conducted. The free parameter is employed. The convergence criterion for the relative residual of the restarted GMRES solver is as suggested in Ref. [22]. In Figure 2, the pressure and contours near the cylinder in a free stream with are presented. No pressure-oscillation near the stagnation point is observed. Table 1 documents the drag coefficient and entropy error from a grid refinement study.

(a) | (b) |

(a) | (b) |

(a) | (b) |

(a) | (b) |

Mesh | ||||
---|---|---|---|---|

The convergence histories of the relative residual of pressures are illustrated in Figure 3. It is observed that can decrease to the level of in steady flow simulation at low Mach numbers. The histories of the pseudo-time step size variation are presented in Figure 4. An observation is that will quickly drop as the pseudo-time step size increases to large values. However, as shown in Figure 5, the iteration numbers needed for the restarted GMRES solver to converge (i.e., to meet the convergence criterion ) will significantly increase when , especially when fine meshes are used.

#### Effects of the Global Cut-off Parameter

We have tested three different values of , namely, , and to study the effects of on convergence and accuracy. The convergence criterion for the relative residual of the restarted GMRES solver is set as in these tests. For brevity, only the third-order FR/CPR scheme is considered on the mesh. The results of , and are presented in Table 2. The histories of relative residuals and the iteration numbers needed for each pseudo-time stepping are presented in Figure 6. When , the global cut-off value is very small. Due to the instability induced by the two stagnation points, the linear solver is not able to converge within the maximum iteration number (5000 for this case), and the simulation quickly diverged. Even though when , can drop to a smaller value using less pseudo-time iterations compared to that with , the accuracy does not get better, and more iterations are needed for the linear solver to converge in one pseudo-time step. The predicted and both increase due to the fact that when a larger is used, more numerical dissipation is added to the approximate Riemann solver.

diverged | diverged | diverged | |
---|---|---|---|

(a) | (b) |

### 4.2 Inviscid Flow over the NACA0012 Airfoil

In this section, the inviscid flow of and angle of attack (AOA) over the NACA0012 airfoil is studied to verify the performance of the preconditioned FR/CPR solver when there exists a singular point on the geometry. The profile of the airfoil is expressed as

(44) |

where . The mesh with 5,168 quadrilateral elements is illustrated in Figure 7.

(a) | (b) |

The polynomial degree is refined from to (i.e., from the to order of accuracy). As presented in Table 3, the prediction of becomes more accurate when the degree of the polynomial increases. With preconditioning, instabilities near wall boundaries are suppressed and the flow field in Figure 8 does not show any pressure oscillation near wall boundaries. We have also tested the impact of different criteria , i.e., , and , on the convergence.

In Ref. [22], a 1,792-element mesh is employed and the maximum iteration number of the restarted GMRES solver when is roughly 120. In this study, a relative denser mesh is used. We observe from Figure 9(b) that the maximum iteration number is about with the same as that used in Ref. [22]. According to the convergence comparison of different criteria for the restarted GMRES solver (see Figure 9(a)), when , both the third-order and fifth-order FR/CPR schemes can converge to correct solutions. However, the simulation using the fourth-order FR/CPR scheme will immediately diverge after a few pseudo-time iterations. A smaller for the linear solver can guarantee convergence, but more iterations of the restarted GMRES solver are needed (see Figure 9(b)). Only negligible differences are observed on and among the converged results of different s, which are not presented for brevity. Therefore, a proper choice of is a trade-off between robustness and efficiency.

(a) | (b) |

(a) |

(b) |

### 4.3 Isentropic Vortex Propagation with Prescribed Grid Motion

The isentropic vortex propagation case is employed to verify both the spatial convergence and temporal convergence of the preconditioned FR/CPR methods for unsteady flow simulation in this section. The isentropic vortex propagation case depicts the superposition of an inviscid uniform flow and an irrotational vortex. The vortex can be regarded as a perturbation added onto the uniform flow. The free stream flow is of and the gas constant for this case. The perturbation is defined as [41]

(45) |

where and are parameters that define the vortex strength. is the distance to the center of the vortex . The vortex is propagated in a periodic domain . A prescribed grid motion is employed to validate the moving grid algorithm. The prescribed grid motion is defined as

(46) |

where are the coordinates at , i.e., coordinates of the uniformly divided grid, and , , , , . Periodic boundary conditions are applied to all boundaries.

The time step size is refined from to , where to validate the accuracy of BDF2. A fifth-order FR/CPR scheme and a uniformly divided mesh are employed to ensure that the error is dominated by the temporal discretization. We require the relative residual of the pseudo-time iterations to drop as low as possible within 100 pseudo-time iterations. The of the restarted GMRES solver is with the maximum iteration number 100. We only simulate this case until when the mesh deformation is maximized to demonstrate the convergence of BDF2 with GCL. An illustration of the pressure contour and the mesh deformation can be found in Figure 10. As observed from Figure 11(a), the optimal convergence rate of BDF2 is achieved for both and .

A grid refinement study is also conducted on a uniformly-divided mesh set with , , and elements, respectively. The time step size is to guarantee that the error is dominated by the spatial discretization. The error of pressure and velocity are presented in Figure 11(b). We observe that for both and , the numerical orders of a degree solution construction are close to with and without dynamic meshes. Note that the optimal convergence rate is . The observed order reduction is reasonable since the Rusanov approximate Riemann solver is employed to calculate the common inviscid fluxes [42]. Hence, the preconditioned FR/CPR methods can demonstrate good convergence for unsteady low-Mach-number flow simulations on dynamic meshes.

(a) | (b) |

### 4.4 Laminar Flow over a Plunging NACA0012 Airfoil

An order of accuracy refinement study is firstly carried out for a steady case, i.e., the viscous flow of , and AOA over the stationary NACA0012 airfoil. The mesh is the same as that in Section 4.2 (see Figure 7). The drag coefficients of the third-, fourth-, and fifth-order FR/CPR schemes are , and , respectively. Therefore, we adopt the fourth-order FR/CPR scheme to simulate laminar flows over a plunging NACA0012 airfoil. The plunging motion is defined as , where is the plunge amplitude, and is the plunge frequency. The reduced frequency is defined as , and dimensionless height , where is the chord length of the airfoil. In this study, , [8]. The whole mesh is considered as an oscillating rigid body. At the far field, , and the approximate Riemann solver is used to enforce boundary conditions.

The of the restarted GMRES solver is with the maximum iteration number 100. For the sake of efficiency, we would require the drop by certain orders instead of as low as possible. In this case, three criteria for of the pseudo-time marching have been tested, i.e., , , and , to examine if the force predictions are sensitive to the convergence criterion of . The time step size is set to , and is set to . All simulations end at . In Figure 12, the vortex shedding process during one typical plunging period is presented. The thrust coefficient and lift coefficient of the tenth period are displayed in Figure 13, and the time-averaged thrust coefficient , root mean square of the lift coefficient , and maximum lift coefficient are documented in Table 4. From Figure 13 and Table 4, we observe that a moderate tolerance () for is sufficient to maintain the accuracy of force prediction. However, a smaller criterion means more iterations.

(a) | (b) |

(c) | (d) |

(e) | (f) |

(g) | (h) |

(a) | (b) |

0.0589 | 2.352 | 3.685 | |

0.0589 | 2.351 | 3.684 | |

0.0589 | 2.351 | 3.684 |

#### Effects of the Mach Number on Force Generation

We further conducted a numerical study to test how the Mach number affects force prediction. Different Mach numbers, i.e., , , and , are studied here. A smaller time step size is employed, i.e., , to ensure that the time integration errors are negligible. is adopted as the convergence criterion for pseudo-time iterations. The and of the tenth period are illustrated in Figure 14. , and are documented in Table 5.

(a) | (b) |

0.0578 | 2.340 | 3.685 | |

0.0578 | 2.341 | 3.692 | |

0.0539 | 2.361 | 3.774 | |

0.0370 | 2.369 | 3.608 | |

Incompressible [8] | 0.0557 | 2.348 | - |

It is observed that when is less than , the compressibility is almost negligible. The force histories for and coincide with each other. When is larger than , phase shift is observed in the force histories. In general, when is increased to a certain level (e.g., in this case), the time-averaged thrust coefficient will decrease when increases; at the same time, the root mean square of the lift coefficient and the maximum lift coefficient will vary slightly. For example, when , the reduction of is over compared to that at ;