Comparison of Two Mesh Smoothing Techniques for Unstructured Grids

In this paper two mesh smoothing techniques are applied to unstructured grids. These are based on the solution of elliptic equations, Laplace and Winslow. In the first case, the equations are solved by a barycentric averaging procedure, and a finite volume scheme is used for the second operator. The comparison of these two techniques are based on two quality measures, shape factor and minimum angle which both are local criteria. Also, a global quality smoothness criterion is introduced and used to assess global smoothing characteristics of these methods.


Introduction
Grid smoothness is a post-processing procedure designed to improve the mesh quality. This can be done with various techniques when the nodal coordinates are modified as the solution of an operator such as Laplace, Winslow smoothing [1], or the use of functionals [2]. These have been successfully used for structured meshes. However, elliptic smoothing has been challenging for unstructured grid generation due to non-conservation form of these equations, as well as the lack of an implied unique computational domain. Knupp has shown unstructured Winslow mesh smoothing on unstructured quadrilateral meshes using a locally defined computational domain [3]. Karman [4] and Sahasrabudhe [5] successfully applied unstructured smoothing by introducing local optimized computational domains called virtual control volumes. Masters [6] then extended this to stretched viscous grids. Arabi et al. [7] addressed two major problems concerning the application of elliptic smoothing methods on unstructured meshes. The first is the non-conservation aspect, and the second addresses the mapping procedure by employing an explicit mapping. However, the method is not unique because the physical space is mapped to a square which can be arbitrary for a general unstructured mesh. The solution of the Winslow equations on an unstructured mesh was achieved by introducing a 9-point finite difference stencil around each node of the mesh to correctly handle the mixed derivatives.
Following the works of Winslow [1] and Sahasrabudhe [5], in the current study local computational space along with a finite volume scheme introduced by [7] is used to manage cross derivatives using an averaging process instead of augmented cells.

Elliptic Smoothing
Elliptic smoothing consists in solving a partial differential equation, where the dependent variables are the coordinates in physical space, (x, y), in terms of independent variables in computational space (ξ, η). Two such techniques will be investigated and compared. Laplacian smoothing is the simplest and consists in solving a Laplace operator for (x, y) In this work, this will be carried out simply by using a barycentric averaging procedure of the coordinates connected to each node of the mesh.
Where n is the number of nodes around the specific node i of the mesh. This yields smooth meshes for most regular geometries, but invalid meshes can result for configurations with concave corners where large changes of curvature occur. The smoothing procedure can be improved with the Winslow equations [1] in computational space (ξ, η) given as, L(y) = g 11 y ξξ − 2g 12 y ξη + g 22 y ηη = 0 where In [9] and [10] it is mentioned that while the Winslow operator guarantees continuum global mapping, truncation errors can lead to folded meshes. In such instances, additional control is needed to adapt the mesh around the boundaries to insure the validity of the results, especially around discontinuous parts of the physical boundary. In [10], a method based on Taylor series expansion to solve this operator on equilateral triangles where all the angles are equal to π/3 has been proposed. Another method is presented in [4], based on generating a virtual control volume in the physical domain locally around each node, as a local computational space with the same number of neighboring nodes as in physical space.
Ref. [11] compared two techniques for the solution of the functional operator using finite volume and finite difference methods. The latter based on inserting a 9-point finite difference stencil presented several near sharp corners. In the current work, Winslow equations are solved using the local computational space method [4] combined with averaging method presented by [11] for evaluation of cross-derivatives.

Computational Domain
Elliptic equations enforce a smoothness condition between the spacing in the computational domain (ξ, η) and the physical (x, y) as a result of the min/max principle. The key to optimizing a mesh through the use of elliptic equations lies in using a regular computational domain. In structured meshes, the computational domain is the same topology as the physical space, and ideal spacing i.e δx = δy, thus effectively enforcing an equal spacing around each node. For an unstructured mesh, there is no such a universal computational domain that matches every unstructured topology; therefore, for each unstructured mesh a computational domain must be devised to match the topology of the physical mesh. A solution to this problem was proposed by Sahasrabudhe with the introduction of computational domains that are only defined locally [12]. For a stencil of elements surrounding a single node, an ideal mesh can be defined by equally spacing the connected points on a unit circle. These stencils, called virtual control volumes, are locally defined and can be used to drive the solution of the smoothing equations. Thus formulated, smoothing by the use of elliptic operators is a distinct boundary value problem with dirichlet boundary conditions for each node.

