2D and 3D Potential Flow Simulation around NACA 0012 Airfoil with Ground E ﬀ ect

. The purpose of this paper is to develop a Boundary Element Based Method (BEM) for determining the steady potential about two and three dimensional airfoil. The numerical investigation of NACA 0012 airfoil with using Boundary Element Method is utilized. Two different physical problems of the NACA 0012 airfoil are examined: potential flow around airfoil in an unbounded fluid and potential flow prediction with ground effect. Computation of potential flow around the airfoil is investigated by the mixed constant strength source and constant strength dipole based panel method. Boundary Element Code is written in FORTRAN. To check the accuracy of the 2D boundary element based code, the validation studies are carried out by comparing the present results obtained for the NACA 0012 airfoil from the XFoil and other published simulation results. 3D results are also evaluating with the available experimental and other numerical simulation results. The numerical outcomes are examined in terms of pressure distribution and lifting force on the foil.


Introduction
It is known that for an airfoil flying close to ground clearance can produce a considerable lift growth as compared to airfoil flying out of the ground effect. In order to investigate the ground effect, it is essential to first calculate the aerodynamic forces on a wing. Dionne and Lee [1] experimentally investigated ground effect on aerodynamics of airfoil. In this study lift increment was found to increase with reducing ground distance.
Numerical simulations are developing a widespread technique to assess the aerodynamic and hydrodynamic performance in preliminary design stage. The model experiments are still very valuable, hence time and cost encourage the use of Computational Fluid Dynamics (CFD). With the improvements in computer technology, scientists have started to study flow around the airfoil with solving Navier-Stokes equations. By the help of computers, it has been able to model the flow around airfoil bodies by using suitable turbulence model. RANS based numerical approaches are widely used in airfoil design. Winslow et al. [2] investigated airfoil characteristics at low Reynolds number. In this study, the computational fluid dynamics tool used is a Reynoldsaveraged Navier-Stokes solver with a Spalart-Allmaras turbulence model. Rubel et al., [3], Sohaib Obeid et al. [4] have carried out some numerical analyses in order to decide airfoil performance. Alsarheed and Sedaghat [5] computationally investigated ground effects using both the panel and finite volume methods. Their study show that the results of both panel and finite volume numerical methods are consistent are in agreement with the experimental results. The flow is incompressible, irrotational and inviscid. Some CFD methods based on potential flow theory because the Reynolds Averaged Navier-Stokes Equations (RANSE) is both computer power and time consuming.
The prediction of pressure distribution of an airfoil is important for the airfoil's performance. Some methods can be found in the literature, but the most reliable ones are developed by Hess and Smith [6]. They used the panel method in solving the problem of lifting and non-lifting bodies by distribution of constant sources and vortices on them. Several references that have been printed such as Morino and Kuo [7], Kerwin [8] and the books by (Anderson [9], Moran [10], Newmann [11], Katz and Plotkin [12], Katsikadelis [13]), so these will not be examined any further here.
Over the last few decades, Potential based numerical approaches are widely used in aerodynamic design. Tarafder et al. [14] and Tarafder et al. [15] studied on the computation of potential flow around 2D and 3D NACA 0012 airfoil by the source based panel methods respectively. They evaluated their results with the analytical and experimental data. Chaves [16] analysed 2D airfoil by using the panel method. His work showed that trends of lift coefficient and pressure distribution results follow the experimental and CFD results. Farinha et al. [17] examined the application of panel method, which is able to solve steady, two-dimensional, potential flow around airfoils including the ground effect. Their study emphasized that at small angles of attack, the findings obtained from this simple approach are qualitatively similar to experimental results.
Potential theory not only used for aerodynamic design but also used for hydrodynamic calculations. Theoretical analyses on estimation of wave pattern generated by hydrofoils and ships have been carried out for along time. As a pioneering study, Dawson [18] suggested an efficient linearized method, which gives convincing results. Nowadays, Dawson's method successfully using for simulation of free surface flow. Tarafder and Suzuki [19] modeled a potential-based boundary element method for solving a nonlinear free-surface flow problem for a Wigley catamaran moving with a uniform speed in deep water. Their study successfully presents a potential-based boundary element method for predicting the hydrodynamic characteristics of a catamaran hull. Xie and Vassalos [20] develop a potential-based panel method for determining the steady potential flow about three dimensional hydrofoil under free surface. An IBEM [iterative boundary element (panel) method] for the solution of hydrofoil transfering under a free surface was explained in detail in Bal and Kinnas [21]. Their iterative method is extended to apply to the surface piercing bodies and some numerical results of the method are given in Bal [22], Bal [23].
The main purpose of this study is to generate the boundary element based solution code. 2D and 3D simulations for various angle of attack in order to be able to evaluate results from the Xfoil and the other published results from literature. The model was solved with a range of different attack of angles from 0 o to 10 o . 3D simulations were carried out different span length. In the following sections, the computational results of aerodynamic coefficient and surface pressure distribution for NACA 0012 are presented. The boundary element code written in FORTRAN, the program is able to find the lift force and pressure distribution. The ability of the panel methods in solving the potential flow problem around 2-D and 3-D lifting airfoils is experienced in this study.

