STUDY OF MESH QUALITY IMPROVEMENT FOR CFD ANALYSIS OF AN AIRFOIL

ABSTRACT:  Airfoils generate lift in engineering applications such as for airplanes, wind turbines, automotive spoilers, etc. For accurate CFD analysis of airfoils, the quality of the mesh is of paramount importance, especially when dealing with turbulent flows commonly encountered in real life applications. Currently there are different tools that are available to improve the quality of the mesh required for CFD studies. This paper describes a study to assess the significant of the quality of the mesh on CFD analyses of NACA 23012 airfoil by using selected open source tools. The turbulence is modeled using the well-known k-Ï‰ Shear Stress Transport model. For validation, results have been compared with experimental datasets which were obtained from â€œTAG Stuttgart #1â€ tunnel. 
ABSTRAK: Sayap pesawat dapat menghasilkan daya angkat dalam aplikasi kejuruteraan seperti kapal terbang, turbin angin, spoiler automotif, dan sebagainya. Kualiti pada jaringan adalah amat penting bagi mendapatkan analisa CFD yang tepat pada sayap pesawat, terutamanya apabila berhadapan situasi aliran turbulen sebenar. Pada masa ini terdapat pelbagai perisian bagi meningkatkan mutu jaringan dalam kajian CFD. Kertas kerja ini membentangkan satu kajian bagi menilai kepentingan kualiti jaringan pada analisis CFD bagi sayap pesawat NACA 23012 dengan menggunakan sumber terpilih perisian terbuka. Model turbulen dibangunkan mengguna pakai model k-Ï‰ Shear Stress Transport (SST) yang terkenal. Bagi pengesahan, keputusan uji kaji telah dibandingkan dengan set data yang diperoleh dari terowong "TAG Stuttgart #1â€."


INTRODUCTION
Airfoils are widely used for different engineering applications like airplanes, wind turbines, automotive spoilers, etc. Simulating the aerodynamic performance of airfoils at different angles of attack (AOA) is essential for most engineering applications which rely on the lift (L) to drag (D) ratio. With the improvement in computing hardware and parallel algorithms, researchers and industry are increasingly relying on CFD simulations.
The airfoil NACA 23012 [1] has been chosen in the present study with airfoil chord length c = 1m. NACA 23012 is interesting for some CFD analyses due to is high aerodynamic performance at low flying velocities. It develops high maximum lift coefficient (CLmax) and high stall angle (αstall). According to reference [1] , its geometric characteristics are as follows:  Thickness: 12%  Maximum thickness position: 30%  Maximum chamber: 1.83%  Maximum chamber position: 13% In the present work, the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations are numerically approximated by using the open source finite volume solver OpenFOAM [2]. In the present research, a parametric CFD study of the flow around a NACA airfoil have been conducted and the significance of mesh quality on prediction of lift and drag coefficients have been discussed.

MESH QUALITY
In CFD simulations, the accuracy of the solutions are highly dependent upon on the quality of the mesh. There are several mesh quality metrics, some of the most critical parameters are orthogonality, skewness, smoothness (also termed as uniformity, growth factor, growth rate or change in cell size) and aspect ratio. It should be noted that the present study is done for two dimensional situations only.
According to reference [3] -"skewness has an adverse effect on the accuracy of interpolation on the face, non-orthogonality increases the error of the approximation of the surface-normal gradient, and nonuniformity reduces the order of the approximation of the surface-normal gradient to first order". Interested readers can consult references [4]- [7] for more information about mesh quality and influence these parameters for successful CFD analysis. Brief explanation regarding each parameters that affect the mesh quality is described below as stated in reference [8]. Fig. 1 shows a simple representation of mesh non-orthogonality. By definition, mesh non-orthogonality is the angular deviation of the vector S (located at face center f) from vector d connecting the two adjacent cell centers P and N. Mesh non-orthogonality mainly affects the diffusive terms as it adds diffusion to the solution [8].

Skewness
According to reference [8], mesh skewness is the deviation of vector d that connects the two adjacent cells for example P and N from the face center, f as shown in Fig. 2. It mainly adds diffusion to the solution and affects convective terms of the solution [8].

Smoothness
Smoothness (also known as expansion rate or growth rate) is the ratio of transition in size between cells. Further explanation is as illustrated in Fig. 3. Large transition ratios between cells add diffusion to the solution [8].

Aspect Ratio
Mesh aspect ratio is the ratio between ∆x and ∆y as shown in Fig. 4. High AR will smear the gradients, therefore it adds numerical diffusion to the solution [8].

METHODOLOGY
To demostrate the importance of mesh qualiy on the solution accuracy, we conducted a parametric study using open source tools. Different stages of the study and the tools used are shown in Fig. 5.

Geometry Preparation
The initial stage in a CFD study is preparation of the geometry. To prepare the geometry of the computational domand aaround a NACA 23012 airfoil, we used GNU Octave [9] and a script was written following the guidelines given in reference [10]. To generate the geometry, this script used the following inputs: Among these parameters,"Bump" and "N" vector are critical for generating good quality meshes around airfoils. A brief description about these two parameters are given below as decribed in [10].
1. Bump: defines the degree of grid concentration around leading and trailing edge of the airfoil. When bump=1 the cell length is uniformly distributed along the airfoil contour. When bump>1 the grid are concentrated along the midchord of the airfoil. If bump is between 0 and 1, the grid will be concentrated around the edges. 6. Specified yPlus (y + ) : The value of y + mesh variable specified at the beginning of the analysis. y + value of 50 is use for the analyses.
In this study, geometry files for generating computational domains around a NACA 23012 airfoil have been prepared for two different angles of attack which are 0 o and 10 o . Fig. 6 shows the computational domains of the analyses.

