Assessing Validity of Mesh Refinement Sequences, with Application to DPW-III Meshes Carl Ollivier-Gooch Advanced Numerical Simulation Laboratory The University of British Columbia Elucidating the relationship between mesh quality and solution quality is a particularly thorny problem, but one that is receiving increasing attention, especially in light of the challenges with accurate solutions for the AIAA Drag Prediction Workshop problems which have attributed, at least in part, to mesh quality. This paper focuses on the assessment of whether two meshes are sufficiently geometrically similar to be considered as members of the same family for purposes of mesh refinement studies. The local length scale, aspect ratio, and cell orientation for each cell in each of a pair of meshes is computed. Data is projected from the coarse mesh to the fine mesh for comparison between the two meshes. Based on comparisons of meshes carefully generated to be nearly ideal mesh refinement pairs, guidelines are given for how well one can expect two meshes in a refinement pair to correspond to each other in terms of cell size, aspect ratio, and orientation. Examples are given from Drag Prediction Workshop meshes.

I. Introduction Among other results of the Third AIAA Drag Prediction Workshop (DPW-III), the panel concluded that issues of mesh quality and mesh refinement are among the chief outstanding difficulties in accurately computing the aerodynamic flows of interest in these workshops.9 Considering the relatively mature state of scientific computing on unstructured meshes, it is surprising how little we know about the effects of mesh quality on solution quality. We know theoretically that commonly used schemes will asymptotically converge to correct solutions to partial differential equations we wish to solve, but the details are the effects of, for example, geometric mesh element quality or differences in the size of adjacent elements are currently poorly understood. There has been some work in the area of the effect of mesh quality on interpolation error for finite element methods,7, 8 but to the best of the author’s knowledge, there is no analogous work for finite volume methods. Also, Knupp3 described some useful solution-based quality metrics, though these have seen little application in practice. The overall problem of elucidating the interaction between mesh and solution is a large and complex one. This paper seeks to address a small but important corner of that problem: the problem of determining whether a sequence of meshes truly constitutes a mesh refinement sequence for purposes of a mesh convergence study. Our goal here is, given a sequence of meshes, to determine whether those meshes are geometrically similar enough and have a uniform enough ratio of length scale from mesh to mesh to be considered plausible for use in a mesh refinement study. For purposes of this study, two meshes will be said to form an ideal refinement pair if their cell sizes are uniformly proportional across the entire domain and if cell aspect ratio and orientation is the same in anisotropic regions of the mesh. Realistically, ideal refinement pairs will not exist; we will also explore what bounds are reasonably attainable for variation from the ideal. We will not, in this context, consider classical mesh quality measures, which presumably the generators and users of the mesh will routinely evaluate using readily available tools. Also, we will not consider any solution related tasks; the test we propose here is a strictly a priori test of the meshes themselves. To determine whether a series of meshes constitutes a valid mesh refinement sequence, we seek to compare cell size, shape, and orientation between pairs of meshes. The size, shape, and orientation of a cell in an unstructured mesh can be computed in more than one way; one viable approach is to compute the moments of the cells, and solve an eigenvalue problem to find the size, aspect ratio, and principal directions of the cell (see Section II). Once these quantities are known for all cells in a mesh, comparison between meshes can be accomplished by projecting data from one mesh to another. To validate the methodology and to explore the outcomes that can be expected both for carefully constructed refinement pairs and for meshes which have been deliberately generated to have clear flaws

1 of 14 American Institute of Aeronautics and Astronautics

in mesh refinement, a number of simple examples are presented in Section III. Results evaluating whether sample published meshes used in DPW-III are in fact reasonable mesh refinement sequences are given in Section IV.

II.

Methodology

Cell size Si is by far the easiest property that we wish to compute; size formulae for triangles, quadrilaterals, tetrahedra, prisms, and √ hexahedra are well-known and need not be repeated here. We will compare cell length scale, which we define as d Si in d dimensions. II.A.

Measuring Cell Anisotropy, and Orientation

Cell anisotropy and orientation can be assessed from the normalized second moments of the cell. These moments can be written as, for example, Ixx

=

Ixz

=

1 x2 dS Si i Z 1 xz dS Si i Z

where the integral is over area in two dimensions and over volume in three dimensions. In practice, these moments can also be evaluated by applying Gauss’s Theorem and integrating around cells: Ixx

=

Ixz

=

1 x3 nx dl Si ∂ i 3 Z x2 z 1 nx dl Si ∂ i 2 Z

where nx is the x-component of the normal to the cell boundary, and dl indicates integration around the boundary4 (either a line integral in two dimensions or a surface integral in three dimensions). In two dimensions the moment of inertia tensor of the cell can be written as # " Ixx Ixy M= Ixy Iyy This form differs from the standard form of the moment of inertia tensor only in that each term has been divided by the size of the cell. We can determine the cell anisotropy either by finding the eigenvalues and eigenvectors of the moment of inertia tensor or by applying a Mohr’s circle analysis to find the principal moments and their directions. Here, we will illustrate the latter. From Mohr’s circle analysis, found in any introductory text mechanics, we know that ron solid 2 Ixx −Iyy I +I 2 . The minimum and the inertia tensor represents a circle with its center at xx 2 yy , 0 and a radius of + Ixy 2

maximum principal moments are, therefore

Imax

=

Imin

=

Ixx + Iyy +

q 2 (Ixx − Iyy ) + 4Ixy

q2 2 Ixx + Iyy − (Ixx − Iyy ) + 4Ixy 2

I −I

The principal axes are rotated by an angle of 12 arctan Ixy / xx 2 yy ratio of maximum and minimum moments: r Imax AR = Imin

. Finally, the aspect ratio is the square root of the

As expected, this definition gives AR = l/h for rectangle of length l and height h. In three dimensions, there is no easy analogy to Moore’s circle, and so instead we solve the 3 × 3 eigenvalue problem numerically. In this case, we compute two aspect ratios, based on the ratios of maximum and median eigenvalues and of the median and minimum eigenvalues. For a typical 3D viscous mesh, these correspond to the cross-stream to 2 of 14 American Institute of Aeronautics and Astronautics

streamwise aspect ratio and the streamwise to wall normal aspect ratio, respectively. Orientation comparisons in three dimensions are most easily done by retaining the unit vectors in the directions of the median and minimum eigenvalues — that is, the stream wise and wall normal directions — and determining the differences between these directions for a pair of meshes by computing dot products. II.B.

Practical Considerations

The description in the preceding section is deceptively simple. In practice, because unstructured meshes are intrinsically non-smooth at the fine level, data must be adjusted somewhat to allow a reasonable comparison. With respect to cell size, averaging is applied to construct an approximate average cell size over a region somewhat larger than an individual cell. For each cell, consider the neighborhood of that cell that would be used in a quadratic reconstruction of solution.4 In practice this results in a neighborhood with about 10 cells in two dimensions and 15 to 18 cells in three dimensions. A weighted average over the stencil is used, with weights computed based on inverse square of distance of a cell centroid from the centroid of the central cell. As Figure 1 illustrates, this technique improves smoothness of measured cell size, which will make comparison of sizes between meshes much easier to do reliably.

(a) Unsmoothed length scale

(b) Smoothed length scale

Figure 1: Effect of weighted averaging on measured cell size. When measuring anisotropy, it is important to keep in mind that even an ideal unstructured mesh will have cell aspect ratios that differ somewhat from one. For instance, in two dimensions, typical high quality mesh generators still produce triangles with √angles as small as 30 degrees. This can lead to aspect ratios as large as 3 (for a 30-30120 triangle) or even 3 + 2 3 3 (for a 30-75-75 triangle), even though such triangles appear not to cause significant difficulties for flow solvers. In three dimensions, the situation is even more complex. Experimentally, we find that isotropic meshes that have good geometric quality by other measures typically have aspect ratios as high as about 2.5. Therefore, in this paper all cells with aspect ratios less than 2.5 will be treated as isotropic — that is, as having unit aspect ratio. In three dimensions, this rule is applied independently to the two measured aspect ratios. Finally, when the mesh is nearly isotropic — that is, when the aspect ratio falls below the isotropic threshold previously imposed — orientation is clearly not an issue. Therefore, when either mesh meets this threshold locally, the orientation difference will arbitrarily be set to zero. II.C.