Geometry and conditions
In this study, the NACA 0012, the well-documented airfoil was used. For the four digit numbers, the first number describes the maximum camber as a percentage of chord length. The second number is the position is the maximum tenths of the chord and third and fourth number corresponds to the maximum thickness of the foil as a percent of the chord length. This airfoil is symmetrical as both first and second number is zero, and maximum thickness of 12% of the chord length. The thickness of the NACA 0012 airfoil can be calculated from a well-known equation (Abbott and Von Doenhoff [24]) and geometry of the airfoil is given in Fig. 1. In order to validate the present simulation process, the operating conditions are mirrored to equalled the operating conditions of NASA Langley Research Center validation cases (Abbott and Von Doenhoff [24]).

Mathematical formulations
It is accepted that the fluid is incompressible, inviscid and that the flow is irrotational. A Cartesian coordinate system is placed on foil surface the components of the free stream velocity U in the x-z frame of reference. The angle of attack α is defined as the angle between the free stream velocity and the x-axis. The total velocity potential is separated into sum of a free stream potential and a perturbation potential (Bal [25], Lee et al. [26]) Where φ is the perturbation potential. The total velocity potential ϕ should satisfy the Laplace's equation in the fluid domain Ω.
(2) The domain Ω is bounded by the airfoil surface S A , wake surface S W and an outer control surface S ∞ enclosing the body and the wake surfaces. The definition of notations, airfoil and the fluid domain is depicted in Fig. 2. The problem can be constructed by specifying the boundary conditions as follows: Kinematic body boundary condition: Where ⃗ is the unit normal vector, which points out of the fluid domain. Radiation condition at infinity: Kutta condition at the trailing edge of a foil: Assumption of the wake surface: streamwise velocity must be continuous across the surface: subscript U and L are the upper and lower surfaces of the wake, and V and P are the velocity and pressure respectively. is a unit vector in the direction of the mean velocity. The wake surface has zero thickness and the pressure jump across S W is zero, while there is a jump in the potential:

Numerical implementation
The velocity potential ϕ can be written an integral equation on the airfoil surface by applying Green's theorem: where S A and S W are the boundaries of the airfoil surface and wake surface respectively. G is the Green function given as = 1 4 in 3D and = 1 2 ln in 2D, r is the location vector of the singularities and the field point.
To construct a numerical solution, the surface S A divided into N panels and the discretized form of Eq. (9) is expressed on the airfoil surface: The values of and are constant on each element, so they can be moved outside the integral. Eq. (10) may be written as:  Fig. 4 shows the surface panels of the NACA 0012 airfoil and span length (S=10).
These elements estimate the actual geometry by a constant panel while the unknown quantity supposed to be constant. Both G ij and H ij is calculated analytically. A more detailed discussion of this point is provided in Katz and Plotkin (2001). Eq. (13) may be written as: is the Kronecker delta, which is defined as = 0 for ≠ and = 1 for = . Determination of the influence coefficient the perturbation potentials at the each of the collocation point will result in an × influence matrix, with N+1 unknown (where the wake potential Δ is the ( + 1)-th unknown). The additional equation is provided by using the Kutta condition: The above set of algebraic equations is solved by using Gaussian elimination method.

