# On the Kutta Condition in Potential Flow over Airfoil

Authors: Farzad Mohebbi; Mathieu Sellier
Journal: Journal of Aerodynamics
 (2014)

---

## 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.


---

## Body

## 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 O-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 (Figure1). The flow is governed by the Laplace’s equation ∇2ψ=0 (ψ is the stream function). The boundary conditions are as shown in Figure 1.Figure 1
Boundary conditions at infinity and on the airfoil surface (no penetration).


### 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 V∞, then the components of the flow velocity can be written as
(1)u=∂ψ∂y=V∞cos⁡α,(2)v=-∂ψ∂x=V∞sinα,
where u and v are components of velocity vector V; that is, V=ui+vj  (i and j are the unit vectors in x and y directions, resp.).


### 2.2. Condition on the Airfoil Surface

For inviscid flow, flow cannot penetrate the airfoil surface. Thus the velocity vector must betangent to the surface. This wall boundary condition can be expressed by
(3)∂ψ∂s=0orψ=constant,
where s is tangent to the surface. In the problem of the flow over an airfoil, if the free stream velocity V∞ 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 coefficientCp is defined as
(4)Cp≡p-p∞(1/2)ρ∞V∞2=1-(VV∞)2.
At standard sea level conditions,
(5)ρ∞=1.23kg/m3,p∞=1.01×105N/m2.


## 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 V∞, then the components of the flow velocity can be written as
(1)u=∂ψ∂y=V∞cos⁡α,(2)v=-∂ψ∂x=V∞sinα,
where u and v are components of velocity vector V; that is, V=ui+vj  (i and j are the unit vectors in x and y directions, resp.).


## 2.2. Condition on the Airfoil Surface

For inviscid flow, flow cannot penetrate the airfoil surface. Thus the velocity vector must betangent to the surface. This wall boundary condition can be expressed by
(3)∂ψ∂s=0orψ=constant,
where s is tangent to the surface. In the problem of the flow over an airfoil, if the free stream velocity V∞ 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 coefficientCp is defined as
(4)Cp≡p-p∞(1/2)ρ∞V∞2=1-(VV∞)2.
At standard sea level conditions,
(5)ρ∞=1.23kg/m3,p∞=1.01×105N/m2.


## 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 x and y 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:
(6)ξxx+ξyy=P(ξ,η),ηxx+ηyy=Q(ξ,η),
where ξ and η are the computational coordinates corresponding to x and y in the physical coordinate, respectively. P and Q aregrid 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 x and y in terms of grid points ξi (i∈[1,M]) and ηj (j∈[1,N]), the following relations may be used:
(7)αxξξ-2βxξη+γxηη=-J2(P(ξ,η)xξ+Q(ξ,η)xη),αyξξ-2βyξη+γyηη=-J2(P(ξ,η)yξ+Q(ξ,η)yη),
where
(8)α=xη2+yη2,β=xξxη+yξyη,γ=xξ2+yξ2J=xξyη-xηyξ(Jacobianoftransformation).
The solution of the above equations (using the finite difference method to discretize the terms) gives x and y coordinates (in the physical domain) of coordinate (i,j) in the computational domain.Figure 2
The physical domain.O-type scheme and discretization of the boundaries.Figure 3
The computational domain showing the discretization of the physical domain boundaries.TheO-type elliptic grid generation is employed here which results in a smooth and orthogonal grid over the airfoil surface. The O-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, M and N=2N1+2N2+N3+6 are the number of nodes in the ξ and η directions, respectively. The resulting O-type grid scheme over an airfoil for the case N2=N1 and N3=2N1-1 or N=6N1+5 is shown in Figure 4.O-type grid (elliptic) around an airfoil. The figure illustrates orthogonality and smoothness of the gridlines especially near airfoil surface.
(a)
Close-up view ofO-type grid around the airfoil
(b)
Magnified view of grid around the leading edge
(c)
Magnified view of grid around the trailing edgeThe 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 (O-type) method from folding.