Comparison Between Meshes

The final technical issue to be considered here is how to compare size, aspect ratio, and orientation between meshes. When comparing cell sizes, our goal is to determine whether all cells in one mesh are proportionately smaller than all cells in the other mesh. The constant of proportionality depends on the relative number of cells in the mesh. Specifically, the cell sizes in each mesh should be proportional to the domain size divided by the number of cells. Because we are comparing cell characteristic length, for meshes covering effectively identical domains, we expect the length scale to behave as follows: r N2 L1 = 2D: L2 N1 r L1 N2 3D: = 3 L2 N1 3 of 14 American Institute of Aeronautics and Astronautics

where Ni is the number of cells in mesh i and Li is the local length scale for mesh i. To confirm this, we first project the length scale from mesh 1 to mesh 2. We do this by performing an approximate numerical integration of the length scale on mesh 1 over the cells of mesh 2. This integration is problematic for high-aspect ratio cells and curved boundaries, where the fine and coarse domains are significantly different, making the correspondence between fine and coarse mesh control volumes for anisotropic meshs especially bad near the wall. To prevent numerical problems for this situation, we evaluate the numerical integration using injection rather than by reconstructing the solution on one mesh and integrating more exactly. While the projection we use is inexact, the “more exact” approach of reconstructing and integrating has its own problems, because the mesh quantities we are transfering are discontinuous and the two meshes are not nested. A more accurate and robust approach to data transfer is definitely desirable. After the projection is completed, we compute a refinement size quality Qr , defined as: r L1 d N1 Qr ≡ (1) L2 N2 where d is the number of space dimensions; this quantity should be identically one in theory. For aspect ratio, we expect the same value in each location on pairs of geometrically similar meshes. In this case, following a projection of aspect ratio from one mesh to the other, we simply take the ratio of aspect ratios for the two meshes, again expecting the results to be exactly one for ideal meshes. Finally, for orientation, we expect that anisotropic cells will be oriented at the same angle. We can compute the angular deviation between orientations based on the dot product or cross product of these vectors.

III. Assessing Mesh Pairs In this section, we will examine the behavior of the mesh pair analysis techniques described in Section II. We will begin by examining pairs of meshes that have been carefully generated to be nearly ideal mesh refinement pairs, for both isotropic and anisotropic meshes. Then we will examine the effects of particular deliberate pathologies in refinement pairs to characterize the signature of these pathologies in the current methodology. Together, these two steps will provide a basis for interpretation of results when refinement pairs of unknown quality are examined in Section IV. III.A.

Near-Ideal Pairs of Isotropic Meshes

(a) Coarse mesh size

(b) Fine mesh size

(c) Refinement ratio (Equation~1)

Figure 2: Isotropic Mesh Refinement Pair in a Square with a Clipped Corner We begin by examining a refinement pair for a mesh in a square with a small corner (1% of the side length) cut off. These two meshes were both generated using the GRUMMP mesh generator,2, 5, 6 which uses a Delaunay insertion algorithm. GRUMMP automatically computes a geometric length scale based on the boundary curvature and proximity to one boundary to another. Also, users can specify that the final cell size be some fraction of this geometric length scale. In this case, the coarse mesh in the refinement pair was created with a mesh length scale one-fourth as large as the geometric length scale, while the finer mesh was created with a mesh length scale one-eighth as large as the geometric length scale. The course mesh has 1021 cells, while the fine mesh has 4045 cells; these meshes are shown in Figure 2. These meshes are isotropic, and as expected our analysis of aspect ratio correctly identified this: the largest measured aspect ratio in either mesh was 2.34. Also, of course, we expect that the cell length scales will be almost exactly proportional to one another. Figure 2 shows the results of length-scale analysis for this pair of

4 of 14 American Institute of Aeronautics and Astronautics

meshes. Figure 2(a) shows the length scale for the coarse mesh, and Figure 2(b), the length scale for the fine mesh; these are shown with the same logarithmic color scale. As expected, there is a difference in length scale at the same location in the mesh of about a factor of two. Finally, Figure 2(c) shows the refinement size quality Qr ; all values fall between 0.75 and 1.21 for this pair of meshes; this implies a length scale ratio between 1.5 and 2.4 everywhere in the mesh. While an ideal mesh refinement pair would have Qr ≡ 1 everywhere in the domain, the range achieved here is acceptable. As we will see, this range is consistent for isotropic mesh refinement pairs that are generated with reasonable care.

(a) Coarse mesh size

(b) Fine mesh size

(c) Refinement ratio (Equation~1)

(d) Histogram of refinement size quality Qr

Figure 3: Isotropic Mesh Refinement Pair for the NACA 0012 Airfoil As an example of more practical interest, let us consider a refinement pair for the NACA 0012 airfoil. Again, the coarse mesh (Figure 3(a)) has a length scale that is locally equal to one-fourth of the domain geometric length scale, while the fine mesh (Figure 3(b)) has a length scale half as large. As these medium views of the mesh suggest, the local length scale differs by about a factor of two between the two meshes. Figure 3(c) shows the refinement size quality Qr ; once again, this quantity is near one throughout the domain, with random excursions to values as large as about 1 1/3 and as low as about three quarters. This good behavior of the refinement ratio is confirmed by the histogram in Figure 3(d).a This is consistent with the previous example. Also consistent with the previous example is that the mesh is correctly identified as isotropic everywhere, with aspect ratio never exceeding 2.44. As a final example of an ideal isotropic refinement pair, we consider two meshes carefully generated in a clipped cube to be a refinement pair. The coarse mesh has 1206 vertices (5564 tetrahedra), while the fine mesh has 7972 vertices (41826 tetrahedra). Cutaways of these meshes are shown in Figure 4(a) and (b). The refinement size quality Qr , shown in Figure 4(c) with a histogram in Figure 4(d), again has a tightly restricted range, with no coarse cells more than about 16% too large or about 42% too small. Just as for the two-dimensional cases, the refinement size quality clusters tightly around 1. a Note that, as a side-effect of the way piecewise constant cell data is exported to Paraview and the histograms are created, each vertex of each cell is counted in the histogram, giving a count that is three times the number of cells in each bin for triangular meshes and four times the number of cells for quadrilateral meshes. The tetrahedral meshes appearing later in the paper are exported differently, with each vertex included only once in the output file, so the vertex counts in those histograms are accurate.

5 of 14 American Institute of Aeronautics and Astronautics

(a) Coarse mesh (5564 cells). Colored by cell length scale.

(b) Fine mesh (41826 cells). Colored by cell length scale.

(c) Refinement size quality Qr .

(d) Histogram of refinement size quality.

Figure 4: Isotropic mesh refinement pair for a cube with one corner clipped off.

6 of 14 American Institute of Aeronautics and Astronautics

III.B.

Near-Ideal Pairs of Anisotropic Meshes

(a) Coarse mesh

(b) Fine mesh

(c) Histogram of refinement ratio

(d) Histogram of relative aspect ratio

Figure 5: Anisotropic Mesh Refinement Pair for the NACA 0012 Airfoil We now turn our attention to an anisotropic mesh refinement pair. For this example, the meshes were generated using an elliptic mesh generator, producing O-meshes around a 0012 airfoil. In each case, the domain had a radius of five chords. the coarse mesh was 50 × 50, with a normal spacing at the wall of 10−3 ; as a consequence, the stretching ratio prescribed along the branch cut in the mesh was 1.137. The fine mesh is 100 × 100, with a normal spacing at the wall of 5 · 10−4; in this case, the stretching ratio was 1.066. These meshes are shown in Figure 5(a) and (b). As for the isotropic case, these meshes should represent as nearly ideal a scenario as it is practical to achieve without generating strictly nested meshes. Figure 5(c) shows the refinement size quality for this refinement pair. As expected, this quantity is quite close to one everywhere; only a handful of cells fall outside the range (0.75, 1.33), and these may be numerical artifacts induced by differences in the shape of the domain for the fine versus the coarse mesh. We also expect that aspect ratios will be nearly the same at the same location in the domain. As Figure 5(d) shows, this is also true. Finally, we expect that the orientation of anisotropic cells will be well aligned between the two meshes. When orientation difference is measured by the dot product of the principal directions on the two meshes, fewer than 1% of cells have a dot product smaller than 0.98. Of these cells, most are almost exactly isotropic and hence the control volumes have nearly equal eigenvalues and no dominant principal direction. The remainder appear, once again, to be boundary artifacts. III.C.