Numerical results and discussion
After perturbation potential on hydrofoil is calculated, the local external tangential velocity on each collocation point can be computed by differentiating the velocity potential along the tangential direction as follows Katz and Plotkin (2001): Pressure coefficient can be obtained by using Bernoulli Equation: , U Airfoils have different geometry and dimensions. So, non-dimensional coefficients (lift coeffic NACA 0012 S/c=5 ients) were considerate to evaluate the advantages and disadvantages of airfoils. The non-dime nsional lift and drag coefficients were given as below: Where L is the lift force and D is the drag force. The contribution of element to the aerodynamic loads ∆ : Where ∆ is the panel area. The total aerodynamic loads are then obtained by adding the contributions of the individual panels.

Unbounded fluid
We examine the dependency of the solution to the number of panels by comparing force results with increasing numbers of panels.  constant as the number of panels increase. It can be seen also from this figure non-uniform spacing with 34 panels arrangement gives accurate results Simulation for several angle of attack were performed in order to able to compare the results from the present simulations and numerical results from the XFoil potential based code and the validate them with existing experimental data from reliable sources. Table 1 shows the computed lift coefficient with using both present study and XFoil results for 2D NACA 0012. The comparison present simulation with XFoil is very good. It is accepted that incoming velocity is unity, because the potential solution in coefficient form is independent of scale. The lift coefficients increased linearly with the increasing angles of attack. Present BEM model had a good agreement with XFoil results angles of attacks from 0 o to 10 o . The lift coefficients are shown for various angle of attack, computed with present study, and compared with the includes numerical results from the XFoil are given Fig. 7.  Fig. 7. From the given figure, it is observed that some distinctions are observed with comparison of experimental results of Abbott & von Doenhoff [14]. For the attack of angle lower than 10° agreement is quite impressive with the present numerical simulation and other numerical computations and experimental results. Flow is attached to the airfoil through this regime. For the attack of angle that are greater than 10° some peculiarities can be seen and viscous effects begin to show up as a reduction in lift with increasing angle of attack. Figure 7 shows that potential theory cannot estimate stall angle. At angle of attack roughly 14-17°, stall began to develop. Variations of computed lift coefficient C L with angle of attack for 3D NACA 0012 for different span wise length are depicted in Fig.8. The triangle in the figure is the numerical results of Tarafder et al. [15]. Some distinctions are seen from this figure. Our explanation for this discrepancy is that the mesh structure is insufficient in the study of Tarafder et al. They used 18 chord wise panels and 5 span wise panels for discretization of geometry.
As the angle of attack is increased, greater values of C P are viewed on the lower surfaces and the region of the high pressure extends towards to the trailing edge. It can be seen from Fig. 9 and

