A High-Order Discontinuous Galerkin Method for Solving Preconditioned Euler Equations

: A high-order discontinuous Galerkin (DG) method is presented for solving the preconditioned Euler equations with an explicit or implicit time marching scheme. A detailed description is given of a practical implementation of a precondition matrix of the type of Weiss and Smith and of the DG spatial discretization scheme employed, with particular emphasis on the artiﬁcial viscosity-based shock capturing techniques. The curved boundary treatment is proposed through adopting a NURBS surface equipped with a radial basis function interpolation to propagate the boundary displacement to the interior of the mesh. The resulting methods are veriﬁed by simulating ﬂows over two-dimensional airfoils, such as symmetric NACA0012 or asymmetric RAE2822, and over three-dimensional bodies, such as an academic hemispherical headform or aerodynamic ONERA M6 wing. Numerical results show that the present method functions for both transonic and nearly incompressible ﬂow simulations, and the proposed treatment of curved boundaries, play an important role in improving the accuracy of the obtained solutions, which are in good agreement with available experimental data or other numerical solutions reported in literature.


Introduction
Under the dominance of second order accurate methods in computational fluid dynamics (CFD) for real applications, the discontinuous Galerkin (DG) method started becoming increasingly popular in recent years due to its ability to achieve higher order spatial accuracy rigorously, even on an unstructured grid [1][2][3][4]. The compressible DG method is used in a wide spectrum of flow regimes [5][6][7][8][9][10]. In order to be able to capture shocks existing in the flow fields, several key ingredients, such as slope limiter [5], weighted essentially non-oscillatory (WENO) based limiter [6,7], and artificial viscosity [8,9] are usually imposed. The application of the slope limiter requires the identification of the elements lying in the shock region, and also requires for the reduction of the interpolating polynomial to be first order in those identified elements, which may destroy the high-order accuracy in the shocked region. In order to maintain the required accuracy, a WENO-based limiter can be selected at the cost of high computational time consuming. In contrast, the artificial viscosity-based approach is reported to be efficient in view of a relatively low cost [8] and to be flexible in terms of independence on the approximation order with different types of mesh elements [9]. Thus, in the present work, artificial viscosity is intended to be imposed on the preconditioned DG method, capable of compressible shocked flow simulations with explicit or implicit time integration schemes.
In general, the schemes with implicit time integration have the advantage of the time steps not being restricted by the severely strict stability limits in comparison to the explicit counterpart [2,10], and hence so-called implicit acceleration can often be achieved in terms of the convergence rate. In fact, successful efforts and applications were reported using

Governing Equations
The three-dimensional compressible Euler equations governing inviscid flows can be expressed in a dimensionless conservative form [20] as follows: where the unknown vector Q of conserved variables and the Cartesian components of the flux function F = (f 1 , f 2 , f 3 ) are given by where ρ, p, and E are the density, pressure, and total energy per unit mass, respectively, and u, v, and w are the Cartesian components of the velocity vector V. This system can be closed up by adding the following equation of state: which is valid for perfect gas, Γ = C p C v is the ratio of specific heats and taken as Γ = 1.4 for air.
The Equation (1) in conservative variables are reported to be ill-conditioned for the flows at low Mach numbers [21]. In order to single out the propagation of acoustic waves [22], the governing equations in primitive variables q = (p, u, v, w, T) T are considered here. By transforming the dependent variables in Equation (1) from conserved quantities to such primitive variables, Equation (1) can be written as where P is the transformation matrix here, H = E + p/ρ is total enthalpy per unit mass of the fluid. For inviscid flows of interest, the derivatives that appear in Equation (5) can be directly determined by the equation of state, p = 1 Γ ρT, which reads: where c is the acoustic velocity. The derivative ρ p that multiplies the pressure time derivative in the continuity equation controls the speed of propagation of acoustic waves [23] in the system (4). This is consistent with the notion of infinite pressure wave speeds in an incompressible flow due to ρ p = 0 for constant density flows. The disparity of wave speeds of the system (4), mainly between acoustic waves and entropy waves, becomes very large for low Mach number flows, and can be ameliorated through adopting techniques such as preconditioning [10,23]. Following the well-known preconditioning technique of Weiss and Smith [23], the matrix P in Equation (4) is modified through replacing the derivative ρ p , with one proportional to the inverse of the local velocity squared to control the disparity of the wave Appl. Sci. 2022, 12, 7040 4 of 22 speeds of the system (4). The resultant matrix, termed as "preconditioning matrix Γ", can be expressed as follows: where Θ is computed by Here U r is a reference velocity defined as follows: where ε is a small number defined by the user for preventing singularities at stagnation points. In the present paper, ε = O(M) is adopted as chosen in Reference [10]. In so doing, the governing equations of interest become the preconditioned Euler equations, which can be expressed as follows: In a strict sense, Equation (9) is not conservative for the time-dependent flows due to the fact that the time accuracy of the equations is destroyed by preconditioning. However, solutions in the steady state remain unchanged [10,23]. Thus, it is not a problem to employ Equation (9) for present steady calculations.
The eigenvalues of the preconditioned system (9) in the direction of the unit vector n = n x , n y , n z , which can be given by where all the eigenvalues that indicate wave speeds of the resulting system are modified to be the same magnitude order (see Reference [23] for detailed expressions of eigenvalues).