Mesh Pairs with Known Mismatches

Taken together, the examples in the previous subsections give a reasonable estimate of how well we expect the size, aspect ratio, and orientation of cells in a mesh refinement pair to correspond to each other. Specifically, it appears to be reasonable to expect that a good refinement pair will have: • refinement size quality Qr between mesh;

3 4

and

4 3

for nearly all cells, regardless of dimension or (an)isotropy of the

7 of 14 American Institute of Aeronautics and Astronautics

• ratio of aspect ratios in anisotropic mesh regions nearly 1 with the possible exception of noise near the boundary; • and difference in orientation between meshes that is comparable to the change in orientation between adjacent cells in the same mesh. With this background, we will now examine several alleged refinement pairs for which we expect, on the basis of how the meshes were generated, that the meshes are not in fact a legitimate refinement pair. For these pairs, we expect that our comparisons of cell size, aspect ratio, and orientation will allow us to identify that a problem exists.

(a) Coarse mesh

(b) Fine mesh

(c) Cells with poor refinement size quality

(d) Histogram of size ratio

Figure 6: Example of a bad isotropic “refinement pair.” For our first example, we will consider an isotropic “refinement pair”. The coarse mesh in this case is the same as the coarse mesh from the 0012 refinement pair shown earlier. The fine mesh has approximately the same number of cells as the previous fine mesh, but with coarser resolution compared with geometric length scale and slower cell size grading. These meshes are shown in figure 6(a) and (b). There is certainly nothing in the visual appearance of these meshes to suggest that they’re bad meshes, and in fact they are not by standard geometric quality measures; for instance, each has a minimum angle in excess of 30 degrees. However, as Figure 6(c) shows, the length scales for these meshes are not uniformly proportional; the cells shown are only those whose refinement size quality Qr falls outside the range 43 , 43 . The histogram in Figure 6(d) confirms that the fine mesh in this pair has many cells that are too small compared with coarse mesh cells at the same location and a few cells that are much too large; the latter are clustered at the leading and trailing edges. As a second example, we will examine an anisotropic “refinement pair”. The coarse mesh in this case is the anisotropic 50 × 50 mesh of the previous section. The fine mesh is again 100 × 100, but the wall normal spacing is the same as the coarse mesh — 0.001 chords; this is not an atypical choice in practice, given the restrictions imposed on near-wall spacing by the requirements of turbulence models. In this case, the stretching ratio is 1.057. Because of √ the way that these meshes were generated, we expect that the refinement size quality will be about 2 near the wall — because the normal spacing on the “fine” mesh is about twice as large as it should be — and that the ratio of aspect ratios will be correspondingly be about 12 near the wall. Figure 7 supports both of these notions. Not surprisingly,

8 of 14 American Institute of Aeronautics and Astronautics

(a) Refinement size quality Qr

(b) Aspect ratio comparison

Figure 7: Example of a bad anisotropic “refinement pair..”

given that these meshes were generated elliptically, these effects damp out as we move away from the boundary. The striped effect in these figures is an artifact of using injection to transfer data from mesh to mesh.

IV. Application to DPW-III Meshes In Section III, we validated our analysis approach using several mesh pairs with known characteristics. We will now apply this approach to analyze a pair of meshes from the DPW-III repository.1 The two meshes considered are the coarse and medium Raytheon wing-body meshes, labeled in the repository as DLR-F6_WB_FX2B. These particular meshes were chosen as an easy starting point among the DPW-III mesh data, as these are all-tetrahedral meshes, and very coarse for RANS simulations: the coarse mesh has 544,000 vertices (3.14 million tetrahedra), while the fine mesh has 1.08 million vertices (6.28 million tetrahedra). The surface of the airplane and a closeup of the trailing edge/body junction for each mesh is shown in Figure 8. There are several things worth noting about these meshes based on visual inspection. First, the surface meshes are isotropic nearly everywhere (the leading and trailing edges are the exceptions); correspondingly, we expect that one of the two aspect ratios computed — specifically, the ratio of largest to second-largest eigenvalues for tetrahedra in the boundary layer — will be approximately one. Second, careful examination shows that each mesh has the same number of control volumes across thickness of the blunt trailing edge. This should be reflected in the statistics for cell size and anisotropy in that region. First, let us examine the refinement size quality for this mesh pair, shown in the same two views in Figure 9(a) and (b). For much of the surface of the wing — though not the trailing edge — the coarse and fine meshes have very nearly the right ratio of cell volumes (Qr ≈ 1). However, on most of the fuselage, the coarse mesh cells are too small relative to the fine mesh cells, by as much as a factor of about 3.5 in characteristic length, or about 40 in volume. The histogram in Figure 9(c) shows clearly that, for the mesh as a whole, there are a large number of cells in the coarse mesh that are somewhat too small (the most common refinement size quality is about 0.85). Similarly to the deliberately poor isotropic “refinement pair” in Section III.C, this mesh pair also features a distribution of size quality Qr with a relatively long “tail” above the ideal value of 1 compared with the relatively sharp dropoff below 1. Because these meshes are intended for RANS simulations, the tetrahedra near the boundary are highly anisotropic in the surface-normal direction. The highest computed anistropy is about 1000:1 for the coarse mesh and about 4000:1 for the fine mesh; note that these values are based on principal moments of the control volumes, and will differ from estimates based on edge lengths. Also, it is important to bear in mind that the vast majority (some 70%) of all tetrahedra in these meshes have overall aspect ratios (based on maximum and minimum principal moments) less than 10. Given this, the ratio of aspect ratios will more or less automatically be tightly clustered around one; we are interested in the behavior of the outliers in the distribution. Figure 10(a) shows the general distribution of the ratio of aspect ratios on the surface of the mesh, while Figure 10(b) shows a close-up at the junction of the body and the wing trailing edge; note that these two plots have vastly different scales. The histogram in Figure 10(c) shows clearly that nearly all the cells in the coarse mesh have a ratio of anistropy within 20% of the ideal value of 1 compared with the same location in the fine mesh. There is a very small number of extreme outliers with ratio of anisotropy exceeding 5.0 — less than 0.6% of all cells fall outside the range of the histogram shown — but nearly all of the cells adjacent to the body appear from Figure 10(a) to fall into this category. This suggests very strongly that the same sort of numerical artifacts that arise in the two-dimensional anistropic case are likely to be the cause of some or all of the problems with this measure near the surface, especially since the wing surface and flatter parts of the body — just above the wing attachment, the flat section for tail attachment, and the belly below the wing — have sharply lower reported values of ratio of

9 of 14 American Institute of Aeronautics and Astronautics

(a) Coarse mesh

(b) Fine mesh

(c) Coarse mesh, trailing edge-body junction

(d) Fine mesh, trailing edge-body junction

Figure 8: Surface meshes for the wing body, colored by cell volume.

(a) Whole aircraft

(b) Closeup at the trailing edge-body junction

(c) Histogram of refinement size quality

Figure 9: Refinement size quality Qr for the Raytheon wing-body meshes.

10 of 14 American Institute of Aeronautics and Astronautics

anisotropy than most of the rest of the body surface.

(a) Whole aircraft

(b) Closeup at the trailing edge-body junction

(c) Histogram of aspect ratio comparison. Approximately 0.6% of cells are omitted from this histogram because of too high a ratio.