Meshing
The geometrical data file generated by the script is then used to generate mesh using an open-source software, GMSH [11]. The mesh files from GMSH can be converted into the OpenFOAM mesh format using the "gmshToFoam" command in the OpenFOAM environment. Tables 1 and 2 show selected input parameters which can alter the quality of the meshes. Numerous meshes with different qualities have been investigated and only the selected ones are presented in these tables. After preparing the meshes, the quality of each mesh have been assessed using the "checkMesh" command in OpenFOAM and the key quality related parameters are listed in Tables 1 and 2.  To determine the quality of the mesh, the output parameters which are maximum aspect ratio (AR), maximum non-orthogonality and maximum skewness of mesh have been considered. For an ideal mesh, the maximum AR should be closer to 1 and the maximum skewness should be as low as possible. Also, the maximum mesh non-orthogonality should be less than 70%.
The high quality mesh case (mesh ID: M14 and U14) and the low quality mesh case (Mesh ID: M22 and U22) were simulated and the resulting aerodynamic coefficients were compared to experimental data [12].

CFD Solutions
In this study, we used the open-source numerical library OpenFOAM version 3.0.1 [2], which is a collocated unstructured finite volume solver. The CFD cases were run in a LINUX (Ubuntu 14.04 LTS) platform with 4 processor, maximum speed of 2.59GHz and a total of 12GB of RAM memory. For turbulence modelling, we used the well validated k-ω SST turbulence model of Menter with wall functions [13]. The flow is specified to be incompressible subsonic flow with turbulence intensity (TI) is set to be 0.02% with Reynolds Number (Re) of 1x10 6 and Mach Number (M 2 ) of 0.04. The cases is solved using the second order unbounded gradient discretization scheme. Each case described in tables 1 and 2, is solved using a method that is at least second order accurate in space. The cases were run until the quantities of interest did not show any oscillations.

Post-processing
The results obtained with OpenFOAM were analyzed using the post-processing software ParaView [14]. Also, quantities of interests (such as the residuals, lift and drag coefficient), were plotted using GNUplot [15].
To show the wake topology, we used the line-integral convulation technique (LIC), implemented in ParaView [14]. This technique, generates high density streamlines plots in two-dimensional view [16]. The LIC technique to show wake topology in this study is more preferable compared to the traditional streamlines technique due to its visualization feature of condense lines which shows the development of wake structure much clearer.

RESULTS AND DISCUSSION
In the present study, several meshes were investigated at angle of attack, AOA, of 0° and 10 o and the best mesh (in terms of mesh quality and its suitability to resolve the physics involve) was selected for subsequent detailed analysis at different AOA. To perform the steady flow simulations, the simpleFoam solver (which uses the SIMPLE algorithm for pressure-velocity coupling) was used. The results obtained from OpenFOAM are compared with experimental results adopted from experimental study at Stuttgart Tunnel #1 in reference [12]. The experimental was carried out in a wind tunnel at six Reynold Numbers ranging from 700,000 to 3,000,000 using a NACA 23012 at TI of 0.02%.
To setup the case in OpenFOAM, proper turbulence related parameters were calculated for a turbulence intensity of 0.02% in order to match the experimental results. Table 3 and 4 show the aerodynamic coefficient of the analysed case studies.  M14 and U14 are the studies using good quality mesh for AOA=0 o and AOA=10 o cases respectively as mentioned before. It can be seen that the obtained results closely matched the experimental results. It can also be seen that although the drag coefficient doesn't show much effect, lift coefficient shows significant difference. In the present study, a wall function has been used to resolve the boundary layer around the NACA airfoil. However, modeling drag forces for flow over airfoils in CFD analysis is quite challenging and better results are expected from transition and turbulence models.  Even with the same computational capability and flow solution, M14 and U14 cases took lesser run time compared to both M22 and U22 case. This shows that with low quality mesh, it is likely to obtain acceptable results but at a cost of increased computation time. Keep in mind that non-orthogonality, skewness and non-uniformity only cause problems in the regions with high gradient variation. Therefore, in this case, mesh quality does not affect drag value much for all 4 cases. Fig. 9 and 10 show the post-processing result from both case analysis of AOA=0 o and AOA=10 o using the Line-integral Convolution (LiC) technique of visualization. It can be observed that there is no wake structure developed from both of the cases.

CONCLUSION
In this article, the effect of mesh quality on the lift and drag coefficients was studied. From the results obtained, it is concluded that the mesh quality plays an important role on the reliability, accuracy and convergence speed of the solution. It was observed that nonorthogonality and skewness has a strong influence on the solution outcome, and these values need to be keep as low as possible, which is difficult, especially for industrial meshes. From the study, we can see that lift and drag coefficient value closely matched the experimental value when using a good quality mesh. Although, drag value for low quality mesh shows an acceptable agreement with the experimental result, lift value differs significantly. It may be due to the cells with high mesh non-orthogonality and high skewness is at an area where does not affect CD.