The Effects of Grid Accuracy on Flow Simulations: A Numerical Assessment

: High-quality, accurate grid generation is a critical challenge in the computational simulation of ﬂuid ﬂows around complex geometries. In particular, the accuracy of the grids is an effective factor in order to achieve a successful numerical simulation. In the current study, we present a series of systematic numerical simulations for ﬂuid ﬂows around a NACA 0012 airfoil using different computational grid generation techniques, including the standard second-order, fourth-order compact, and Theodorsen transformation approaches, to assess the effects of grid accuracy on the ﬂow solutions. The ﬂow solvers are based on the second- and fourth-order schemes for spatial discretizations and Beam-Warming linearization method for time advancement. The obtained grids, as well as the metrics and the corresponding numerical ﬂow solution for each grid generation technique, are compared and studied in detail. It is demonstrated that the quality and orthogonality of the grids is improved by using the fourth-order compact scheme. Moreover, the numerical assessment showed that the accuracy and the quality of the grids directly inﬂuence the numerical ﬂow solutions. Finally, the higher-order accurate ﬂow solvers are found to be more sensitive to the accuracy of the generated grid.


Introduction
One of the most critical subjects in computational fluid dynamics (CFD) simulation and modeling of flows around complex geometries is to generate high-quality grids on the flow region. The quality of the generated grid can strongly affect the accuracy of the numerical solution of a partial differential equation system [1][2][3][4]. In the past few decades, it has been established that the generation of high-quality grids is imperative for effective computations, and that can be equally important as improving the accuracy of the flow solver. Here, the grid quality means that the grid should be as smooth and orthogonal as possible in order to reduce the truncation error and carry out successful computational processes.
The structured and unstructured mesh generation techniques have been widely applied to a variety of flow solvers. Grid generation requires computational mapping from physical to computational space and geometry modeling. One of the most popular approaches for generating structured grids is solving a set of partial differential equations. The commonly used methods are elliptic and hyperbolic systems. Elliptic systems apply an appropriate control function to achieve grid smoothness and orthogonality. Another common method is conformal mapping containing an orthogonal transformation that yields nearly orthogonal grids in the physical space [5].
In addition to high-quality grids, it is necessary to employ high-order computational techniques to achieve high-accurate flow solutions. High-order accurate methods have been widely used for the solution of practical problems in fluid mechanics [6][7][8][9][10][11]. In Direct Numerical Simulation (DNS) and Large-Eddy Simulation (LES) of turbulent fluid flows and computational aeroacoustics, applying compact schemes is prevalent because of their spectral-like resolution, short size stencil, low dispersion and diffusion errors, and high-order accuracy [12][13][14][15][16][17][18]. It is, therefore, essential to carefully investigate the accuracy of these solutions.
Developing high-order unstructured mesh techniques is of particular interest, especially in the automatic mesh generation for commercial software. Furthermore, it is essential to develop high-order unstructured grids for high-order CFD solvers. Karman et al. [19] developed a three-dimensional mesh curving method based on an optimized mesh smoothing scheme, where they generated high-quality high-order meshes by transforming the smoothed nodes to the high-order meshes. Engvall and Evans [20] also investigated the generation of exact unstructured meshes from a geometrical point of view and proposed a surface reconstruction methodology to recover the exact computer-aided design geometries. Using the classical Winslow equations, Fortunato and Persson [21] introduced a new methodology for producing high-order unstructured curved meshes. In comparison with the optimization and elasticity analogies-based techniques, their methodology reduced the computational cost while preserving the accuracy of the curved elements. Additionally, Zhao et al. [22] developed a grid generation approach based on an improved radius basic function technique in order to generate a high-order curvilinear mesh. The proposed methodology was demonstrated to improve the accuracy of numerical solutions of curve-boundary problems.
The computational simulations of fluid mechanics problems are largely influenced by the accuracy and quality of the generated grids. Indeed, there exists an important link between the solution accuracy, mesh size, and mesh type [23]. Wang et al. [24] derived a fourth-order compact finite-difference scheme, which allowed unequal mesh-sizes to be used in all three dimensions for solving the Poisson's equation. Their results verified that the proposed method delivered formal fourth-order accuracy of numerical solutions. In another study, Shirayama [25] evaluated the effects of grid quality on the accuracy and convergence of flow solutions and improved the grid quality by applying an adaptive gridding technique based on weight functions and model equations. To achieve high-quality Cartesian grids, a new simplified ghost-point approach was introduced by Farooq and Müller [26] for embedded solid boundaries. They specifically improved the accuracy of the Cartesian grids for compressible Euler equations to simulate transonic flow over a circular bump. In a recent study, Karcz and Kacperski [27] further concluded that the numerical results of fluid flow in the agitated vessel depend significantly on the grid density and geometrical mesh quality.
The grid orthogonality and smoothness are key factors in high-quality grid generation techniques. A number of studies in the literature have investigated the impacts of orthogonality of the grids on the predicted results and solution accuracy for a wide range of diverse CFD applications (e.g., [28][29][30][31][32][33][34]). Matsuno [2] developed a new orthogonal hyperbolic grid technique based on a high-order upwind approach and showed that the orthogonality and robustness of the generated grids result in highly accurate flow solutions. Moreover, it was concluded that low-order accurate methods are too dissipative to achieve orthogonal grids. Sankaranarayanan and Spaulding [28] also investigated the effects of the grid orthogonality and resolution on the solution accuracy of shallow water equations. The error analysis showed that the truncation error of first-and second-order derivative terms are reduced for orthogonal grids. Besides, it was shown by Bagade et al. [29] that the performance and accuracy of transient flow simulations strongly depend on the orthogonality of the grid. Bagade et al. [29] further performed a grid sensitivity analysis to demonstrate that the grid characteristics and properties are effective parameters in high accuracy flow computations. Finally, in a comprehensive computational study, Akcelik et al. [30] investigated the effects of boundary type, grid quality control function, and grid properties on their proposed orthogonal grid generation methodology. This grid generation technique was shown to provide an acceptable balance between controlling the grid orthogonality and aspect ratio.
It is well established that a numerical solution of fluid flow strongly depends on the quality of the generated grid. Yet, there exist no reliable criteria for measuring grid quality and its effects on the flow solution accuracy. For the structured grids, the transformation metrics and the corresponding Jacobian appear in the transformed governing equations in computational space. In the present work, therefore, we aim to systematically investigate the effects of grid accuracy on the numerical solution of flow problems with complex geometries. In that regard, a high-order compact finite-difference scheme is employed to assess the accuracy and quality of the grid. To this end, the accuracy of the geometric metrics of the grid produced by the second-order method is first compared with those obtained by the fourth-order compact method. The orthogonal grids produced by Theodorsen transformation are then compared with those generated through the second-and fourth-order schemes. Finally, to assess the effects of grid accuracy on flow solutions, the numerical simulation of compressible flows around NACA 0012 with different grid generation techniques are implemented and investigated in detail. The rest of this paper is organized as follows. In Section 2, the governing equations for compressible flows are derived. The numerical technique, including spacial discretization, temporal advancement, and spatial filtering, is then explained in Section 3. Section 4 provides the results of the high-order flow solver for several test cases. The elliptic grid generation and conformal mapping procedures are explained in Sections 5 and 6, respectively. Section 7 compares the different grid generation techniques. Finally, the effect of grid accuracy on the flow solutions is studied in Section 8.