Figure 10: Ratio of aspect ratios for the Raytheon wing-body meshes. One fringe benefit of the ratio of anisotropy test is that some of the outliers it identifies are actually genuine outliers: sliver tetrahedra embedded between (subdivided) prisms in the interior of the mesh. An example of two such tetrahedra in the very near wake, just downstream of the blunt trailing edge, is shown in Figure 11. The surface in this figure is colored by cell size ratio (scale not shown), while the two slivers in the foreground are colored by ratio of anisotropy, as indicated on the scale in the figure. The red sliver here is aligned with the flow, and is the most extreme tetrahedron in the coarse mesh by this measure. The blue tetrahedron is aligned spanwise and roughly parallel to the upper and lower surfaces of the wing. The red tetrahedron, but probably not the blue one, could also have been identified by the misalignment between its minimum principal direction (normal to its flattest “plane”) and the minimum principal directions of its neighbors, which are approximately vertical in the figure. Finally, we turn our attention to cell orientation. Because the geometry for this problem is comparatively simple, and the surface mesh is almost entirely isotropic, there is only one significant direction of anisotropy to consider in the mesh. We expect that, for both meshes, this direction will be normal to the wall. Figure 12 shows a comparison of cell orientation between meshes, measured by the absolute value of the dot product of cell orientations on the coarse and fine meshes. The line of apparent poor alignment high on the side of the fuselage is an artifact, caused by normals switching from pointing inward to outward.b To ensure that the results of this test are well-behaved throughout the boundary layer, a domain cut is included in Figure 12(b). As expected, the boundary layer meshes are similarly aligned in both meshes. Also, the isotropic far field meshes show no correlation in mesh orientation, a result which is also expected. b Heuristically, unit normal vectors have their direction adjusted so that the largest component is positive rather than negative, rather than trying to determine which choice points away from the body. Had this domain been for the right half of the airplane, this line would have been at 45◦ below horizontal on the fuselage rather than 45◦ above.

11 of 14 American Institute of Aeronautics and Astronautics

Figure 11: Slivers identified based on the ratio of anisotropy

(a) Whole aircraft

(b) Wing-body junction, including a field cut.

Figure 12: Comparison of orientation for the Raytheon wing-body meshes.

12 of 14 American Institute of Aeronautics and Astronautics

V. Summary This paper has described a new approach for analyzing whether two meshes are geometrically similar and therefore appropriate for use in a mesh refinement study. The size, anisotropy, and orientation of cells are calculated based on the moments of the cells. These quantities are projected from one mesh to another to enable comparison. The refinement size quality — the ratio of actual to desired cell size ratio — has been defined (see Equation 1). Anisotropy is compared by a simple ratio; in three dimensions, there are two anisotropies (spanwise-streamwise and streamwise-normal), each of which is compared independently. Finally, the directions of anisotropy are compared using the magnitude of the dot product. In each case, the ideal value of the computed comparisons is one. We have validated the methods used to compare meshes by applying them to pairs of meshes that are known to be valid. For both two- and three-dimensional isotropic meshes, the new approach has been shown to give excellent results for cell size. For two-dimensional anisotropic meshes, the results for aspect ratio and orientation are equally good. As an additional validation step, tests were done for mesh pairs that were generated deliberately to not be geometrically similar. For the isotropic pair, this non-similarity was exhibited as a shift in the distribution of refinement size quality for the meshes. For the anisotropic pair, there were distinctive changes in the computed aspect ratio and cell size, consistent with the parameters used to generate the meshes. Finally, we have examined a pair of meshes from the DPW-III repository, the two coarsest Raytheon wing-body meshes. The results of the analysis show that the refinement size quality is not ideally distributed, with an excess in the coarse mesh both of cells that are slightly smaller and larger than expected at the expense of the ideal size. The aspect ratio and orientation results for this mesh pair both appear to be good, with the former pointing up a possible issue in the mesh comparison approach.

VI. Remaining Challenges The approach described in this paper appears to be quite promising. Given its stage of maturity, there are — not surprisingly — a number of technical issues that remain to be resolved. These are listed, roughly, in order of increasing difficulty. 1. Support for vertex-centered meshes and meshes with mixed cell types. These are simple extensions of the current cell-centered tetrahedral code. Routines are needed for computing the moments of these other commonly-used control volume types. 2. Mesh-to-mesh transfer accuracy for highly anisotropic meshes near curved boundaries. The main difficulty here is that the fine and coarse meshes cover slightly different subdomains, making the obvious transfer operators — injection, interpolation, and reconstruction and integration, for instance — all suspect for one reason or another. This problem is known to affect cell size and anisotropy for two-dimensional anisotropic meshes, and there is reasonably strong circumstantial evidence from the wing-body case that three-dimensional anisotropic meshes suffer from the same problem. 3. Mesh-to-mesh transfer efficiency. Projection of data, by whatever method, requires the ability to identify which control volume contains a given point. This search requires some sort of geometric tree, but the query times dominate total CPU time. A more efficient tree structure than the standard alternating digital tree is apparently required. 4. Parallelization. The present code is tightly memory bound in three dimensions. At an absolute minimum, each control volume must store ten moments (size, centroid location, and six second moments), each vertex must store its coordinates (not the same as the centroid, even for vertex-centered control volumes), and each cell must know its vertices. In addition, a total of nine quantities per control volume must be stored to completely characterize the size, aspect ratios, and principal directions for a control volume. This bare minimum amounts to 948 bytes per vertex (172 bytes per cell) for cell-centered tetrahedral meshes, and 198 bytes per vertex for vertex-centered tetrahedral meshes. In practice, with the addition of an averaging stencil (which can do double duty for reconstruction if and when needed), the total in the current implementation is somewhat more than twice that, without accounting for memory used to enable efficient geometric search while projecting data from one mesh to another. All told, the present implementation required only slightly less than one kilobyte per cell for the wing-body case shown here. While memory footprint could doubtless be reduced, the total for meshes in the range of 40-80 million vertices, as anticipated for DPW-IV, will still exceed 100 GB. The challenge in

13 of 14 American Institute of Aeronautics and Astronautics

parallelizing this process lies in mesh-to-mesh transfer, which will scale well only if partitioning for both meshes is done with the same geometric bounds on partition shape. 5. Interpretation of results. In the end, this new analysis technique will only prove useful if we as a community can develop, through collective experience, criteria for when a mesh pair really is geometrically similar enough to be usable. This experience can only come by analyzing many mesh pairs and correlating their performance in mesh refinement studies with their compute quality metrics. An important but non-trivial step in this process will be developing metrics that summarize the overall geometric similarity of two meshes based on the control volume by control volume measures discussed here.

Acknowledgments This work was supported by the Canadian Natural Sciences and Engineering Research Council under Grant OPG0194467.

References 1 AIAA Applied Aerodynamics Technical Committee. Third AIAA CFD Drag Prediction Workshop. http://aaac.larc.nasa.gov/ tsab/cfdlarc/aiaa-dpw/Workshop3/workshop3.html, June 2006. 2 Charles Boivin and Carl F. Ollivier-Gooch. Guaranteed-quality triangular mesh generation for domains with curved boundaries. International Journal for Numerical Methods in Engineering, 55(10):1185–1213, December 2002. 3 Patrick Knupp. Remarks on mesh quality. In Forty-sixth AIAA Aerospace Sciences Meeting, 2008. 4 Carl Ollivier-Gooch, Amir Nejat, and Christopher Michalak. On obtaining high-order finite-volume solutions to the Euler equations on unstructured meshes. In Proceedings of the Eighteenth AIAA Computational Fluid Dynamics Conference, 2007. AIAA-2007-4464. 5 Carl F. Ollivier-Gooch. GRUMMP — Generation and Refinement of Unstructured, Mixed-element Meshes in Parallel. http://tetra.mech.ubc.ca/GRUMMP, 1998–2009. 6 Carl F. Ollivier-Gooch and Charles Boivin. Guaranteed-quality simplicial mesh generation with cell size and grading control. Engineering with Computers, 17(3):269–286, 2001. 7 Shmuel Rippa. Long and thin triangles can be good for linear interpolation. SIAM Journal on Numerical Analysis, 29(1):257–270, February 1992. 8 Jonathan Shewchuk. What is a good linear element? Interpolation, conditioning, and quality measures. In Proceedings of the Eleventh International Meshing Roundtable, 2002. 9 John C. Vassberg, Edward N. Tinoco, Mori Mani, Olaf P Brodersen, Bernhard Eisfeld, Richard A Wahls, Joseph H Morrison, Tom Zickuhr, Kelly R Laflin, and Dimitri J Mavriplis. Summary of the Third AIAA CFD Drag Prediction Workshop. Forty-fifth AIAA Aerospace Sciences Meeting, 2007.