Spatial Discretization
In this work, the discontinuous Galerkin finite element method is considered to be implemented on each element with polynomial functions of a maximum degree equal to p, which results in the method being (p + 1)th order accurate. Let d be the number of space dimensions, and the number of degrees of freedom per equation per element is Assume that the computational domain Ω ⊂ R d is subdivided into a collection of non-overlapping finite elements {Ω i } and define the following space for both solution and test function: where V m p is the space of vector-valued polynomials of a degree, at most p, and m is the dimension of the conservative state vector (note: m = 4 for present 2D case and m = 5 for 3D).
Multiplying Equation (9) by a test function W h ∈ V p h and integrating by parts on each element Ω i , we obtain the weak formulation of the DG method: where q h is the approximate solution at element Ω i and can be represented by basis functions φ = φ j ∈ V p h as follows: where the unknowns q i,j (t) are the solution coefficients of the element to be solved. According to Equation (12), we have the following system of N p equations for each Ω i Thus, the unknowns q i,j (t) can be determined through solving system (14). Keep in mind that the flux terms are not uniquely defined at element interface due to discontinuous function approximation. In order to treat this issue properly, the flux term F(q h ) · n, appearing in the second term of Equation (14), is usually required to be replaced by a numerical flux H q + h , q − h , n , where q + h is the neighboring element interface state, q − h the internal interface state, and n the direction normal to the interface. It is reported that a Roe-type flux function [24] satisfies the flux consistency conditions [25] associated with ensuring the formal accuracy of the scheme, and hence is adopted in the present paper. As mentioned above, the wave speeds of the governing equations were changed by preconditioning, and the related dissipative part of the numerical fluxes H q + h , q − h , n are supposed to be modified accordingly.
Following the work of Reference [21], the resultant numerical flux can be written as It can be noted that the effect of preconditioning is considered in the dissipative part of the last term of the Equation (15). σ ∈ (0, 1] is a user-tuned factor that controls the amount of the dissipation [26].

Shock Capturing
In order to simulate high-speed shocked flows, a kind of shock capturing technique [8] is adopted through adding a small amount of artificial viscosity to the preconditioned governing Equation (9) in order to deal with possible strong shocks, which reads: where ∇q is the gradient of the primitive variable vector and is computed by the BR2 scheme [27] in the present work, and µ is the artificial viscosity coefficient. Keep in mind that the artificial viscosity introduced is intended to suppress non-physical oscillations in the vicinities of the strong shocks without affecting the solution accuracy in smooth regions.
The smoothness of solution is measured by variable, which could be pressure or density. In the present paper, pressure is used. Following the work of Reference [8], the resultant µ i associated with Ω i is computed by where µ 0 , s 0 , and κ are the control parameters that are usually chosen empirically by the users. It can be noted that the coefficients for each component in expansion (13) decay very quickly in smooth regions, while the strength of the discontinuity in shock regions will slow down the rate of such decay [8]. In other words, the elemental artificial viscosity coefficient µ i is close to zero in smooth regions and increases in shock regions. By adopting such an artificial viscosity coefficient, the shock capturing technique can restrain the spurious numerical oscillation in shock regions, causing almost no degrading of solution accuracy in the smooth regions. However, the above piecewise constant artificial viscosities are discontinuous at the element interfaces, which may still cause small numerical oscillations in the vicinities of the strong shocks [28]. The artificial viscosities, made to be continuous, are usually preferred. In the present work, such continuity is achieved through adopting a post-smoothing process, realized in two steps: Step 1. Transferring the elemental viscosity to the vertex of the element. The viscosity v k at vertex k is determined by where n is the number of the adjacent elements with weights Step 2. Computing the viscosity at each integral point of the element through a linear interpolation.
It can be noted that a kind of smoothing is realized through a weighted average process in Step 1.
In so doing, the semi-discrete weak form of the Equation (16) can be obtained from Equation (14) through adding the corresponding artificial viscosity terms, which reads: Following the treatment of BR2 scheme [27], the additional flux term associated with artificial viscosity in Equation (18) is determined by The quantities at the boundary are determined through enforcing the related boundary conditions.

