A Feature-Based Grid Adaptation Method with GradientBased Smoothing Weiyang Lin 1 and Steve L. Karman 2 University of Tennessee at Chattanooga, Chattanooga, TN 37403 An investigation on a feature-based grid adaptation method with gradient-based smoothing is presented. The method uses sub-division and deletion to refine and coarsen mesh points according to the statistics of gradients. Then the optimization-based smoothing is used to obtain a high-quality mesh. The adaptive optimization is available as an option to the smoothing so it is capable of distributing the grid points with respect to the gradients.

T

I. Introduction to the Methods

HE primary interest of this report is the feature-based adaptation methods, as opposed to the so-called errorbased or adjoint-based adaptations. The overall algorithm may be described as an iterative procedure of adding and deleting points so that the places of high-gradients may be sufficiently resolved. A statistic technique is used to decide if an edge is to be refined or coarsened:  e  p AF =∇φ ⋅  e (1) e In Eq. (1), AF is the statistic variable which has dependencies on the gradient vector ∇φ to adapt to, and the edge  vector e . The following edge-based criteria is used AFR AFavg + CRσ = (2)

AFC AFavg − CC σ =

(3)

where σ is the standard deviation (STD) of AF ne

σ=

∑ ( AF − AF ) i =1

i

avg

2

(4)

nn and CR and CC are user-defined coefficients. With a sufficiently large CR, no mesh points will be added, and with a sufficiently small CR, the mesh would be uniformly refined. Similar situation happens in the choices of CC. A triangular sub-division is processed after the edges being marked as either case. An example is shown in Fig. 1. In this circle-case, the triangles are sub-divided according to where the high-gradients are. The refinement coefficient is chosen to be CR = 1.0. The deletion has to be accompanied with a re-connection of mesh points. For simplicity the Trimesh algorithm is used along after each deletion. The Trimesh algorithm can also be employed to optimize the connectivity of the mesh after each refinement pass. A smoothing procedure is usually called after each adaptation pass to improve the mesh quality. In this report, Figure 1. An adapted grid of the circle case the optimization-based smoothing is used. The traditional with sub-division only smoothing, optimization-based adaptation (OBA) and another 1 2

Ph.D. Student, Graduate School of Computational Engineering, 701 East M.L. King Boulevard, Student Member. Professor, Graduate School of Computational Engineering, 701 East M.L. King Boulevard, Associate Fellow. 1 American Institute of Aeronautics and Astronautics

point-movement method are described in the following sections. A. Isotropic Optimization-Based Smoothing The basic algorithm of optimization-based smoothing uses a unique weight matrix for all triangles to adapt to. The weight matrix is usually defined as an equilateral triangle, and the condition number of each triangle in the mesh may be defined by AW −1 WA−1 (5) CN = 3 The cost function may be defined by 1 − J , if J ≤ 0.0  cost =  (6) 1 1 − CN ,if J > 0.0 The sensitivity derivatives of the cost function can thus be calculated for each triangle, and each node can be moved towards the opposite direction of the weighted gradient ne  ( sn )e ⋅ cost e ∑  pn = − e =1 ne (7) cost ∑ e e =1

B. Optimization-Based Adaptation The OBA method uses local gradients to adjust the edge lengths to get local weight matrices, thus the OBA  e method tries to smooth the mesh to be anisotropic if the gradient field is  not zero [1]; otherwise the OBA method is identical to the basic method e as is described in section A. After the reduction of edge lengths, a user defined coefficient hmin is used in an iterative process to make sure that e* each reference triangle (weight matrix) is valid.  Each edge length is reduced from 1.0 to be less-than or equal-to 1.0, e⊥* depending on the direction and the magnitude of the gradient vector, as well as a coefficient Ac (may be calculated based on hmin, see Appendix   B). The edge vector e and the gradient vector ∇f is used to define the Figure 2. The depiction of the  algorithm to get the desired edge target edge vector e * .  vector e *       * ∇f ⋅ e  e 1  ∇f e = (8)  +    2   1 − e ∇f ⋅ e  1 + Ac ∇f ⋅ e  ∇f   In Eq. (8), the normalized dot-product is used to define the direction and maximum magnitude of the vector used to   get the resulting edge vector. If, for example, the unnormalized dot-product ∇f ⋅ e is sufficiently large so that

(

1   1 + Ac ∇f ⋅ e

(

)

2

)

→0

(9)

   then the resulting edge vector e⊥* is perpendicular to the gradient vector. Usually the angle between e and e *   (solid-line in Fig. 2) is less than the angle between e and e⊥* (dashed-line in Fig. 2). If the gradient is perpendicular   to the edge vector (i.e. ∇f ⋅ e =0 ), then no reduction will be done. To ensure valid reference triangles, a coefficient hmin is used so that the sum of each two edges shall be no less than the other edge plus hmin: d1 + d 2 ≥ d 0 + hmin (10) If Eq. 10 is violated, edge d1 and d2 will be adjusted to d +h d1* = 0 min d1 (11) d1 + d 2 2 American Institute of Aeronautics and Astronautics

