Computational fluid dynamics (CFD) is essentially a series of approximations; the role of the practitioner is managing the magnitude of the errors in these approximations. In CFD, one of the most basic approximations lies in the type of cells used for the simulation. This article will address some of the potential pitfalls encountered with a generic cell-centered finite volume scheme.

Along those lines, the choices made in a CFD workflow will ultimately determine both the applicability and the credibility of the solution. CFD practitioners have no shortage of decisions to make in regard to their simulations: solution algorithms, turbulence models, convergence schemes, etc. Even before these parameters are considered, the choice of solver must be made from an abundance of possibilities that may include different physics and use varying degrees of complexity to handle, mitigate, or reduce numerical errors.

Step back even further, before CFD even enters the picture, and there is the choice of the discretization scheme. Although the decision of which type of computational mesh to use is perhaps one of the least “physical” approximations that occurs – meaning for a given resolution, the cell type or quality does not alter the physics being considered – the choice can drastically alter the results in potentially meaningful ways. Decisions regarding the computational grid require both knowledge of your solver – What type of cells does it support? How susceptible is it to grid-related errors? – and control over your computational mesh.

Ideally, grids would consist of only flow-aligned, orthogonal hexahedral cells, but perhaps the geometry under consideration is simply too complex to consider such an approach. Further, a grid that is insufficiently refined in areas of high gradients can adversely affect the solution by under-representing the shear present in the flow. In particular, poorly resolved wakes or separation regions can even have a serious impact on integrated quantities like lift and drag. Understanding the consequences of the choices regarding the computational mesh is essential for efficiently generating high-fidelity simulation results.

Assessing Grid Induced Error
It is not a trivial task to assess the error introduced by the grid for an arbitrary case because of the complexity of the dynamics inherent in solutions to the Navier-Stokes equations. For this reason, it was decided to try to isolate individual grid effects with simpler test cases. That is, efforts were made to evaluate the error introduced by the grid on a term-by-term basis and remove any convolution arising from multiple terms (advective, diffusive, sources, etc.). The convective derivative in the Navier-Stokes equations, (v ⋅ ∇) v, seems a fitting candidate to begin with when examining the effects of the grid because it is in many ways the source of the richness of the dynamics (due to its non-linear nature).

However, this term from the Navier-Stokes equations could not be considered in the absence of a viscous dissipation term since a true convective derivative operating on a velocity field would eventually give rise to turbulence in the presence of a persistent gradient. For this reason, a passive scalar advection case was considered. Such a scheme contains the same mathematics in terms of computational errors, but the operator affects a scalar field independent of the background flow. For these tests, an advection-only solver implemented in OpenFOAM® was employed. OpenFOAM is a collection of finite volume algorithms that can be flexibly pieced together to take into account various physics. Although other solvers are certain to have different (and perhaps more sophisticated) algorithms to handle the types of errors encountered in these studies, the intent of this work is to demonstrate baseline discrepancies between grid types for a given solver, with the understanding that more sophisticated algorithms may improve the quality of the solution, but the trends will remain.

Finite volume codes work on the principle of balancing fluxes within differential control volumes (cells), and numerical errors are introduced when calculating the flux between adjacent cells. Node-centered codes are also subject to deficiencies in grid quality, but for illustrative purposes, these examples will be based on cell-centered formulations. A computed solution is susceptible to errors arising from many grid quality short-comings including non-orthogonality or large jumps in cell volume. However, for a purely advective test case, cell skewness is the critical metric to consider. This is due primarily to interpolation errors introduced when projecting cell-centered values to cell faces when computing advective fluxes. The source of this error can be visualized in Figure 1, where the face center is offset a distance δ from the cell centers. In the presence of a vertical gradient, the computed flux through that face will be underrepresented due to this displacement.

Figure 1: In the presence of a vertical gradient, the offset δ between the face center and the line connecting cell centers (cell skewness) can affect the accuracy of the computed flux for the shared face, introducing numerical errors.

As suggested, this type of error predominantly affects advective terms such as the convective derivative in the Navier-Stokes equation. This term drives many of the non-linearities present in turbulent flow, and any misrepresentation can adversely affect the simulation accuracy, particularly for separated flows with off-body gradients. The skewness in question can be quantified in Pointwise using the Centroid Skewness metric. Ideally, values of the Centroid Skewness will be as close to zero as possible, as non-zero skew will result in some small accumulation of error. As an aside, the Equiangle Skewness metric is computed quite differently in Pointwise, and tends to inhibit convergence when excessively large for some solvers. Values for Equiangle Skewness less than 0.7 are suggested, although some codes can handle cells with higher values.

Note that in the absence of gradients, cell quality has little effect on the solution, since any interpolation does not introduce appreciable error. However, in regions of sharp gradients, the net effect of this type of cell quality deficiency is most broadly described as “numerical diffusion,” as the numerical artifacts produce artificial diffusion blurring the sharp gradients present in the flow.

Advection of a Thermal Gradient in 2D
In order to demonstrate these effects, first consider a simple 2D case in which a thermal profile is advected across a simulation domain as illustrated in Figure 2.