Boundary Conditions
It can be learned from Equation (18) that the enforcement of high-order accurate boundary conditions is straightforward and is mainly involved with the numerical fluxes H q b , q − h , n , which are weakly prescribed at the integral points of the boundary. The state of the primitive variables at each integral point of the boundary, T , is determined by boundary conditions, such as far field and solid wall boundary conditions.

Subsonic Far Field Boundary Condition
The far field boundary conditions are usually implemented to be nonreflecting. It is well known that the boundary state q b determined by Riemann invariants is nonreflecting for traditional Euler equations. However, relatively complicated characteristic boundary conditions [29] should be used for the present preconditioned Euler equations due to the fact that the wave speeds of the governing equations were changed by preconditioning. In the present paper, an alternative simplified treatment is adopted following the work of Turkel [29].
In order to make it clear, take the far field subsonic flow as an example; the boundary state at each integral point can be assigned simply by where superscript '∞' denotes the freestream variables, and superscript '−' the internal variables. Similar treatment can be taken for the subsonic outflow situation of interest.

Solid Wall Boundary Condition
In the case of inviscid flows, the fluid slips over the wall surface. In other words, the normal component of the velocity vanishes at the solid wall boundary. Such a slip wall condition employed is based on following boundary state: With the boundary variables determined by Equation (20), the numerical flux of Equation (15) at the solid wall becomes F q b · n. This means that the computation of fluxes on the wall boundary are in the same manner as the one in non-preconditioned DG schemes.

Treatment of Curved Boundaries
As pointed out in the Reference [10], an appropriate representation of the possibly curved boundary geometry is essential for preserving numerical accuracy. Here, some details concerning the present treatment of the curved boundary are described.

Geometry Representation
The integrals that appear in Equation (18) are intended to be calculated in reference elements [30] using the Gaussian quadrature rules. The physical coordinate x = (x, y, z) associated with the reference coordinate ξ = (ξ, η, ζ) can be expressed as where n d is the total number of the control points x i , and ψ i (ξ) are the corresponding shape functions that are fixed in the reference element κ. As illustrated in Figure 1, with the mapping relationship (21), the linear reference element can be used for representing simple straight-edged elements exactly, while quadratic reference elements can be applied for approximating curve-edged elements only. For enhancing the accuracy of high-order approximation, special treatment, so-called high-order re-mesh, is further required [9] and addressed in the following.
shape functions that are fixed in the reference element .
As illustrated in Figure 1, with the mapping relationship (21), the linear element can be used for representing simple straight-edged elements exact quadratic reference elements can be applied for approximating curve-edged only. For enhancing the accuracy of high-order approximation, special tr so-called high-order re-mesh, is further required [9] and addressed in the follow