Governing Equations
The Navier-Stokes equations for compressible flows can be expressed in a general curvilinear frame of reference by ∂Q ∂t where t is the time, and ξ, η, and ζ are coordinate axes in the streamwise, spanwise, and vertical directions, respectively. Moreover, in the above equation, Q represents the conserved variable vector, F, G, and H denote the convective flux vectors, F υ , G υ , and H υ imply the viscous flux terms, and Re is the reference Reynolds number. The vector of conserved variables is defined as in which the subscript i denotes the corresponding quantity at node i and h refers to the grid spacing difference. For points close to the boundaries, we also considered high-order forward and backward compact schemes to preserve the tridiagonal form of the equations. The following formulations, therefore, were used to treat the physical boundary points at nodes i = 1 and i = I N , respectively: This fourth-order accurate scheme near the boundaries allowed us to retain the accuracy of the numerical solution over the entire computational domain.

Temporal Discretization
The implicit linearization method for the finite-difference scheme proposed by Beam and Warming [35] is applied for the time advancement. The implicit algorithm is further extended to three dimensions. Using Equation (1) and the trapezoidal formula (see [35]), we can express the time-marching equation as in which Q n+1 is the updated value of Q n at the iteration n + 1, and ∆t is the time step size. A local Taylor expansion about Q n results in a linear relationship between the convective flux vector F n+1 and the conserved variable vector Q n+1 in a way that, where A is the Jacobian matrix and A n = ∂F n /∂Q n . Moreover, the Alternating-Direction-Implicit (ADI) method [36] and the fractional step algorithm [37] were next applied to Equation (9) to take advantage of efficient multi-dimensional implicit algorithms. By applying similar linearization to other flux terms, Equation (9) can be reduced to where I is the identity matrix, and B n = ∂G n /∂Q n and C n = ∂H n /∂Q n are Jacobian matrices. Similarly, Using the procedure explained above, we were able to reduce the complexity of the problem from an immense matrix inversion (see Equation 9) to a series of tridiagonal matrix inversion problems for which efficient solution algorithms exist. To solve Equation (11), three one-dimensional sweeps are employed, in which ∆Q * * and ∆Q * are mediator variables, ∆Q n = Q n+1 − Q n , and RHS refers the right-hand side of Equation (11). Equation (13a)-(13c) are sweeps in the ξ, η, and ζ directions, respectively. By substituting Equation (6) with first-order derivatives into Equation (13a)-(13c), we obtain the following system of equations: Here, the updated value of the variable vector Q n+1 is obtained through calculating ∆Q n = Q n+1 − Q n from Equation (14c). The above set of equations describes a system of tridiagonal matrix inversion problem solved by Thomas algorithm.

