COMPARISON OF VARIOUS SPRING ANALOGY RELATED MESH DEFORMATION TECHNIQUES IN TWO-DIMENSIONAL AIRFOIL DESIGN OPTIMIZATION

During the last few decades, CFD (Computational Fluid Dynamics) has developed greatly and has become a more reliable tool for the conceptual phase of aircraft design. This tool is generally combined with an optimization algorithm. In the optimization phase, the need for regenerating the computational mesh might become cumbersome, especially when the number of design parameters is high. For this reason, several mesh generation and deformation techniques have been developed in the past decades. One of the most widely used techniques is the Spring Analogy. There are numerous spring analogy related techniques reported in the literature: linear spring analogy, torsional spring analogy, semitorsional spring analogy, and ball vertex spring analogy. This paper gives the explanation of linear spring analogy method and angle inclusion in the spring analogy method. In the latter case, two di ̈erent solution methods are proposed. The best feasible method will later be used for two-dimensional (2D) Airfoil Design Optimization with objective function being to minimize sectional drag for a required lift coe©cient at di ̈erent speeds. Design variables used in the optimization include camber and thickness distribution of the airfoil. SU2 CFD is chosen as the §ow solver during the optimization procedure. The optimization is done by using Phoenix ModelCenter Optimization Tool.


INTRODUCTION
The presence of dynamic mesh concept has greatly evolved in the past decades.The idea to deform the mesh has been greatly developed in order to perform some computations easier, especially the ones dealing with aeroelasticity computations, airfoil oscillations, or even shape optimizations.In all those purposes, the mesh has to be updated according to the updated boundary domain.Lin et al. [1] categorized some methods that can be applied to update the mesh: remeshing, mesh deformation, and combination of both remeshing and mesh deformation.Compared to mesh deformation, the concept of remeshing is more expensive, since new mesh needs to be updated for each updated boundary domain case, especially for a complex boundary problem.On the other hand, the mesh deformation technique, the computational mesh, is updated to the new displaced boundary without changing the mesh connectivity [2].
Mesh deformation technique reported in the literature, in general, can be clas-si¦ed into several methods such as: partial di¨erential equation (PDE) methods, physical analogies methods, and algebraic methods [2].In the PDE methods, Laplacian or bi-harmonic operators are generally used.The physical analogies methods [3] consider each edge of the mesh to behave as a spring which has its own sti¨ness value.On the other hand, the algebraic methods compute the movement of grid nodes as a function of boundary nodes which has no actually attached physical meaning.Sample of algebraic methods that have been reported in the literature include: Inverse Distance Weighting Interpolation [4], Delaunay Mapping [5], and Radial Basis Function Interpolation [6].Among these mesh deformation techniques, the physical spring analogy method is the most commonly used.
The mesh deformation technique with linear spring analogy has evolved since ¦rstly introduced by Batina [3] for the 2D unstructured mesh.Some improvement in the methods have been greatly developed such as: the torsional spring [7], semitorsional spring [8], and ball-vertex spring analogy [9].The same authors have also conducted another research about improvement in spring analogy technique by using ball-center spring analogy [10].
In this paper, a detailed comparison between a spring analogy method and a modi¦ed spring analogy method by taking angle into account is explained.Furthermore, in the case where edge angle is considered in the formulation, two di¨erent solution methods are presented.These two methods are later compared in terms of their accuracy and e¨ectiveness.
From the developed methods used in various spring analogy methods considered, the best viable method is used for an airfoil shape optimization.The airfoil shape optimization is accompanied by several basic design parameterizations comprised of camber, thickness, and their combination.In the optimization procedure, the objective function is to minimize the sectional drag of an airfoil by specifying the required airfoil£s lift coe©cient.The optimization is conducted by the aid of Gradient Optimizer of Phoenix Model Center.

SPRING ANALOGY MESH DEFORMATION
Based on the spring analogy mesh categorization by Blom [11], there are two major categorizations: vertex spring method and segment spring method.In the vertex method, the nodal position is considered, while in the segment method, the nodal displacement is considered.The ¦rst spring analogy developed by Batina [3] is considered as the segment method.The improvement in the traditional spring analogy considered in here is to include the angle made by each spring edge.This enhancement is similar to the approach that was also suggested by Burg [12].However, a di¨erent solution technique is explained here.