High-Order Re-Mesh
It is recognized that the order of the DG method is limited by the ord boundary representation [25]. Thus, it is necessary to have a proper treatme curved boundary. Here, a practical approach is proposed through high-order rein view of easy implementation and curved-mesh quality.
The approach is NURBS surface-based, and a NURBS surface is enhance RBF interpolation in the high-order re-mesh procedure to have an appropriate t of curved boundaries. Specifically, NURBS-based surface representations a generated based on the initial surface elements of the curved boundaries (often

High-Order Re-Mesh
It is recognized that the order of the DG method is limited by the order of the boundary representation [25]. Thus, it is necessary to have a proper treatment of the curved boundary. Here, a practical approach is proposed through high-order re-meshing in view of easy implementation and curved-mesh quality.
The approach is NURBS surface-based, and a NURBS surface is enhanced with a RBF interpolation in the high-order re-mesh procedure to have an appropriate treatment of curved boundaries. Specifically, NURBS-based surface representations are firstly generated based on the initial surface elements of the curved boundaries (often available in AutoCAD shapes) (see Figure 2a), and then used to create high-order-required control points on the curved surfaces or edges (see Figure 2b-d).
A key to form an appropriate set of high-order-required control points lies in adding new control points, which can be easily operated from the initial surface elements. Take the second order approximation as an example; the midpoints of the vertexes of the initial surface elements are usually selected to be control points, as illustrated by the red points in Figure 2c. Such newly added control points are certainly required to be mapped to the curved NURBS surface generated, which can be done by solving the problem associated with the rule of shortest displacement [11], which reads: where S(a, b) denotes a point on NURBS surface and x new is a newly added control point. Keep in mind that with the NURBS surface generated, the gradient at the point S(a, b) can be directly computed. Thus, the above problem can be efficiently solved by gradient-based optimization, such as the Gauss-Newton algorithm, to obtain the required new control points (see red points in Figure 2d). in AutoCAD shapes) (see Figure 2a), and then used to create high-order-required contro points on the curved surfaces or edges (see Figure 2b-d). A key to form an appropriate set of high-order-required control points lies in adding new control points, which can be easily operated from the initial surface elements. Take the second order approximation as an example; the midpoints of the vertexes of the ini tial surface elements are usually selected to be control points, as illustrated by the red points in Figure 2c. Such newly added control points are certainly required to be mapped to the curved NURBS surface generated, which can be done by solving the problem as sociated with the rule of shortest displacement [11], which reads: where ( , ) denotes a point on NURBS surface and is a newly added contro point. Keep in mind that with the NURBS surface generated, the gradient at the poin ( , ) can be directly computed. Thus, the above problem can be efficiently solved by gradient-based optimization, such as the Gauss-Newton algorithm, to obtain the re quired new control points (see red points in Figure 2d). However, as pointed out in References [9,31], such kind mapping may cause low quality or even invalid elements in the vicinity of the curved boundary with high cur vature (see intersected surfaces existed at the vertex B of the element in Figure 3b). Such a problem is reported to be ameliorated by improving the quality of the mesh near the boundary through adopting optimization processes [32]. An alternative approach used widely is to propagate the boundary movement to the interior through mesh-deforming techniques, including the RBF interpolation method [31], the linear elasticity method However, as pointed out in References [9,31], such kind mapping may cause low quality or even invalid elements in the vicinity of the curved boundary with high curvature (see intersected surfaces existed at the vertex B of the element in Figure 3b). Such a problem is reported to be ameliorated by improving the quality of the mesh near the boundary through adopting optimization processes [32]. An alternative approach used widely is to propagate the boundary movement to the interior through mesh-deforming techniques, including the RBF interpolation method [31], the linear elasticity method [9,33], etc. Among them, the RBF interpolation method has the ability to achieve high-quality meshes at a relatively low computational costs [31], and hence, whether the invalid elements exist or not, the RBF interpolation is applied as the last step of the present approach to retain the element quality (see Figure 3c). In order to assure the quality of the resulting mesh, an additional step of quality check is usually required. But it is learned from the test cases presented (see Section 5) that it is not necessary after having the RBF interpolation. sent approach to retain the element quality (see Figure 3c). In order to assure the quality of the resulting mesh, an additional step of quality check is usually required. But it is learned from the test cases presented (see Section 5) that it is not necessary after having the RBF interpolation. Applying the RBF interpolation method, the displacement ( ) of the point in the computational domain can be expressed as where , = 1,2, . . . , are the positions of the boundary control points, ( − ) is the radius basis function for the displacement field , and is the corresponding unknown coefficient in that basis. With the given boundary displacements, ( ), = 1,2, . . . , , the unknown coefficients can be determined by solving the following system of equations: Once the coefficients are determined, the interior displacement of any point in the domain can be calculated by Equation (22). Concerning further details of the RBF interpolation, we refer to the work of Reference [31].

Time Stepping Schemes
By assembling all the elemental contributions together, the semi-discrete form of Equation (18) governing the evolution in time of the discrete solution can be written as where is the block diagonal matrix in which each block corresponds to a specific element, is the vector of all the solution coefficients to be solved, and is the vector of the residual. Both explicit and implicit time stepping schemes are used in the present paper for discretizing the semi-discrete Equation (24) in time. Specifically, an implicit LU-SGS scheme will be mainly addressed after having a brief description of the explicit Runge-Kutta scheme.

Explicit Runge-Kutta Scheme
The explicit three-stage third-order TVD Runge-Kutta scheme [5] is used to advance the solution of Equation (24) from time to time + ∆ , which reads: Applying the RBF interpolation method, the displacement s(x) of the point x in the computational domain can be expressed as where x bj , j = 1, 2, . . . , n b are the positions of the boundary control points, ψ x − x bj is the j th radius basis function for the displacement field s, and a j is the corresponding unknown coefficient in that basis. With the given boundary displacements, s(x bi ), i = 1, 2, . . . , n b , the unknown coefficients a j can be determined by solving the following system of n b equations: Once the coefficients a j are determined, the interior displacement of any point in the domain can be calculated by Equation (22). Concerning further details of the RBF interpolation, we refer to the work of Reference [31].