Spatial Filtering
The nonlinear convective terms in the Navier-Stokes equations are the major source of instabilities and irregularities. The numerical instabilities result in a continuous and unrestricted generation of high-frequency wave modes. Therefore, incorporating a high-order low-pass filter into the numerical methodology is essential to handle such instabilities. Here, we implemented a high-order non-dispersive compact spatial filtering technique [38,39] as part of the numerical approach. The implicit filter applied to the solution vector can be expressed as where φ i is the filtered function of φ i at grid point i, β f is the coefficient of implicit terms , b n is the coefficient of explicit terms for different orders of filter, and 2 × N is the order of the filter. For a sixth-order filter, we have, where β f = 0.4. The filtering is carried out on the conserved variable vector obtained from Equation (14) after each time step and estimated using Thomas algorithm, a typical tridiagonal matrix solver. Furthermore, the sixth-order forward and backward filtering technique was applied for near boundary points, The coefficients of the right-hand side of the above equations, a n , are derived using Taylor and Fourier-series analyses and can be found in, for example, Gaitonde and Visbal [40]. The values at the boundaries, nodes i = 1 and I N , are determined explicitly through prescribed boundary conditions. Thus, the tridiagonal structure of the filtering was maintained entirely.

Flow Solution Results
In this section, we first carry out a series of numerical simulations to verify and validate the current high-order compact flow solver against experimental results for fluid flow over a flat plate and a NACA 0012 airfoil. The accuracy of the numerical method and hence the flow solver is also assessed for the flow inside a Shubinn nozzle.