Basic Spring Analogy Method
In the spring analogy segment method, the idea is to consider each edge of the grid as a spring, which has a spring sti¨ness.Mathematically, the force exerted on this spring is computed as The sti¨ness of the spring is assumed to be inversely proportional to the length of edge ij.The following equation shows the mathematical formulation of the sti¨ness of the spring: .
The updated coordinates of the mesh nodes are computed by imposing force equilibrium for both directions on each node.As a result, the force exerted on node i by neighbor nodes j can be computed as The equilibrium can be achieved by equating the right-hand sides of Eqs.(1) to zero.Based on that, the updated nodal coordinates can be found iteratively as

Angle Inclusion in Spring Analogy Method
The previous mentioned technique lacks control for overlapped nodes since the method does not include the angle into the formulation.In this section, a detailed explanation regarding the inclusion of angle in the spring analogy formulation is brie §y explained.By taking the angle into consideration, interaction between abscissa and ordinate of the node will be presented.This interaction is modeled in the same manner like formulation used for Finite Element Analysis [13] for a beam element.In the formulation, each edge is modeled as a beam element which has its own sti¨ness matrix.This sti¨ness matrix can be computed as a function of angle made by the edge and is shown below: The solution for this improvement method can be divided into two di¨erent categories: direct and indirect methods.In the direct method, the above formulation is solved iteratively for each node located on the mesh domain.On the other hand, the solution by indirect method is solved by constructing a big square matrix corresponding to all unknown displacements in the grid.

Direct method
In the direct method, the idea is to solve the unknown displacements, namely, -x and -y for each nodes iteratively.This can be achieved by summing all forces exerted on node i from all the neighbor nodes j based on Eq. ( 2).Mathemat-ically, the force exerted on node i by only one neighbor node j can be computed as In order to ¦nd the resultant forces exerted on node i, all the spring forces coming from neighbor nodes j are summed up.The total forces exerted on node i are shown below: The unknown displacements -x i and -y i are solved iteratively by imposing Dirichlet£s boundary condition on the moving nodes.These unknowns are solved by equating Eqs.(3) to zero.The ¦nal expression for the solution method used is shown below: This equation can be solved by using Cramer£s rule.The solution is said to be converged if there exists no more variation in the nodal displacements.

Indirect method
Unlike the previous method on which the solution is achieved by solving the displacement for each node iteratively, the indirect method proposes a method by solving the unknown displacement for whole nodes iteratively.This idea is very similar to the solution of Finite Element Analysis [13].
As mentioned earlier, each edge of the mesh has its own local sti¨ness matrix shown in Eq. ( 2).These local matrices are later assembled into a global sti¨ness matrix based on the global node numbering used in the solution.The node numbering is made in such a way that the nodes corresponding to the moving boundary are numbered ¦rst and the nodes corresponding to the active nodes (movable nodes) are numbered later.Based on this numbering, the global sti¨ness matrix is partitioned.This global sti¨ness matrix, K, can be written and partitioned as where Q b de¦ne the boundary displacements and Q a de¦ne the active displacements.K bb and K aa are square partitioned matrices based on node numbering de¦ned earlier.The solutions for active displacements are computed based on the last equation of Eqs.(4).In the observed case, the solution is computed based on the Conjugate-Gradient Method.