## 4. Solution Approach

SinceV∞ and α are known, the stream function ψ at any point in the physical domain can be obtained from (1) and (2) as follows:
(9)ψi=ψj+(yi-yj)V∞cos⁡α,ψi=ψj-(xi-xj)V∞sinα,
where subscripts i and j 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 x and y. By using the chain rule, it can be concluded that
(10)∂ψ∂x=∂ψ∂ξ∂ξ∂x+∂ψ∂η∂η∂x=∂ψ∂ξξx+∂ψ∂ηηx,∂ψ∂y=∂ψ∂ξ∂ξ∂y+∂ψ∂η∂η∂y=∂ψ∂ξξy+∂ψ∂ηηy.
By interchanging x and ξ, and y and η, the following relationships can also be derived:
(11)∂ψ∂ξ=∂ψ∂x∂x∂ξ+∂ψ∂y∂y∂ξ=∂ψ∂xxξ+∂ψ∂yyξ,∂ψ∂η=∂ψ∂x∂x∂η+∂ψ∂y∂y∂η=∂ψ∂xxη+∂ψ∂yyη.
By solving (11) for ∂ψ/∂x and ∂ψ/∂y, we finally obtain
(12)∂ψ∂x=1J(yη∂ψ∂ξ-yξ∂ψ∂η),(13)∂ψ∂y=1J(-xη∂ψ∂ξ+xξ∂ψ∂η),
where J=xξyη-xηyξ 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(x,y),
(14)∇2ψ=∂2ψ∂x2+∂2ψ∂y2=0.After transformation, in the computational domain(ξ,η),
(15)∇2ψ=1J2(αψξξ-2βψξη+γψηη)+[(∇2ξ)ψξ+(∇2η)ψη],
where
(16)α=xη2+yη2,β=xξxη+yξyη,γ=xξ2+yξ2,J=xξyη-xηyξ(Jacobianoftransformation)
and ∇2ξ=P and ∇2η=Q are control functions which may be assumed to be zero in both the grid generation and the flow solver sections (P=Q=0). These assumptions lead to the following equation to solve the above Laplace’s equation and obtain ψ at every grid point of the physical domain:
(17)((xη2+yη2)ψξξ-2(xξxη+yξyη)ψξη+(xξ2+yξ2)ψηη)=0.
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).Different possible shapes of the trailing edge and their relation to the Kutta condition.
(a)
(b)Suppose that the velocities along the top surface and bottom surface areV1 and V2, 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 (V1=V2=0). For the cusped trailing edge (Figure 5(b)), having two velocities in the same directions at point a shows that both V1 and V2 can be finite. However, the pressure at point a is unique and Bernoulli equation states that [2]
(18)pa+12ρV12=pa+12ρV22
or
(19)V1=V2.
In order to obtain relationships for the Kutta condition in terms of stream function ψ, consider the finite-angle trailing edge in the O-type grid scheme shown in Figure 6.Figure 6
Grid notation of the trailing edge.From (1), we have
(20)u=ψy.
From the transformation relationship (see (13)),
(21)ψy=1J[-(xη)(ψξ)+(xξ)(ψη)].
If V1 and VN are the velocities of the grid points (1,1) and (1,N), respectively, the Kutta condition V1=V2=0 gives
(22)V1=VN=0⟹u1=uN=0,1J[-(xη)(ψξ)+(xξ)(ψη)]|1=1J[-(xη)(ψξ)+(xξ)(ψη)]|N=0-xηψξ+xξψη|1=0.
By discretizing (22) in the computational domain, we get
(23)xηψξ=xξψη,[(x1,2-x1,1)][(ψ2,1-ψ1,1)]=[(x2,1-x1,1)][(ψ1,2-ψ1,1)],ψ1,1=ψ2,1(x1,2-x1,1)-ψ1,2(x2,1-x1,1)x1,2-x2,1.
By considering the wall boundary condition (ψ1,1=ψ1,2), we can simplify (23) to get
(24)ψ1,1=ψ2,1.
Since the grid points (1,1) and (1,N) are the same points in the physical domain, we have
(25)ψ1,1=ψ1,N=ψ2,1.
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 Figure7.Figure 7
Cusped trailing edge and the associated grid notation.Since for the cusped trailing edge both vectorsV1 and VN are equal in the magnitude and direction, from the Kutta condition for the cusped trailing edge (V1=VN) we can write
(26)V1=VN⟹u1=uN,[1J(-xηψξ+xξψη)]1=[1J(-xηψξ+xξψη)]N.
But
(27)xξ|1=x2,1-x1,1,xξ|N=x2,N-x1,N.
Since x2,N=x2,1 and x1,N=x1,1 we have
(28)xξ|1=xξ|N.
In similar approach, we have
(29)yξ|1=yξ|N.
Furthermore, ψξ|1=ψ2,1-ψ1,1 and ψξ|N=ψ2,N-ψ1,N. Since ψ1,1=ψ1,N and ψ2,1=ψ2,N,
(30)ψξ|1=ψξ|N.
Moreover, ψη|1=ψ1,2-ψ1,1 and ψη|N=ψ1,N-ψ1,N-1. Since ψ1,1=ψ1,2=ψ1,N-1=ψ1,N (wall boundary condition), we obtain
(31)ψη|1=ψη|N=0.
By substituting (28) through (31) into (26), we have
(32)1(x2,1-x1,1)(y1,2-y1,1)-(x1,2-x1,1)(y2,1-y1,1)×[-(x1,2-x1,1)(ψ2,1-ψ1,1)+0]=1(x2,1-x1,1)(y1,N-y1,N-1)-(x1,N-x1,N-1)(y2,1-y1,1)×[-(x1,N-x1,N-1)(ψ2,1-ψ1,1)+0].
By solving (32) for ψ1,1(using software Maple), we get
(33)ψ1,1=ψ2,1.
In addition, ψ1,N=ψ1,1=ψ2,1. Equation (33) is the required expression for the cusped trailing angle.Figures8 and 9 show the stream function ψ for both the finite-angle (NACA 0012 airfoil with angle of attack of α=40∘ and a free stream velocity of V∞=70m/s) and the cusped (NACA 64012 with angle of attack of α=40∘ and a free stream velocity of V∞=70m/s) trailing edge, respectively.Figure 8
Stream function for a finite-angle trailing edge. The figure shows the Kutta condition at the trailing edge.Figure 9
Stream function for a cusped trailing edge. The figure shows the Kutta condition at the trailing edge.


