High Order Mesh Generation at Pointwise

---

Introduction

Finite-element 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 H-refinement) or by increasing the polynomial degree of the basis function assumed for the solution within each element (termed P-refinement). P-refinement 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 no-slip 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 high-order (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. Linear-elastic (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 two-dimensional 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.

This article describes a mesh curving method in three dimensions that uses an extension of an optimization-based mesh smoothing scheme to produce valid curved-edge, high-order meshes. The mesh smoothing method uses a cost function to enforce desired element shapes of linear sub-elements. Once valid sub-elements are obtained, the smoothed node locations are transferred to the high-order mesh.

Mesh Curving

There are three steps in the process for curving meshes.

1. The first step is to elevate linear elements of an input mesh to the prescribed polynomial degree such as quadratic (polynomial degree 2 or P2), cubic (polynomial degree 3 or P3), and quartic (polynomial degree 4 or P4) by inserting new nodes on element edges and interiors.
2. The second step is to perturb the inserted edge nodes on boundary elements to the shape of the geometry model and define a perturbation vector field for all interior nodes.
3. The third step is to smooth the HO mesh to improve the mesh quality.

Degree Elevation of Linear Meshes

The process begins with a linear mesh (P1) generated in Pointwise. Cells in this P1 mesh include:

• tetrahedra with 4 nodes,
• pyramids with 5 nodes,
• prisms with 6 nodes, and
• hexahedra with 8 nodes.

Degree elevation to quadratic (P2) proceeds as follows.

• Insert a mid-edge node on each edge.
• Insert a mid-face node on each quadrilateral face.
• Insert an interior node in each hexahedral element.

The inserted nodes in a P2 element are connected in a way that allows linear sub-elements (used directly by the WCN smoothing in the third step) to be constructed. The P2 tetrahedron, for instance, is comprised of four sub-tetrahedra at each corner and four additional tetrahedra in the inner region formed by the six mid-edge nodes. The inner tetrahedra could be formed by drawing a diagonal edge between two opposing mid-edge nodes. Sub-elements of HO pyramids are comprised on linear pyramids and linear tetrahedra. Sub-elements 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.

• Insert two mid-edge nodes on each edge.
• Insert a mid-face node on each triangular face.
• Insert four mid-face nodes on each quadrilateral face.
• Insert one interior node in each pyramid element.
• Insert two interior nodes in each prism element.
• Insert eight interior nodes in each hexahedral element.

Degree elevation to P4 and higher and construction of the sub-elements proceeds in a similar manner.

Define a Perturbation Field

The 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 inverse-distance weighted smoothing method is used to spread the surface perturbations into the field. A node-to-node connectivity, based on the linear sub-elements 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 optimization-based smoothing of the actual node locations. The number of perturbation smoothing passes is user defined, typically ten to hundreds of passes.

Optimization-Based Mesh Improvement of a HO Mesh

Optimization-based 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 sub-element 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 sub-elements. The shape of those sub-elements is valid prior to projecting the new surface nodes to the true boundaries. These sub-elements and their unperturbed shapes will be used to curve edges of the HO mesh.

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 mid-edge node to the true geometry without use of sub-elements, creating an invalid element. The right side of the figure shows the sub-elements used to elevate the element to P2, depicted with dashed blue lines. After elevating the surface mid-edge node to the geometry and applying the smoothing a valid HO element results. The sub-elements can then be discarded.

The sub-elements must be monitored for validity and quality. For extremely skewed, unperturbed HO elements the sub-elements can become even more skewed. Therefore, the shapes are evaluated for corner inversion and for included angles in excess of a user-defined 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 user-specified parameters are as follows.

• Relaxation parameter Ω are 0.1 to 0.5.
• Maximum number of iterations is typically 500 –1000.
• Cost threshold ranges from 0.9 to 0.999.

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.

Results

Sphere

A 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.

Drag Prediction Workshop III, Wing-Only

The 3^rd AIAA Drag Prediction Workshop included a wing-only 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.

ROBIN Fuselage

Rotorcraft 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 three-dimensionality 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 flow-field 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.

BANC III, Landing Gear Geometry

The 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.

Drag Prediction Workshop VI, NASA CRM Wing-Body

The sixth AIAA CFD Drag Prediction workshop featured the NASA Common Research Model (CRM) in a wing-body configuration. The coarse resolution, linear, unstructured tetrahedral mesh was used as the basis for P2 (Figure 13) and P3 (Figure 14) HO meshes.

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 Work

Weighted condition number smoothing has been shown to be an effective technique for blending the curvature of high-order 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 high-order 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.

1. Geometry models used in the original research were faceted. The method has now be extended to use NURBS geometry.
2. The inverse distance weighted smoothing method used to propagate the boundary perturbations into the interior has been replaced with a simple method that directly transfers the boundary perturbation to the nearest interior neighboring nodes.

Learn More About High Order Mesh Generation Using Pointwise

If you would like to Improve CFD solution accuracy by generating boundary conforming, high-order meshes by using Pointwise, please watch, “High-Order Mesh Generation Using Pointwise.”

Acknowledgment

This research was originally published as “High-Order Mesh Curving Using WCN Optimization” by Steve L. Karman, J. Taylor Erwin, Ryan S. Glasby, and Douglas L. Stefanski, AIAA-2016-3178.

All references to prior work cited herein may be found in the AIAA paper.