Yiin WANG,Ning QIN,Ning ZHAO,*

aCollege of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China

bDepartment of Mechanical Engineering,The University of Sheffield,Sheffield S1 3JD,United Kingdom

KEYWORDS Delaunay graph mapping;Dynamic mesh;Grid deformation;Mesh;Moving mesh

Abstract The fluid–structure interaction and aerodynamic shape optimization usually involve the moving or deforming boundaries,thus the dynamic mesh techniques are the key techniques to cope with such deformation.A novel dynamic mesh method was developed based on the Delaunay graph in this paper.According to the Delaunay graph,the mesh points were divided into groups.In each group,a factor ranging from 0 to 1 was calculated based on the area/volume ratio.By introducing a proper function for this factor,this method can control the mesh quality with high efficiency.Several test cases were compared with other dynamic mesh methods regarding mesh quality and CPU time,such as radial basis function method and Delaunay graph mapping method.

1.Introduction

The dynamic mesh techniques are wieldy used in the fluid–structure interaction and aerodynamic shape optimization which involve the moving or deforming boundaries.For these numerical simulations,two aspects of the dynamic mesh method are very critical,one is the mesh quality after deforming,and the other is the efficiency.From the viewpoint of the mesh quality,physical analogy method,1–4such as the spring analogy approach5and elastic analogy,6,7can well maintain the mesh quality after the deformation.However,these methods need to solve certain equations,which normally require large computational cost.Particularly,for large shape changing,some of these methods may cause invalid cells,therefore,Farhat et al.8introduced the torsional springs to prevent the mesh from becoming invalid.Recently,some interpolation methods which can largely preserve the mesh quality near the boundaries were developed,e.g.Radial Basis Function(RBF)method9–11and Inverse Distance Weighting(IDW)interpolation method.12,13The RBF method,by using the displacement of the boundary nodes,constructs a matrix to interpolate the interior nodes from the original position to the new position.This method needs to solve two(for 2D)or three(for 3D)nbby nbmatrixes,where nbis the total number of boundary nodes.As the mesh size increases,the computational cost can be significant.To reduce the computational cost,Rendall and Allen14proposed an approximate RBF method with the data reduction algorithms.Though it improves the efficiency,some surface errors of mismatch may be introduced to the boundary.Later on,they developed a surface correction step to fix the problem.15The IDW interpolation method,developed by Witteveen and Bijl,12,13simply uses the inverse distance as a weight for the displacement of each boundary to calculate the new position of the interior nodes.This method is faster than the original RBF method;however,it still needs to loop all the boundary nodes for each interior node,resulting in a large computational cost for large 3D problems.Luke et al.16proposed an explicit interpolation method based on the IDW method,which shows similar mesh quality as the RBF method,but with a relatively faster speed.Recently,Zhou and Li17developed a 2D dynamic mesh method based on disk relaxation.Later on,they managed to further improve the method and successfully applied it to the 3D cases.18

For the fluid–structure interaction problem,the mesh needs to be updated at each time step;while for the aerodynamic shape optimization,each new configuration needs a renewed mesh.Hence,due to the frequent use of the dynamic mesh,the efficiency of the dynamic method becomes very important for these problems.From the viewpoint of efficiency,the interpolation method,such as transfinite interpolation19–21or Delaunay Graph Mapping(DGM)method,22can quickly deform the mesh.The transfinite interpolation,however,can only be used on simple structured mesh.The Delaunay Graph Mapping(DGM)method can handle most of the deformation with quality mesh,23–25but it is difficult to deal with the rotation motion.Considering its high efficiency,for large 3D problems with small deformation,it is a better choice than the RBF and IDW methods.In general,most of the current methods either generate high quality mesh with large computational cost,or generate lower quality mesh with high efficiency.Recently,a new method successfully combines the RBF and DGM methods with both advantages,26which is called the DGRBF.It is based on the Delaunay graph to divide the mesh nodes into groups,then uses the RBF method to interpolate the mesh nodes to its new position.As a result,the large matrix becomes a series of 3 by 3(for 2D)or 4 by 4(for 3D)small matrixes;therefore the computational cost is substantially decreased.Later on,the IDW function was implemented into the DGM method which shows similar mesh quality as the DGRBF but slightly better efficiency.

Fig.1 Contours of f with different r0s