Figure 2: This illustration shows the test case consisting of a rectangular domain and imposed thermal gradient that is to be advected across the simulation volume.

Since physical diffusion was deactivated by construction, the sharp gradient should be advected across the domain without any change to the shape of the initial temperature profile in the absence of any artificial (numerical) diffusion. By examining the simulation result with the expected value, it is possible to quantify the error arising from numerical diffusion.

To begin, consider the error from a purely quadrilateral grid shown in Figure 3. The small error present in the solution is simply the baseline interpolation error resulting from the finite discretization size. That is, the error present in the solution is wholly due to grid resolution limitations, not a result of poor cell quality.

Figure 3: The baseline error in the solution computed on an orthogonal grid with quadrilateral elements shows errors due entirely to the limit in resolution.

These results can then be compared to the computed solutions from two unstructured grids with triangular elements. The first, generated with a Delaunay algorithm, has a distribution of cells with non-zero skewness, as shown in Figure 4, and the effect of skewness on the accuracy of the solution is clearly seen. A heatmap of the grid's centroid skewness is shown in Figure 5.

Figure 4: The error in the solution computed on a Delaunay triangular grid is considerably higher than the error from the structured grid.

Figure 5: Examining the heatmap of centroid skewness illustrates the source of the errors.

The second unstructured grid, generated using an Advancing Front algorithm, contains equilateral triangles with zero skewness in the majority of the computational domain as shown in Figure 7. By comparing the solution accuracy in Figure 6 with the heatmap of cell skewness in Figure 7, error is only introduced to the solution in areas where opposing fronts collide, resulting in local non-zero skewness.

Figure 6: The error in the solution computed on an Advancing Front triangular grid is considerably less than that computed from the Delaunay mesh, but still notably higher than the structured grid.

Figure 7: The heatmap of centroid skewness for the Advancing Front mesh clearly shows the origin of the numerical errors.

Figure 8 demonstrates the skewness that arises from the irregularity inherent in the Delaunay algorithm.

Figure 8: Close examination of Delaunay (a) and Advancing Front (b) triangular meshes highlights the source of the non-zero skewness for Delaunay meshes.

It should be noted that non-zero skewness is likely an unavoidable consequence of modeling any geometry of useful purpose, but this study is aimed at providing an incentive to researchers wishing to minimize grid-related errors. Intuitive foresight into the anticipated flow field should direct an experienced CFD practitioner to isolate critical regions of interest and maximize the quality of the cells in that area.

Advection in 3D
These simplistic two-dimensional examples can be extended to three dimensions, as shown in Figures 9-11. In Figure 9, the advection of a Gaussian thermal profile is shown on a hexahedral grid with 100³ cells (101³ nodes). As in the structured case shown in Figure 3, the resulting error is due solely to the limit in resolution, demonstrated in Figure 10.

Figure 9: The cubic domain with imposed Gaussian temperature profile extended the simple problem to three dimensions.

Figure 10: As in 2D, the error in the solution computed on a grid consisting of 1.0M orthogonal hexahedral cells was due entirely to the limit in resolution.

Figure 11 shows the results from a tetrahedral grid with over 6.5 million cells. Despite the significantly higher cell count, the maximum error present in the solution is an order of magnitude higher than that in the hexahedral grid because of the effects of numerical diffusion.

Figure 11: The error in the solution computed on a grid consisting of 6.5M tetrahedral cells showed a considerable cumulative error as before in two-dimensions, despite a significantly larger number of elements.

These results are not to say that unstructured grids are to be avoided. As mentioned at the onset, these results are from an isolated test case specifically designed to highlight the effects of the grid on the solution. Further, the errors quoted in the figures, while quantitative in nature, are still qualitative measures suited only for direct comparisons. In practice, modern CFD solvers can produce quite robust solutions on both hexahedral and tetrahedral grids (depending on the solver, of course). Further, the simplistic algorithm used in computing these results lacked many existing techniques for mitigating effects of poor cell quality that may be present in your solver. Lastly, the errors being discussed are generally diminished with increasing grid resolution. Therefore, if the mesh is sufficiently refined, these errors are inconsequential.

An Ongoing Study
To reiterate, the purpose of this work was to clearly demonstrate the effects of cell type and quality and provide motivation and guidance for obtaining more accurate solution results via improved grid quality. This work is ongoing and will eventually be expanded to the other operators found in the Navier-Stokes equations. Nothing can replace a dedicated grid independence study for ultimate rigor, but understanding the basic interplay between the grid and the solution and the consequences of the implicit approximations being made should help you construct better quality grids for any application.

References
The preliminary work in two dimensions was performed in preparation for the 9th Annual OpenFOAM Workshop earlier this year. More details regarding the three dimensional effects will be presented at the 67th annual American Physical Society Division of Fluid Dynamics conference on 23-25 November.

1. “How grid quality affects solution accuracy,” J. Rhoads and T. Carrigan, 9th Annual OpenFOAM Workshop.
2. “How grid quality affects solution accuracy,” J. Rhoads, 67th Annual American Physical Society Division of Fluid Dynamics Meeting.