Finite Volume Scheme
The Winslow operators are in non-conservative form and the three coefficients, g11, g12 and g22 are functions of gradients of the dependent variables in the computational space. In this work the coefficients are frozen during a number of sub iterations and updated for the next iteration. Using a linearization procedure, the various terms of Eqns. 3 can be integrated separately over a control volume defined around each point of the mesh in computational space. The integration Figure 1. Local mapping from physical to computational space path for the application of Green's theorem is formed by joining the centroid of each triangular element to the midpoints of its sides, as shown by the dashed lines in Fig. 1. The cell edges divide each triangular element into three equal areas, and, collectively, these areas form nonoverlapping contiguous control volumes associated with a node in the mesh. The hashed region in Fig. 1 indicates a control volume with a centroid node which is the storage location of all dependent variables. Linearizing Eqns. 4, this results in the following integral form Applying the divergence theorem to the second order derivative terms, x ξξ and x ηη , gives where the components of function F (x ξ , 0). Similarly, for the cross derivative terms, for instance x ξη we have While mathematically x ξη = x ηξ for continuous functions, numerically this integration procedure yield different results depending on the order of the integration. Ref. [4] has proposed an rigorous procedure resulting in a complex scheme using equilateral cells. To simplify the algorithm, [7] has proposed to compute the term Q in two manners. Using Q 1 = (0, x η ) and Q 2 = (x ξ , 0), gives Q as an average of these two Integrating over the control volume and applying the divergence theorem for each dependent variable, for example ξ, gives The term on the RHS integral represents the net flux through the surface of the volume and, for the Winslow operator, can be evaluated as where in counter clockwise direction, the lengths of the sides of each control volume are calculated by Ref. [7] concluded that for the cross derivative terms, applying the Green's theorem only for one component on each control volume side around node (ξ i , η i ) yields a degenerated final mesh in most cases of geometries with severe boundary curvature variations. In other words, for arbitrary deformations in the (x, y) plane, the values of the calculated fluxes in (ξ i , η i ) are dominated by the values from the cross derivatives terms. Moreover, taking only one component of the cross derivative term after applying the Green's theorem, deflects the final mesh in one direction. In [4] these terms were calculated using a set of augmented cells attached to the control volumes. Therefore, the shape of virtual control volumes are different for the cross derivative terms. In [11], it was found that taking an average by splitting the cross derivatives in two components and applying the Green's theorem in both directions avoids the generation of folded cells in the physical domain. In addition, this averaging procedure maintains the symmetry of the final mesh when the deformation of the physical boundary is symmetric. Finally, the discrete operators for L(x) and L(y) are solved using an SOR-type iterative procedure.

Mesh Quality and Mesh Smoothness
Several criteria can be used to measure the quality of a mesh such as minimum angle, maximum angle, minimum edge, maximum edge, minimum Jacobian, shape factor. We have used minimum angle and shape factor as mesh quality criteria. The minimum angle criterion means that elements with small angles are considered to be of a worse quality than ones with larger angles. The shape factor criterion measures the likeness of an element to a reference equilateral triangle as where A i is the area of the triangle, and l i (i = 1, 2, 3) are the lengths of the triangle's edges. However, in order to devise a mesh mapping appropriate for arbitrary boundary shapes in the physical domain, in addition to positive Jacobian for all cells, other measurable criteria must also be considered. In contrast to the traditional definition of mesh quality, which usually considers individual criteria of each element, smoothness of a mesh has a global definition. Thus, these two distinct measures of mesh quality and mesh smoothness may be contradictory for some cases. Indeed, a smoother mesh does not necessarily imply better mesh quality as we are going to show in this study. The smoothness criterion was introduced by [11] and will be used to compare meshes resulting from the two smoothing techniques investigated in this study. The mesh smoothness is quantified for each cell as where SR i represents the smoothness ration, A i is the area of cell i and the denominator represents the maximum area of its adjacent cells. An ideal values for SR i is as close as possible to one. The Smoothness Factor (SF ) of a mesh, was defined in [11] as follows, where N e is the total number of elements in the mesh. The range of values for this factor is 0 < SF ≤ 1, and hence, the greater SF , the smoother the mesh.

Results and Discussion
The two smoothing techniques were implemented and investigated in the context of the preceding framework for their effectiveness. Numerous test cases for different geometries and grid sizes were performed. Two representative examples, a sharp spike and a slit inside a circle, will be used to illustrate the results. The effect of both smoothing techniques on an initial raw mesh, generated using a frontal unstructured grid procedure, is shown in Fig.2. These qualitative results are quantified in Figs.3 to 5 for the shape measure, minimal angle and smoothness ratio, respectively.
As can be seen from Fig.3, from shape factor viewpoint, the barycentric method gives a more satisfactory distribution. This is also the case for the minimum angle criterion (Fig.4), as the smallest angles for the barycentric method is around 40 degree, while for the Winslow method they are around 25 degree. However, from smoothness point of view, it is clear that the Winslow operator gives better results than the barycentric method. Table (1) shows that the smoothness factor for Winslow's method is higher than that resulting from the barycentric method. The comparisons were repeated for a finer grid to show the independence of the results from grid size ( Fig.7 and Table (1)).
Another test was carried out to investigate the behavior of each method starting from an invalid mesh shown in Fig.6(a). Applying the two smoothing techniques shows that both are capable to yield a valid mesh. However, the grid obtained by the Winslow is not identical to    The final test case is a slit with several sharp corners inside a circle (Fig.8). Results show better mesh quality for the barycentric method but smoothness remains better for the Winslow method.

Conclusion
Two elliptic mesh smoothing methods, barycentric and Winslow, have been compared for 2D unstructured grids. An averaging method along with the local computational space method was employed to introduce a new discretization scheme for cross derivative terms of Winslow equations. Both methods show reasonable results concerning grid smoothness and grid quality criteria. From mesh quality point of view, barycentric method almost always gives equal or better results than the Winslow equations. However, the Winslow operator always shows better results from smoothness view point. Our investigation shows that final results based on the Winslow operators is dependent on the initial mesh. Then, barycentric technique can be employed as a pre-smoother and/or initial guess for the Winslow smoothing method. Moreover, since the Winslow operators respect the positive Jacobian for all the cells in the physical domain, they can be used for many applications such as geometry optimization or fluid-structure interaction problems as well. These two smoothing methods can be extended to three dimensional cases.