Time Stepping Schemes
By assembling all the elemental contributions together, the semi-discrete form of Equation (18) governing the evolution in time of the discrete solution can be written as where M is the block diagonal matrix in which each block corresponds to a specific element, q is the vector of all the solution coefficients to be solved, and R is the vector of the residual. Both explicit and implicit time stepping schemes are used in the present paper for discretizing the semi-discrete Equation (24) in time. Specifically, an implicit LU-SGS scheme will be mainly addressed after having a brief description of the explicit Runge-Kutta scheme.
In addition to the preconditioning mentioned, the local time stepping is also used to improve the convergence speed to steady state solutions. In specific, the local time step ∆t i is computed for each element Ω i by where p is the approximate order as mentioned, and CFL is the Courant number. In general, CFL could be 1 for 2D and less than 1 for 3D applications of the DG method.

Implicit LU-SGS Scheme
The implicit backward Euler time discretization of Equation (24) can be written as where ∆ q n = q n+1 − q n , in which the superscripts n and n + 1 denote the current and the next time steps, respectively. ∂R n ∂ q is the Jacobian matrix of the DG spatial discretization, which is computed analytically (except for the dissipative part of the numerical flux using divided differencing [34]) in this paper.
Recall that M = diag(M 1 , M 2 , . . . , M N c ), and let C i be a set of the indexes of the neighboring connected elements of Ω i , define A = M ∆t − ∂R n ∂ q as the global system matrix, and the Equation (27) can then be written as in which The resultant system of Equation (28) encapsulates the implicit iteration schemes, and it can be solved iteratively to converge to a steady state. The traditional LU-SGS scheme consists of a forward and a backward sweeping through all the computational elements in a sequential order [16], which can be written as Backward sweep :

Results
In this section, the performance of the present DG method is investigated by simulating a set of typical two-and three-dimensional flow problems. The test cases include flows over a symmetric NACA0012 airfoil, an asymmetric RAE2822 airfoil, a hemispherical headform, as well as more practical flows over the ONERA M6 wing. All the simulations are carried out with the present implicit LU-SGS algorithm in view of having implicit accelerations. The explicit scheme is used only for the cases with comparisons between implicit and explicit schemes.