### 4.2. Velocity Calculation

There are three sections where the velocity must be known:(1)
the outer boundaries (four sidesCD, DE, EF, and FC of the rectangle shown in Figure 2),(2)
the airfoil surface (AH 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, x-component of the velocity vector (u) on all the outer boundaries is equal to V∞cos⁡α and y-component of the velocity vector (v) on all the outer boundaries is equal to V∞sinα. For the inside of the physical domain and the airfoil surface, we can use (12) and (13) as follows:
(34)ui,j=∂ψ∂y|i,j=1J[-(xη)i,j(ψξ)i,j+(xξ)i,j(ψη)i,j],vi,j=-∂ψ∂x|i,j=-1J[(yη)i,j(ψξ)i,j-(yξ)i,j(ψη)i,j].
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
(35)Vi,j=ui,j2+vi,j2.
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,
(36)pi,j=12ρ(V∞2-Vi,j2)+p∞.


### 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 ϕ, (∂ϕ/∂n=0), or the stream function ψ, (∂ψ/∂s=0), where n and s 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
(37)∂ϕ∂nairfoilsurface=-1Jα(αϕξ-βϕη)=0,
where J, α, and β are defined in (16). The solution of the above equation for the airfoil surface using the finite difference method gives the value for ϕ1,j (j=1,…,N). From the definition of the velocity potential,
(38)v=ϕy.
In a similar way to the derivation for Kutta condition in terms of the stream function given in (21) to (23), we get
(39)ϕy=1J[-(xη)(ϕξ)+(xξ)(ϕη)],V1=VN=0⟹v1=vN=0.
And finally
(40)ϕ1,1=ϕ2,1(x1,2-x1,1)-ϕ1,2(x2,1-x1,1)x1,2-x2,1.
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.


## 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).Different possible shapes of the trailing edge and their relation to the Kutta condition.
(a)
(b)Suppose that the velocities along the top surface and bottom surface areV1 and V2, 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 (V1=V2=0). For the cusped trailing edge (Figure 5(b)), having two velocities in the same directions at point a shows that both V1 and V2 can be finite. However, the pressure at point a is unique and Bernoulli equation states that [2]
(18)pa+12ρV12=pa+12ρV22
or
(19)V1=V2.
In order to obtain relationships for the Kutta condition in terms of stream function ψ, consider the finite-angle trailing edge in the O-type grid scheme shown in Figure 6.Figure 6
Grid notation of the trailing edge.From (1), we have
(20)u=ψy.
From the transformation relationship (see (13)),
(21)ψy=1J[-(xη)(ψξ)+(xξ)(ψη)].
If V1 and VN are the velocities of the grid points (1,1) and (1,N), respectively, the Kutta condition V1=V2=0 gives
(22)V1=VN=0⟹u1=uN=0,1J[-(xη)(ψξ)+(xξ)(ψη)]|1=1J[-(xη)(ψξ)+(xξ)(ψη)]|N=0-xηψξ+xξψη|1=0.
By discretizing (22) in the computational domain, we get
(23)xηψξ=xξψη,[(x1,2-x1,1)][(ψ2,1-ψ1,1)]=[(x2,1-x1,1)][(ψ1,2-ψ1,1)],ψ1,1=ψ2,1(x1,2-x1,1)-ψ1,2(x2,1-x1,1)x1,2-x2,1.
By considering the wall boundary condition (ψ1,1=ψ1,2), we can simplify (23) to get
(24)ψ1,1=ψ2,1.
Since the grid points (1,1) and (1,N) are the same points in the physical domain, we have
(25)ψ1,1=ψ1,N=ψ2,1.
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 Figure7.Figure 7
Cusped trailing edge and the associated grid notation.Since for the cusped trailing edge both vectorsV1 and VN are equal in the magnitude and direction, from the Kutta condition for the cusped trailing edge (V1=VN) we can write
(26)V1=VN⟹u1=uN,[1J(-xηψξ+xξψη)]1=[1J(-xηψξ+xξψη)]N.
But
(27)xξ|1=x2,1-x1,1,xξ|N=x2,N-x1,N.
Since x2,N=x2,1 and x1,N=x1,1 we have
(28)xξ|1=xξ|N.
In similar approach, we have
(29)yξ|1=yξ|N.
Furthermore, ψξ|1=ψ2,1-ψ1,1 and ψξ|N=ψ2,N-ψ1,N. Since ψ1,1=ψ1,N and ψ2,1=ψ2,N,
(30)ψξ|1=ψξ|N.
Moreover, ψη|1=ψ1,2-ψ1,1 and ψη|N=ψ1,N-ψ1,N-1. Since ψ1,1=ψ1,2=ψ1,N-1=ψ1,N (wall boundary condition), we obtain
(31)ψη|1=ψη|N=0.
By substituting (28) through (31) into (26), we have
(32)1(x2,1-x1,1)(y1,2-y1,1)-(x1,2-x1,1)(y2,1-y1,1)×[-(x1,2-x1,1)(ψ2,1-ψ1,1)+0]=1(x2,1-x1,1)(y1,N-y1,N-1)-(x1,N-x1,N-1)(y2,1-y1,1)×[-(x1,N-x1,N-1)(ψ2,1-ψ1,1)+0].
By solving (32) for ψ1,1(using software Maple), we get
(33)ψ1,1=ψ2,1.
In addition, ψ1,N=ψ1,1=ψ2,1. Equation (33) is the required expression for the cusped trailing angle.Figures8 and 9 show the stream function ψ for both the finite-angle (NACA 0012 airfoil with angle of attack of α=40∘ and a free stream velocity of V∞=70m/s) and the cusped (NACA 64012 with angle of attack of α=40∘ and a free stream velocity of V∞=70m/s) trailing edge, respectively.Figure 8
Stream function for a finite-angle trailing edge. The figure shows the Kutta condition at the trailing edge.Figure 9
Stream function for a cusped trailing edge. The figure shows the Kutta condition at the trailing edge.


## 4.2. Velocity Calculation

There are three sections where the velocity must be known:(1)
the outer boundaries (four sidesCD, DE, EF, and FC of the rectangle shown in Figure 2),(2)
the airfoil surface (AH 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, x-component of the velocity vector (u) on all the outer boundaries is equal to V∞cos⁡α and y-component of the velocity vector (v) on all the outer boundaries is equal to V∞sinα. For the inside of the physical domain and the airfoil surface, we can use (12) and (13) as follows:
(34)ui,j=∂ψ∂y|i,j=1J[-(xη)i,j(ψξ)i,j+(xξ)i,j(ψη)i,j],vi,j=-∂ψ∂x|i,j=-1J[(yη)i,j(ψξ)i,j-(yξ)i,j(ψη)i,j].
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
(35)Vi,j=ui,j2+vi,j2.
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,
(36)pi,j=12ρ(V∞2-Vi,j2)+p∞.


## 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 ϕ, (∂ϕ/∂n=0), or the stream function ψ, (∂ψ/∂s=0), where n and s 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
(37)∂ϕ∂nairfoilsurface=-1Jα(αϕξ-βϕη)=0,
where J, α, and β are defined in (16). The solution of the above equation for the airfoil surface using the finite difference method gives the value for ϕ1,j (j=1,…,N). From the definition of the velocity potential,
(38)v=ϕy.
In a similar way to the derivation for Kutta condition in terms of the stream function given in (21) to (23), we get
(39)ϕy=1J[-(xη)(ϕξ)+(xξ)(ϕη)],V1=VN=0⟹v1=vN=0.
And finally
(40)ϕ1,1=ϕ2,1(x1,2-x1,1)-ϕ1,2(x2,1-x1,1)x1,2-x2,1.
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) is10-8.


### 5.2. Trailing Edge with Finite-Angle

Validation Case  1. The pressure coefficient distribution (Cp) over a NACA 0012 airfoil at an angle of attack of α=9∘ is plotted. The results are compared with the results from [2]. The O-type grid size used in the computation is 155×155. The computation time is 53 seconds (see Figure 10).Figure 10
Comparison between the results from [2] and the results from our method for validation case 1. The figure shows an excellent agreement between the results.Validation Case  2. The pressure coefficient distribution (Cp) over a NACA 0024 airfoil at an angle of attack of α=0∘ is plotted. The results are compared with the results from [9]. The O-type grid size used in the computation is 155×155. The computation time is 41 seconds (see Figure 11).Figure 11
Comparison between the results from [9] and the results from our method for validation case 2. The figure shows an excellent agreement between the results.Validation Case  3. The pressure coefficient distribution (Cp) over a NACA 4414 airfoil at an angle of attack of α=2∘ is plotted. The results are compared with the results from the software Xfoil [10]. The O-type grid size used in the computation is 155×155. The computation time is 51 seconds (see Figure 12).Figure 12
Comparison between the results from [10] and the results from our method for validation case 3. The figure shows an excellent agreement between the results.Validation Case  4. The pressure coefficient distribution (Cp) over a NACA 4412 airfoil at an angle of attack of α=10∘ is plotted. The results are compared with the results in [5]. The O-type grid size used in the computation is 155×155. The computation time is 55 seconds (see Figure 13).Figure 13
Comparison between the results from [5] and the results from our method for validation case 4. The figure shows an excellent agreement between the results.


### 5.3. Cusped Trailing Edge

Validation Case  1. The pressure coefficient distribution (Cp) over a NACA 64012 airfoil at an angle of attack of α=6∘ is plotted. The results are compared with the results from the software XFLR5 [11]. The O-type grid size used in the computation is 245×245. The computation time is 4 minutes and 15 seconds (see Figures 14, 15, and 16).Figure 14
Pressure distribution over the airfoil (NACA 64012) used in validation case 1 (cusped trailing edge).Figure 15
Pressure coefficient distribution over a NACA 64012 airfoil at a6∘ angle of attack obtained by XFLR5.Figure 16
Comparison between the results from the software XFLR5 and the results from our method for validation case 1 (cusped trailing edge). The figure shows an excellent agreement between the results.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 forCp 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,
(41)VT.E.=0,Cp=1-VT.E.2V∞2=1-0=1.
For the cusped trailing edge, VT.E.≠0. Thus the value of Cp at T.E. is not equal to 1  (Cp≠1), as shown in Figure 16.


## 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) is10-8.


## 5.2. Trailing Edge with Finite-Angle

Validation Case  1. The pressure coefficient distribution (Cp) over a NACA 0012 airfoil at an angle of attack of α=9∘ is plotted. The results are compared with the results from [2]. The O-type grid size used in the computation is 155×155. The computation time is 53 seconds (see Figure 10).Figure 10
Comparison between the results from [2] and the results from our method for validation case 1. The figure shows an excellent agreement between the results.Validation Case  2. The pressure coefficient distribution (Cp) over a NACA 0024 airfoil at an angle of attack of α=0∘ is plotted. The results are compared with the results from [9]. The O-type grid size used in the computation is 155×155. The computation time is 41 seconds (see Figure 11).Figure 11
Comparison between the results from [9] and the results from our method for validation case 2. The figure shows an excellent agreement between the results.Validation Case  3. The pressure coefficient distribution (Cp) over a NACA 4414 airfoil at an angle of attack of α=2∘ is plotted. The results are compared with the results from the software Xfoil [10]. The O-type grid size used in the computation is 155×155. The computation time is 51 seconds (see Figure 12).Figure 12
Comparison between the results from [10] and the results from our method for validation case 3. The figure shows an excellent agreement between the results.Validation Case  4. The pressure coefficient distribution (Cp) over a NACA 4412 airfoil at an angle of attack of α=10∘ is plotted. The results are compared with the results in [5]. The O-type grid size used in the computation is 155×155. The computation time is 55 seconds (see Figure 13).Figure 13
Comparison between the results from [5] and the results from our method for validation case 4. The figure shows an excellent agreement between the results.


## 5.3. Cusped Trailing Edge

Validation Case  1. The pressure coefficient distribution (Cp) over a NACA 64012 airfoil at an angle of attack of α=6∘ is plotted. The results are compared with the results from the software XFLR5 [11]. The O-type grid size used in the computation is 245×245. The computation time is 4 minutes and 15 seconds (see Figures 14, 15, and 16).Figure 14
Pressure distribution over the airfoil (NACA 64012) used in validation case 1 (cusped trailing edge).Figure 15
Pressure coefficient distribution over a NACA 64012 airfoil at a6∘ angle of attack obtained by XFLR5.Figure 16
Comparison between the results from the software XFLR5 and the results from our method for validation case 1 (cusped trailing edge). The figure shows an excellent agreement between the results.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 forCp 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,
(41)VT.E.=0,Cp=1-VT.E.2V∞2=1-0=1.
For the cusped trailing edge, VT.E.≠0. Thus the value of Cp at T.E. is not equal to 1  (Cp≠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 (O-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.



---

*Source: 676912-2014-04-01.xml*