d 2* =

d 0 + hmin d2 d1 + d 2

(12)

C. On-the-Fly Instead of using weight matrices to optimize the entire mesh, one could move each mesh points towards the direction of gradient directly (On-the-Fly, OTF). Inspired by the algorithm described in Section B, the vector used to reduce the normalized edge vector in Eq. (8) may be used to define the direction for point movement       ∇f ⋅ e  1  ∇f = s (13)    2   q 1 − ∇f ⋅ e  1 + Ac ∇f ⋅ e  ∇f   Similar to the p used in Eq. (1) and hmin in Eq. (10), q is used in Eq. (13) to control the desired strength of clustering and convergence rate (typically 0 ≤ q < 1.25). Again, no movement will be contributed if the gradient is perpendicular to the edge vector.

(

)

II. Results of Adaptation In this section, the adaptation methods are applied to two cases. The first case is the circle case, where adding or deleting mesh points is disabled for the purpose of testing OBA and OTF. The second case is a case coupled with a CFD solver with 2nd-order accuracy in space.

Figure 3. The original grid of the circle case

Figure 4. An adapted grid of the circle case with OBA only

Figure 5. An adapted grid of the circle case with OBA and Trimesh

Figure 6. An adapted grid of the circle case with OTF and Trimesh

A. Circle Case without Sub-Division or Deletion 3 American Institute of Aeronautics and Astronautics

The results of the circle case with fixed number of points are shown in Figs. 3-6. Starting from the original unmodified mesh (Fig. 3), three different ways to distribute the mesh points are conducted for 5 adaptation passes. The way of using OBA only results in a mesh distribution (Fig. 4) not as satisfactory as the one with additional connectivity optimization by Trimesh (Fig. 5); that is, the clustering of points does not align with the high-gradients if Trimesh is not used along. The result obtained by moving mesh points on-the-fly shows the alignment of clustering with favor to the highgradients (Fig. 6). However, at the same time, the triangles near the clustering are stretched in such a way that mesh quality is not guaranteed. B. Fishtail Case (NACA0012, Mach = 0.95, α = 0°) The case with Mach number 0.95 and 0 angle-of-attack on NACA0012 typically has a flow-field presenting a shape like fishtail, which consists of two oblique shocks and a normal shock near the trailing-edge of the airfoil. The reported location of the aft shock at y = 0 from 3.32 to 3.35 chord-length away from the trailing edge [2]. To demonstrate the importance of the exponent p in Eq. (1), the adaptation is first done with p = 1.0. The results after 8 adaptation refinement passes are shown in Figs. 7-9, where the merit variables are chosen to be pressure, velocity magnitude and Mach number, respectively. The statistic coefficients are CR = 1.0 and CC = 1.0e4, and the isotropic smoothing is conducted after each adaptation pass. It is evident that the grids are refined at the discontinuities.

Figure 7. An adapted grid of the airfoil case with p = 1.0 and isotropic smoothing (Pressure)

Figure 8. An adapted grid of the airfoil case with p = 1.0 and isotropic smoothing (Velocity)

Figure 9. An adapted grid of the airfoil case with p = 1.0 and isotropic smoothing (Mach)

Figure 10. An adapted grid of the airfoil case with p = 2.0 and isotropic smoothing (Pressure)

4 American Institute of Aeronautics and Astronautics

Figure 11. An adapted grid of the airfoil case with p = 2.0 and adaptive smoothing (Pressure)

Figure 12. Flow-field of the case with p = 2.0 and adaptive smoothing (Pressure)

Figure 13. Close-up view of the case with p = 2.0 and isotropic smoothing (Pressure)

Figure 14. Close-up view of the case with p = 2.0 and adaptive smoothing (Pressure)

Figure 15. An adapted grid of the airfoil case with p = 2.0 and isotropic smoothing after 4 passes (Pressure)

Figure 16. An adapted grid of the airfoil case with p = 2.0 and adaptive smoothing after 4 passes (Pressure)

5 American Institute of Aeronautics and Astronautics

Figure 17. An adapted grid of the case with p = 2.0 and OTF (Pressure)

Figure 18. Close-up view of the case with p = 2.0 and OTF (Pressure)

