Note : The formulation for vortex panel method is taken from the book “ Foundations of Aerodynamics: Bases of
Aerodynamic Design” by Arnold M. Kuethe and Chuen-Yen Chow
Vortex Panel Method
The vortex panel method is a discretized computational method for vortex sheet. Vortex panel method
is one of the panel methods which are used for solving incompressible, potential and lifting flows over
arbitrary bodies like airfoil. The mathematical model of vortex panel method along with the matlab code
is described in this section.
The airfoil is represented by a closed polygon of vortex panels. In this method, the circulation density on
each panel varies linearly from one corner to the other and is continuous across the corner, as indicated
in Fig 1. The n number of planes are assumed to be planar and are named in clockwise direction, starting
from the trailing edge and ending at the trailing edge.
Fig. 1 Vortex panels with linear varying strength distributed over the surface of an airfoil
The following program creates nodes on airfoil
nacaseries = input('Enter the 4-digit naca series = ');
c = input('Enter the chord length = ');
n = input('Enter the number of nodes = ');
a = input('Enter the angle of attack in degree = ');
s = num2str(nacaseries);
% creating nodes on airfoil
if numel(s)==2
s1 = str2double(s(1));
s2 = str2double(s(2));
m=0;p=0;t=(10*s1+s2)*0.01;
else
s1 = str2double(s(1));
s2 = str2double(s(2));
s3 = str2double(s(3));
s4 = str2double(s(4));
m = s1*0.01; p = s2*0.1 ; t = (10*s3+s4)*0.01;
end
for i= n:-1:1
theta = (i-n)*2*pi/n;
x = 0.5*c*(1+cos(theta));
if(x/c)<p
yc(n+1-i) = m*c/p^2*(2*p*(x/c)-(x/c)^2);
dydx(n+1-i) = (2*m/p^2)* (p-x/c);
beta(n+1-i) = atan(dydx(n+1-i));
else
yc(n+1-i) = m*c/(1-p)^2 * ((1-2*p)+2*p*(x/c)-(x/c)^2);
dydx(n+1-i) = (2*m/(1-p)^2)* (p-x/c);
beta(n+1-i) = atan(dydx(n+1-i));
end
yt=5*t*c*(0.2969*sqrt(x/c)-0.1260*(x/c)...
-0.3516*(x/c)^2+0.2843*(x/c)^3-0.1036*(x/c)^4);
if(i<(0.5*n+1))
xa(n+1-i)=x - yt*sin(beta(n+1-i));
ya(n+1-i)=yc(n+1-i)+yt*cos(beta(n+1-i));
else
xa(n+1-i)=x + yt*sin(beta(n+1-i));
ya(n+1-i)=yc(n+1-i)-yt*cos(beta(n+1-i));
end
end
xa(n+1)= c;
ya(n+1) = 0;
yc(n+1) = 0; % trailing edge
The following program creates control points on each panel
%This loop calculates the location of the control points
for i = 1:n
xmid(i) = (xa(i)+xa(i+1))/2;
ymid(i) = (ya(i)+ya(i+1))/2;
Sj(i) = sqrt((xa(i+1)-xa(i))^2+(ya(i+1)-ya(i))^2);%array of panel lengths
phi(i) = atan2((ya(i+1)-ya(i)),(xa(i+1)-xa(i)));
rhs(i)= sin(phi(i)-alpha);% RHS of equation 4
end
In the presence of a uniform flow V , at an angle of attack and n vortex panels, the velocity
potential at the ith control point ( xi , yi ) due to free stream and all other panels is
n ( s j ) −1 yi − y j
( xi , yi ) = V ( xi cos + yi sin ) − tan ds j
2 x −x (1)
j =1 j i j
Where,
sj
( s j ) = j + ( j +1 − j ) (2)
Sj
Applying boundary condition of impermeability, we get
( xi , yi ) = 0; i = 1, 2,..., n (3)
ni
Carrying out the involved differentiation and integration, we obtain
(C 'j + Cn 2ij 'j +1 ) = sin(i − );
n
n1ij i = 1, 2,..., n (4)
j =1
Here, ' = / 2 V is a dimensionless circulation density , i is the orientation angle of the ith panel
measured from the x axis to the panel surface , and the coefficients are expressed as follows:
Cn1ij = 0.5DF + CG − Cn 2ij
Cn 2ij = D + 0.5QF / S j − ( AC + DE )G / S j
The constants are defined as :
A = − ( xi − X j ) cos j − ( yi − Y j ) sin j
B = ( xi − X j ) + ( yi − Y j )
2 2
C = sin(i − j )
D = cos(i − j )
E = ( xi − X j ) sin j − ( yi − Y j ) cos j
S 2j + 2 AS j
F = ln 1 +
B
ES j
G = tan −1
B + AS
j
P = ( xi − X j ) sin(i − 2 j ) + ( yi − Y j ) cos(i − 2 j )
Q = ( xi − X j ) cos(i − 2 j ) − ( yi − Y j ) sin(i − 2 j )
(X − X j ) + (Y j +1 − Y j )
2 2
Sj = j +1
,
For i = j , the coefficients become:
Cn1ii = −1 and Cn1ii = 1
% This loop calculates the coefficients to find the matrix Aij for equation
% 6.Variables are as given in report.
for i = 1:n
for j = 1:n
if i==j
cn1(i,j) = -1;
cn2(i,j) = 1;
else
A=-(xmid(i)-xa(j))*cos(phi(j))-(ymid(i)-ya(j))*sin(phi(j));
B = (xmid(i)-xa(j))^2+(ymid(i)-ya(j))^2;
C = sin(phi(i)-phi(j));
D = cos(phi(i)-phi(j));
E = (xmid(i)-xa(j))*sin(phi(j))-(ymid(i)-ya(j))*cos(phi(j));
F = log(1+(Sj(j)^2+2*A*Sj(j))/B);
G = atan2((E*Sj(j)),(B+A*Sj(j)));
P = (xmid(i)-xa(j))*sin(phi(i)-2*phi(j))+(ymid(i)-ya(j))*cos(phi(i)-
2*phi(j));
Q = (xmid(i)-xa(j))*cos(phi(i)-2*phi(j))-(ymid(i)-ya(j))*sin(phi(i)-
2*phi(j));
cn2(i,j)=D+0.5*Q*F/Sj(j)-(A*C+D*E)*G/Sj(j);
cn1(i,j)=0.5*D*F+C*G-cn2(i,j);
end
end
end
which describe the self-induced normal velocity at the ith control point.
To ensure a smooth flow at the trailing edge ,the kutta condition is applied as follows:
(kutta condition demands that the strength of the vorticity at the trailing edge be zero )
1' + n' +1 = 0 (5)
There are altogether n+1 equations after combining equations (4) and (5) which are sufficient to solve
n+1 unknown j values. Rewriting this system of simultaneous equations in a more convenient form,
'
we get
n +1
A
j =1
nij 'j = RHSi ; i = 1, 2,..., n + 1 (6)
in which , for i n + 1 :
Ani1 = Cn1i1
Anij = Cn1ij + Cn 2ij−1 ; j = 2,3,..., n
Anin+1 = Cn2 in
RHSi = sin(i − )
And for i = n + 1 :
Ani1 = Anin+1 = 1
Anii = 0; j = 2,3,..., n
RHSi = 0
Except for i = n + 1 , Anij is called the normal-velocity influence coefficients representing the influences
of j on the normal velocity at the ith control point.
'
%This loop calculates the coefficients Anij and Atij used in eq. 6.
for i=1:n
an(i,1)=cn1(i,1);
an(i,n+1)=cn2(i,n);
for j=2:n
an(i,j)=cn1(i,j)+cn2(i,(j-1));
end
end
%kutta condition
an(n+1,1)=1;
an(n+1,n+1)=1;
rhs(n+1)=0;
for j=2:n
an(n+1,j)=0;
end
% calculating circulation density by solving linear equations given by
% eq. 6
g = an\rhs';
After solving the above equations we obtain the values of circulation densities ' at each panel. We use
these values to compute the velocity and pressure at the control points. At such point, the velocity has
only tangential component at the panel surface because of the vanishing of the normal component
there. If t i designate the unit tangential vector on the ith panel, the local dimensionless velocity
defined as Vi = ( ti ) / V can be computed as
Vi = cos(i − ) + ( Ct1ij 'j + Ct 2ij 'j +1 ) ;
n
i = 1, 2,..., n (7)
j =1
in which
Ct1ij = 0.5CF − DG − Ct 2ij
Ct 2ij = C + 0.5PF / S j + ( AD − CE )G / S j
Ct1ij = Ct 2ij =
2
% This loop calculates the coefficients to find the matrix Aij for equation
% 7.Variables are as given in report.
for i = 1:n
for j = 1:n
if i==j
ct1(i,j) = pi/2;
ct2(i,j) = pi/2;
else
A=-(xmid(i)-xa(j))*cos(phi(j))-(ymid(i)-ya(j))*sin(phi(j));
B = (xmid(i)-xa(j))^2+(ymid(i)-ya(j))^2;
C = sin(phi(i)-phi(j));
D = cos(phi(i)-phi(j));
E = (xmid(i)-xa(j))*sin(phi(j))-(ymid(i)-ya(j))*cos(phi(j));
F = log(1+(Sj(j)^2+2*A*Sj(j))/B);
G = atan2((E*Sj(j)),(B+A*Sj(j)));
P = (xmid(i)-xa(j))*sin(phi(i)-2*phi(j))+(ymid(i)-ya(j))*cos(phi(i)-
2*phi(j));
Q = (xmid(i)-xa(j))*cos(phi(i)-2*phi(j))-(ymid(i)-ya(j))*sin(phi(i)-
2*phi(j));
ct2(i,j)=C+0.5*P*F/Sj(j)+(A*D-C*E)*G/Sj(j);
ct1(i,j)=0.5*C*F-D*G-ct2(i,j);
end
end
end
Equation (7) can be further written as
n +1
Vi = cos(i − ) + Atij 'j ; i = 1, 2,..., n (8)
j =1
where tangential-velocity influence coefficients are defined as follows:
Ati1 = Ct1i1
Atij = Ct1ij + Ct 2ij−1 ; j = 2,3,..., n
Atin+1 = Ct2 in
%This loop calculates the coefficients Anij and Atij used in
%eq. 8
for i=1:n
an(i,1)=cn1(i,1);
an(i,n+1)=cn2(i,n);
at(i,1)=ct1(i,1);
at(i,n+1)=ct2(i,n);
for j=2:n
an(i,j)=cn1(i,j)+cn2(i,(j-1));
at(i,j)=ct1(i,j)+ct2(i,(j-1));
end
end
The pressure coefficient at the ith control point is
C pi = 1 − Vi 2 (9)
% calculating tangential velocity and coefficent of pressure
for i=1:n
sum=0;
for j=1:n+1;
sum=sum+at(i,j)*g(j);
end
v(i) = (cos(phi(i)-alpha)+sum);
cp(i) = 1-v(i)*v(i);
end
The matlab code to compute coefficient of pressure using linearly varying strength vortex panel method
is given below.
clc
clear all
close all
nacaseries = input('Enter the 4-digit naca series = ');
c = input('Enter the chord length = ');
n = input('Enter the number of nodes = ');
a = input('Enter the angle of attack in degree = ');
s = num2str(nacaseries);
% creating nodes on airfoil
s1 = str2double(s(1));
s2 = str2double(s(2));
s3 = str2double(s(3:4));
m = s1*0.01; p = s2*0.1 ; t = s3*0.01;
end
for i= n:-1:1
theta = (i-n)*2*pi/n;
x = 0.5*c*(1+cos(theta));
if(x/c)<p
yc(n+1-i) = m*c/p^2*(2*p*(x/c)-(x/c)^2);
dydx(n+1-i) = (2*m/p^2)* (p-x/c);
beta(n+1-i) = atan(dydx(n+1-i));
else
yc(n+1-i) = m*c/(1-p)^2 * ((1-2*p)+2*p*(x/c)-(x/c)^2);
dydx(n+1-i) = (2*m/(1-p)^2)* (p-x/c);
beta(n+1-i) = atan(dydx(n+1-i));
end
yt=5*t*c*(0.2969*sqrt(x/c)-0.1260*(x/c)...
-0.3516*(x/c)^2+0.2843*(x/c)^3-0.1036*(x/c)^4);
if(i<(0.5*n+1))
xa(n+1-i)=x - yt*sin(beta(n+1-i));
ya(n+1-i)=yc(n+1-i)+yt*cos(beta(n+1-i));
else
xa(n+1-i)=x + yt*sin(beta(n+1-i));
ya(n+1-i)=yc(n+1-i)-yt*cos(beta(n+1-i));
end
end
xa(n+1)= c;
ya(n+1) = 0;
yc(n+1) = 0; % trailing edge
alpha = a*pi/180;
%This loop calculates the location of the control points
for i = 1:n
xmid(i) = (xa(i)+xa(i+1))/2;
ymid(i) = (ya(i)+ya(i+1))/2;
Sj(i) = sqrt((xa(i+1)-xa(i))^2+(ya(i+1)-ya(i))^2);%array of panel lengths
phi(i) = atan2((ya(i+1)-ya(i)),(xa(i+1)-xa(i)));
rhs(i)= sin(phi(i)-alpha);% RHS of equation 4
end
% This loop calculates the coefficients to find the matrix Aij for equation
% 6.Variables are as given in report.
for i = 1:n
for j = 1:n
if i==j
cn1(i,j) = -1;
cn2(i,j) = 1;
ct1(i,j) = pi/2;
ct2(i,j) = pi/2;
else
A=-(xmid(i)-xa(j))*cos(phi(j))-(ymid(i)-ya(j))*sin(phi(j));
B = (xmid(i)-xa(j))^2+(ymid(i)-ya(j))^2;
C = sin(phi(i)-phi(j));
D = cos(phi(i)-phi(j));
E = (xmid(i)-xa(j))*sin(phi(j))-(ymid(i)-ya(j))*cos(phi(j));
F = log(1+(Sj(j)^2+2*A*Sj(j))/B);
G = atan2((E*Sj(j)),(B+A*Sj(j)));
P = (xmid(i)-xa(j))*sin(phi(i)-2*phi(j))+(ymid(i)-ya(j))*cos(phi(i)-
2*phi(j));
Q = (xmid(i)-xa(j))*cos(phi(i)-2*phi(j))-(ymid(i)-ya(j))*sin(phi(i)-
2*phi(j));
cn2(i,j)=D+0.5*Q*F/Sj(j)-(A*C+D*E)*G/Sj(j);
cn1(i,j)=0.5*D*F+C*G-cn2(i,j);
ct2(i,j)=C+0.5*P*F/Sj(j)+(A*D-C*E)*G/Sj(j);
ct1(i,j)=0.5*C*F-D*G-ct2(i,j);
end
end
end
%This loop calculates the coefficients Anij and Atij used in eq. 6 and
%eq. 8 respectively
for i=1:n
an(i,1)=cn1(i,1);
an(i,n+1)=cn2(i,n);
at(i,1)=ct1(i,1);
at(i,n+1)=ct2(i,n);
for j=2:n
an(i,j)=cn1(i,j)+cn2(i,(j-1));
at(i,j)=ct1(i,j)+ct2(i,(j-1));
end
end
%kutta condition
an(n+1,1)=1;
an(n+1,n+1)=1;
rhs(n+1)=0;
for j=2:n
an(n+1,j)=0;
end
% calculating circulation density by solving linear equations given by
% eq. 6
g = an\rhs';
% calculating tangential velocity and coefficent of pressure
for i=1:n
sum=0;
for j=1:n+1;
sum=sum+at(i,j)*g(j);
end
v(i) = (cos(phi(i)-alpha)+sum);
cp(i) = 1-v(i)*v(i);
end
% compute coefficients of lift and drag
cy = 0.0;
cx = 0.0;
for i = 1:n
dx = xa(i+1)-xa(i);
dy = ya(i+1)-ya(i);
cy = cy-cp(i)*dx;
cx = cx+cp(i)*dy;
end
% print cl and cd on the screen
cl = cy*cos(alpha)-cx*sin(alpha)
cd = cy*sin(alpha)+cx*cos(alpha)
plot(xmid(1:n/2)/c,cp(1:n/2),'-*r')
set(gca,'Ydir','reverse')
hold on
plot(xmid(n/2+1:n)/c,cp(n/2+1:n),'-*b')
hold on
splot(xa,ya,'-k')
xlabel('x/c')
ylabel('cp')
title('cp vs x/c')
Fig.2 Pressure distribution over the surface of NACA0012 airfoil at 90 AOA