14 of 14 American Institute of Aeronautics and Astronautics

I. Introduction Among other results of the Third AIAA Drag Prediction Workshop (DPW-III), the panel concluded that issues of mesh quality and mesh refinement are among the chief outstanding difficulties in accurately computing the aerodynamic flows of interest in these workshops.9 Considering the relatively mature state of scientific computing on unstructured meshes, it is surprising how little we know about the effects of mesh quality on solution quality. We know theoretically that commonly used schemes will asymptotically converge to correct solutions to partial differential equations we wish to solve, but the details are the effects of, for example, geometric mesh element quality or differences in the size of adjacent elements are currently poorly understood. There has been some work in the area of the effect of mesh quality on interpolation error for finite element methods,7, 8 but to the best of the author’s knowledge, there is no analogous work for finite volume methods. Also, Knupp3 described some useful solution-based quality metrics, though these have seen little application in practice. The overall problem of elucidating the interaction between mesh and solution is a large and complex one. This paper seeks to address a small but important corner of that problem: the problem of determining whether a sequence of meshes truly constitutes a mesh refinement sequence for purposes of a mesh convergence study. Our goal here is, given a sequence of meshes, to determine whether those meshes are geometrically similar enough and have a uniform enough ratio of length scale from mesh to mesh to be considered plausible for use in a mesh refinement study. For purposes of this study, two meshes will be said to form an ideal refinement pair if their cell sizes are uniformly proportional across the entire domain and if cell aspect ratio and orientation is the same in anisotropic regions of the mesh. Realistically, ideal refinement pairs will not exist; we will also explore what bounds are reasonably attainable for variation from the ideal. We will not, in this context, consider classical mesh quality measures, which presumably the generators and users of the mesh will routinely evaluate using readily available tools. Also, we will not consider any solution related tasks; the test we propose here is a strictly a priori test of the meshes themselves. To determine whether a series of meshes constitutes a valid mesh refinement sequence, we seek to compare cell size, shape, and orientation between pairs of meshes. The size, shape, and orientation of a cell in an unstructured mesh can be computed in more than one way; one viable approach is to compute the moments of the cells, and solve an eigenvalue problem to find the size, aspect ratio, and principal directions of the cell (see Section II). Once these quantities are known for all cells in a mesh, comparison between meshes can be accomplished by projecting data from one mesh to another. To validate the methodology and to explore the outcomes that can be expected both for carefully constructed refinement pairs and for meshes which have been deliberately generated to have clear flaws

1 of 14 American Institute of Aeronautics and Astronautics

in mesh refinement, a number of simple examples are presented in Section III. Results evaluating whether sample published meshes used in DPW-III are in fact reasonable mesh refinement sequences are given in Section IV.

II.

Methodology

Cell size Si is by far the easiest property that we wish to compute; size formulae for triangles, quadrilaterals, tetrahedra, prisms, and √ hexahedra are well-known and need not be repeated here. We will compare cell length scale, which we define as d Si in d dimensions. II.A.

Measuring Cell Anisotropy, and Orientation

Cell anisotropy and orientation can be assessed from the normalized second moments of the cell. These moments can be written as, for example, Ixx

=

Ixz

=

1 x2 dS Si i Z 1 xz dS Si i Z

where the integral is over area in two dimensions and over volume in three dimensions. In practice, these moments can also be evaluated by applying Gauss’s Theorem and integrating around cells: Ixx

=

Ixz

=

1 x3 nx dl Si ∂ i 3 Z x2 z 1 nx dl Si ∂ i 2 Z

where nx is the x-component of the normal to the cell boundary, and dl indicates integration around the boundary4 (either a line integral in two dimensions or a surface integral in three dimensions). In two dimensions the moment of inertia tensor of the cell can be written as # " Ixx Ixy M= Ixy Iyy This form differs from the standard form of the moment of inertia tensor only in that each term has been divided by the size of the cell. We can determine the cell anisotropy either by finding the eigenvalues and eigenvectors of the moment of inertia tensor or by applying a Mohr’s circle analysis to find the principal moments and their directions. Here, we will illustrate the latter. From Mohr’s circle analysis, found in any introductory text mechanics, we know that ron solid 2 Ixx −Iyy I +I 2 . The minimum and the inertia tensor represents a circle with its center at xx 2 yy , 0 and a radius of + Ixy 2

maximum principal moments are, therefore

Imax

=

Imin

=

Ixx + Iyy +

q 2 (Ixx − Iyy ) + 4Ixy

q2 2 Ixx + Iyy − (Ixx − Iyy ) + 4Ixy 2

I −I

The principal axes are rotated by an angle of 12 arctan Ixy / xx 2 yy ratio of maximum and minimum moments: r Imax AR = Imin

. Finally, the aspect ratio is the square root of the

As expected, this definition gives AR = l/h for rectangle of length l and height h. In three dimensions, there is no easy analogy to Moore’s circle, and so instead we solve the 3 × 3 eigenvalue problem numerically. In this case, we compute two aspect ratios, based on the ratios of maximum and median eigenvalues and of the median and minimum eigenvalues. For a typical 3D viscous mesh, these correspond to the cross-stream to 2 of 14 American Institute of Aeronautics and Astronautics

streamwise aspect ratio and the streamwise to wall normal aspect ratio, respectively. Orientation comparisons in three dimensions are most easily done by retaining the unit vectors in the directions of the median and minimum eigenvalues — that is, the stream wise and wall normal directions — and determining the differences between these directions for a pair of meshes by computing dot products. II.B.

Practical Considerations

The description in the preceding section is deceptively simple. In practice, because unstructured meshes are intrinsically non-smooth at the fine level, data must be adjusted somewhat to allow a reasonable comparison. With respect to cell size, averaging is applied to construct an approximate average cell size over a region somewhat larger than an individual cell. For each cell, consider the neighborhood of that cell that would be used in a quadratic reconstruction of solution.4 In practice this results in a neighborhood with about 10 cells in two dimensions and 15 to 18 cells in three dimensions. A weighted average over the stencil is used, with weights computed based on inverse square of distance of a cell centroid from the centroid of the central cell. As Figure 1 illustrates, this technique improves smoothness of measured cell size, which will make comparison of sizes between meshes much easier to do reliably.

(a) Unsmoothed length scale

(b) Smoothed length scale

Figure 1: Effect of weighted averaging on measured cell size. When measuring anisotropy, it is important to keep in mind that even an ideal unstructured mesh will have cell aspect ratios that differ somewhat from one. For instance, in two dimensions, typical high quality mesh generators still produce triangles with √angles as small as 30 degrees. This can lead to aspect ratios as large as 3 (for a 30-30120 triangle) or even 3 + 2 3 3 (for a 30-75-75 triangle), even though such triangles appear not to cause significant difficulties for flow solvers. In three dimensions, the situation is even more complex. Experimentally, we find that isotropic meshes that have good geometric quality by other measures typically have aspect ratios as high as about 2.5. Therefore, in this paper all cells with aspect ratios less than 2.5 will be treated as isotropic — that is, as having unit aspect ratio. In three dimensions, this rule is applied independently to the two measured aspect ratios. Finally, when the mesh is nearly isotropic — that is, when the aspect ratio falls below the isotropic threshold previously imposed — orientation is clearly not an issue. Therefore, when either mesh meets this threshold locally, the orientation difference will arbitrarily be set to zero. II.C.

Comparison Between Meshes