Numerical Validation
The two-dimensional laminar viscous flow over a flat plate with upstream conditions of Ma = 0.2 and Re = 200 × 10 3 is solved using the high-order compact method. The profiles of horizontal velocity and friction coefficient are shown and compared with the classical Blasius solution for the flat-plate problem in Figure 1. The velocity profile is plotted against the non-dimensional position variable η * ; in the Blasius solution, η * = y √ u ∞ /νx is a non-dimensional quantity that combines both x and y positions. The skin fiction coefficient, C f = 2τ w /ρu ∞ where τ w is the wall shear stress, is also plotted along the length of the plate in Figure 1b. As it can be observed from Figure 1, there exist an excellent agreement between the current numerical solution and those obtained from the Blasius boundary layer solution for flow over a flat plate. The rise of the skin friction coefficient at the trailing edge of the flat plate is explained by the triple-deck theory that is not considered in the Blasius solution (see [41] for more details). Further, we examined the fourth-order compact flow solver by simulating the flow field around a NACA 0012 airfoil with a Mach number of 0.85 and a Reynolds number of 500. The profiles of the pressure and skin friction coefficients are presented in Figure 2. These profiles are well-matched with those of GAMM-Workshop [42]. Besides, we considered an inviscid supersonic flow around the NACA 0012 airfoil with Ma = 1.2 to validate the correct behavior of the solver for supersonic inviscid flows as well. Figure 3 presents the pressure coefficient distribution over the airfoil, where the numerical results are compared with those obtained by Arias et al. [43]. Again, there is an excellent agreement between the two. It should be noted here that the grid around the NACA 0012 airfoil is generated using the elliptic method in x and y directions where N x = 129 and N y = 37. In the z direction the grid is uniform with N z = 37. In this study, the flow results of the NACA 0012 airfoil are presented in mid-plane in the z direction. The inflow and outflow boundary conditions are applied at the inlet and outlet, respectively. The wall boundary condition as well as adiabatic condition is used for the surface of the airfoil.

Accuracy Study
In order to assess the accuracy of the current compact finite-difference solver, we simulated the one-dimensional flow inside the Shubbin nozzle (see [44]) and compared the results with its corresponding exact solution. The cross-sectional area of a Shubbin nozzle is defined as The inflow conditions are described as ρ in = 0.05008261, p in = 0.2712900, u in = 1.0991840, Ma = 1.262214, and for the outflow p out = 0.5156000. The pressure changes inside the Nozzle is demonstrated in Figure 4a. The analytical solution to the problem is also indicated in this figure, and it can be observed that the result of the present compact method agrees quite well with the exact solution.  Further, the L2-norm error is calculated to estimate the accuracy order of the solution. The L2-norm error can be defined as where p and p e are the values of the pressure for the compact method solver and the exact solution, respectively. Figure 4b shows the L2-norm error, and it can be observed that the compact method presented in this study provides a fourth-order accurate solution.

Elliptic Method
The elliptic grid generation method is based on solving a pair of Poisson equations in the computational domain with the desired grid points on the boundary of the physical domain, where the interior points can be determined by the solution of the following equations [45]: Here, P and Q are the control functions. By interchanging the roles of dependent (ξ, η, ζ) and independent variables (x, y, z), one can arrive at where and J is the transformed Jacobian defined by J = x ξ y η − x η y ξ . In order to control the grid attributes in the boundaries, grid spacing, and the orientation of coordinate lines, the control functions P and Q are defined as [46] where a L , a U , b L , and b U are positive constants that control the orthogonality and stretching grid nodes near the boundaries. Moreover, the functions P L , Q L , P U , and Q U are defined near the lower and upper boundaries (see Sorenson [46] for more detail). Furthermore, to impose the orthogonality condition near the boundaries, the following equation is employed: ∇ξ · ∇η = |∇ξ| · |∇η| cos θ.
After applying some simplifications, Equation (24) is reduced to where x η and y η can be expressed for θ = π/2 as in which Furthermore, in order to avoid the numerical instability, especially when the control functions P and Q impose large values, the following upwind scheme is utilized for the first-order spatial derivatives: For the second-order derivatives, the fourth-order compact scheme is used [12], Equations (28) and (29) are also applied in the η direction. As a result, there are five coupled equations at each point, namely, where subscripts i and j indicates the ξ-and η-direction, respectively. The above equations can be written in the following matrix form: where and The coefficients a 1 to a 6 are reported in Table 1. The above system of equations is then solved using the Jacobi point iteration method. A similar approach is also applied in the y-direction. Again, it is worth mentioning here that the grid is uniform in the z-direction.  (32) and (33).

