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

TL;DR

This paper introduces a fully implicit high-order preconditioned flux reconstruction method for low-Mach-number flow simulations with dynamic meshes, ensuring accuracy and efficiency in unsteady flow analysis.

Contribution

It develops a novel implicit high-order preconditioned FR/CPR method combined with dual-time stepping and GCL enforcement for low-Mach flows with moving meshes.

Findings

Demonstrates good convergence and accuracy in benchmark tests.

Shows significant impact of weak compressibility on thrust generation.

Validates the method's effectiveness for unsteady flow simulations.

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 resulted 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.

Tables

Table 1. Table 1: The grid refinement study for the inviscid flow over a circular cylinder at M ​ a ∞ = 10 − 3 𝑀 subscript 𝑎 superscript 10 3 Ma_{\infty}=10^{-3} with t ​ o ​ l r ​ e ​ s = 10 − 6 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑠 superscript 10 6 tol_{res}=10^{-6} .

	$3^{rd}$		$4^{th}$
Mesh	$C_d$	$E_s$	$C_d$	$E_s$
$20\times 24$	$5.2720\times {10}^{-4}$	$7.0302\times {10}^{-12}$	$9.3369\times {10}^{-5}$	$2.6179\times {10}^{-13}$
$40\times 48$	$1.7897\times {10}^{-4}$	$6.9203\times {10}^{-13}$	$6.7059\times {10}^{-6}$	$2.0882\times {10}^{-14}$
$80\times 96$	$2.8318\times {10}^{-5}$	$7.3010\times {10}^{-14}$	$4.5020\times {10}^{-7}$	$3.2653\times {10}^{-15}$

Table 2. Table 2: C l subscript 𝐶 𝑙 C_{l} , C d subscript 𝐶 𝑑 C_{d} and E s subscript 𝐸 𝑠 E_{s} for different κ ′ superscript 𝜅 ′ \kappa^{\prime} s when the third-order FR/CPR scheme is used to simulate the inviscid flow over a circular cylinder at M ​ a ∞ = 10 − 3 𝑀 subscript 𝑎 superscript 10 3 Ma_{\infty}=10^{-3} on the 80 × 96 80 96 80\times 96 mesh with t ​ o ​ l r ​ e ​ s = 10 − 9 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑠 superscript 10 9 tol_{res}=10^{-9} .

$\kappa$	$C_d$	$C_l$	$E_s$
$\sqrt{0.2}$	diverged	diverged	diverged
$1$	$2.8319\times {10}^{-5}$	$6.6839\times {10}^{-7}$	$7.3013\times {10}^{-14}$
$\sqrt{5}$	$5.1692\times {10}^{-5}$	$1.2015\times {10}^{-6}$	$1.3620\times {10}^{-13}$

Table 3. Table 3: The order of accuracy refinement study for the inviscid flow over the NACA0012 airfoil at M ​ a ∞ = 10 − 3 𝑀 subscript 𝑎 superscript 10 3 Ma_{\infty}=10^{-3} with t ​ o ​ l r ​ e ​ s = 10 − 6 𝑡 𝑜 subscript 𝑙 𝑟 𝑒 𝑠 superscript 10 6 tol_{res}=10^{-6} .

	$3^{rd}$	$4^{th}$	$5^{th}$
$C_d$	$7.5664\times {10}^{-6}$	$2.7816\times {10}^{-6}$	$2.1640\times {10}^{-6}$
$E_s$	$2.4161\times {10}^{-12}$	$4.0361\times {10}^{-13}$	$1.8262\times {10}^{-13}$

Table 4. Table 4: Study on the effects of the convergence criterion t ​ o ​ l p ​ s ​ e ​ u ​ d ​ o 𝑡 𝑜 subscript 𝑙 𝑝 𝑠 𝑒 𝑢 𝑑 𝑜 tol_{pseudo} on force prediction of a plunging NACA0012 airfoil, Δ ​ t = T / 100 Δ 𝑡 𝑇 100 \Delta t=T/100 , and M ​ a ∞ = 10 − 3 𝑀 subscript 𝑎 superscript 10 3 Ma_{\infty}=10^{-3} .

$to{l}_{pseudo}$	$\overline{C_t}$	$C_{l,rms}$	$C_{l,max}$
${10}^{-2}$	0.0589	2.352	3.685
${10}^{-4}$	0.0589	2.351	3.684
${10}^{-6}$	0.0589	2.351	3.684

Table 5. Table 5: Study on the effects of the Mach number on force prediction of a plunging NACA0012 airfoil, Δ ​ t = T / 1000 Δ 𝑡 𝑇 1000 \Delta t=T/1000 .

$Ma$	$\overline{C_t}$	$C_{l,rms}$	$C_{l,max}$
${10}^{-3}$	0.0578	2.340	3.685
${10}^{-2}$	0.0578	2.341	3.692
${10}^{-1}$	0.0539	2.361	3.774
$0.2$	0.0370	2.369	3.608
Incompressible [8]	0.0557	2.348	-

Equations

$$\frac{\partial {\mathbf{q}}^c}{\partial t}+\frac{\partial \mathbf{f}}{\partial x}+\frac{\partial \mathbf{g}}{\partial y}=0,$$

$$\{\begin{array}{l}t=\bar{t}/(L_{ref}/u_{ref}),x=\bar{x}/L_{ref},\\ \rho =\bar{\rho }/{\rho }_{ref},\mathbf{v}=\bar{\mathbf{v}}/u_{ref},p=\bar{p}/{\rho }_{ref}{u}_{ref}^2,T=\bar{T}/T_{ref},\\ \mu =\bar{\mu }/{\mu }_{ref},R=1/(\gamma M{a}_{ref}^2),C_p=1/((\gamma -1)M{a}_{ref}^2).\end{array}$$

$${\mathbf{f}}_{inv}=\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ (E+p)u\end{array}\right),{\mathbf{g}}_{inv}=\left(\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ (E+p)v\end{array}\right),$$