The final technical issue to be considered here is how to compare size, aspect ratio, and orientation between meshes. When comparing cell sizes, our goal is to determine whether all cells in one mesh are proportionately smaller than all cells in the other mesh. The constant of proportionality depends on the relative number of cells in the mesh. Specifically, the cell sizes in each mesh should be proportional to the domain size divided by the number of cells. Because we are comparing cell characteristic length, for meshes covering effectively identical domains, we expect the length scale to behave as follows: r N2 L1 = 2D: L2 N1 r L1 N2 3D: = 3 L2 N1 3 of 14 American Institute of Aeronautics and Astronautics

where Ni is the number of cells in mesh i and Li is the local length scale for mesh i. To confirm this, we first project the length scale from mesh 1 to mesh 2. We do this by performing an approximate numerical integration of the length scale on mesh 1 over the cells of mesh 2. This integration is problematic for high-aspect ratio cells and curved boundaries, where the fine and coarse domains are significantly different, making the correspondence between fine and coarse mesh control volumes for anisotropic meshs especially bad near the wall. To prevent numerical problems for this situation, we evaluate the numerical integration using injection rather than by reconstructing the solution on one mesh and integrating more exactly. While the projection we use is inexact, the “more exact” approach of reconstructing and integrating has its own problems, because the mesh quantities we are transfering are discontinuous and the two meshes are not nested. A more accurate and robust approach to data transfer is definitely desirable. After the projection is completed, we compute a refinement size quality Qr , defined as: r L1 d N1 Qr ≡ (1) L2 N2 where d is the number of space dimensions; this quantity should be identically one in theory. For aspect ratio, we expect the same value in each location on pairs of geometrically similar meshes. In this case, following a projection of aspect ratio from one mesh to the other, we simply take the ratio of aspect ratios for the two meshes, again expecting the results to be exactly one for ideal meshes. Finally, for orientation, we expect that anisotropic cells will be oriented at the same angle. We can compute the angular deviation between orientations based on the dot product or cross product of these vectors.

III. Assessing Mesh Pairs In this section, we will examine the behavior of the mesh pair analysis techniques described in Section II. We will begin by examining pairs of meshes that have been carefully generated to be nearly ideal mesh refinement pairs, for both isotropic and anisotropic meshes. Then we will examine the effects of particular deliberate pathologies in refinement pairs to characterize the signature of these pathologies in the current methodology. Together, these two steps will provide a basis for interpretation of results when refinement pairs of unknown quality are examined in Section IV. III.A.

Near-Ideal Pairs of Isotropic Meshes

(a) Coarse mesh size

(b) Fine mesh size

(c) Refinement ratio (Equation~1)

Figure 2: Isotropic Mesh Refinement Pair in a Square with a Clipped Corner We begin by examining a refinement pair for a mesh in a square with a small corner (1% of the side length) cut off. These two meshes were both generated using the GRUMMP mesh generator,2, 5, 6 which uses a Delaunay insertion algorithm. GRUMMP automatically computes a geometric length scale based on the boundary curvature and proximity to one boundary to another. Also, users can specify that the final cell size be some fraction of this geometric length scale. In this case, the coarse mesh in the refinement pair was created with a mesh length scale one-fourth as large as the geometric length scale, while the finer mesh was created with a mesh length scale one-eighth as large as the geometric length scale. The course mesh has 1021 cells, while the fine mesh has 4045 cells; these meshes are shown in Figure 2. These meshes are isotropic, and as expected our analysis of aspect ratio correctly identified this: the largest measured aspect ratio in either mesh was 2.34. Also, of course, we expect that the cell length scales will be almost exactly proportional to one another. Figure 2 shows the results of length-scale analysis for this pair of

4 of 14 American Institute of Aeronautics and Astronautics

meshes. Figure 2(a) shows the length scale for the coarse mesh, and Figure 2(b), the length scale for the fine mesh; these are shown with the same logarithmic color scale. As expected, there is a difference in length scale at the same location in the mesh of about a factor of two. Finally, Figure 2(c) shows the refinement size quality Qr ; all values fall between 0.75 and 1.21 for this pair of meshes; this implies a length scale ratio between 1.5 and 2.4 everywhere in the mesh. While an ideal mesh refinement pair would have Qr ≡ 1 everywhere in the domain, the range achieved here is acceptable. As we will see, this range is consistent for isotropic mesh refinement pairs that are generated with reasonable care.

(a) Coarse mesh size

(b) Fine mesh size

(c) Refinement ratio (Equation~1)

(d) Histogram of refinement size quality Qr

Figure 3: Isotropic Mesh Refinement Pair for the NACA 0012 Airfoil As an example of more practical interest, let us consider a refinement pair for the NACA 0012 airfoil. Again, the coarse mesh (Figure 3(a)) has a length scale that is locally equal to one-fourth of the domain geometric length scale, while the fine mesh (Figure 3(b)) has a length scale half as large. As these medium views of the mesh suggest, the local length scale differs by about a factor of two between the two meshes. Figure 3(c) shows the refinement size quality Qr ; once again, this quantity is near one throughout the domain, with random excursions to values as large as about 1 1/3 and as low as about three quarters. This good behavior of the refinement ratio is confirmed by the histogram in Figure 3(d).a This is consistent with the previous example. Also consistent with the previous example is that the mesh is correctly identified as isotropic everywhere, with aspect ratio never exceeding 2.44. As a final example of an ideal isotropic refinement pair, we consider two meshes carefully generated in a clipped cube to be a refinement pair. The coarse mesh has 1206 vertices (5564 tetrahedra), while the fine mesh has 7972 vertices (41826 tetrahedra). Cutaways of these meshes are shown in Figure 4(a) and (b). The refinement size quality Qr , shown in Figure 4(c) with a histogram in Figure 4(d), again has a tightly restricted range, with no coarse cells more than about 16% too large or about 42% too small. Just as for the two-dimensional cases, the refinement size quality clusters tightly around 1. a Note that, as a side-effect of the way piecewise constant cell data is exported to Paraview and the histograms are created, each vertex of each cell is counted in the histogram, giving a count that is three times the number of cells in each bin for triangular meshes and four times the number of cells for quadrilateral meshes. The tetrahedral meshes appearing later in the paper are exported differently, with each vertex included only once in the output file, so the vertex counts in those histograms are accurate.

5 of 14 American Institute of Aeronautics and Astronautics

(a) Coarse mesh (5564 cells). Colored by cell length scale.

(b) Fine mesh (41826 cells). Colored by cell length scale.

(c) Refinement size quality Qr .

(d) Histogram of refinement size quality.

Figure 4: Isotropic mesh refinement pair for a cube with one corner clipped off.

6 of 14 American Institute of Aeronautics and Astronautics

III.B.

Near-Ideal Pairs of Anisotropic Meshes

(a) Coarse mesh

(b) Fine mesh

(c) Histogram of refinement ratio

(d) Histogram of relative aspect ratio

Figure 5: Anisotropic Mesh Refinement Pair for the NACA 0012 Airfoil We now turn our attention to an anisotropic mesh refinement pair. For this example, the meshes were generated using an elliptic mesh generator, producing O-meshes around a 0012 airfoil. In each case, the domain had a radius of five chords. the coarse mesh was 50 × 50, with a normal spacing at the wall of 10−3 ; as a consequence, the stretching ratio prescribed along the branch cut in the mesh was 1.137. The fine mesh is 100 × 100, with a normal spacing at the wall of 5 · 10−4; in this case, the stretching ratio was 1.066. These meshes are shown in Figure 5(a) and (b). As for the isotropic case, these meshes should represent as nearly ideal a scenario as it is practical to achieve without generating strictly nested meshes. Figure 5(c) shows the refinement size quality for this refinement pair. As expected, this quantity is quite close to one everywhere; only a handful of cells fall outside the range (0.75, 1.33), and these may be numerical artifacts induced by differences in the shape of the domain for the fine versus the coarse mesh. We also expect that aspect ratios will be nearly the same at the same location in the domain. As Figure 5(d) shows, this is also true. Finally, we expect that the orientation of anisotropic cells will be well aligned between the two meshes. When orientation difference is measured by the dot product of the principal directions on the two meshes, fewer than 1% of cells have a dot product smaller than 0.98. Of these cells, most are almost exactly isotropic and hence the control volumes have nearly equal eigenvalues and no dominant principal direction. The remainder appear, once again, to be boundary artifacts. III.C.