Conformal Mapping
In the calculation of flow around airfoils, it is often useful to map the exterior flow field of an airfoil conformally onto the interior of a circle. A direct, point-by-point transformation of an airfoil onto a circle can be performed in two distinct steps using the method suggested by Theodorsen [47,48]. The first step is to transform the airfoil onto a quasi-circle by Joukowski transformation, which can be written as where ζ and z are the complex variables associated with the airfoil and quasi-circle planes, respectively. The second step is to map the area outside of the quasi-circle onto the interior of the circle using the following transformation: in which z 0 is the center of the base circle of the radius r 0 . This transformation not only shifts the origin of the coordinate system to the center of the quasi-circle but transforms the near-circle into an almost perfect circle. Transforming from a quasi-circle to a perfect circle with an unknown radius imposes the following boundary conditions: where r 0 = e ψ 0 is the radius of the objective circle. Further, ψ 0 , A n /r n 0 , and B n /r n 0 can be calculated using Fourier transformation relations. The subscript Ω implies that the equations are valid at the boundary points. Moreover, the quasi circle and perfect circle planes are defined as z = e ψ+iθ and z = e ψ 0 +iφ , respectively.
The conformal mapping procedure is schematically illustrated in Figure 5, where the airfoil section with an arbitrary shape is first transformed onto a quasi-circle using the Joukowski transformation (see Equation (34)). The quasi-circle is then transformed into a circle using Equation (35); the areas outside the closed curves in the two later planes are conformal representations of the area outside the airfoil in the physical plane.

Grid Generation Results
To study the impact of different grid generation methods, we first considered the grid generated in the area enclosed by two semicircles using standard second-order and fourth-order compact schemes, as shown in Figure 6. Here, the term standard second-order scheme indicates that the second-order finite-difference method is applied to produce grid points by solving Equation (20). Figure 7a,b also represent the deviation of the coordinate lines generated by the second-and fourth-order finite-difference schemes from the orthogonal lines in the polar coordinates close to the lower and upper boundaries. It is observed that the error corresponding to the fourth-order compact method is significantly less than the second-order scheme by almost, on average, 90%.   Next, the O-and C-type generated grids around a NACA 0012 airfoil are respectively represented in Figure 8a,b. To generate these two types of meshes, we used three different schemes, including the fourth-order compact scheme, the Theodorsen's scheme, and the standard second-order scheme. The set of equations expressed in Equation (21) were solved to obtain the elliptic O-and C-type grids using both second-order and compact finite-difference schemes. It can be observed from Figure 8 that the high-order compact method generates high-quality grids, which are as orthogonal as Theodorsen grids. Further, Figures 9 and 10 represent a comparison of the transformation metrics ξ x and ξ y that are obtained by different methods and located at x = 0.5 close to the airfoil profile for both O-and C-type grids. The observed differences in the transformation metrics are of significant importance since the coordinate metrics appear in the transformed governing equations (see Equation (1)), and consequently, they can strongly affect the numerical solutions. The metrics calculated using the standard second-order and fourth-order compact schemes for the C-type grid are shown in Figure 9. In addition, the metrics of Theodorsen's grid are compared with those of second-order and compact methods in Figure 10 for the O-type grid. Overall, the transformation metrics ξ x and ξ y of the corresponding compact method are well-matched with those obtained by Theodorsen's approach indicating the capability of the presented high-order compact technique for accurate grid generation.
x x y 19 19.  In order to assess the quality of the grids, the smoothness and orthogonality properties can be considered as follows: in which I o and I s indicates the orthogonality and smoothness of the grid, respectively. The lower value of I o and I s denotes a higher quality of the generated grid. The computed values of I o and I s are reported in Table 2 for the standard second-order, fourth-order compact, and Theodorsen approaches. Here, the orthogonality factor I o for the compact method is improved by approximately 5% in comparison with the second-order scheme. Such improvements in the orthogonality factor are expected to be more substantial for grid points near the airfoil surface. Moreover, it is observed that the smoothness factor I s for the elliptic grids is generally less than the Theodorsen grid. This is due to the fact that the grid points are solved through a system of coupled equations (Equation (20)) for the compact and second-order elliptic grid generation.

