Abstract
This paper proposes a novel method to implement the Kutta condition in irrotational, inviscid, incompressible flow (potential flow) over an airfoil. In contrast to common practice, this method is not based on the panel method. It is based on a finite difference scheme formulated on a boundary-fitted grid using an O-type elliptic grid generation technique. The proposed algorithm uses a novel and fast procedure to implement the Kutta condition by calculating the stream function over the airfoil surface through the derived expression for the airfoils with both finite trailing edge angle and cusped trailing edge. The results obtained show the excellent agreement with the results from analytical and panel methods thereby confirming the accuracy and correctness of the proposed method.
1. Introduction
The advent of high speed digital computers has revolutionized the numerical treatment of fluid dynamics problems. Numerical methods, nowadays, have become a routine tool to investigate fluid flows over the bodies such as airfoil. Amongst such fluid flows, incompressible potential flows are of crucial importance in studying the low-speed aerodynamics problems. The limitations associated with the exact (analytical) solutions with complex variables methods (conformal mapping) motivated fluid dynamicists to develop numerical techniques to solve incompressible potential flow problems (the Laplace’s equation) over an airfoil. Since the late 1960s, the panel methods have become the standard aerodynamic tools to numerically treat such flows [1]. Panel methods are applicable to any fluid-dynamic problem governed by Laplace’s equation. In these methods, the airfoil surface is divided into piecewise straight line segments or panels and singularities such as source, doublet, and vortex of unknown strength are distributed on each panel. Panel method used for the simulation of an incompressible potential flow past an airfoil is concerned with the vortex panel strength and circulation quantities and the evaluation of such quantities results in the calculation of the velocity distribution over the airfoil surface and hence the determination of the pressure coefficients. These methods have been extensively investigated in the aerodynamics literature [2–6], so these will not be discussed any further here. The interested reader can refer to the above references for further information. However, dealing with the panels and their attributes is numerically much more complex than the method proposed in this paper and of high programming effort. The Kutta condition should be introduced into the computational loop in order to solve the derived system of equations for the vortex panel strengths. In this paper, we propose a novel method to numerically solve the incompressible potential flow over an airfoil which is exempt from considering the quantities such as the vortex panel strength and circulation. This method takes advantage of an -type elliptic grid generation technique to generate the grid over the flow domain and approximate the flow field quantities such as stream function, velocity, and pressure at the grid points. The Kutta condition is implemented into the computational loop by an exact derived expression. An expression is derived for the finite-angle and cusped trailing edges. Finally, the obtained results from the proposed method are compared to those from the standard literature (both analytical and numerical) through several test cases.
2. Governing Equation for Irrotational, Incompressible Flow: Laplace Equation
Consider the irrotational, incompressible flow over an airfoil (Figure 1). The flow is governed by the Laplace’s equation ( is the stream function). The boundary conditions are as shown in Figure 1.
2.1. Conditions at Infinity
Far away from the airfoil surface (toward infinity), in all directions, the flow approaches the uniform free stream conditions. If the angle of attack (AOA) is and the free stream velocity is , then the components of the flow velocity can be written as where and are components of velocity vector ; that is, ( and are the unit vectors in and directions, resp.).
2.2. Condition on the Airfoil Surface
For inviscid flow, flow cannot penetrate the airfoil surface. Thus the velocity vector must be tangent to the surface. This wall boundary condition can be expressed by where is tangent to the surface. In the problem of the flow over an airfoil, if the free stream velocity and the angle of attack are known, from the boundary conditions at infinity (see (1) and (2)) and the wall boundary condition (see (3)) one can compute the stream function at any point of the physical domain (flow region). Then, by knowing , one can compute the velocity of all points in the physical domain. Since, for an incompressible flow, the pressure coefficient is a function of velocity only, one can obtain the pressure at any point in the flow region, as will be shown.
2.3. Pressure Coefficient
The pressure coefficient is defined as At standard sea level conditions,
3. Grid Generation
We have now presented all relations needed to obtain the pressure distribution in an incompressible, irrotational, inviscid flow over an airfoil. To calculate the pressure at any point in the flow region, a grid should be generated over the region. The elliptic grid generation proposed by Thompson et al. [7] is based on solving a system of elliptic partial differential equations to distribute nodes in the interior of the physical domain by mapping the irregular physical domain from the and physical plane (Figure 2) onto the and computational plane (Figure 3), which is a regular region. It is based on solving the Poisson equations as follows: where and are the computational coordinates corresponding to and in the physical coordinate, respectively. and are grid control functions which control the density of grids towards a specified coordinate line or about a specific grid point. To find an explicit relation for and in terms of grid points () and (), the following relations may be used: where The solution of the above equations (using the finite difference method to discretize the terms) gives and coordinates (in the physical domain) of coordinate in the computational domain.
The -type elliptic grid generation is employed here which results in a smooth and orthogonal grid over the airfoil surface. The -type elliptic grid generation technique has the advantage that the grid around the airfoil is orthogonal. The discretization of the physical domain and the corresponding computational domain are shown in Figures 2 and 3, respectively. In the computational domain, and are the number of nodes in the and directions, respectively. The resulting -type grid scheme over an airfoil for the case and or is shown in Figure 4.
| (a) Close-up view of -type grid around the airfoil |
(b) Magnified view of grid around the leading edge
(c) Magnified view of grid around the trailing edge
The initial guess for the elliptic grid generation is performed using the Transfinite Interpolation (TFI) method. Since the TFI method is an algebraic technique and does not need much time to generate grids over the physical domain, it will be an appropriate initial guess for the elliptic grid generation method and accelerate the convergence time for the elliptic grid generation method. Another advantage of using the TFI method as an initial guess is that it prevents the grids generated by elliptic (-type) method from folding.
4. Solution Approach
Since and are known, the stream function at any point in the physical domain can be obtained from (1) and (2) as follows: where subscripts and refer to any two arbitrary grid points at the physical domain boundaries. Equations (9) are applied to vertical and horizontal boundaries of the physical domain, respectively. By knowing the values of stream function on boundaries of the physical domain as well as on the airfoil surface (from wall boundary condition), we can obtain the values of over the physical domain. Since we deal with Laplace’s equation, it is necessary to find relationships for the transformation of the first and second derivatives of the field variable with respect to the position variables and . By using the chain rule, it can be concluded that By interchanging and , and and , the following relationships can also be derived: By solving (11) for and , we finally obtain where is Jacobian of the transformation. To transform terms in the Laplace equation, the second order derivatives are needed. Therefore, one has the following.
In the physical domain ,
After transformation, in the computational domain , where and and are control functions which may be assumed to be zero in both the grid generation and the flow solver sections (). These assumptions lead to the following equation to solve the above Laplace’s equation and obtain at every grid point of the physical domain: To solve the above equation, the finite difference method may be conveniently used.
4.1. Kutta Condition
The Kutta condition states that the flow leaves the sharp trailing edge of an airfoil smoothly [8]. To apply the Kutta condition in our calculation, we need to consider two possible configurations of the trailing edge. The trailing edge can have a finite-angle or can be cusped (Figure 5).
(a)
(b)
Suppose that the velocities along the top surface and bottom surface are and , respectively. For a finite-angle trailing edge, having two finite velocities in two different directions at the same point is physically impossible (Figure 5(a)) and, therefore, the only possibility is that both velocities should be zero (). For the cusped trailing edge (Figure 5(b)), having two velocities in the same directions at point shows that both and can be finite. However, the pressure at point is unique and Bernoulli equation states that [2] or In order to obtain relationships for the Kutta condition in terms of stream function , consider the finite-angle trailing edge in the -type grid scheme shown in Figure 6.
From (1), we have From the transformation relationship (see (13)), If and are the velocities of the grid points and , respectively, the Kutta condition gives By discretizing (22) in the computational domain, we get By considering the wall boundary condition (), we can simplify (23) to get Since the grid points and are the same points in the physical domain, we have This value is constant on the airfoil surface due to the wall boundary condition.
The derivation of an equation for the cusped trailing edge is more complicated. Consider the cusped trailing edge and the associated grid notation shown in Figure 7.
Since for the cusped trailing edge both vectors and are equal in the magnitude and direction, from the Kutta condition for the cusped trailing edge () we can write But Since and we have In similar approach, we have Furthermore, and . Since and , Moreover, and . Since (wall boundary condition), we obtain By substituting (28) through (31) into (26), we have By solving (32) for (using software Maple), we get In addition, . Equation (33) is the required expression for the cusped trailing angle.
Figures 8 and 9 show the stream function for both the finite-angle (NACA 0012 airfoil with angle of attack of and a free stream velocity of ) and the cusped (NACA 64012 with angle of attack of and a free stream velocity of ) trailing edge, respectively.
4.2. Velocity Calculation
There are three sections where the velocity must be known:(1)the outer boundaries (four sides , , , and of the rectangle shown in Figure 2),(2)the airfoil surface ( in Figure 2),(3)the inside of the physical domain.The velocity values on the outer boundaries are known from the conditions at infinity (using (1) and (2)). In other words, -component of the velocity vector () on all the outer boundaries is equal to and -component of the velocity vector () on all the outer boundaries is equal to . For the inside of the physical domain and the airfoil surface, we can use (12) and (13) as follows: The central and one-sided difference schemes are used for the inside of the physical domain and the airfoil surface, respectively. After obtaining the components of the velocity vector, the total velocity can be computed by As stated before, for an incompressible flow, the pressure coefficient can be expressed in terms of velocity only. Thus (4) can be used to determine the pressure of any grid point in the domain. Therefore,
4.3. Kutta Condition in Terms of the Velocity Potential
The proposed method can be easily developed in terms of the velocity potential . The wall boundary condition may be expressed in terms of either the velocity potential , , or the stream function , , where and are the unit vector normal to the airfoil surface and the distance along the body (airfoil) surface, respectively. Using the transformation relationships for mapping the physical domain onto the computational one, we can write where , , and are defined in (16). The solution of the above equation for the airfoil surface using the finite difference method gives the value for (). From the definition of the velocity potential, In a similar way to the derivation for Kutta condition in terms of the stream function given in (21) to (23), we get And finally By including (37) and (40) into the solution loops, we can find the velocity potential over the domain. The above procedure also can be extended to the three dimension case.
5. Results
5.1. Validation of the Results for the Pressure Distribution
The results obtained here are compared with the results from using the panel method. The results are obtained by a Fortran compiler (PGI) and computations are run on a PC with Intel Pentium Dual 1.73 and 1 G RAM. The tolerance used in the iterative loops (the mesh generation and the stream function) is .
5.2. Trailing Edge with Finite-Angle
Validation Case 1. The pressure coefficient distribution () over a NACA 0012 airfoil at an angle of attack of is plotted. The results are compared with the results from [2]. The -type grid size used in the computation is . The computation time is 53 seconds (see Figure 10).
Validation Case 2. The pressure coefficient distribution () over a NACA 0024 airfoil at an angle of attack of is plotted. The results are compared with the results from [9]. The -type grid size used in the computation is . The computation time is 41 seconds (see Figure 11).
Validation Case 3. The pressure coefficient distribution () over a NACA 4414 airfoil at an angle of attack of is plotted. The results are compared with the results from the software Xfoil [10]. The -type grid size used in the computation is . The computation time is 51 seconds (see Figure 12).
Validation Case 4. The pressure coefficient distribution () over a NACA 4412 airfoil at an angle of attack of is plotted. The results are compared with the results in [5]. The -type grid size used in the computation is . The computation time is 55 seconds (see Figure 13).
5.3. Cusped Trailing Edge
Validation Case 1. The pressure coefficient distribution () over a NACA 64012 airfoil at an angle of attack of is plotted. The results are compared with the results from the software XFLR5 [11]. The -type grid size used in the computation is . The computation time is 4 minutes and 15 seconds (see Figures 14, 15, and 16).
Excellent agreement can be obtained by comparing the results from our method and the ones from the panel method given in validation cases for both finite-angle and the cusped trailing edges. As shown in the validation cases results, the maximum value for is exactly equal to 1. The pressure coefficient at the trailing edge (T.E.) is equal to unity because the velocity is zero at this stagnation point. Accordingly, For the cusped trailing edge, . Thus the value of at T.E. is not equal to 1 (), as shown in Figure 16.
6. Conclusion
This paper presents a novel method to implement the Kutta condition in the numerical solution of two-dimensional incompressible potential flow over an airfoil. The proposed method is based on solving the Laplace’s equation for the stream function at each grid point generated by the elliptic grid generation technique (-type). Therefore, it is exempt from considering the panels and the quantities such as the vortex panel strength and circulation used in the panel method. It applies for both finite-angle and cusped trailing edges. A novel and very easy to implement expression for the stream function for the finite-angle and the cusped trailing edges is derived. The accurate results obtained for both cases show the correctness and accuracy of the numerical scheme.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.