Mesh Pairs with Known Mismatches

Taken together, the examples in the previous subsections give a reasonable estimate of how well we expect the size, aspect ratio, and orientation of cells in a mesh refinement pair to correspond to each other. Specifically, it appears to be reasonable to expect that a good refinement pair will have: • refinement size quality Qr between mesh;

3 4

and

4 3

for nearly all cells, regardless of dimension or (an)isotropy of the

7 of 14 American Institute of Aeronautics and Astronautics

• ratio of aspect ratios in anisotropic mesh regions nearly 1 with the possible exception of noise near the boundary; • and difference in orientation between meshes that is comparable to the change in orientation between adjacent cells in the same mesh. With this background, we will now examine several alleged refinement pairs for which we expect, on the basis of how the meshes were generated, that the meshes are not in fact a legitimate refinement pair. For these pairs, we expect that our comparisons of cell size, aspect ratio, and orientation will allow us to identify that a problem exists.

(a) Coarse mesh

(b) Fine mesh

(c) Cells with poor refinement size quality

(d) Histogram of size ratio

Figure 6: Example of a bad isotropic “refinement pair.” For our first example, we will consider an isotropic “refinement pair”. The coarse mesh in this case is the same as the coarse mesh from the 0012 refinement pair shown earlier. The fine mesh has approximately the same number of cells as the previous fine mesh, but with coarser resolution compared with geometric length scale and slower cell size grading. These meshes are shown in figure 6(a) and (b). There is certainly nothing in the visual appearance of these meshes to suggest that they’re bad meshes, and in fact they are not by standard geometric quality measures; for instance, each has a minimum angle in excess of 30 degrees. However, as Figure 6(c) shows, the length scales for these meshes are not uniformly proportional; the cells shown are only those whose refinement size quality Qr falls outside the range 43 , 43 . The histogram in Figure 6(d) confirms that the fine mesh in this pair has many cells that are too small compared with coarse mesh cells at the same location and a few cells that are much too large; the latter are clustered at the leading and trailing edges. As a second example, we will examine an anisotropic “refinement pair”. The coarse mesh in this case is the anisotropic 50 × 50 mesh of the previous section. The fine mesh is again 100 × 100, but the wall normal spacing is the same as the coarse mesh — 0.001 chords; this is not an atypical choice in practice, given the restrictions imposed on near-wall spacing by the requirements of turbulence models. In this case, the stretching ratio is 1.057. Because of √ the way that these meshes were generated, we expect that the refinement size quality will be about 2 near the wall — because the normal spacing on the “fine” mesh is about twice as large as it should be — and that the ratio of aspect ratios will be correspondingly be about 12 near the wall. Figure 7 supports both of these notions. Not surprisingly,

8 of 14 American Institute of Aeronautics and Astronautics

(a) Refinement size quality Qr

(b) Aspect ratio comparison

Figure 7: Example of a bad anisotropic “refinement pair..”

given that these meshes were generated elliptically, these effects damp out as we move away from the boundary. The striped effect in these figures is an artifact of using injection to transfer data from mesh to mesh.

IV. Application to DPW-III Meshes In Section III, we validated our analysis approach using several mesh pairs with known characteristics. We will now apply this approach to analyze a pair of meshes from the DPW-III repository.1 The two meshes considered are the coarse and medium Raytheon wing-body meshes, labeled in the repository as DLR-F6_WB_FX2B. These particular meshes were chosen as an easy starting point among the DPW-III mesh data, as these are all-tetrahedral meshes, and very coarse for RANS simulations: the coarse mesh has 544,000 vertices (3.14 million tetrahedra), while the fine mesh has 1.08 million vertices (6.28 million tetrahedra). The surface of the airplane and a closeup of the trailing edge/body junction for each mesh is shown in Figure 8. There are several things worth noting about these meshes based on visual inspection. First, the surface meshes are isotropic nearly everywhere (the leading and trailing edges are the exceptions); correspondingly, we expect that one of the two aspect ratios computed — specifically, the ratio of largest to second-largest eigenvalues for tetrahedra in the boundary layer — will be approximately one. Second, careful examination shows that each mesh has the same number of control volumes across thickness of the blunt trailing edge. This should be reflected in the statistics for cell size and anisotropy in that region. First, let us examine the refinement size quality for this mesh pair, shown in the same two views in Figure 9(a) and (b). For much of the surface of the wing — though not the trailing edge — the coarse and fine meshes have very nearly the right ratio of cell volumes (Qr ≈ 1). However, on most of the fuselage, the coarse mesh cells are too small relative to the fine mesh cells, by as much as a factor of about 3.5 in characteristic length, or about 40 in volume. The histogram in Figure 9(c) shows clearly that, for the mesh as a whole, there are a large number of cells in the coarse mesh that are somewhat too small (the most common refinement size quality is about 0.85). Similarly to the deliberately poor isotropic “refinement pair” in Section III.C, this mesh pair also features a distribution of size quality Qr with a relatively long “tail” above the ideal value of 1 compared with the relatively sharp dropoff below 1. Because these meshes are intended for RANS simulations, the tetrahedra near the boundary are highly anisotropic in the surface-normal direction. The highest computed anistropy is about 1000:1 for the coarse mesh and about 4000:1 for the fine mesh; note that these values are based on principal moments of the control volumes, and will differ from estimates based on edge lengths. Also, it is important to bear in mind that the vast majority (some 70%) of all tetrahedra in these meshes have overall aspect ratios (based on maximum and minimum principal moments) less than 10. Given this, the ratio of aspect ratios will more or less automatically be tightly clustered around one; we are interested in the behavior of the outliers in the distribution. Figure 10(a) shows the general distribution of the ratio of aspect ratios on the surface of the mesh, while Figure 10(b) shows a close-up at the junction of the body and the wing trailing edge; note that these two plots have vastly different scales. The histogram in Figure 10(c) shows clearly that nearly all the cells in the coarse mesh have a ratio of anistropy within 20% of the ideal value of 1 compared with the same location in the fine mesh. There is a very small number of extreme outliers with ratio of anisotropy exceeding 5.0 — less than 0.6% of all cells fall outside the range of the histogram shown — but nearly all of the cells adjacent to the body appear from Figure 10(a) to fall into this category. This suggests very strongly that the same sort of numerical artifacts that arise in the two-dimensional anistropic case are likely to be the cause of some or all of the problems with this measure near the surface, especially since the wing surface and flatter parts of the body — just above the wing attachment, the flat section for tail attachment, and the belly below the wing — have sharply lower reported values of ratio of

9 of 14 American Institute of Aeronautics and Astronautics

(a) Coarse mesh

(b) Fine mesh

(c) Coarse mesh, trailing edge-body junction

(d) Fine mesh, trailing edge-body junction

Figure 8: Surface meshes for the wing body, colored by cell volume.

(a) Whole aircraft

(b) Closeup at the trailing edge-body junction

(c) Histogram of refinement size quality

Figure 9: Refinement size quality Qr for the Raytheon wing-body meshes.

10 of 14 American Institute of Aeronautics and Astronautics

anisotropy than most of the rest of the body surface.

(a) Whole aircraft

(b) Closeup at the trailing edge-body junction

(c) Histogram of aspect ratio comparison. Approximately 0.6% of cells are omitted from this histogram because of too high a ratio.