The Effect of the Grid Accuracy on Flow Solutions
In this section, we investigate the effects of grid accuracy on the numerical solutions of fluid flow around a NACA 0012 airfoil. The high-order compact flow solver was implemented in which the grid was generated through three different methods, including the standard second-order method, the compact scheme, and the Theodorsen technique. The distributions of pressure coefficient C p = 2(p − p ∞ )/ρu 2 ∞ and skin friction coefficient C f = 2τ w /ρu 2 ∞ on the airfoil surface are respectively represented in Figure 11a,b for the subsonic, laminar flow (see, for example, [49][50][51]) with Ma = 0.5 and Re = 10 × 10 3 . The fourth-order compact grid generation method shows improvements in both pressure and skin friction coefficients compared to the second-order method and presents an excellent agreement with those results obtained from the Theodorsen grid generation technique. In addition to the subsonic flow, the supersonic flow over a NACA 0012 airfoil was also calculated. Figure 12 exhibits the pressure coefficient distributions for the inviscid supersonic flow with a Mach number of 1.2. Again, the zoomed-in views of this plot clearly demonstrate that the compact method agrees well with the Theodorsen method compared to the results of the second-order method. Although these figures indicate that there exists possibly an insignificant difference (approximately 5%, on average) between the results obtained from the second-order and the compact schemes, it is expected that the flow solutions are more affected by the grid accuracy for more complex geometries and turbulent flows [52,53]. As a final case, we examined the subsonic flow over the NACA 0012 airfoil but with Ma = 0.85 and Re = 500. Here, however, we employed the second-order flow solver in addition to the high-order compact solver to emphasize the importance of high-order solvers. The grid generation techniques remained the same as previous, i.e., second-order, fourth-order, and Theodorsen methods. Figures 13  and 14 present the results of flow solutions for different flow solvers and grids. In general, the results of the high-order grid implemented with the high-order flow solver demonstrates a distinct difference compared to the results of the second-order grid coupled with the second-order solver. Moreover, the high-order compact finite-difference flow solver is more sensitive to the grid accuracy than the second-order solver. It can be argued that the truncation errors of the second-order flow solutions eliminate the effect of the grid accuracy. It can be inferred from Figure 13 that the pressure coefficient might be likely more susceptible to the accuracy of the flow solver than the skin friction coefficient. Finally, as shown in Figure 14, both horizontal and vertical velocity profiles present slight variations by either changing the accuracy of the grid generation technique or the flow solver.
x/c   4th order solver, 4th order grid 4th order solver, 2nd order grid 2nd order solver, 4th order grid 2nd order solver, 2nd order grid (b) Figure 14. Comparison of the (a) horizontal and (b) vertical velocity components obtained using second-and fourth-order flow solvers with different generated grids for flow around a NACA 0012 airfoil with Ma = 0.85 and Re = 500.

Conclusions
In the present work, we presented a systematic numerical investigation to examine the influences of grid accuracy on the numerical solution of fluid flow problems with complex-shaped geometries. To this end, the flow field around a NACA 0012 airfoil was simulated using a high-order flow solver with various grid generation techniques. The grids around the NACA 0012 airfoil were generated using the standard second-order method and the fourth-order compact finite-difference scheme based on the elliptic approach. Moreover, following the Theodorsen approach, a conformal mapping is further employed to generate an orthogonal grid around the NACA 0012 airfoil. Here, a clear difference was observed between the grid metrics of the second-and fourth-order finite-difference schemes. Besides, compared to the Theodorsen's orthogonal grid, the orthogonality of the grids generated by the fourth-order compact scheme is noticeably better than that of the second-order approach.
The effects of the grid accuracy on flow solutions were then evaluated for the flow fields around a NACA 0012 airfoil for which the computational grids were generated using the second-order, fourth-order compact, and Theodorsen schemes. The flow solvers were based on the standard second-order and fourth-order compact finite-difference methods. The results demonstrated that the accuracy of flow solutions are directly linked to the accuracy of the generated grids. In that regard, the grids generated using the fourth-order compact scheme showed improved results in comparison with grids generated by the second-order method. In particular, the pressure and skin-friction coefficients obtained from the numerical simulations with grids produced by the high-order compact method were well-matched with those obtained from Theodorsen conformal mapping. Furthermore, it was found that the fourth-order compact flow solver is more sensitive to the grid accuracy than the second-order one.

Conflicts of Interest:
The authors declare no conflict of interest.