$${\mathbf{f}}_v=\frac{1}{R{e}_{ref}}\left(\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{xy}\\ u{\tau }_{xx}+v{\tau }_{xy}+k\frac{\partial T}{\partial x}\end{array}\right),{\mathbf{g}}_v=\frac{1}{R{e}_{ref}}\left(\begin{array}{c}0\\ {\tau }_{xy}\\ {\tau }_{yy}\\ u{\tau }_{xy}+v{\tau }_{yy}+k\frac{\partial T}{\partial y}\end{array}\right).$$

$$\{\begin{array}{l}{\tau }_{xx}=\mu (\frac{4}{3}\frac{\partial u}{\partial x}-\frac{2}{3}\frac{\partial v}{\partial y}),\\ {\tau }_{yy}=\mu (\frac{4}{3}\frac{\partial v}{\partial y}-\frac{2}{3}\frac{\partial u}{\partial x}),\\ {\tau }_{xy}={\tau }_{yx}=\mu (\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}).\end{array}$$

$$\frac{\partial {\mathbf{q}}^c}{\partial {\mathbf{q}}^p}\frac{\partial {\mathbf{q}}^p}{\partial \tau }+\frac{\partial \mathbf{f}}{\partial x}+\frac{\partial \mathbf{g}}{\partial y}=0,$$

$$\mathbf{M}=\frac{\partial {\mathbf{q}}^c}{\partial {\mathbf{q}}^p}=\left(\begin{array}{cccc}{\rho }_p& 0& 0& {\rho }_T\\ {\rho }_{p}u& \rho & 0& {\rho }_{T}u\\ {\rho }_{p}v& 0& \rho & {\rho }_{T}v\\ {\rho }_{p}h-1& \rho u& \rho v& {\rho }_{T}h+\rho C_p\end{array}\right).$$

$$\rho_{p}=\frac{\partial\rho}{\partial p}\Bigr{|}_{T}=\frac{1}{RT},\rho_{T}=\frac{\partial\rho}{\partial T}\Bigr{|}_{p}=-\frac{\rho}{T}.$$

$$\Gamma =\left(\begin{array}{cccc}\Theta & 0& 0& {\rho }_T\\ \Theta u& \rho & 0& {\rho }_{T}u\\ \Theta v& 0& \rho & {\rho }_{T}v\\ \Theta h-1& \rho u& \rho v& {\rho }_{T}h+\rho C_p\end{array}\right).$$

$$\varTheta=\Bigl{(}\frac{1}{U^{2}_{r}}-\frac{\rho_{T}}{\rho C_{p}}\Bigr{)},$$

$$\Gamma \frac{\partial {\mathbf{q}}^p}{\partial \tau }+\frac{\partial \mathbf{f}}{\partial x}+\frac{\partial \mathbf{g}}{\partial y}=0.$$

