Best Practices for Unstructured Grid Shock-Fitting Peter L. McCloud ⇤ ERC-Inc, Houston, TX 77058 Unstructured grid solvers have well-known issues predicting surface heat fluxes when strong shocks are present. Various e↵orts have been made to address the underlying nu- mericalissuesthatcausetheerroneouspredictions. Thepresentworkaddressessomeofthe shortcomings of unstructured grid solvers, not by addressing the numerics, but by apply- ing structured grid best practices to unstructured grids. A methodology for robust shock detection and shock-fitting is outlined and applied to production-relevant cases. Results achieved by using the Loci-CHEM Computational Fluid Dynamics solver are provided. I. Introduction Traditionally, Computational Fluid Dynamics (CFD) solvers have used point-matched, structured grids foranalyzinghypersoniccases. Twoofthesesolvers,LAURA1,2 andDPLR3–5 havelongbeenconsideredthe standardforhypersonicapplications. Therequirementtohavepoint-matchedstructuredgridscanmakegrid generation tedious or nearly impossible when the geometry becomes complex. To address grid generation complexity, e↵orts are being made to create hypersonic CFD codes that utilize either structured overset grids6,7 or unstructured grids. Unstructured grids o↵er the promise of easier grid generation processes, but haveinherentnumericalcomplicationswhenstrongshocksarepresent. Satisfactoryinviscidfluxformulations to handle strong shocks are currently an ongoing area of research. Two unstructured grid solvers currently being developed as the next generation of hypersonic solvers are US3D8 and FUN3D.9 While US3D has shown great promise for unstructured grids, it should be noted that the preferred grid topology for US3D is to use purely hexahedrals to take advantage of the higher order methods. US3D also requires that a cell connectivity map be created. The present work sidesteps the underlying numerical issues and addresses the issues of strong shocks on unstructured grids by adopting structured grid solver best practices. This requires creating a more complex grid topology, but the flexibility of unstructured grid generation can easily incorporate the increased com- plexity. Amethodologyofperformingunstructuredshock-fittingonproduction-relevantcaseswasdeveloped, tested and verified against test data. Whiletheframeworkforperformingtheunstructuredshockadaptionwasbuilttobesolveragnostic,and could even work with structured grid solvers, the results shown were computed using the Loci-CHEM10,11 CFD solver. Loci-CHEM is an unstructured, finite-rate chemistry solver developed at Mississippi State University (MSU) by Dr. Ed Luke. II. Structured Shock-Fitting Methodology Structured grids for hypersonic problems are typically point-matched and the coordinates are labeled i, j,andk,withthekdirectionbeingtheo↵-bodydirection. Onceasolutionisobtainedontheinitialgrid,the shock-fittingprocessisstarted. TheuserspecifiesapercentageofthefreestreamMachnumber. ForDPLR, this is typically 95%. The adaption algorithm analyzes the k-lines and determines where the Mach number reachesthespecifiedvalue,whichrepresentstheshocklocation. Itshouldbenotedthattheidentifiedpoints on all of the k-lines essentially represent a Mach iso-surface, which is utilized in the unstructured process. Once the shock surface is identified, the locations are typically smoothed relative to each other. Finally, the k-lines are then modified such that the k point lies just outside of the shock, with a minimal margin for max maximum computational e�ciency. ⇤AerothermodynamicsEngineer,JETS-EG3NASAJSC,LifetimeMember 1of14 AmericanInstituteofAeronauticsandAstronautics An example of the structured shock-fitting methodology on an axisymmetric heatshield forebody at a specified angle of attack is provided in figure 1. Figure 1a is the center-line Mach iso-contour for the initial unadapted grid, while the center-line Mach iso-contour for the final shock-adapted grid can be seen in figure 1b. Together, these two figures demonstrate how the k-lines are modified such that the endpoints lie just beyond the bow shock. (a) Center-lineMachiso-contourfortheinitialunadaptedgrid. (b) Center-line Mach iso-contour for the final shock-adapted grid. Figure 1. Before and after shock-fitting results on an axisymmetric heatshield forebody at Mach 10 using structured grids. The drawback to this methodology is that it does not allow for the existence of internal shocks or the presence of plumes in the flow field. As the algorithm follows the k-lines and encounters a location where the Mach number meets the criteria, even if not the bow shock, it will treat that feature as the bow shock and the process will fail. III. Unstructured Shock-Fitting Methodology The overall goal of the shock-fitting process is to create a grid with element faces that are aligned to the shock, decreasing the numerical noise. The present work concentrated on developing a shock-fitting methodology that leveraged structured grid best practices and made the process solver agnostic wherever possible in order for it to be applicable to a wide range of tools. The shock-fitting process developed for the present work is outlined in figure 2. The process is composed of an iterative loop composed of several sub-steps with the first loop starting from the initial solution. Thefirststepoftheloopappliesashockdetectionalgorithmtobuildasurfacethatrepresentstheshock. Next, the surface is smoothed to better align the cell faces to the shock. To maximize the cell quality at the shock, the shock surface is then re-meshed. With the shock surface finalized, a prismatic layer is then extruded on both sides of the surface. The volume grid is then re-meshed using the prismatic layer as a new outer boundary. The final step of the loop is to solve for the updated solution. The loop is repeated until the solution achieves su�cient accuracy. Implementation of the unstructured shock-fitting methodology is part of the mesh tools tool suite de- veloped at the NASA Johnson Space Center (JSC) Aerosciences Branch. The tool suite is a collection of Python scripts and modules that utilize the Visualization Tool Kit (VTK).12 The VTK libraries provide e�cient data structures and algorithms that maximize the performance of the Python scripts and utilize a pipeline process where filters are applied to data sets to obtain the desired results. 2of14 AmericanInstituteofAeronauticsandAstronautics Figure 2. Shock-fitting process. III.A. Initial Solution As with structured grids, the first step in the process is to obtain an initial grid and solution. The present work utilized Altair’s Hypermesh13 to build a quad-dominant surface mesh and AFLR314,15 to generate the volume mesh. The initial volume mesh consists of a prismatic layer grown from the viscous walls and the remaining volume is filled with tetrahedra. Obtaining the initial solution can sometimes be problematic due to carbuncles forming at the bow shock location. Typically, the best results are obtained by averaging the solution over of a number of iterations to account for the shock movement due to carbuncles. III.B. Shock Detection One of the main challenges with the unstructured shock-fitting is identifying the shock location e�ciently without the connectivity information inherent in structured grids. As mentioned previously, the structured shock-fitting process is essentially identifying a Mach iso-contour. To start the unstructured shock-fitting process, the solution is read in and the Mach iso-contour is then found by simply applying the vtkContour- Filter. Foratypicalsmoothbodyhypersoniccase, theresultingMachiso-contourrepresentsthebowshock. Anexampleoftheunstructuredshockdetectionprocessperformedontheaxisymmetricheatshieldforebody is found in figure 3. The initial center-line Mach iso-contour is seen in figure 3a. The resulting shock surface using the 95% freestream Mach number criteria is shown in figure 3b. (a) Center-lineMachiso-contourfortheinitialunadapted,un- (b) Center-line Mach iso-contour with a Mach iso-surface set structuredgrid to95%ofthefreestreamMachnumber Figure 3. Shock detection results for an axisymmetric heatshield forebody at Mach 10 using unstructured grids. 3of14 AmericanInstituteofAeronauticsandAstronautics This simplified approach to identifying a hypersonic bow shock breaks down when internal shocks or plumes are present in the flow field. However, additional processing of the iso-contour can provide the isolated surface that represents the bow shock. The connectivity of the iso-contour is first analyzed using the vtkPolyDataConnectivityFilter. This allows the iso-contour to be broken into separate surfaces, each associated with the various flow phenomena (bow shock, internal shock and plume). Identification of the bowshockisachievedbyexaminingtheboundingboxforeachofthesurfacesthatmakeuptheiso-contour. For hypersonic conditions, the bow shock will encompass all of the other flow features, so the bounds of the bow shock will match the bounds of the iso-contour as a whole. This approach is particularly e�cient since the bounding box information is automatically computed when the VTK data structures are created. Figure 4 is an example of a complex case with internal flow features. This is the Orion Exploration Flight Test (EFT)-1 launch configuration. The Launch Abort Vehicle (LAV) that surrounds the Orion capsuleduringascenthaslargecavitiesandprotrusionsthatcreatenumerousinternalshocks. Thiscasewas chosen because it’s currently not possible to perform shock adaption on this geometry using either DPLR or LAURA. Figure 4a shows the center-line Mach iso-contour results for the initial solution. In figure 4b, the initial step of the shock detection process can be seen, as well as the internal flow features. The final results of the shock detection process where the internal flow features have been removed using the outlined process are apparent in figure 4c. (a) Center-lineMachiso-contour. (b) Center-line Mach iso-contour with (c) Center-line Mach iso-contour with a Mach iso-surface set to 95% of the a Mach iso-surface set to 95% of the freestream Mach number and internal freesteam Mach number and internal featurespresent. featuresremoved. Figure 4. Iso-contour results for a complex case. III.C. Shock Smoothing The overall goal of the shock-fitting process is to create a grid with element faces that are aligned to the shock, decreasing the numerical noise. The raw shock surface identified in the prior step can be quite noisy, particularly when the shock lies within a region of tetrahedrons, as seen in figure 3. This noise is eliminated by smoothing the raw shock surface so that it is aligned to the physical shock. The smoothing algorithms used for structured shock-fitting depend on the structured nature of the grids and aren’t directly applicable to unstructured surfaces. Therefore, di↵erent algorithms had to be identified for the present work. There are numerous algorithms16–18 available for smoothing unstructured grids which were studied. III.C.1. Algorithms Two of the algorithms studied are based on a di↵usion process based on Eq. (1), where X represents the polyhedral surface mesh, � is a scale factor and L(X) represents the Laplacian. @X =�L(X) (1) @t Gaussian smoothing was one of the algorithms studied that is based on this di↵usion process. The Laplacian in Eq. (1) is approximated by Eq. (2), where m is the number of neighbors, N (i) are the first 1 4of14 AmericanInstituteofAeronauticsandAstronautics ring neighbors and x represents the vertex positions. While this algorithm reduces the noise present in the mesh, the side e↵ect is that the mesh will shrink with each iteration. 1 L(X)= x x (2) j i m � j2XN1(i) The second di↵usion-based algorithm studied was Taubin16 smoothing. This algorithm consists of two successive Gaussian iterations, the first with a positive factor � and the second with a negative scale factor µ, greater in magnitude than the first scale factor (0<�< µ). Together, the two passes act as a low pass � filter, eliminating large changes in curvature while preserving the smaller changes in curvature that define the general mesh shape, minimizing the mesh shrinkage. The downside to this algorithm is that numerous iterations are often required and the e↵ects tend to plateau quickly. Meancurvatureflow17 wasanotheralgorithmstudied. Aspartofthecurvatureflowfamilyofalgorithms, mean curvature flow works by moving each vertex along the point normals by a speed equal to the surface curvature , as shown in Eq. (3). This family of algorithms perform better at eliminating low frequency noise than the Gaussian and Taubin algorithms. The downside to using curvature flow is that small time steps must be chosen to prevent the algorithm from becoming unstable. @x i =  n (3) i i @t � There are many definitions of curvature, but for the present work, the mean curvature ¯, defined as the average of the principal curvatures,  and  (see Eq. (4)), was used exclusively. This was due to mean 1 2 curvature being easily computed in the VTK framework by applying the vtkCurvatures filter. 1 ¯ = ( + ) (4) 1 2 2 The last algorithm studied was the Two-Step smoothing algorithm.18 The first step averages the face normals, while the second step adjusts the vertex positions to fit the averaged normals. This process is iterated, typically 30-50 times. III.C.2. Implementation The Gaussian, Taubin and mean curvature flow algorithms were all implemented inside of the mesh tools suite,inanexplicitmanner. Duetoitscomplexity,theTwo-Stepsmoothingalgorithmwasnotimplemented internally. Instead, since the algorithm has already been implemented as part of Meshlab,19 the mesh tools suite calls Meshlab externally when the algorithm is required. III.C.3. Testing All of the above algorithms were tested to identify a robust method of smoothing the shock surface. During eachtest,themeancurvatureoftheshocksurfacewascomputedtoevaluatethee↵ectivenessofthesmoothing process. Thesurfacecurvatureisanexcellentcriteriaforevaluatingtheshocksmoothinge↵ectivenessbecause the enthalpy change across the shock is related to the shock angle. If there is noise in the surface curvature, then the enthalpy downstream of the shock exhibits a proportional amount of noise, leading to erroneous heat flux predictions at the walls. Testingonavarietyofcasesfoundthatacombinationofalgorithmsappliedinpasseswasthemostrobust approach. No single algorithm was able to provide a su�ciently smoothed surface. Gaussian smoothing had theunwantedsidee↵ectofshrinkingthemesh,butwasgoodatremovingparticularlynoisyregions. Taubin smoothing was perhaps the most conservative algorithm, as it did not shrink the surface, but its e↵ects tended to plateau as successive iterations were applied. The mean curvature flow algorithm provided the smoothest surface with the fewest number of iterations and without any shrinkage, but became unstable easily. The Two-Step smoothing algorithm was the superior method for cleaning particularly noisy surfaces, but had the tendency to introduce noise in smoother surfaces. 5of14 AmericanInstituteofAeronauticsandAstronautics (a) Rawiso-contourforthefirstshock-fittingstep. (b) Smoothediso-contourforthefirstshock-fittingstep. (c) Rawiso-contourforthefinalshock-fittingstep. (d) Smoothediso-contourforthefinalshock-fittingstep. Figure 5. Mean curvature results for the axisymmetric heatshield forebody. III.C.4. Final Smoothing Process After the shock detection process, the filter vtkCurvatures is applied to compute the mean curvature. Based on the amount of noise present, one of two possible routines will be used. For noisy surfaces, such as when the shock lies in a field of tetrahedrons, the best combination was found to be a series of four passes, with each pass consisting of 20 initial Taubin iterations, followed by three iterations of the Two-Step smoothing algorithmandfiveGaussiansmoothingiterations. Forsmoothersurfaces,thebestcombinationwasfoundto be successive passes with 25 Taubin iterations followed by 80 mean curvature smoothing iterations and five iterations of Gaussian smoothing. Both the Taubin and Gaussian algorithms stabilize the mean curvature iterations. After each pass, the mean curvature is re-computed and the scripts evaluate whether the most recent pass su�ciently improved the result. If no significant improvements in the surface quality are identified, the resulting surface is passed onto the next step in the process. If the quality is still improving significantly, 6of14 AmericanInstituteofAeronauticsandAstronautics the process is continued up to a maximum of 30 passes. Results of the final smoothing process for the axisymmetric heatshield forebody can be seen in figure 5. For the initial shock-fitting step where the shock lies in the tetrahedral region, the surface has significant noise, apparent by large maximum and minimum values of mean curvature in figure 5a. The smoothing process is then appliedand the results in figure 5b show that thenoise in the mean curvature is significantly reduced. Also evident in the smoothed surface is a carbuncle in the stagnation region. This indicates that the smoothing process preserves lower frequency features while eliminating the higher frequency noise. In figure 5c, the raw iso-contour results are shown for the final shock-fitting step. There are still some regions of high and low mean curvature values, where the shock steps between the prismatic layers, but overall, the shock surface is much improved from the initial solution. Additionally, it can be seen that the large carbuncle has been removed during the intermediate shock-fitting steps. The final result of the shock fitting process in figure 5d is now significantly smoother, which will lead to better predictions of the surface heat fluxes. While some noise is still present in the surface, the ratio of the noise to the underlying surface curvature is much less. III.D. Re-Meshing the Shock Surface The smoothing process creates an unstructured surface that has minimal noise in the mean curvature but doesn’t take into account the element sizing distribution. The smoothing process tends to create face sizes thatleadtopoorvolumeelementquality. Togetthebestresultsfromtheunstructuredshock-fittingprocess, the smoothed shock surface must be re-meshed to achieve the highest volume element quality possible. To accomplish this, the smoothed shock surface is passed externally to Gmsh,20 an open source grid generation software. Gmsh then re-meshes the shock surface using a user-defined spacing distribution. Occasionally during the re-meshing process, the surface normals can become flipped. Therefore, after Gmsh has finished, the final shock surface grid is checked to ensure that the surface normals are oriented correctly. Before and after results of the re-meshing process can be seen in figure 6. Figure 6a is the shock surface after the smoothing process where the faces have various sizes. The results of the re-meshing can be seen in figure 6b, where the final face sizes are distributed much more evenly. (a) Beforere-meshingthesurface. (b) Afterre-meshingthesurface. Figure 6. Surface mesh face areas before and after re-meshing. 7of14 AmericanInstituteofAeronauticsandAstronautics III.E. Extruding a Prismatic Layer around the Shock Surface Creating elements at the shock surface that have faces aligned with the shock has been demonstrated by Bonfiglioli, et al.21 to improve the accuracy of the solution. For the present work, it was found that going a step further and creating layers of prismatic elements aligned with the shock provided the best results. The prismatic layer is formed by marching the re-meshed shock surface both upstream and downstream for several layers. The upstream layers are necessary to provide margin for the shock as it settles into the final location throughout successive shock-fitting loops. The number of layers required can depend on the problem, but for the present work 10 to 30 layers were typically used. The prismatic layers were grown using the Pointwise22 grid generation package. The mesh tools scripts write out the appropriate Pointwise glyph script and run Pointwise externally. After Pointwise completes theextrusionprocess,itoutputsboththeresultingvolumegridandthedownstreamsurfaceoftheprismatic shock layer. This downstream surface will be referred to as the shock interface surface for the rest of the paper. A slice though the final volume grid, including the cells around the shock created by the extrusion process, can be seen in figure 7b. III.F. Rebuilding the Volume Grid Thefinalstepoftheunstructuredshock-fittingprocessistorebuildthevolumegrid. Anewouterboundary is created using the shock interface surface. The shock interface has an open boundary on the downstream sideandanexitsurfaceiscreatedtofilltheopeningandcompletethewater-tightouterboundary. Depending onthetypeofproblem,theexitsurfacewilleithersimplyfilltheholeintheshockinterfacesurfaceorbridge the space between the viscous walls and the shock interface surface. The completed outer boundary is then merged with the original wall surfaces to complete the boundary surfaces that define the volume. A new volume grid for the region downstream of the shock interface surface is then built in amanner similar tothe initial volume grid. Lastly, this volume grid is merged with the prismatic layer around the shock to make the final updated volume grid. IV. Results The unstructured shock-fitting methodology outlined was applied to a variety of cases in a production capacity. The mesh tools scripts are not integral to any solver and can be run alongside a solver to complete the process in an automated manner. For the present work however, all of the results presented were performed with the Loci-CHEM solver. IV.A. Axisymmetric Heatshield Forebody One test case used to validate the unstructured shock-fitting process was the axisymmetric heatshield fore- body. This class of problem is solved more e�ciently using a structured grid solver, but serves as a basic test of the shock-fitting process. Validation of the process was performed by making heat flux predictions using the Loci-CHEM solver and comparing the predictions to test data obtained from a series of shock tunnel tests. The data source chosen was the Orion 126CH23 test run at the Calspan-University at Bu↵alo Research Center (CUBRC) Lens I shock tunnel. The primary purpose of the test was to define aerothermal environments for the Orion Exploration Mission-1 compression pad design, but there were a number of thin film sensors placed on the vehicle center-line, away from the compression pads, that provide excellent data for the axisymmetric heatshield forebody. All the 126CH shock tunnel tests were performed with a Mach 10 freestream condition. A number of runs were simulated using the Loci-CHEM solver, but only two representative cases will be shown. Run 9, a low Reynolds number case with purely laminar flow on the forebody and run 8, a high Reynolds number case with turbulent flow over most of the forebody. IV.A.1. Run 9 - Laminar The first case was run using laminar conditions, with perfect air and an iso-thermal wall boundary. After obtaining the initial solution, the shock-fitting process was repeated for three loops. Figure 7 shows the center-line Mach iso-contour results for the initial solution and the final shock-fit solution. For the initial 8of14 AmericanInstituteofAeronauticsandAstronautics solution in figure 7a, the shock lies in a region of tetrahedrons. The noise in the shock surface is apparent in the upper portion of the figure. Additionally, the initial grid is ine�cient due to the large number of cells upstream of the shock. The final results in figure 7b show that the shock adaption process has created a prismatic layer of cells, well-aligned with the shock, and that the shock surface is smoother. Another benefit to the shock-fitting process is that only a small portion of the cells lie outside the shock, making the calculations more e�cient. (a) Center-line Mach iso-contour for the initial unstructured (b) Center-line Mach iso-contour for the final shock-adapted grid. unstructuredgrid. Figure7. Beforeandaftershock-fittingresultsonanaxisymmetricheatshieldforebodyatMach10usingunstructured grids. The axisymmetric heatshield forebody laminar heat flux results for each shock-fitting iteration can be seen in figure 8. The results for the initial solution in figure 8a show such a high level of noise in the heat flux that the results are completely unusable. A dramatic improvement in the heat flux predictions can be seen in the first shock adaption in figure 8b, due to the addition of the prismatic layer around the shock. Subsequent shock-fitting steps in figures 8c-d show further improvement in the heat flux predictions. It can also be seen in figure 8 that the shock-fitting process outlined in the present work has some room for improvement. The predicted heat flux isn’t completely free of noise, particularly around the stagnation point. Closer inspection of the results suggests this is due to a residual amount of noise present in the shock surface, which can be seen in figure 5d. Further improvements in the shock smoothing process should improve the heat flux predictions. The predicted center-line heat flux using the shock-fitting process is compared to the run 9 test data and the center-line heat flux predicted by DPLR in figure 9. As expected, the DPLR predictions match 9of14 AmericanInstituteofAeronauticsandAstronautics (a) Initialgrid (b) Shockfit1 (c) Shockfit2 (d) Shockfit3 Figure 8. Surface heat flux results for the axisymmetric heatshield forebody at Mach 10 using unstructured grids. 10of14 AmericanInstituteofAeronauticsandAstronautics