Ground effect
In this part our numerical results are compared with experimental simulation for symmetrical 3D NACA 0012 geometries. To provide certainty in the current results, a validation has been completed against printed results for similar studies. Experimental measurements were carried out by N. Moore et al. [27] in Southampton University. The main dimension of test case is given Table  2. In Boundary element method, the ground effect simulated with the method of images, which consist on placing a mirrored image of the airfoils with the wall as a symmetry axis, as shown in Figure 11. Figure 12 shows the measured and computed lift coefficient values for given speeds for the NACA 0012 airfoil model. For ground effect, the comparison with the experimental measures of N. Moore et al. [27] is very good for attack of angles between 1 and 6. Both experimental and numerical results reveal that for attack of angle between 3° and 6°, lift coefficient decrease with the growth in the altitude. This is due to the greater pressure between airfoil and ground. The BEM computations generally predicted higher lift coefficient values. The lift curve slope gradients (∂C L /∂α) are summarized in Table 3, showing comparison with the experimental data Moore et al. [27] in Figure 13. Comparisons between experimental and numerical results reveal that for altitude (h/c) between 0.1 and 0.3, curve slope gradient differences (∂C L /∂α) are increased. The obtained differences are believed due to the potential theories deficiency. Figure 14 illustrates the comparisons of lift performance of the NACA 0012 section between h/c 0.3 to 0.1. An remarkable observation is reported by Moore et al. When the NACA 0012 is at h/c=0.1 large downward forces produced when AoA reduces flow the three degrees. Potential theory cannot simulate this forces but good agreement the other altitudes. Differences between our numerical results and experiment Moore et al. [27] are reduced with the increase of ground clearance. Figures 16 shows the simulation results of pressure coefficient at angles of attack 1°, 6°, 10° and h/c=0.5 with BEM method. The lower surface of the aerofoil obtained a higher pressure coefficient compared to upper surface. As the angle of attack increased the lower surface pressure distribution was much higher than the pressure distribution of the upper surface. That was awaited result. Figure 16 depicts ground effect on the aerofoil pressure distribution. It can be seen that the effect of the ground on aerodynamic performance is significant when the airfoil is located near the ground surface.

Conclusion
The paper deals with the source and dipole based panel methods for computing the potential flow around 2D and 3D NACA 0012 airfoil moving with uniform speed in an unbounded fluid. This paper shows the importance of employing BEM method technique for the solution of flow around 2D and 3D bodies. The following conclusions can be obtained from the present study: • Author's BEM code is a low order panel method, panels do not have to represent exactly boundaries, • This code does not taking into account viscosity effect, • Results cannot estimate airfoil stall angle, • The given simulation is limited to both single panel structure and panel number, • Present 2D BEM simulations are in good agreement with XFoil results, • Author's BEM code fairly simulates the flow around the wing in ground effect, • Wing in ground effect with clear wing becomes lift coefficient greater than in free stream flow, • The BEM method has been showed to be computational efficient in analyzing the airfoil geometry, • Although the accurate estimation of the lift coefficient cannot be awaited with based on potential theory, numerical solutions are generally good agreement with experimental results, • The method is restricted to the evaluation of lift coefficient, • Potential theory can simulate the suction effect on the lower surface at small ground clearance, • This paper shows validation of BEM method for only one model NACA 0012 symmetrical hydrofoil, • Our BEM code is a relatively basic program, with restrictions. As compared with engineering CFD codes. In this study viscosity effect is not taken into account, • The method needs less computer power and small analyzing times, • Potential theory is very practical method compared to RANSE techniques.
It may be concluded that the proposed method can evaluate lift coefficients with sufficient accuracy by using potential flow theory and can be used as a design tool for airfoils. For the future work, examination of the effect of panel structure and different panel numbers is intended. Moreover, it is planned to investigate the usefulness of the proposed method on different airfoil types.  Discretization of 3D NACA 0012 airfoil by 34x10 panels (S=5) Figure 5 Change of drag coe cient with number of panels Comparison between experimental data from Abbott & von Doenhoff [24] and Xfoil and present simulation of lift coe cient for 2D airfoil Figure 8 Comparison between potential ow simulation data from Tarafder et al. [15] and present study of lift coe cient for 3D airfoil Figure 9 Variation of surface pressure coe cient for ow past 2D NACA 0012 airfoil, α=4°, α=8° Variation of CL with altitude for NACA 0012 airfoil   Variation of surface pressure coe cient for ow past 3D NACA 0012 airfoil, α=1°, α=6°, α=10° and h/c=0.5

Figure 16
Variation of surface pressure coe cient for ow past 3D NACA 0012 airfoil, h/c=0.1, h/c=0.3, h/c=1 and α=1° Figure 17 Pressure distrubution on the NACA 0012 airfoil surface for α=1 deg at different values of h/c at y/c=0 section