$$\{\begin{array}{l}{\lambda }_1={\lambda }_2=u_n\\ {\lambda }_3=u_n^{\prime }+c^{\prime }\\ {\lambda }_4=u_n^{\prime }-c^{\prime },\end{array}$$

$$\{\begin{array}{l}{u}_n=\mathbf{v}\cdot \mathbf{n}\\ u_n^{\prime }=u_n(1-\alpha )\\ c^{\prime }=\sqrt{{\alpha }^2{u}_n^2+U_r^2}\\ \alpha =(1-\beta U_r^2)/2\\ \beta =\left({\rho }_p+\frac{{\rho }_T}{\rho C_p}\right).\end{array}$$

$$ϵ=\min (1,\max (\kappa M{a}_{\infty },Ma))$$

$$\Gamma \frac{\partial {\mathbf{q}}^p}{\partial \tau }+\frac{\partial {\mathbf{q}}^c}{\partial t}+\frac{\partial \mathbf{f}}{\partial x}+\frac{\partial \mathbf{g}}{\partial y}=0,$$

$$\Gamma |\mathbf{J}|\frac{\partial {\mathbf{q}}^p}{\partial \tau }+\frac{\partial |\mathbf{J}|{\mathbf{q}}^c}{\partial t^{\ast }}+\frac{\partial \mathbf{F}}{\partial \xi }+\frac{\partial \mathbf{G}}{\partial \eta }=0,$$

$$\{\begin{array}{l}\mathbf{F}=\left|\mathbf{J}\right|\left({\mathbf{q}}^c{\xi }_t+\mathbf{f}{\xi }_x+\mathbf{g}{\xi }_y\right)\\ \mathbf{G}=\left|\mathbf{J}\right|\left({\mathbf{q}}^c{\eta }_t+\mathbf{f}{\eta }_x+\mathbf{g}{\eta }_y\right).\end{array}$$

$$\{\begin{array}{l}\frac{\partial }{\partial \xi }(|\mathbf{J}|{\xi }_x)+\frac{\partial }{\partial \eta }(|\mathbf{J}|{\eta }_x)=0\\ \frac{\partial }{\partial \xi }(|\mathbf{J}|{\xi }_y)+\frac{\partial }{\partial \eta }(|\mathbf{J}|{\eta }_y)=0\\ \frac{\partial |\mathbf{J}|}{\partial t^{\ast }}+\frac{\partial }{\partial \xi }(|\mathbf{J}|{\xi }_t)+\frac{\partial }{\partial \eta }(|\mathbf{J}|{\eta }_t)=0.\end{array}$$

$$\mathbf{J}=\frac{\partial (x,y,t)}{\partial (\xi ,\eta ,t^{\ast })}=\left(\begin{array}{ccc}{x}_{\xi }& x_{\eta }& x_{t^{\ast }}\\ y_{\xi }& y_{\eta }& y_{t^{\ast }}\\ 0& 0& 1\end{array}\right)$$

$$|\mathbf{J}|=x_{\xi }{y}_{\eta }-x_{\eta }{y}_{\xi },\{\begin{array}{l}{\xi }_t=-{\mathbf{v}}_g\cdot \nabla \xi ,\\ {\eta }_t=-{\mathbf{v}}_g\cdot \nabla \eta .\end{array}$$

$$\Gamma |\mathbf{J}|\frac{\partial {\mathbf{q}}^p}{\partial \tau }+|\mathbf{J}|\frac{\partial {\mathbf{q}}^c}{\partial t^{\ast }}-{\mathbf{q}}^c\left(\frac{\partial |\mathbf{J}|{\xi }_t}{\partial \xi }+\frac{\partial |\mathbf{J}|{\eta }_t}{\partial \eta }\right)+\frac{\partial \mathbf{F}}{\partial \xi }+\frac{\partial \mathbf{G}}{\partial \eta }=0,$$

$$\Gamma \frac{\partial {\mathbf{q}}^p}{\partial \tau }+\frac{\partial {\mathbf{q}}^c}{\partial t^{\ast }}+\frac{\partial \mathbf{f}}{\partial x}+\frac{\partial \mathbf{g}}{\partial y}-u_g\frac{\partial {\mathbf{q}}^c}{\partial x}-v_g\frac{\partial {\mathbf{q}}^c}{\partial y}=0.$$

$${\mathbf{A}}_{\mathbf{n},mov}=n_x\frac{\partial {\mathbf{f}}_{inv}}{\partial {\mathbf{q}}^p}+n_y\frac{\partial {\mathbf{g}}_{inv}}{\partial {\mathbf{q}}^p}-u_g\frac{\partial {\mathbf{q}}^c}{\partial {\mathbf{q}}^p}{n}_x-v_g\frac{\partial {\mathbf{q}}^c}{\partial {\mathbf{q}}^p}{n}_y.$$

$$u_n=(\mathbf{v}-{\mathbf{v}}_g)\cdot \mathbf{n},$$

$$ϵ=\min (1,\max (\kappa M{a}_{\infty },M{a}_{mov})),$$

$$\Gamma \frac{\partial {\mathbf{q}}^p}{\partial \tau }$$

$$+\frac{1}{\left|\mathbf{J}\right|}\left(\frac{\partial {\mathbf{F}}^{cor}}{\partial \xi }+\frac{\partial {\mathbf{G}}^{cor}}{\partial \eta }\right)+\frac{\partial {\mathbf{q}}^c}{\partial t^{\ast }}=0.$$

$$\{\begin{array}{l}\begin{array}{rl}{\mathbf{F}}^{cor}(\xi ,\eta )& =\left({\mathbf{F}}_L^{num}(\eta )-{\mathbf{F}}^{loc}({\xi }_L,\eta )\right)g_L^c(\xi )\\ & +\left({\mathbf{F}}_R^{num}(\eta )-{\mathbf{F}}^{loc}({\xi }_R,\eta )\right)g_R^c(\xi ),\end{array}\\ \begin{array}{rl}{\mathbf{G}}^{cor}(\xi ,\eta )& =\left({\mathbf{G}}_L^{num}(\xi )-{\mathbf{G}}^{loc}(\xi ,{\eta }_L)\right)g_L^c(\eta )\\ & +\left({\mathbf{G}}_R^{num}(\xi )-{\mathbf{G}}^{loc}(\xi ,{\eta }_R)\right)g_R^c(\eta ),\end{array}\end{array}$$

$$\{\begin{array}{l}{\mathbf{F}}^{num}=|\mathbf{J}||\nabla \xi |{\mathbf{f}}_{\mathbf{n}}^{com}\mathrm{sign}(\mathbf{n}\cdot \nabla \xi ),\\ {\mathbf{G}}^{num}=|\mathbf{J}||\nabla \eta |{\mathbf{f}}_{\mathbf{n}}^{com}\mathrm{sign}(\mathbf{n}\cdot \nabla \eta ),\end{array}$$

$${\mathbf{f}}_{\mathbf{n},inv}=\mathbf{f}{n}_x+\mathbf{g}{n}_y-\left({\mathbf{v}}_g\cdot \mathbf{n}\right){\mathbf{q}}^c.$$

$${\mathbf{f}}_{\mathbf{n},inv}^{com}=\frac{1}{2}\left({\mathbf{f}}_{\mathbf{n},inv}^{+}+{\mathbf{f}}_{\mathbf{n},inv}^{-}\right)-\frac{1}{2}{\left|\lambda \right|}_{max}\tilde{\Gamma }\left({\mathbf{q}}^{p,+}-{\mathbf{q}}^{p,-}\right).$$

Full text

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

Lai Wang  and Meilin Yu Graduate student. Email: bx58858@umbc.eduAssistant Professor. Corresponding Author. Email: mlyu@umbc.edu Department of Mechanical Engineering,

University of Maryland, Baltimore County, Baltimore, MD 21250

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 resulted 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 [3], Turkel adopted $(p,\boldsymbol{v},dS)$, where $p$ is the pressure, $\boldsymbol{v}$ is the velocity vector, and $dS$ is the change of entropy, as working variables in a series of work [9, 11, 12, 13]. Choi and Merkle  [14], and Weiss and Smith [15] employed $(p,\boldsymbol{v},T)$ as working variables. Herein, $T$ is the temperature. Van Leer et al. [16, 17] developed a symmetric Van Leer–Lee–Roe preconditioner using $(dp/(\rho c),d\boldsymbol{v},dS)$ as working variables, where $\rho$ is the fluid density, and $c$ is the speed of sound. We note that different sets of working variables have been assessed by Turkel et al. [12]. Hauke and Hughes in Ref. [18] pointed out that the conservative incompressible formulation is well defined only for the entropy variables and the primitive variables. As reported in Refs. [14, 15, 19, 20, 21], 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 [15] is adopted. 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. [22, 23, 21, 24, 25]. In addition to local preconditioning methods, we note that researchers have also proposed to directly modify the approximate Riemann solvers [26, 27, 28, 29] to balance numerical dissipation on different characteristic waves. A comprehensive comparison of several approximate Riemann solvers at low Mach numbers in the context of high-order flux reconstruction methods can be found in Ref. [30].

Contributions. 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. [31, 32]. In a previous work [33], 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 Section 2, the non-dimensionalized compressible Navier–Stokes equations are presented first ( Section 2.1). The preconditioning methodology for steady flow problems is then explained in Section 2.2. After that, preconditioning for unsteady flow problems on moving grids with implicit GCL enforcement is elaborated in Section 2.3. Numerical methods are introduced in Section 3. A brief review of the FR/CPR method is given in Section 3.1, and Section 3.2 formulates the concrete discrete forms of time marching methods. Boundary conditions for preconditioned systems are then discussed in Section 3.3. In Section 4, numerical results from both steady and unsteady simulations using the preconditioned FR/CPR methods are presented. The present study is summarized in Section 5.

2 Low-Mach-Number Preconditioning

2.1 Governing Equations

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

[TABLE]

with the following non-dimensionalizations [34],

[TABLE]

Herein, $(\bar{*})$ is the variable in the physical domain, $(*_{ref})$ is the reference variable, $\boldsymbol{q}^{c}=(\rho,\rho u,\rho v,E)^{\top}$ are the non-dimensionalized conservative variables, $\rho$ is the non-dimensionalized density, $\boldsymbol{v}=(u,v)^{\top}$ is the non-dimensionalized velocity vector, $p=\rho RT$ is the non-dimensionalized pressure, $T$ is the non-dimensionalized temperature, $E=\frac{p}{\gamma-1}+\frac{1}{2}\rho(u^{2}+v^{2})$ is the non-dimensionalized total energy, the heat capacity ratio $\gamma=1.4$, and $Ma$ is the Mach number. The fluxes $\boldsymbol{f}$ and $\boldsymbol{g}$ consist of inviscid and viscous parts, i.e., $\boldsymbol{f}=\boldsymbol{f}_{inv}-\boldsymbol{f}_{vis}$ and $\boldsymbol{g}=\boldsymbol{g}_{inv}-\boldsymbol{g}_{vis}$. The components in $\boldsymbol{f}_{inv}$, $\boldsymbol{g}_{inv}$, $\boldsymbol{f}_{vis}$, $\boldsymbol{g}_{vis}$ can be found below,

[TABLE]

[TABLE]

Herein, the non-dimensionalized thermal conductivity $k=\frac{\mu C_{p}}{Pr}$, the Prandtl number $Pr=0.72$ and $\mu$ is the non-dimensionalized dynamic viscosity, which is treated as a constant in this work. $Re_{ref}=\rho_{ref}u_{ref}L_{ref}/\mu_{ref}$ is the Reynolds number of the reference flow at far fields. Non-dimensionalized viscous stresses are given by

[TABLE]

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

[TABLE]

where $\tau$ is the pseudo-time, $\boldsymbol{q}^{p}=(p,u,v,T)^{\top}$ are the non-dimensionalized primitive variables, $\boldsymbol{M}=\frac{\partial\boldsymbol{q}^{c}}{\partial\boldsymbol{q}^{p}}$ is the Jacobian matrix of the two sets of variables, i.e., the conservative and primitive variables, which can be expressed as

[TABLE]

Herein, $h$ is the non-dimensionalized specific total enthalpy defined as $h=C_{p}T+\frac{1}{2}(u^{2}+v^{2})$.

Recall that the state equation $p=\rho RT$. Thus the derivatives in Eq. (7) can be calculated as [21]

[TABLE]

In the preconditioning method adopted here, $\boldsymbol{M}$ is replaced with a matrix $\boldsymbol{\varGamma}$ called preconditioning matrix, which can cluster the eigenvalues of the hyperbolic system. The preconditioning matrix is given by

[TABLE]

Note that $\rho_{p}=\frac{\partial\rho}{\partial p}$ in $\boldsymbol{M}$ is substituted with $\varTheta$ in the preconditioning matrix $\boldsymbol{\varGamma}$. The definition of $\varTheta$ follows that in Ref. [15]

[TABLE]

where $U_{r}=\epsilon c$, $c$ is the local speed of sound, and $\epsilon$ is a free parameter related to the local Mach number $Ma$, and the global cut-off Mach number which is usually chosen as the free stream Mach number $Ma_{\infty}$.

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

[TABLE]

The preconditioned Jacobian matrix related to inviscid fluxes in the normal direction of any surface becomes $\boldsymbol{\varGamma}^{-1}\boldsymbol{A}_{\boldsymbol{n}}$, where $\boldsymbol{A}_{\boldsymbol{n}}=n_{x}\frac{\partial\boldsymbol{f}_{inv}}{\partial\boldsymbol{q}^{p}}+n_{y}\frac{\partial\boldsymbol{g}_{inv}}{\partial\boldsymbol{q}^{p}}$, and $\boldsymbol{n}=(n_{x},n_{y})^{\top}$ is the unit normal vector of the surface. The eigenvalues of the preconditioned Jacobian matrix in the normal direction are [15]

[TABLE]

where

[TABLE]

For an ideal gas, $\beta=(\gamma RT)^{-1}=1/c^{2}$. At low speed, when $U_{r}$ approaches zero, $\alpha$ will approach $\frac{1}{2}$. All the eigenvalues will then have the same magnitude as $u_{n}$. Thus, the stiffness of the compressible Navier–Stokes equations is significantly decreased.

In the present study, when no moving grids are involved, the parameter $\epsilon$ in $U_{r}=\epsilon c$ is defined as

[TABLE]

where $\kappa$ is a free parameter. Typically, $\kappa=1$. Turkel suggested that $\kappa$ can be up to $\sqrt{3}\sim\sqrt{5}$ in Ref. [13]. The global cut-off parameter $\kappa Ma_{\infty}$ [19, 13] is employed to prevent robustness deterioration instabilities near stagnation points.

2.3 Preconditioning for Unsteady Flow Problems with Dynamic Meshes

The preconditioning matrix $\boldsymbol{\varGamma}$ in Eq. (6) destroys the temporal accuracy of Eq. (1). The dual-time stepping method [35] formulated as

[TABLE]

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, in order to facilitate the incorporation of GCL, one can transfer Eq. (15) from the physical domain $(t,x,y)$ into the computational domain $(t^{*},\xi,\eta)$ as

[TABLE]

where

[TABLE]

The GCL of the transformation can be formulated as

[TABLE]

In this study, $t^{*}=t$, thus,

[TABLE]

and

[TABLE]

Herein, $\boldsymbol{v}_{g}=(u_{g},v_{g})^{\top}$ is the grid velocity vector. Since grid velocities are only related to the physical time $t$ or $t^{*}$, when marching along the pseudo-time $\tau$, GCL can be enforced analytically to eliminate the error discrepancy between BDF2 and the spatial FR/CPR schemes following the approach in Refs. [8, 31, 32]. Specifically, on substituting the last equality in Eq. (18) into Eq. (16), Eq. (16) can be written as

[TABLE]

One can further substitute Eq. (17) into the above equation and transfer it back onto the physical domain using the first two formulae in Eq. (18), which reads

[TABLE]

Therefore, the GCL is analytically enforced. The error are only related to the numerical discretizations. Interested readers are referred to our previous work [8] for a more detailed derivation. There the governing equations for incompressible flows on dynamic meshes are solved with FR/CPR via the artificial compressibility approach.

When dynamic meshes are involved, the preconditioned Jacobian matrix in the normal direction of any surface for the convection terms becomes $\boldsymbol{\varGamma}^{-1}\boldsymbol{A}_{\boldsymbol{n},mov}$, where

[TABLE]

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

[TABLE]

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

[TABLE]

where $Ma_{mov}$ is the Mach number calculated from $\boldsymbol{v}-\boldsymbol{v}_{g}$.

3 Numerical Methods

3.1 Spatial Discretization Method

For completeness, a brief review of the FR/CPR method [37, 38, 39, 40] 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

[TABLE]

For quadrilateral elements, $F^{cor}$ and $G^{cor}$ are given by

[TABLE]

where $\boldsymbol{F}^{num}$ and $\boldsymbol{G}^{num}$ 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. $g_{L}^{c},g_{R}^{c}$ 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, $g_{L}^{c}$ and $g_{R}^{c}$ 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 $\boldsymbol{F}^{num}$ and $\boldsymbol{G}^{num}$ can be expresses as

[TABLE]

where $\boldsymbol{f}^{com}_{\boldsymbol{n}}$ are the common fluxes in the normal direction $\boldsymbol{n}$ of a physical surface. On element interfaces, the inviscid fluxes along $\boldsymbol{n}$ can be written as

[TABLE]

Therefore, $\boldsymbol{f}^{com}_{\boldsymbol{n},inv}$ can be formulated as, taking the Rusanov approximate Riemann solver as an example,

[TABLE]

where $\boldsymbol{\widetilde{\varGamma}}$ is the local preconditioning matrix calculated by averaging primitive variables on elements interfaces, and superscripts ‘$+$’ and ‘$-$’ denote the right and left side of the interface, respectively. $\lambda$ stands for the eigenvalue of the Jacobian matrix $\boldsymbol{\varGamma}^{-1}\boldsymbol{A}_{\boldsymbol{n}}$ or $\boldsymbol{\varGamma}^{-1}\boldsymbol{A}_{\boldsymbol{n},mov}$.

The preconditioning does not affect the viscous terms and $\boldsymbol{f}_{n,vis}^{com}=\boldsymbol{f}_{vis}(\boldsymbol{q}^{p,+},\nabla\boldsymbol{q}^{p,+},\boldsymbol{q}^{p,-},\nabla\boldsymbol{q}^{p,-})$. To calculate the common viscous fluxes, we need to define common $\boldsymbol{q}^{p,com}$ and common $\nabla\boldsymbol{q}^{p,com}$ on element interfaces. Simply taking average of the primitive variables yields

[TABLE]

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

[TABLE]

where $\boldsymbol{r}^{+}$ and $\boldsymbol{r}^{-}$ are the corrections to the gradients on the interfaces. See Ref. [42] for more information.

3.2 Time Marching Method

For steady flow problems, define $\boldsymbol{R}=-\left(\frac{\partial\boldsymbol{f}}{\partial x}+\frac{\partial\boldsymbol{g}}{\partial y}\right)$ and substitute it into Eq. (11) to obtain

[TABLE]

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

[TABLE]

where $\Delta\left(\boldsymbol{q}^{p}\right)^{m}=\left(\boldsymbol{q}^{p}\right)^{m+1}-\left(\boldsymbol{q}^{p}\right)^{m}$. A widely used approach to update $\Delta\tau$ is the switched evolution relaxation (SER) originally developed by Mulder and van Leer [43], which reads

[TABLE]

where $Res_{m}/Res_{0}$ is the relative residual of pressure at the pseudo-time step $m$, and $r_{ser}$ is a free parameter to control the adaptation rate of $\Delta\tau$. If not specifically mentioned, $\Delta\tau_{min}=\Delta\tau_{0}=0.01$, $\Delta\tau_{max}=10^{20}$, and $r_{ser}=2$ are employed for all steady problems. We note that with $r_{ser}=2$ the SER employed in this study is a relatively aggressive approach to update $\Delta\tau$, which could lead to divergence of the pseudo transient continuation. As will be demonstrated in Sections 4.1 and 4.2, two approaches are recommended to facilitate convergence: in the first one, a tight (or small) tolerance of the linear solver, i.e. restarted GMRES, is used to ensure its convergence at each pseudo-time step; in the other, the maximum pseudo-time step size $\Delta\tau_{max}$ is reduced. We note that Ceze and Fidkowski have proposed more sophisticated strategies for updating $\Delta\tau$ in Refs. [44, 45]. They found that a constraint to avoid non-physical state is important to improve the robustness of the pseudo transient continuation. Interested readers are referred to their work; further discussion on this topic is beyond the scope of this paper.

For unsteady flow problems with the dual-time stepping method, define $\boldsymbol{R}=-\left(\frac{\partial\boldsymbol{f}}{\partial x}+\frac{\partial\boldsymbol{g}}{\partial y}-u_{g}\frac{\partial\boldsymbol{q}^{c}}{\partial x}-v_{g}\frac{\partial\boldsymbol{q}^{c}}{\partial y}\right)$. On plugging it into Eq. (22), we obtain

[TABLE]

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

[TABLE]

where $(\boldsymbol{q}^{p})^{m=0}=(\boldsymbol{q}^{p})^{n}$ and when $m\to\infty$, $(\boldsymbol{q}^{p})^{m+1}\to(\boldsymbol{q}^{p})^{n+1}$. To achieve a second-order temporal accuracy, BDF2 is used to discretize the derivative with respect to $t^{*}$. Thus, Eq. (37) is discretized as

[TABLE]

Note that $\Delta t=\Delta t^{*}$. After linearizing $\boldsymbol{R}^{m+1}$ with respect to $\boldsymbol{R}^{m}$, Eq. (38) can be further written as

[TABLE]

where $\Delta(\boldsymbol{q}^{p})^{m}=(\boldsymbol{q}^{p})^{m+1}-(\boldsymbol{q}^{p})^{m}$. For unsteady flow problems tested here, $\Delta\tau$ is always fixed as $\Delta\tau=10^{20}$ for pseudo-time iterations.

In the present study, the full matrix $\frac{\partial\boldsymbol{R}}{\partial\boldsymbol{q}^{p}}$ is computed numerically using the first-order finite difference approach. Note that if only block diagonal terms of the Jacobian matrix were supplied, the convergence would be dramatically deteriorated [21]. The restarted GMRES solver in the PETSc library [46] is used to solve Eq. (34) and Eq. (39).

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. [12, 21]. 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 $\boldsymbol{q}^{p,b}=\boldsymbol{q}^{p,\infty}$ 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, $P^{4}$ 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 $L_{2}$ error is used for the accuracy study. Specifically, for any solution variables $\phi$ on domain $\Omega$, the $L_{2}$ error is defined as

[TABLE]

Exclusively, the entropy error $E_{s}$ is defined as

[TABLE]

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 $Res_{m}/Res_{0}$, denoted as $tol_{pseudo}$; the other is for the restarted GMRES linear solver, denoted as $tol_{res}$. For steady flow simulation using a pseudo transient continuation (see Eq. (33)), the goal is to decrease the steady relative residual $Res_{m}/Res_{0}$ as much as possible with a reasonable criterion for $tol_{res}$; for unsteady flow simulation (see Eq. (36)), the goal is to decrease the unsteady relative residual $Res_{m}/Res_{0}$ to a satisfactory level with a reasonable criterion for $tol_{res}$.

The drag coefficient $C_{d}$, lift coefficient $C_{l}$ and thrust coefficient $C_{t}$ are defined as

[TABLE]

where $A$ is the area that the force acts on, and $U_{\infty}$ is the free stream velocity. The pressure is normalized as

[TABLE]

for better visualization, where $p_{min}$ and $p_{max}$ 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 [47]. The dimension of the Krylov subspace is 150 for steady flow problems and 30 for unsteady ones. The Jacobian matrix $\frac{\partial\boldsymbol{R}}{\partial\boldsymbol{q}^{p}}$ 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 $200\sqrt{2}$. 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 $40\times 48$ and $80\times 96$ meshes, respectively.

Grid refinements for the third- ($P^{2}$) and fourth-order ($P^{3}$) FR/CPR schemes have been conducted. The free parameter $\kappa=1$ is employed. The convergence criterion for the relative residual of the restarted GMRES solver is $tol_{res}={10^{-6}}$ as suggested in Ref. [21]. In Figure 2, the pressure and $Ma$ contours near the cylinder in a free stream with $Ma_{\infty}=10^{-3}$ are presented. No pressure oscillation near the stagnation point is observed. Table 1 documents the drag coefficient $C_{d}$ and entropy error $E_{s}$ from a grid refinement study.

The convergence histories of the relative residual $Res_{m}/Res_{0}$ of pressures are illustrated in Figure 3. It is observed that $Res_{m}/Res_{0}$ can decrease to the level of $10^{-9}$ 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 $Res_{m}/Res_{0}$ will quickly drop as the pseudo-time step size increases to large values. This is due to that the Newton’s method will be recovered as $\Delta\tau\to\infty$. However, as shown in Figure 5, the iteration numbers needed for the restarted GMRES solver to converge (i.e., to meet the convergence criterion $tol_{res}={10^{-6}}$) will significantly increase when $\Delta\tau\to\infty$, especially when fine meshes are used. A remedy could be decreasing the maximum value of $\Delta\tau$. As presented in Figure 5(b), when solving the problem with the fourth-order FR method on the $80\times 96$ mesh, $\Delta\tau_{max}=50$ is employed to ensure that the iteration number needed by the GMRES solver to converge is around 400 when $\Delta\tau=\Delta\tau_{max}$. Consequently, more pseudo-time iterations are needed for the pseudo transient procedure to converge as shown in Figure 3(b).

Effects of the Global Cut-off Parameter $\kappa$

We have tested three different values of $\kappa$, namely, $\kappa=\sqrt{0.2}$, $\kappa=1$ and $\kappa=\sqrt{5}$ to study the effects of $\kappa$ on convergence and accuracy. The convergence criterion for the relative residual of the restarted GMRES solver is set as $tol_{res}={10^{-9}}$ in these tests. For brevity, only the third-order FR/CPR scheme is considered on the $80\times 96$ mesh. The results of $C_{d}$, $C_{l}$ and $E_{s}$ 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 $\kappa=\sqrt{0.2}$, 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 $\kappa=\sqrt{5}$, $Res_{m}/Res_{0}$ can drop to a smaller value using less pseudo-time iterations compared to that with $\kappa=1$, the accuracy does not get better, and more iterations are needed for the linear solver to converge in one pseudo-time step. The predicted $C_{d}$ and $E_{s}$ both increase due to the fact that when a larger $\kappa$ is used, the maximum absolute eigenvalue $|\lambda|_{max}$ in Eq. (12) is increased. Therefore, more numerical dissipation is added to the approximate Riemann solver.

4.2 Inviscid Flow over a NACA0012 Airfoil

In this section, the inviscid flow of $Ma=10^{-3}$ and $0^{\circ}$ 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

[TABLE]

where $x\in[0,1]$. The mesh with 5,168 quadrilateral elements is illustrated in Figure 7.

The polynomial degree is refined from $2$ to $4$ (i.e., from the $3^{rd}$ to $5^{th}$ order of accuracy). As presented in Table 3, the prediction of $C_{d}$ 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 $tol_{res}$, i.e., $10^{-3}$, $10^{-6}$ and $10^{-9}$, on the convergence.

In Ref. [21], a 1,792-element mesh is employed and the maximum iteration number of the restarted GMRES solver when $\Delta\tau\to\infty$ 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 $140$ with the same $tol_{res}=10^{-6}$ as that used in Ref. [21]. According to the convergence comparison of different criteria for the restarted GMRES solver (see Figure 9(a)), when $tol_{res}=10^{-3}$, 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. With the SER approach employed in this study, a smaller $tol_{res}$ for the linear solver is suggested for the sake of convergence, but more iterations of the restarted GMRES solver are needed (see Figure 9(b)). Only negligible differences are observed on $C_{d}$ and $E_{s}$ among the converged results of different $tol^{\prime}_{res}$s, which are not presented for brevity. Therefore, a proper choice of $tol_{res}$ is a trade-off between robustness and efficiency.

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 $(\rho,u,v,Ma)=(1.0,\sqrt{2}/2,\sqrt{2}/2,0.05)$ and the gas constant $R=1.0$ for this case. The perturbation is defined as [47]

[TABLE]

where $\phi=\frac{1}{2}$ and $\alpha=5$ are parameters that define the vortex strength. $r=(x-x_{0})^{2}+(y-y_{0})^{2}$ is the distance to the center of the vortex $(x_{0},y_{0})$. The vortex is propagated in a periodic domain $[-10,10]\times[-10,10]$. A prescribed grid motion is employed to validate the moving grid algorithm. The prescribed grid motion is defined as

[TABLE]

where $x_{r},y_{r}$ are the coordinates at $t=0$, i.e., coordinates of the uniformly divided grid, and $a_{x}=1$, $a_{y}=1$, $f_{n}=1.0$, $f_{x}=0.1$, $f_{y}=0.1$. Periodic boundary conditions are applied to all boundaries. Since there is no presence of wall boundaries, the global cut-off parameter $\kappa$ is set as zero for this case.

The time step size $\Delta t$ is refined from $T_{n}/100$ to $T_{n}/800$, where $T_{n}=1/f_{n}$ to validate the accuracy of BDF2. A fifth-order FR/CPR scheme and a $96\times 96$ uniformly divided mesh are employed to ensure that the error is dominated by the temporal discretization. We require the relative residual $Res_{m}/Res_{0}$ of the pseudo-time iterations to drop as low as possible within 100 pseudo-time iterations. The $tol_{res}$ of the restarted GMRES solver is $10^{-6}$ with the maximum iteration number 100. We only simulate this case until $\frac{1}{4}T_{n}$ 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 $p$ and $u$.

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

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 $Re=500$, $Ma=10^{-3}$ and $0^{\circ}$ 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 $0.1723$, $0.1722$ and $0.1722$, 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 $y=Y\sin(2\pi ft)$, where $Y$ is the plunge amplitude, and $f$ is the plunge frequency. The reduced frequency $k$ is defined as $k={2\pi fC}/{U_{\infty}}$, and dimensionless height $h=Y/C$, where $C$ is the chord length of the airfoil. In this study, $k=2$, $h=0.4$ [8]. The whole mesh is considered as an oscillating rigid body. At the far field, $\boldsymbol{q}^{p,b}=\boldsymbol{q}^{p,\infty}$, and the approximate Riemann solver is used to enforce boundary conditions.

The $tol_{res}$ of the restarted GMRES solver is $10^{-6}$ with the maximum iteration number 100. For the sake of efficiency, we would require the $Res_{m}/Res_{0}$ drop by certain orders instead of as low as possible. In this case, three criteria for $Res_{m}/Res_{0}$ of the pseudo-time marching have been tested, i.e., $tol_{pseudo}=10^{-2}$, $10^{-4}$, and $10^{-6}$, to examine if the force predictions are sensitive to the convergence criterion of $Res_{m}/Res_{0}$. The time step size $\Delta t$ is set to $T/100$, and $\kappa$ is set to $\sqrt{5}$. All simulations end at $t_{end}=10T$. In Figure 12, the vortex shedding process during one typical plunging period is presented. The thrust coefficient $C_{t}$ and lift coefficient $C_{l}$ of the tenth period are displayed in Figure 13, and the time-averaged thrust coefficient $\bar{C_{t}}$, root mean square of the lift coefficient $C_{l,rms}$, and maximum lift coefficient $C_{l,max}$ are documented in Table 4. From Figure 13 and Table 4, we observe that a moderate tolerance ($10^{-4}$) for $Res_{m}/Res_{0}$ is sufficient to maintain the accuracy of force prediction. However, a smaller $Res_{m}/Res_{0}$ criterion means more iterations.

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., $10^{-3}$, $10^{-2}$, $10^{-1}$ and $0.2$, are studied here. A smaller time step size is employed, i.e., $\Delta t=T/1000$, to ensure that the time integration errors are negligible. $tol_{pseudo}=10^{-4}$ is adopted as the convergence criterion for pseudo-time iterations. The $C_{l}$ and $C_{t}$ of the tenth period are illustrated in Figure 14. $\bar{C_{t}}$, $C_{l,rms}$ and $C_{l,max}$ are documented in Table 5.

It is observed that when $Ma$ is less than $10^{-2}$, the compressibility is almost negligible. The force histories for $Ma=10^{-3}$ and $Ma=10^{-2}$ coincide with each other. When $Ma$ is larger than $10^{-1}$, phase shift is observed in the force histories. In general, when $Ma$ is increased to a certain level (e.g., $10^{-1}$ in this case), the time-averaged thrust coefficient $\bar{C_{t}}$ will decrease when $Ma$ increases; at the same time, the root mean square of the lift coefficient $C_{l,rms}$ and the maximum lift coefficient $C_{l,max}$ will vary slightly. For example, when $Ma=0.2$, the reduction of $\bar{C_{t}}$ is over $36\%$ compared to that at $Ma=10^{-3}$; $C_{l,rms}$ will increase about $1\%$; $C_{l,max}$ will first increase and then decrease within a small range.

5 Conclusions

In this study, a high-order preconditioned FR/CPR method is developed to solve compressible Navier–Stokes equations at low Mach numbers on moving grids. The dual-time stepping method is used to handle unsteady flows. Specifically, BDF2 is adopted to discretize the physical time derivative term, and the restarted GMRES linear solver with the modified SER pseudo-time step size adaptation strategy is used to converge unsteady residuals at each physical time. For simulations with moving/deforming grids, the GCL is enforced implicitly to eliminate the discrepancy between the temporal and spatial discretization. The preconditioned FR/CPR method has demonstrated good accuracy and convergence through various benchmark tests.

We first study the convergence properties of the preconditioned FR/CPR method for steady inviscid flow simulations. Two cases, namely, low-Mach-number (i.e., $Ma=10^{-3}$ in this study) flows over a cylinder and a NACA0012 airfoil are simulated. We observe that when the SER algorithm is used, the convergence criterion $tol_{res}$ of the restarted GMRES solver should be sufficiently tight (or small) for the sake of convergence. When $\Delta\tau\to\infty$ and SER is aggressive, the number of iterations for the GMRES solver to converge can be significantly increased, especially for higher-order spatial schemes on fine meshes. A remedy could be decreasing the maximum pseudo-time step $\Delta\tau_{max}$. As a trade-off, the convergence speed of pseudo transient continuation will reduce correspondingly. A robust and efficient strategy for pseudo transient continuation is worthy of more research efforts in the future.

The preconditioned FR/CPR method is then used to simulate low-Mach-number unsteady flows. The temporal and spatial orders of accuracy are first verified using isentropic vortex propagation on both static and dynamic meshes. After that, we study the numerical properties of the preconditioned FR/CPR method with low-Mach-number flows over a plunging NACA0012 airfoil. We find that the convergence tolerance $tol_{pseudo}$ for $Res_{m}/Res_{0}$ of the pseudo-time marching does not need to be very small to preserve simulation accuracy. Since less pseudo-time iterations are needed when a larger criterion for $Res_{m}/Res_{0}$ is used, computational cost of unsteady flows can be decreased. We also conduct numerical experiments to test the effects of the Mach number on unsteady force generation. We find that when the Mach number is small enough (e.g., $Ma=10^{-2}$), flow compressibility has negligible effect on force prediction of a plunging airfoil. However, when the Mach number becomes larger (e.g., $Ma=0.2$), thrust generation of a plunging airfoil will change significantly while the change of lift and its root mean square are not apparent.

Acknowledgment

The authors gratefully acknowledge the support of the Office of Naval Research through the award N00014-16-1-2735, and the faculty startup support from the department of mechanical engineering at University of Maryland, Baltimore County (UMBC).