Low-Speed Flows over a NACA0012 Airfoil
In order to verify whether our preconditioned high-order DG method can be used for nearly incompressible flow simulations, flows over a symmetric NACA0012 airfoil with zero angle of attack (α = 0 • ) at four selected Mach numbers, 0.1, 0.01, 0.001, and 0.0001, are first simulated. Without loss of generality, the order of DG method is fixed to be fourth (p = 3) during the simulations, and a fully unstructured mesh having only 3527 triangle elements (see Figure 4) is used in view of adopting a fourth-order accurate DG method.
In order to verify whether our preconditioned high-order DG metho for nearly incompressible flow simulations, flows over a symmetric NA with zero angle of attack ( = 0°) at four selected Mach numbers, 0.1, 0 0.0001, are first simulated. Without loss of generality, the order of DG met be fourth ( = 3) during the simulations, and a fully unstructured mesh ha triangle elements (see Figure 4) is used in view of adopting a fourth-orde method. The behaviors of low speed in convergence (stiffness problem [18,19 numbers can be learned from the convergence histories of three repres numbers, as shown in Figure 5. It seems that the high-order DG method preconditioning) has the ability to cope with the flows at relatively low M (see = 0.1 in Figure 5a), which is different from low order methods difficulties in convergence appear as the Mach number further reduces. S is greatly ameliorated, as illustrated in Figure 5b Figure 6. It can be noted that the developed code, with tioning, fails to obtain the correct physical solution (reaching unphysical so the left-hand side) in view of the chaotic, asymmetric, and poor smoothn tours, while the code with preconditioning can capture smooth contour number (see the right-hand side of Figure 6). The behaviors of low speed in convergence (stiffness problem [18,19]) at low Mach numbers can be learned from the convergence histories of three representative Mach numbers, as shown in Figure 5. It seems that the high-order DG method itself (without preconditioning) has the ability to cope with the flows at relatively low Mach numbers (see Ma = 0.1 in Figure 5a), which is different from low order methods [35]. However, difficulties in convergence appear as the Mach number further reduces. Such a situation is greatly ameliorated, as illustrated in Figure 5b,c, through the help of preconditioning. 0.0001, are first simulated. Without loss of generality, the order of DG method is fixed to be fourth ( = 3) during the simulations, and a fully unstructured mesh having only 3527 triangle elements (see Figure 4) is used in view of adopting a fourth-order accurate DG method. The behaviors of low speed in convergence (stiffness problem [18,19]) at low Mach numbers can be learned from the convergence histories of three representative Mach numbers, as shown in Figure 5. It seems that the high-order DG method itself (without preconditioning) has the ability to cope with the flows at relatively low Mach numbers (see = 0.1 in Figure 5a), which is different from low order methods [35]. However, difficulties in convergence appear as the Mach number further reduces. Such a situation is greatly ameliorated, as illustrated in Figure 5b,c, through the help of preconditioning. The effects of very low speed on the accuracy of the solutions can be learned from the corresponding contours of the Mach number. The typical contours of = 0.0001 are presented in Figure 6. It can be noted that the developed code, without preconditioning, fails to obtain the correct physical solution (reaching unphysical solution [18], on the left-hand side) in view of the chaotic, asymmetric, and poor smoothness of the contours, while the code with preconditioning can capture smooth contours of the Mach number (see the right-hand side of Figure 6).  Figure 6. It can be noted that the developed code, without preconditioning, fails to obtain the correct physical solution (reaching unphysical solution [18], on the lefthand side) in view of the chaotic, asymmetric, and poor smoothness of the contours, while the code with preconditioning can capture smooth contours of the Mach number (see the right-hand side of Figure 6).

Flows over an Asymmetric RAE2822 Airfoil
The capabilities of the developed code equipped with shock capturing, d Section 2.3, are further verified to deal with both transonic compressible an compressible flows over an asymmetric RAE2822 airfoil. The computationa discretized by a fully unstructured mesh having 4386 triangle elements, a Figure 7. A set of simulations with different flow conditions was carried o ferent accurate orders. Here, numerical results of third-order ( = 2) accurate are presented in Figure 8 for the nearly incompressible flow ( = 0.01 and and in Figure 9 for the transonic shocked flow ( = 0.729 and = 2.31 noted from Figure 8 that the obtained surface pressure coefficients (see Figu good agreement with other researchers' results [35,36] and the Mach conto tributed smoothly (see Figure 8b). Additionally included in Figure 9a are mental data [37] and other researchers' results [38]. Quite a good agreement noted in view of the strength or location of the captured shock. Frankly, su solution is based on the strength and the affected region of the artificial v mentioned in the References [8,28], the region functioned with artificial vi narrow as possible for the given mesh. The results can be further improved local refinement of the mesh near the shock, but this is not considered here.

Flows over an Asymmetric RAE2822 Airfoil
The capabilities of the developed code equipped with shock capturing, described in Section 2.3, are further verified to deal with both transonic compressible and nearly incompressible flows over an asymmetric RAE2822 airfoil. The computational domain is discretized by a fully unstructured mesh having 4386 triangle elements, as shown in Figure 7. A set of simulations with different flow conditions was carried out with different accurate orders. Here, numerical results of third-order (p = 2) accurate simulations are presented in Figure 8 for the nearly incompressible flow ( Ma = 0.01 and α = 1.89 • ), and in Figure 9 for the transonic shocked flow (Ma = 0.729 and α = 2.31 • ). It can be noted from Figure 8 that the obtained surface pressure coefficients (see Figure 8a) are in good agreement with other researchers' results [35,36] and the Mach contours are distributed smoothly (see Figure 8b). Additionally included in Figure 9a are the experimental data [37] and other researchers' results [38]. Quite a good agreement can also be noted in view of the strength or location of the captured shock. Frankly, such shocked solution is based on the strength and the affected region of the artificial viscosity. As mentioned in the References [8,28], the region functioned with artificial viscosity is as narrow as possible for the given mesh. The results can be further improved when using local refinement of the mesh near the shock, but this is not considered here.
(a) (b) Figure 6. Contours of the Mach number for flows over NACA0012 airfoil. (a) W tioning; (b) with preconditioning.

Flows over an Asymmetric RAE2822 Airfoil
The capabilities of the developed code equipped with shock capturin Section 2.3, are further verified to deal with both transonic compressible compressible flows over an asymmetric RAE2822 airfoil. The computatio discretized by a fully unstructured mesh having 4386 triangle elements Figure 7. A set of simulations with different flow conditions was carried ferent accurate orders. Here, numerical results of third-order ( = 2) accura are presented in Figure 8 for the nearly incompressible flow ( = 0.01 a and in Figure 9 for the transonic shocked flow ( = 0.729 and = 2. noted from Figure 8 that the obtained surface pressure coefficients (see Fi good agreement with other researchers' results [35,36] and the Mach co tributed smoothly (see Figure 8b). Additionally included in Figure 9a a mental data [37] and other researchers' results [38]. Quite a good agreem noted in view of the strength or location of the captured shock. Frankly solution is based on the strength and the affected region of the artificia mentioned in the References [8,28], the region functioned with artificial narrow as possible for the given mesh. The results can be further improv local refinement of the mesh near the shock, but this is not considered here  It should be mentioned that based on the preconditioned Euler equation, the shock capturing is implemented and only functions when shock appears, and hence the per formance of the developed code in convergence is still kept for the nearly incompressible flow. Actually, both explicit and implicit time marching schemes described in Section 4 are used for the code developed. The implicit speed up of convergence histories is re markable for the above nearly incompressible flow, as illustrated in Figure 10, in view o being consuming in both iteration and time.  It should be mentioned that based on the preconditioned Euler equation, the shock capturing is implemented and only functions when shock appears, and hence the per formance of the developed code in convergence is still kept for the nearly incompressible flow. Actually, both explicit and implicit time marching schemes described in Section 4 are used for the code developed. The implicit speed up of convergence histories is re markable for the above nearly incompressible flow, as illustrated in Figure 10  It should be mentioned that based on the preconditioned Euler equation, the shock capturing is implemented and only functions when shock appears, and hence the performance of the developed code in convergence is still kept for the nearly incompressible flow. Actually, both explicit and implicit time marching schemes described in Section 4 are used for the code developed. The implicit speed up of convergence histories is remarkable for the above nearly incompressible flow, as illustrated in Figure 10, in view of being consuming in both iteration and time.
formance of the developed code in convergence is still kept for the nearly inco flow. Actually, both explicit and implicit time marching schemes described i are used for the code developed. The implicit speed up of convergence his markable for the above nearly incompressible flow, as illustrated in Figure 10 being consuming in both iteration and time.

Flows over a Hemispherical Headform
It is pointed out that the above relatively accurate results of the flows over the airfoils are obtained with the help of the curved boundary treatment presented in Section 3. In order to have a view of the effect of curved boundary treatment on the numerical solution, a typical flow over a hemispherical headform [35][36][37][38][39] at the low Mach number 0.01 with zero angle of attack is simulated. Following the computational domain defined in the Reference [39], a relatively coarse mesh, having 27,976 tetrahedral elements (see Figure 11), is used in the present third-order (p = 2) accurate simulation. As shown in Figure 12a, a meaningful high-order accurate solution can be obtained with the present curved boundary treatment, while it seems that the DG method fails to obtain the correct solution without any curved boundary treatment, judging by the poor smoothness of the Mach contours evidenced in Figure 12b. Such a phenomenon can also be observed in the distribution of corresponding surface pressure coefficients against the ratio of distance along the boundary of the head diameter, as shown in Figure 13. Additionally in Figure 13 are the experimental data [40] and other researchers' results [35,39], which are included in order to have a view of agreement. In particular, the deviation between the numerical results with or without curved boundary treatment can be noted, particularly on the surface with relatively coarse meshes, which indicates that an appropriate curved boundary treatment is most helpful for seeking high-order accurate numerical solutions.

Flows over a Hemispherical Headform
It is pointed out that the above relatively accurate results of the flow foils are obtained with the help of the curved boundary treatment present In order to have a view of the effect of curved boundary treatment on the lution, a typical flow over a hemispherical headform [35][36][37][38][39] at the low 0.01 with zero angle of attack is simulated. Following the computational d in the Reference [39], a relatively coarse mesh, having 27,976 tetrahedra Figure 11), is used in the present third-order ( = 2) accurate simulation Figure 12a, a meaningful high-order accurate solution can be obtained w curved boundary treatment, while it seems that the DG method fails to ob solution without any curved boundary treatment, judging by the poor sm Mach contours evidenced in Figure 12b. Such a phenomenon can also be o distribution of corresponding surface pressure coefficients against the ra along the boundary of the head diameter, as shown in Figure 13. Additio 13 are the experimental data [40] and other researchers' results [35,39], cluded in order to have a view of agreement. In particular, the deviatio numerical results with or without curved boundary treatment can be note on the surface with relatively coarse meshes, which indicates that an appr boundary treatment is most helpful for seeking high-order accurate nume    ppl. Sci. 2022, 12, x FOR PEER REVIEW Figure 13. Surface pressure coefficients over hemispherical headform (Cao [35], Kaz periment [40]).

Flows over the ONERA M6 Wing
Finally, flows over an aerodynamic ONERA M6 wing [41][42][43] are validate the developed DG code in relatively more practical simulations. N pressible flow at Mach number = 0.01 and typical transonic shocked 0.84 with the same angle of attack ( = 3.06°) are selected and tested to s ties of coping with both nearly incompressible and compressible flows for developed. A relatively coarse mesh, having only 62,897 tetrahedron elem ated and used in the present third-order ( = 2) accurate flow simulation The corresponding meshes on the surface of the wing and in the symme presented in Figure 14.  (Cao [35], Kazem [39], Experiment [40]).

Flows over the ONERA M6 Wing
Finally, flows over an aerodynamic ONERA M6 wing [41][42][43] are considered to validate the developed DG code in relatively more practical simulations. Nearly incompressible flow at Mach number Ma = 0.01 and typical transonic shocked flow at Ma = 0.84 with the same angle of attack (α = 3.06 • ) are selected and tested to show the abilities of coping with both nearly incompressible and compressible flows for the methods developed. A relatively coarse mesh, having only 62,897 tetrahedron elements, is generated and used in the present third-order (p = 2) accurate flow simulations mentioned. The corresponding meshes on the surface of the wing and in the symmetric plane are presented in Figure 14.
ties of coping with both nearly incompressible and compressible flows fo developed. A relatively coarse mesh, having only 62,897 tetrahedron elem ated and used in the present third-order ( = 2) accurate flow simulatio The corresponding meshes on the surface of the wing and in the symm presented in Figure 14. In Figure 15, the obtained contours of the pressure coefficient on the wing are first presented for the mentioned transonic case, together with the pressure coefficient in the symmetric plane, in order to have a global vie flow fields. The corresponding distributions of the surface pressure co typical sections along the wingspan are also presented in Figure 16. It ca that the obtained results at the sections are all in a reasonable agreemen merical results [41] and experimental data [42] in view of the strength or captured shock. In Figure 15, the obtained contours of the pressure coefficient on the surface of the wing are first presented for the mentioned transonic case, together with the isolines of the pressure coefficient in the symmetric plane, in order to have a global view of captured flow fields. The corresponding distributions of the surface pressure coefficient at six typical sections along the wingspan are also presented in Figure 16. It can be observed that the obtained results at the sections are all in a reasonable agreement with the numerical results [41] and experimental data [42] in view of the strength or location of the captured shock.
Appl. Sci. 2022, 12, x FOR PEER REVIEW    Following the Reference [43], the obtained distributions of the surface pressure coefficient at the specific wingspan section η = 0.80 are then presented (see Figure 17) for the case of the nearly incompressible flow mentioned. It can be noted that the results captured are in good agreement with other numerical results [43]. The corresponding contours of the surface pressure coefficient of the wing, together with the field isolines of the pressure coefficient in the symmetric plane, are also presented in Figure 18 to have a view of flow fields. As shown in Figure 19, a significant improvement is again achieved in the convergence rate when preconditioning is adopted. The successful flow simulations of the above tested Mach numbers indicate that the DG methods presented function for both nearly incompressible and transonic shocked flows. captured are in good agreement with other numerical results [43]. Th contours of the surface pressure coefficient of the wing, together with th the pressure coefficient in the symmetric plane, are also presented in Fig view of flow fields. As shown in Figure 19, a significant improvement i in the convergence rate when preconditioning is adopted. The success tions of the above tested Mach numbers indicate that the DG methods pr for both nearly incompressible and transonic shocked flows.    contours of the surface pressure coefficient of the wing, together with th the pressure coefficient in the symmetric plane, are also presented in F view of flow fields. As shown in Figure 19, a significant improvement in the convergence rate when preconditioning is adopted. The succes tions of the above tested Mach numbers indicate that the DG methods p for both nearly incompressible and transonic shocked flows.    view of flow fields. As shown in Figure 19, a significant improvement i in the convergence rate when preconditioning is adopted. The success tions of the above tested Mach numbers indicate that the DG methods pr for both nearly incompressible and transonic shocked flows.

Conclusions Remarks
A high-order DG method for solving the preconditioned Euler equations with explicit or implicit time marching scheme is presented. A detailed description is given of a practical implementation of a precondition matrix of the type of Weiss and Smith, as well as of the DG spatial discretization scheme employed, with particular emphasis on the artificial viscosity-based shock capturing techniques and the proposed treatment of curved boundaries through adopting a NURBS surface equipped with RBF interpolation in order to propagate the boundary displacement to the interior of the mesh. The method is verified by considering a series of two-and three-dimensional test cases, including flows around a NACA0012 airfoil, a RAE2822 airfoil, a hemispherical headform, and an aerodynamic M6 wing. Numerical tests show that the present DG method functions for both transonic and nearly incompressible flow simulations. Numerical results also show that when a freestream Mach number reduces to an incompressible limit, the preconditioning adopted becomes necessary to accelerate convergence toward a steady state in the process of obtaining relatively accurate solutions. The deviation between the numerical results with or without curved boundary treatment confirms that an appropriate curved boundary treatment is helpful for seeking high-order accurate numerical solutions under the situation of relatively coarse meshes, particularly on the surface and with high curvature. The present DG method is implemented with unstructured meshes, and hence has the ability to cope with practical flows over more complex geometries. However, the computational cost grows rapidly with an increasing order of approximation. Thus, extension for further speedup, such as parallelization, can be considered. In addition, the governing equations can be easily updated with physical viscous terms to treat more practical viscous flows, which will be addressed in future research.