



Share This Article 
High Order Mesh Generation at PointwiseIntroductionFiniteelement methods (FEM) for computational fluid dynamics (CFD) simulations are increasing in popularity for many reasons, one of which is the relative ease with which more accurate solutions can be obtained. This increased accuracy can be obtained by adding points to the mesh (termed Hrefinement) or by increasing the polynomial degree of the basis function assumed for the solution within each element (termed Prefinement). Prefinement requires elevating the polynomial degree of linear meshes. Degree elevation is not usually a problem except for elements near highly curved geometric boundaries. The high aspect ratio mesh cells near noslip boundaries can become inverted or suffer other quality problems when their linear edges are updated to curved (higher polynomial order) shapes for boundary conformance. Therefore, some type of mesh smoothing or improvement is often required before these meshes can be used in highorder (HO) FEM simulations. (Note: For the purpose of characterizing an HO mesh, its order is equal to its polynomial degree plus one. Therefore, a linear mesh has degree 1 and order 2, a quadratic mesh has degree 2 and order 3, etc. This definition is important because the terms order and degree are often used synonymously.) Curving of HO meshes is typically accomplished using deformation techniques that employ smoothing or interpolation of the boundary deformation into the mesh's interior. Linearelastic (LE) smoothing treats the mesh elements as a collection of deformable viscoelastic solids with stiffness of the elements imposed through specification of Young's Modulus and Poisson's Ratio. Interpolation of the boundary deformation onto the interior can employ simple weighted interpolation using nearest surface points, radial basis functions (RBF) to define interior point motion, mean value coordinates (MVC), or other bases. Previous mesh curving research performed at Pointwise explored the LE, RBF, and MVC methods in two dimensions. Each method was successful in most of the cases attempted. One particularly challenging case was the wavy wall illustrated in Figure 1 that was meshed with HO triangles of degree 4. A twodimensional version of the weighted condition number (WCN) smoothing approach described in this article produced the result shown in the lower right of Figure 1 (the white mesh). Unlike the other three methods, there is no edge crossing in the WCN mesh. Figure 1: Examples of interpolation schemes applied to an HO mesh adjacent to a wavy wall. Clockwise from upper left: linearelastic, mean value coordinate, weighted condition number, and radial basis function smoothing. The green and magenta lines are the linear and fourth degree element edges, respectively. Dots are additional HO nodes in each cell’s interior. The lack of curvature in the WCN mesh (lower right) is an artifact of the software used to draw the mesh. Its edges are smoothly curved. This article describes a mesh curving method in three dimensions that uses an extension of an optimizationbased mesh smoothing scheme to produce valid curvededge, highorder meshes. The mesh smoothing method uses a cost function to enforce desired element shapes of linear subelements. Once valid subelements are obtained, the smoothed node locations are transferred to the highorder mesh. Mesh CurvingThere are three steps in the process for curving meshes.
Degree Elevation of Linear MeshesThe process begins with a linear mesh (P1) generated in Pointwise. Cells in this P1 mesh include:
Degree elevation to quadratic (P2) proceeds as follows.
The inserted nodes in a P2 element are connected in a way that allows linear subelements (used directly by the WCN smoothing in the third step) to be constructed. The P2 tetrahedron, for instance, is comprised of four subtetrahedra at each corner and four additional tetrahedra in the inner region formed by the six midedge nodes. The inner tetrahedra could be formed by drawing a diagonal edge between two opposing midedge nodes. Subelements of HO pyramids are comprised on linear pyramids and linear tetrahedra. Subelements of HO prisms and HO hexahedra are comprised solely of linear prisms and linear hexahedra, respectively. Cubic (P3) elements are created in a similar fashion.
Degree elevation to P4 and higher and construction of the subelements proceeds in a similar manner. Define a Perturbation FieldThe second step of the mesh curving process is to project onto the geometry model the inserted edge and face nodes from boundary elements. The projection distances (i.e. the perturbations) on the boundaries are used to define a perturbation field for all the interior nodes. At this stage of the curving process, a simple inversedistance weighted smoothing method is used to spread the surface perturbations into the field. A nodetonode connectivity, based on the linear subelements of the mesh is used in this perturbation smoothing process. The weighting between nodes is based on the inverse distance between the nodes raised to the fourth power. There is additional weighting imposed when the perturbation vector aligns with the connecting edge. The perturbation equations are shown below. The perturbation vector p→ is the displacement from the original position. The distance between nodes is Δs. A vector between the node n and neighboring node m is given by v→_{nm}. The summation is over the total number of a node's neighbors. A relaxation parameter, Ω, is applied to the perturbation vector update. This perturbation smoothing is an optional step prior to the optimizationbased smoothing of the actual node locations. The number of perturbation smoothing passes is user defined, typically ten to hundreds of passes. OptimizationBased Mesh Improvement of a HO MeshOptimizationbased mesh smoothing seeks to optimize a cost function related to element quality. The use of element condition number as the quality measure was published by Freitag and Knupp and is used elsewhere in Pointwise for smoothing extruded meshes. The fundamental equations for weighted condition number smoothing can be found in “Smooth Extrusion of Boundary Layer Meshes” in the 2016 Q2 issue of The Connector. The smoothing process seeks to reposition nodes to drive each subelement in the HO element toward its ideal shape. The shape of the HO element is still the linear shape. Each enhanced element can be subdivided into a collection of linear subelements. The shape of those subelements is valid prior to projecting the new surface nodes to the true boundaries. These subelements and their unperturbed shapes will be used to curve edges of the HO mesh. Figure 2: Schematic illustration of perturbing the inserted nodes in a P2 HO mesh both without (left side) and with (right side) subelements. In Figure 2 the geometry of the object being meshed is depicted in green and shows a high degree of curvature with respect to the element size. The left side of the image shows the process of converting the element to P2 and elevating the surface midedge node to the true geometry without use of subelements, creating an invalid element. The right side of the figure shows the subelements used to elevate the element to P2, depicted with dashed blue lines. After elevating the surface midedge node to the geometry and applying the smoothing a valid HO element results. The subelements can then be discarded. The subelements must be monitored for validity and quality. For extremely skewed, unperturbed HO elements the subelements can become even more skewed. Therefore, the shapes are evaluated for corner inversion and for included angles in excess of a userdefined threshold (e.g. 175 degrees). The maximum included angle is the maximum of the angle in any corner of any face of the element and the angle between adjacent faces in the element. If any shape exceeds the threshold it is modified towards the ideal shape until the threshold is met. All nodes in the interior of the mesh and nodes on floating boundaries are smoothed using the WCN smoothing method. The movement of each node is defined in the following equations. During each smoothing pass a perturbation vector, p→_{n}, for each node is determined based on a blending of the sensitivity vectors, s→_{n}, of the node with respect to the cost function. The sensitivity vector is the derivative of the X, Y and Z coordinates of node n with respect to the element j cost. The average node perturbation vector is computed by averaging all the sensitivity vectors from that node for all the elements surrounding the node. The final node perturbation vector is then a combination of the average sensitivity vector and vector from the worst cost element for that node. If the minimum element cost, C_{min}, is less than or equal to zero the weight factor, F, is 1 and the minimum sensitivity vector is used. Otherwise a blending of the minimum vector and the average vector is used. The blending exponent, p, is user specified and controls the amount of influence the minimum cost has on the final perturbation vector. Typical values of userspecified parameters are as follows.
For the more difficult cases the imposed boundary perturbations can be specified in a number of fractional steps. During each fractional step the imposed perturbations are scaled by the fractional step factor, such as 0.1 for the first step of 10 total steps. After each fractional step the perturbation field is recomputed to correspond to the full step using the resulting movement from the fractional step. In addition, the perturbation smoothing can be applied prior to each fractional step. This tends to aid in convergence by distributing and smoothing the latest perturbation field prior to the WCN smoothing. ResultsSphereA sphere was meshed with layers of anisotropic cells clustered toward the sphere's surface that transition to an isotropic far field tetrahedral mesh. Meshes of three cell types (only tetrahedra, prisms and tetrahedra, and mixed elements) and four polynomial degrees (linear, quadratic, cubic, and quartic) were generated. Figure 3 shows six of those meshes which were run in the COFFE FEM solver as a proof of concept. A resulting CFD solution for those six meshes is shown in Figure 4. Figure 3: Sphere meshes colored by cell type. Meshes in top row are P1, those in the bottom row are P2. Image made in ParaView. Figure 4: Temperature profiles as computed around the sphere using the COFFE FEM solver. Image made in ParaView. Drag Prediction Workshop III, WingOnlyThe 3^{rd} AIAA Drag Prediction Workshop included a wingonly geometry on which an HO mesh was generated. Pointwise Version 18 was used to generate unstructured surface meshes consisting of a mix of linear triangles and quadrilaterals. For this test case, the number of cells in the original linear mesh was decreased as the polynomial degree increased in order to keep the number of degrees of freedom relatively constant. The volume mesh consists of layers of anisotropic cells near the surface that transition to isotropic tetrahedra in the far field. Figure 5: Closeup views of the DPW 3 wing surface mesh, tip region, leading edge. Clockwise from upper left: P1, P2, P4, and P3. Figure 6: Closeup views of the DPW 3 wing surface mesh, tip region, trailing edge. Clockwise from upper left: P1, P2, P4, and P3. ROBIN FuselageRotorcraft fuselage drag experiments and computational analyses were performed on a generic ROBIN fuselage. This case was meshed with Pointwise and then elevated to a P2 mesh. Two linear meshes and two P2 meshes were created and simulations run with COFFE. The goal was to have approximately the same number of nodes in the P1 and P2 mesh for the each coarse and fine version. The coarse P1 and P2 meshes contained approximately 2.5 million nodes. The fine P1 and P2 meshes contained approximately 7.7 million nodes. These four meshes are shown in Figure 7. An illustration of the threedimensionality of the flowfield is seen in Figure 8, a COFEE solution on the coarse P2 mesh. The effect of the HO mesh on the solution is seen in Figure 9. This latter figure emphasizes the difference in the flowfield separation region between the P1 and P2 meshes. Specifically, the character of the separation region is the same between the two P2 meshes, which stresses the superior mesh convergence properties of the higher order scheme. Figure 9: ROBIN fuselage separation region as computed in COFFE and visualized in ParaView. From top to bottom: P1 fine mesh, P2 coarse mesh, P2 fine mesh. BANC III, Landing Gear GeometryThe 3^{rd} AIAA Workshop on Benchmark Problems for Airframe Noise Computations featured an aircraft nose landing gear configuration. A coarse linear mesh was elevated to P2 and used in the PyFR flow solver to compute a preliminary solution. Figure 12: CFD solution around the BANC III landing gear configuration computed using PyFR and displayed using ParaView. Drag Prediction Workshop VI, NASA CRM WingBodyThe sixth AIAA CFD Drag Prediction workshop featured the NASA Common Research Model (CRM) in a wingbody configuration. The coarse resolution, linear, unstructured tetrahedral mesh was used as the basis for P2 (Figure 13) and P3 (Figure 14) HO meshes. Figure 13: Axial cut through a P2 mesh for the DPW6 CRM wingbody configuration near the wing tip trailing edge. Image made in FieldView. A new surface mesh consisting of mixed elements was generated using Pointwise Version 18 and served as the basis for a P4 mesh (Figure 15). Conclusions, Improvements, and Future WorkWeighted condition number smoothing has been shown to be an effective technique for blending the curvature of highorder mesh elements from the geometric boundaries onto the interior of the mesh such that the resulting meshes can be used as the basis for computations in a highorder CFD flow solver. When combined with the other two steps of the mesh curving process – degree elevation through insertion of additional nodes and perturbation of the curved mesh to the geometric boundary – a useful tool for HO mesh generation results. Since publication of the original AIAA paper on which this article is based, two main improvements have been made to the technique.
Learn More About Pointwise Version 18If you'd like to learn more about Pointwise Version 18 and our plans its continual improvement, the Pointwise User Group Meeting 2016 is the place to be. On 21 September you will have a full day of handson training with V18. On 22 September you will hear presentations from other Pointwise users on their integration of Pointwise into their CFD process and from Pointwise on our product plans. AcknowledgmentThis research was originally published as “HighOrder Mesh Curving Using WCN Optimization” by Steve L. Karman, J. Taylor Erwin, Ryan S. Glasby, and Douglas L. Stefanski, AIAA20163178. All references to prior work cited herein may be found in the AIAA paper. 