The adaptation is then performed with p = 2.0 with isotropic smoothing, and the merit variable is chosen to be pressure (Fig. 10). In comparison with Fig. 7, the places other than the shocks (especially in the upstream of the oblique shocks) are resolved better. The process is repeated with p = 2.0 with adaptive smoothing, where Trimesh is employed afterwards (Figs. 1112). Although the overall mesh quality might be worse than the case shown in Fig. 10, the mesh points lines-up better with the continuities. A close look reveals the difference brought by the OBA technique. Figure 13 shows the close-up view of the adapted mesh with basic smoothing, whereas Fig. 14 shows the one with adaptive smoothing. Another way to demonstrate the effect of the adaptive smoothing is to examine the meshes obtained after 4 iterations (Figs 15-16). It is obvious from the figures that the adaptive smoothing (Fig. 16) helps the clustering of mesh points to the high-gradients. The adaptation is also conducted with p = 2.0 and the OTF method with q = 1.0 (Figs. 17-18). Similar to the case shown in Fig. 6, the mesh points move towards the high-gradients without showing respect to the mesh quality. The recorded x-locations at y = 0 for the aft shock are shown in Table 1. The values of x-locations clearly demonstrate the importance of the exponent p used in adaptation processes. Without sufficient Figure 19. Drag coefficient Cd plotted in comparison with resolution in other appropriate locations reference values than the discontinuities the position of the aft shock may not be captured very well. The numbers of nodes after the adaptation processes are also documented in the third column of Table 1. Table 1. The x-locations at y = 0 for the aft shock and the numbers of nodes Case x-location Number of nodes p = 1.0, isotropic smoothing 4.08 56127 p = 2.0, isotropic smoothing 4.38 61556 p = 2.0, adaptive smoothing 4.39 58039 p = 2.0, OTF 4.38 51347 The drag coefficient Cd is plotted for each case in Fig. 19 in comparison with the value obtained via uniform refinement, as well as a reference value [3]. Note that the value obtained by the flow-solver used in this report 6 American Institute of Aeronautics and Astronautics

through uniform refinement is different than the one reported. The OBA and the OTF methods exhibit similar convergence behavior, but the reason is not clear.

III. Discussion and Future Work In this report, the capability of using the OBA method to adaptively optimize a given mesh is presented. The basic smoothing method uses a unique weight matrix, resulting in more isotropic mesh but may compromise the adapted results. The OBA method smoothes a mesh with respect to to the gradients, leading to a high-quality adapted mesh. In comparison to the OTF method introduced in this report, the OBA is slower. However, the advantage of the OBA method is that, while certain mesh points move towards the high-gradients, the rest of the points move in a way such that the mesh can be optimized. Another advantage is that, the connectivity is remained if the Trimesh is not enabled, though this would introduce other aforementioned problems. There is more work to be done on the OTF method. The OBA method can be applied to time-dependent problems, in which cases the geometric conservation law [4] shall be used to predict the locations of mesh points by considering the grid speeds. More potential algorithms can be considered in the future. Adaptive smoothing can also be achieved by some other techniques such as the Winslow smoothing with virtual control volume (VCV), it may be coupled with the error-based adaptation as well.

Appendix A. User Manual The user-manual of the code developed with this report is shown in Table 2. Table 2. User Manual usage: refine -i mesh -n niter -i input mesh file * -n number of adaptation passes * -o output (options: vtk/gnu/mesh) -f function number (default 0 = flowsolver) advanced options: -e exponent (default 1) -r refinement coefficient (default 1.0) -c coarsen coefficient (default 1.0e4) -d deletion ratio of the maximum (default 0.2) -m modified matrices (default 0, 1 = enable) -h hmin (default 0.1) set to negative values to use the smallest inscribed radius -v variable index to adapt to (default 4) 0 = density 1 = x-momentum 2 = y-momentum 3 = total energy 4 = Mach number 5 = velocity magnitude 6 = pressure -a advanced tri = use Trimesh opt = use optimization fly = use on-the-fly instead of optimization-based-adaptation tips: use big numbers of r/c to disable refine/coarsen The easiest way to run refinement with a flow-solver is to execute the command “refine -i mesh -n iter”. For example, “refine -i naca0012.mesh -n 8” will run the adaptation code on the naca0012.mesh with 8 adaptation 7 American Institute of Aeronautics and Astronautics

passes. The default values are described in Table 1, and they can be changed by adding parameters. For example, to run the case shown in Fig. 11, use “refine -i naca0012.mesh -n 8 -e 2 -r 1.0 -c 1.0e4 -m 1 -h 0.01 -a tri -a opt -v 6”. B. Calculation of Ac The coefficient Ac in Eq. (8) may be calculated automatically by Eq. (14) below 1 h −1 Ac =  min 2 ∇f ⋅ e

)

(

so that

  ∀∇f ⋅ e :

1   1 + Ac ∇f ⋅ e

(

)

2

≥

(14)

max

1   1 + Ac ∇f ⋅ e

(

)

2

≥ hmin

(15)

max

That is, the smallest edge length hmin is not violated for any set of gradients and edges.

References 1. 2. 3. 4.

Karman, S. L. "Adaptive Optimization-Based Smoothing for Tetrahedral Meshes," AIAA SciTech 2015 conference. Vol. 14, 2015. Warren, G. P., Anderson, W. K., Thomas, J. L., and Krist, S. L. "Grid convergence for adaptive methods," AIAA paper Vol. 1592, 1991, p. 1991. Venditti, D. A., and Darmofal, D. L. "Grid adaptation for functional outputs: application to two-dimensional inviscid flows," Journal of Computational Physics Vol. 176, No. 1, 2002, pp. 40-69. Thomas, P., and Lombard, C. "Geometric conservation law and its application to flow computations on moving grids," AIAA journal Vol. 17, No. 10, 1979, pp. 1030-1037.

8 American Institute of Aeronautics and Astronautics