Figure 10: Ratio of aspect ratios for the Raytheon wing-body meshes. One fringe benefit of the ratio of anisotropy test is that some of the outliers it identifies are actually genuine outliers: sliver tetrahedra embedded between (subdivided) prisms in the interior of the mesh. An example of two such tetrahedra in the very near wake, just downstream of the blunt trailing edge, is shown in Figure 11. The surface in this figure is colored by cell size ratio (scale not shown), while the two slivers in the foreground are colored by ratio of anisotropy, as indicated on the scale in the figure. The red sliver here is aligned with the flow, and is the most extreme tetrahedron in the coarse mesh by this measure. The blue tetrahedron is aligned spanwise and roughly parallel to the upper and lower surfaces of the wing. The red tetrahedron, but probably not the blue one, could also have been identified by the misalignment between its minimum principal direction (normal to its flattest “plane”) and the minimum principal directions of its neighbors, which are approximately vertical in the figure. Finally, we turn our attention to cell orientation. Because the geometry for this problem is comparatively simple, and the surface mesh is almost entirely isotropic, there is only one significant direction of anisotropy to consider in the mesh. We expect that, for both meshes, this direction will be normal to the wall. Figure 12 shows a comparison of cell orientation between meshes, measured by the absolute value of the dot product of cell orientations on the coarse and fine meshes. The line of apparent poor alignment high on the side of the fuselage is an artifact, caused by normals switching from pointing inward to outward.b To ensure that the results of this test are well-behaved throughout the boundary layer, a domain cut is included in Figure 12(b). As expected, the boundary layer meshes are similarly aligned in both meshes. Also, the isotropic far field meshes show no correlation in mesh orientation, a result which is also expected. b Heuristically, unit normal vectors have their direction adjusted so that the largest component is positive rather than negative, rather than trying to determine which choice points away from the body. Had this domain been for the right half of the airplane, this line would have been at 45◦ below horizontal on the fuselage rather than 45◦ above.

11 of 14 American Institute of Aeronautics and Astronautics

Figure 11: Slivers identified based on the ratio of anisotropy

(a) Whole aircraft

(b) Wing-body junction, including a field cut.

Figure 12: Comparison of orientation for the Raytheon wing-body meshes.

12 of 14 American Institute of Aeronautics and Astronautics

V. Summary This paper has described a new approach for analyzing whether two meshes are geometrically similar and therefore appropriate for use in a mesh refinement study. The size, anisotropy, and orientation of cells are calculated based on the moments of the cells. These quantities are projected from one mesh to another to enable comparison. The refinement size quality — the ratio of actual to desired cell size ratio — has been defined (see Equation 1). Anisotropy is compared by a simple ratio; in three dimensions, there are two anisotropies (spanwise-streamwise and streamwise-normal), each of which is compared independently. Finally, the directions of anisotropy are compared using the magnitude of the dot product. In each case, the ideal value of the computed comparisons is one. We have validated the methods used to compare meshes by applying them to pairs of meshes that are known to be valid. For both two- and three-dimensional isotropic meshes, the new approach has been shown to give excellent results for cell size. For two-dimensional anisotropic meshes, the results for aspect ratio and orientation are equally good. As an additional validation step, tests were done for mesh pairs that were generated deliberately to not be geometrically similar. For the isotropic pair, this non-similarity was exhibited as a shift in the distribution of refinement size quality for the meshes. For the anisotropic pair, there were distinctive changes in the computed aspect ratio and cell size, consistent with the parameters used to generate the meshes. Finally, we have examined a pair of meshes from the DPW-III repository, the two coarsest Raytheon wing-body meshes. The results of the analysis show that the refinement size quality is not ideally distributed, with an excess in the coarse mesh both of cells that are slightly smaller and larger than expected at the expense of the ideal size. The aspect ratio and orientation results for this mesh pair both appear to be good, with the former pointing up a possible issue in the mesh comparison approach.

VI. Remaining Challenges The approach described in this paper appears to be quite promising. Given its stage of maturity, there are — not surprisingly — a number of technical issues that remain to be resolved. These are listed, roughly, in order of increasing difficulty. 1. Support for vertex-centered meshes and meshes with mixed cell types. These are simple extensions of the current cell-centered tetrahedral code. Routines are needed for computing the moments of these other commonly-used control volume types. 2. Mesh-to-mesh transfer accuracy for highly anisotropic meshes near curved boundaries. The main difficulty here is that the fine and coarse meshes cover slightly different subdomains, making the obvious transfer operators — injection, interpolation, and reconstruction and integration, for instance — all suspect for one reason or another. This problem is known to affect cell size and anisotropy for two-dimensional anisotropic meshes, and there is reasonably strong circumstantial evidence from the wing-body case that three-dimensional anisotropic meshes suffer from the same problem. 3. Mesh-to-mesh transfer efficiency. Projection of data, by whatever method, requires the ability to identify which control volume contains a given point. This search requires some sort of geometric tree, but the query times dominate total CPU time. A more efficient tree structure than the standard alternating digital tree is apparently required. 4. Parallelization. The present code is tightly memory bound in three dimensions. At an absolute minimum, each control volume must store ten moments (size, centroid location, and six second moments), each vertex must store its coordinates (not the same as the centroid, even for vertex-centered control volumes), and each cell must know its vertices. In addition, a total of nine quantities per control volume must be stored to completely characterize the size, aspect ratios, and principal directions for a control volume. This bare minimum amounts to 948 bytes per vertex (172 bytes per cell) for cell-centered tetrahedral meshes, and 198 bytes per vertex for vertex-centered tetrahedral meshes. In practice, with the addition of an averaging stencil (which can do double duty for reconstruction if and when needed), the total in the current implementation is somewhat more than twice that, without accounting for memory used to enable efficient geometric search while projecting data from one mesh to another. All told, the present implementation required only slightly less than one kilobyte per cell for the wing-body case shown here. While memory footprint could doubtless be reduced, the total for meshes in the range of 40-80 million vertices, as anticipated for DPW-IV, will still exceed 100 GB. The challenge in

13 of 14 American Institute of Aeronautics and Astronautics

parallelizing this process lies in mesh-to-mesh transfer, which will scale well only if partitioning for both meshes is done with the same geometric bounds on partition shape. 5. Interpretation of results. In the end, this new analysis technique will only prove useful if we as a community can develop, through collective experience, criteria for when a mesh pair really is geometrically similar enough to be usable. This experience can only come by analyzing many mesh pairs and correlating their performance in mesh refinement studies with their compute quality metrics. An important but non-trivial step in this process will be developing metrics that summarize the overall geometric similarity of two meshes based on the control volume by control volume measures discussed here.

Acknowledgments This work was supported by the Canadian Natural Sciences and Engineering Research Council under Grant OPG0194467.

References 1 AIAA Applied Aerodynamics Technical Committee. Third AIAA CFD Drag Prediction Workshop. http://aaac.larc.nasa.gov/ tsab/cfdlarc/aiaa-dpw/Workshop3/workshop3.html, June 2006. 2 Charles Boivin and Carl F. Ollivier-Gooch. Guaranteed-quality triangular mesh generation for domains with curved boundaries. International Journal for Numerical Methods in Engineering, 55(10):1185–1213, December 2002. 3 Patrick Knupp. Remarks on mesh quality. In Forty-sixth AIAA Aerospace Sciences Meeting, 2008. 4 Carl Ollivier-Gooch, Amir Nejat, and Christopher Michalak. On obtaining high-order finite-volume solutions to the Euler equations on unstructured meshes. In Proceedings of the Eighteenth AIAA Computational Fluid Dynamics Conference, 2007. AIAA-2007-4464. 5 Carl F. Ollivier-Gooch. GRUMMP — Generation and Refinement of Unstructured, Mixed-element Meshes in Parallel. http://tetra.mech.ubc.ca/GRUMMP, 1998–2009. 6 Carl F. Ollivier-Gooch and Charles Boivin. Guaranteed-quality simplicial mesh generation with cell size and grading control. Engineering with Computers, 17(3):269–286, 2001. 7 Shmuel Rippa. Long and thin triangles can be good for linear interpolation. SIAM Journal on Numerical Analysis, 29(1):257–270, February 1992. 8 Jonathan Shewchuk. What is a good linear element? Interpolation, conditioning, and quality measures. In Proceedings of the Eleventh International Meshing Roundtable, 2002. 9 John C. Vassberg, Edward N. Tinoco, Mori Mani, Olaf P Brodersen, Bernhard Eisfeld, Richard A Wahls, Joseph H Morrison, Tom Zickuhr, Kelly R Laflin, and Dimitri J Mavriplis. Summary of the Third AIAA CFD Drag Prediction Workshop. Forty-fifth AIAA Aerospace Sciences Meeting, 2007.

14 of 14 American Institute of Aeronautics and Astronautics