MESH DEFORMATION RESULTS
In this section, the proposed methods de¦ned earlier are tested for a rotating airfoil case.NACA 2412 is chosen as the initial airfoil geometry.There are two di¨erent mesh types considered during the analysis: inviscid mesh and viscous mesh.Figure 1 illustrates the initial baseline mesh for these two mesh types, respectively.The zoom view near the trailing edge is made as well since near that region, a huge displacement is expected to occur and the region where spring analogy might fail.The mesh deformation capability is checked by rotating the airfoil by a ¦xed amount with di¨erent mesh deformation capabilities: basic spring analogy, angle inclusion in spring analogy with direct method, and angle inclusion in spring analogy with indirect method.The deformed mesh using basic spring analogy is depicted in Fig. 2. It can be seen clearly that this method cannot handle a huge displacement for the inviscid case.In the viscous case, this method also fails in terms of not being able to maintain the right angle that viscous cells near the boundary have.In order to overcome this problem, an improvement in terms of angle consideration is applied into the basic spring analogy.The deformed meshes based on this method are depicted in Fig. 3.It can be seen clearly that by considering edge angle into the spring analogy formulation, 50 • can still be applied for inviscid mesh and 25 • can be applied for viscous mesh.Viscous mesh cannot be rotated more than 25 • since the presence of high aspect ratio cell near the boundary becomes a hindrance to the spring analogy formulation to perform a large deformation.However, it can be seen clearly that in the rotated viscous mesh condition, the right angle condition can still be maintained.Although the rotation made by viscous mesh is smaller, the required deformation for shape parameterization is already encompassed by this modi¦ed spring analogy method.
Another thing which can be observed is that the solutions obtained from direct and indirect methods are almost similar to one another.This is already expected that the solution method does not change the main idea of the spring analogy methods.However, in terms of e©ciency, the direct method is much faster compared to the indirect method.The fact that indirect method deals with solving a huge matrix requires a huge memory and consumes much time to get a convergence.On the other hand, solving 2 × 2 matrix in direct methods requires relatively smaller memory and less time.As a result, the angle inclusion spring analogy with direct method is implemented during the airfoil shape optimization.

SHAPE PARAMETERIZATION AND OPTIMIZATION
The best method described in the previous section is later coupled with CFD solver to perform an airfoil optimization.In the optimization, several shape parameterizations are considered: camber variation, thickness variation, and combined camber and thickness variation.In all these shape parameterizations, initial camber line and thickness distribution should be determined.The camber line is estimated as the average between the ordinate of the upper and lower airfoil surfaces with the same abscissa value.In the case where airfoil points are not located on the same abscissa, spline interpolation is utilized.The ¦nal expression for the camber line is On the other hand, the airfoil thickness is estimated as the ordinate di¨erence between the camber line and the airfoil surface.The expression for airfoil thickness is y thick upper = y upper − y camber ; y thick lower = y lower − y camber .Based on this parameterization, the design variables that can be applied are: camber factor which changes the camber variation of the airfoil and thickness factor which changes the thickness variation of the airfoil.The new camber and thickness are calculated by multiplying the initial camber and thickness by camber and thickness factor, respectively.These factors have their own ranges which are tabulated in Table 1.
After computing the new camber and thickness distribution, the airfoil£s ordinate value is updated.The updated value for the airfoil£s ordinate is computed based on y upper final = y camber final + y thick upper ; y lower final = y camber final + y thick lower .
The optimization itself was conducted by using a gradient optimizer provided by Phoenix ModelCenter, OPTLIB Gradient Optimizer.During the optimization, the airfoil drag is considered as the objective function and a speci¦ed lift coe©cient is considered.The lift coe©cient is computed based on the required sectional lift for the airfoil and chord length of the airfoil.In this study, the sectional lift is taken as 61.3125N/m with a chord length of 0.6 m.
There are two di¨erent §ight conditions considered during the optimization scheme, loiter and takeo¨.In each case, di¨erent angle of attack constraint is also imposed.In loiter phase, the angle of attack is constrained to be between 0 • and 7 • .On the other hand, angle of attacks are constrained to be between −3 • and 6 • for takeo¨ §ight condition.The details of the §ow parameters for these two §ight conditions are tabulated in Table 2.
The lift and drag coe©cients of the airfoil are computed based on CFD analysis using SU2 CFD Solver.Viscous solver equipped with SpalartAllmaras is chosen as the §ow solver.The input mesh used for the analysis is based on the deformed mesh obtained from the mesh deformation analysis.The optimization scheme using Phoenix ModelCenter contains three di¨erent modules: input module, mesh deformation module, and CFD solver module.The relations between each module are shown in Fig. 4.