As the RBF and IDW function can improve the mesh quality of the DGM method,are there any other functions with simpler form which can further improve both the mesh quality and efficiency?To achieve this,the functions should meet two critical conditions, firstly the functions can maintain the mesh quality near the wall;secondly the form of the function should be as simple as possible to maintain the high efficiency.In this paper,a new dynamic mesh method based on the above idea is developed.This method is based on the Delaunay graph to divide the mesh nodes into groups,and then it introduces a damp function for each group to control the mesh deforma-tion;therefore the overall mesh quality can be effectively controlled and improved.In addition,the damp function only uses local information;therefore the computational cost is low.Different forms of the functions are compared and discussed.Several test cases are compared with the RBF,IDW,DGM and DGRBF methods,which show higher mesh quality with higher efficiency.

Table 1 Damp functions.

Fig.2 Discontinuity occurs for nodes on the edge of Delaunay triangle.

Fig.3 Original mesh of rotating rectangle.

2.Delaunay graph mapping method with damp function

Fig.4 Comparison of mesh quality with different damp functions and DGM method.

Fig.5 Mesh quality contour(DGMDF by f1).

The general procedure of the Delaunay Graph Mapping with damp Function(DGMF)method is similar to the original DGM method.However,being different from the original one,it utilizes the area ratio to calculate an intermediate variable r rather than the coordinates.By using this variable r,a damp function which controls the displacement of the interior nodes is calculated.The procedures are set as follows:

Fig.6 Mesh quality contour(DGM).

Fig.7 Original mesh of rotating circle.

Fig.8 Comparison of mesh quality for rotating circle.

(1)Generate the Delaunay graph using all the boundary nodes of the original mesh.

(2)Find the parent Delaunay element for each mesh node point so that all the mesh points are grouped by the Delaunay Graph elements.

(3)Calculate r for each node using the area or volume ratio of Delaunay elements(triangle for 2D and tetrahedron for 3D).

(4)Calculate the damp function f（r）for each node.

(5)Calculate the new node position for all the mesh nodes using the damp function f（r）.

At step(3),an intermediate variable r is calculated by the area ratio as follows:where eiis the area or volume ratio(detailed definition can be found in Ref.22),and n=3 for 2D,n=4 for 3D.r0is the user defined constant for the boundary nodes,for all static boundaries r0s≥1,while for all dynamic boundaries r0d=0.By changing r0s,the impact region of f（ r）(the detailed form of the damp function f（r）will be discussed in the next paragraph)can be adjusted(Fig.1).In the figures,the outer boundary is the static boundary,while the inner boundary is the dynamic boundary.The mesh nodes in the red region which are adjacent to the dynamic boundary will move similarly as the dynamic boundary since f is close to 1,thus the mesh quality can be maintained in this region.The mesh nodes in the blue region which is close to the static boundary will almost keep unchanged,because f is close to 0.From the figures,it is clear that by changing the value of r0s,the size of the blue and red region can be adjusted.By increasing r0s,one can compress the red region and enlarge the blue region.

The damp function used in this paper is in the form of

where φ is a function ranging from 0 to 1.In order to maintain high efficiency,computational expensive operations are not considered,such as exponent calculation used in the IDW function.Therefore,the simple polynomial is used in this paper.In Table 1,the functions tested and compared in this paper are listed.These functions are generally the radial basis function or modified radial basis function whose value ranges from 0 to 1;different forms and orders can affect the change of the deformation of the mesh.The comparison of these functions will be discussed in Section 3.1.

The change of the boundary can be mainly classified into three categories,i.e.,rigid body translation,rigid body rotation and shape deformation.Let the displacement of translation be t,the Eulerian angle be θ,and the displacement of shape deformation be d,then the displacement of the interior nodes can be calculated as follows:

where x and y are the original coordinates,and xcandycare coordinates of the center.

Fig.10 Original mesh of twisted bar.

For the 3D problem,let θ = （α，β，γ）be the Eulerian angle,the total displacement of the interior nodes is in the form of

where x0is the original coordinates.Let q(fθ)be the quaternion,and Eq.(5)can be rewritten as

Fig.11 Deformation of the bar as A increases.

Fig.12 Comparison of mesh quality for twisted bar.

Fig.13 Mesh quality contour(A=16).

Normally for the numerical problems involving rigid body motion,the translation and rotation of the moving boundary are known.For deformation,if the moving boundary can be expressed as a known function,d can be directly computed;if not,d can be calculated in the following way.