RESULTS
This section presents the optimization results based on shape parameterization described in the earlier section.The optimization is done for both loiter and takeo¨ §ight conditions de¦ned earlier.In each case, NACA 2412 airfoil is chosen as the initial airfoil shape.By considering the camber e¨ect only in the optimization, the optimum airfoil has a relatively increase in camber by a factor of 1.747.On the other hand, by considering thickness only, the optimum airfoil has a relatively decrease in thickness by a factor of 0.714.Finally, the combined camber and thickness consideration in the optimization scheme leads to an optimum airfoil which has an increase in camber by a factor of 2.485 and decrease in thickness by a factor of 0.6, which is the minimum value.
The drag and angle of attack values corresponding to initial and ¦nal con¦gurations are tabulated in Table 3.It can be seen clearly that the optimization Figure 6 Optimized design for loiter (a) and takeo¨(b) con¦gurations with several shape parameterizations: upper row ¡ optimum airfoil shapes; lower row ¡ Cp distribution for optimum airfoil; 1 ¡ camber only; 2 ¡ thickness only; and 3 ¡ camber and thickness combined which considers both e¨ect of camber and thickness variation leads to a better optimum value among the proposed shape parameterization schemes.Furthermore, the optimum airfoil has a relatively low angle of attack compared to the initial airfoil.The shapes of optimum airfoils with their pressure distribution are shown in Fig. 6a.

Takeo¨Phase Optimization
The required airfoil coe©cient for these §ow parameters is computed as 0.3725.Summary of iteration history for this optimization is shown in Fig. 5b.
The summary of optimization results for takeo¨phase is tabulated in Table 4.It can be clearly seen that the optimization conducted by considering both camber and thickness variation leads to a better design.For a better illustration, the shape and pressure distributions for the optimum airfoils are shown in Fig. 6b.Unlike the loiter phase optimization, where an optimum airfoil has a high camber, the optimum airfoil for takeo¨phase by considering the presence of camber only has a camber reduction.This is due to the fact that the baseline airfoil has a relatively high thickness and the required lift coe©cient is quite low.In fact, the optimum airfoil has a reduction in camber by a factor of 0.8333.In the case where shape modi¦cation is solely based on thickness factor, the optimum airfoil for this §ight phase has a similar trend to loiter phase optimization.The optimum airfoil has a reduction in thickness by a factor of 0.6, which is the minimum allowable value.For the optimization by considering both camber and thickness variation, the optimum airfoil has increasing in camber by a factor of 1.452 and decreasing in thickness by a factor of 0.6.

CONCLUDING REMARKS
This paper has elaborated the idea of mesh deformation technique by using spring analogy and its modi¦cation by considering the angle made by the edge of the mesh.This modi¦ed method can be solved by using either direct solution method which solves unknown displacement iteratively for each node or indirect method by constructing a huge matrix that corresponds to the whole unknown displacements.It is found that the modi¦ed spring analogy can deform inviscid mesh to a higher amount of deformation and viscous mesh to a smaller amount of deformation.Furthermore, the solution achieved by using direct method is better than the indirect method in terms of computation time.
In the airfoil CFD design optimization for two di¨erent §ight conditions, it is found that performing an optimization by considering both camber and thickness variation leads to a better optimum design.In both cases of §ight conditions, it is found that the optimum airfoil has a relatively high camber and low thickness distribution.However, the loiter optimization requires a high camber compared to takeo¨optimization since loiter phase is accompanied by a relatively high lift coe©cient.

Figure 1
Figure 1 Initial baseline inviscid (a) and viscous (b) meshes around NACA 2412: left column ¡ outer view; and right column ¡ zoom view near the trailing edge

Figure 3
Figure 3 Deformed mesh case using angle inclusion in spring analogy with direct (a) and indirect (b) methods: left column ¡ 50 • rotated inviscid mesh airfoil; and right column ¡ 25 • rotated viscous mesh airfoil

Figure 4
Figure 4 Flow chart used in the optimization scheme

Table 1
Range of design variables used in the

Table 2
Summary of §ow parameters used in the optimization

Table 3
Summary of optimization results for loiter phase

Table 4
Summary of optimization results for takeo¨phase