where dithe deformation displacement of the dynamic boundary nodes in a Delaunay triangle/tetrahedron,and ndthe number of the dynamic boundary nodes in the Delaunay triangle/tetrahedron.

Due to the interpolation based on different triangles/tetrahedrons,discontinuity may occur for the DGRBF method.As shown in Fig.2,Node n is on the sharing edge of the Delaunay triangle△ABC and△ACD.Obviously its displacement can be computed by using either triangle.However,for shape deformation,the results can be inconsistent when different triangles are used.For the unstructured meshes,the nodes appearing on an edge of the Delaunay triangle maybe rare;however,for the structured meshes it happens from time to time.For the DGMF method,the discontinuity problem aforementioned is resolved.Taking Node n in Fig.2 as an example and presuming that Nodes B,C and D are dynamic boundary nodes,since Node n is on Edge AC,eband edare equal to 0;thus,for either Delaunay triangle,the displacement for the shape deformation computed from Eqs.(3),(4)and(7)are the same dn=fdc.

3.Results and discussion

3.1.Rotating rectangle

One of the major disadvantages of the original Delaunay mapping method is its difficulty in handling large rotational problems.Hence,as the first test case,the capability for dealing with the large rotational problem was tested.The original mesh is shown in Fig.3.

The rectangle in the middle rotated around its center from 10°to 90°.Therefore,the interior nodes of the mesh can be computed as

It should be noted that all the meshes were deformed directly from the original mesh without any intermediate steps in this paper.Therefore,only one Delaunay graph was needed.The r0s=1 was used for all the test cases in this paper;one can restrain the impact region by increasing this value.The effect of different damp functions is compared in Fig.4,regarding the averaged mesh quality and minimal mesh quality.In this figure,the result of the DGM method is also included.The size-skew metric Ft was used in this paper which assesses the skewness and the size of the mesh.27The essential properties of the size-skew metric are

Ft=1 element has equal angles and the same size as the initial element.

Ft=0 element is degenerated.

The damp function f2,f3and f4show better averaged mesh quality but worse minimal mesh quality,while f5and f6show worse averaged mesh quality with better minimal mesh quality.Function f1is in the middle place among the all.The high order functions(f2,f3and f4)slowly decay in the near boundary region,whereby the mesh quality in this region is well preserved;out of this region,they quickly decrease,resulting in a sudden transition whereby low-quality cells may be generated in this region,which cause a low minimal mesh quality.In contrast,the lower order functions(f4and f5)decrease in a relatively even way,hence the transition is quite smooth,resulting in a high minimal mesh quality.However,the mesh quality in the near boundary region is not as good as the high order functions.It is found in this figure that function f1shows moderate averaged mesh quality and minimal mesh quality,and the mesh quality near the boundary can be largely maintained.Therefore in this paper,this function was used as the damp function for the rest of the test cases.The deformed meshes are shown in Figs.5 and 6.The new method well preserves the mesh quality near the moving boundary;while the DGM method does not.

3.2.Rotating circle

To further test the capability of the new method in handling rotation motion,a rotating circle mesh(in Fig.7)was tested.The radius of the circle is 1 unit and the outer boundary is a square with a 10-unit-side.The mesh quality is compared in Fig.8.In this case,the DGMF method is compared with both the DGM and RBF methods.The DGMF method shows the best mesh quality among the three.As shown in Fig.9,the DGMF method well reserves the mesh quality for 90°rotation,whereby the mesh quality contour further verifies that the DGMF method can well handle the rotation motion.The DGM method,as mentioned earlier,was difficult to handle the rotation problem,therefore for 90°rotation it failed to generate a valid mesh.Based on the displacement of the boundary nodes,the RBF method interpolates the displacement of interior nodes.However,for this case,it does not well preserve the mesh quality near the boundary.From the two test cases,it is found that this new method can well deform the mesh for rigid body rotation problems.

3.3.Twisted bar with rotation

Normally,the DGM method is able to generate the high quality mesh for surface deformation problems such as twisted surface.In this case,the capability in shape deformation problem was tested along with the rotation motion.A bar with 10 unitlong and 1 unit thick was first twisted and then rotated.The original mesh is shown in Fig.10.The shape deformation is based on the bump function,which is

Table 2 Comparison of mesh quality for twisted bar with rotation.

Fig.14 Mesh quality contour of twisted bar with rotation.

The bar was bended to form a horseshoe shape,and as the parameter A increases the deformation becomes more severe in Fig.11.The displacement of the interior node for shape deformation can be computed as

or by Eq.(7).Based on the differences between the calculations of the shape deformation displacement,the former one is called DGMF1,and the latter one which uses Eq.(7)is called DGMF2.The overall mesh qualities of the resultant meshes by different methods are illustrated in Figs.12 and 13,the mesh quality contours are compared for A=16 case.As the A increases,both the averaged quality and minimal mesh quality decrease for all the methods.Especially the minimal mesh quality quickly descends to 0 for the RBF and DGRBF methods,i.e.,the meshes become invalid for A=32 case.The two DGMF methods and the DGM method show similar high averaged mesh quality and minimal mesh quality,for large deformation(A=32),all three methods can still generate valid mesh with proper averaged mesh quality.

Fig.15 Wing-body configurations.

The bar was twisted and then rotated by 45°.The deformed mesh was still based on the original mesh,particularly for the DGM,DGRBF,DGMF1 and DGMF2 methods,and the Delaunay graph was the same as the previous twisted case.The overall mesh quality and mesh quality contours are demonstrated in Table 2 and Fig.14 respectively.Both the RBF and DGM methods fail to preserve the mesh quality near the two ends of the bar.Otherwise,the DGRBF and DGMF1&2 methods well reserve the mesh quality,and show the best averaged and minimal mesh quality among the all.From all the 2D cases,it can be concluded that the DGMF method can well handle the rotation problems,and for the shape deformation,it shows similar quality mesh as the DGM and DGRBF methods do.

3.4.3D test cases

An unstructured mesh with 257909 nodes and 1507462 tetrahedrons based on a wing-body configuration was tested by the proposed DGMF1 method.The initial and deformed configurations are shown in Fig.15,in which the wing was folded by 30°.The deformation is

Therefore the displacement of the interior nodes is

The total CPU time for the DGM method is 0.356 s,while for the DGMF1 method,it is 0.356 s,too.The DGMF1 method shows the same efficiency as the DGM method.Further efficiency comparison will be discussed later.Fig.16 shows the mesh quality contour,where the mesh near the deformed wing is well preserved,and most of the meshes show very good mesh quality.

Fig.16 Mesh quality contour for deformed wing.

Table 3 Comparison of CPU time.

The efficiency of the method was also evaluated by different sizes of the meshes.In Table 3 the CPU time calculated by different methods are compared.All the performances are measured on a single core 2.7 GHz Inteli7 processor.In the comparison four different methods are used,namely the DGM,DGMF,DGRBF and IDW.The RBF method is too time-consuming,so it is not included in this comparison,and the detailed information can be found in Ref.26.The DGM,DGMF and DGRBF methods show quite similar efficiency,which are significantly better than the IDW method.The DGM and DGMF methods are slightly better than the DGRBF method.The CPU time of the three Delaunay graph-based methods increases almost linearly as the total nodes of the mesh increase;while for the IDW method,it increases as the product of increase in the total nodes and total boundary nodes.The required CPU time for the Delaunay graph-based methods can be roughly divided into three parts.The first part is the generation of the Delaunay graph,of which the cost is normally Nbplog(Nbp),where Nbprefers to the number of boundary points.The second part is the identification process which is NsNvpwhere Nsis the number of searches(according to the numerical test,the Nsis normally less than 10 for 3D)26and vp stands for volume points.The computational costs of these two parts are exactly the same for all the Delaunay graph-based methods.The third part is the mapping process in which the cost is Nvpfor all the Delaunay graph-based methods.The cost of the IDW method,however,is NbpNvp.From the perspective of computational complexity,it is clear that the costs of the Delaunay graph-based methods are similar but much lower than that of the IDW method.Thus it can be concluded that the Delaunay graph-based methods are much more efficient than the IDW method.

4.Conclusions

A novel dynamic method based on the Delaunay graph was developed.This method can deform the mesh with both high quality and high efficiency, especially for large rotation problem.A range of large deformation and rotation test cases shows that this method can preserve the mesh quality near the boundary for both structured and unstructured mesh.This method inherits the advantages of the efficiency of the Delaunay graph method while addressing its robustness problem for large rotation. From the large 3D test cases, the required CPU time using this method increases almost linearly as the total node number increases,indicating its suitability for extremely large mesh computations.

Acknowledgements

This work is partially funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions of China,National Natural Science Foundation of China(No.11432007)and Natural Science Foundation of Jiangsu Province of China(No.BK20140805).