Question

CONVERT THE FOLLOWING MATLAB CODE FROM SOURCE PANEL METHOD TO VORTEX PANEL METHOD:

clc;clear all;close all;

Vinf=100; % freestream velocity

R=1; % cylinder radius

n=4; % number of panels

alpha=2; % angle of attack

dtheta=2*pi/n;

theta=pi+pi/n:-dtheta:-pi+pi/n;

X=R*cos(theta);

Y=R*sin(theta);

for i=1:n

% angle of flow with tangent of panel

phi(i)=-alpha+atan2((Y(i+1)-Y(i)),(X(i+1)-X(i)));

% angle of flow with normal of panel

beta(i)=phi(i)+pi/2;

x_mid(i)=(X(i+1)+X(i))/2;

y_mid(i)=(Y(i+1)+Y(i))/2;

S(i)=sqrt((Y(i+1)-Y(i))^2+(X(i+1)-X(i))^2);

end

% Source Panel Method

for j=1:n

neighbors(:,j)=[1:j-1 j+1:n];

xi=x_mid(j);

yi=y_mid(j);

for i=1:n-1

m=neighbors(i,j);

Xj=X(m);

Yj=Y(m);

Xj1=X(m+1);

Yj1=Y(m+1);

A=-(xi-Xj)*cos(phi(m))-(yi-Yj)*sin(phi(m));

B=(xi-Xj)^2+(yi-Yj)^2;

C=sin(phi(j)-phi(m));

D=(yi-Yj)*cos(phi(j))-(xi-Xj)*sin(phi(j));

E=sqrt(B-A^2);

Sj=S(m);

I(j,m)=C/2*log((Sj^2+2*A*Sj+B)/B)+(D-A*C)/E*(atan2((Sj+A),E)-atan2(A,E));

J(j,m)=(D-A*C)/2/E*log((Sj^2+2*A*Sj+B)/B)-C*(atan2((Sj+A),E)-atan2(A,E));

end

F(j,1)=Vinf*cos(beta(j));

end

M=I/2/pi+eye(n)/2;

lambda=-inv(M)*F;

V=Vinf*sin(beta)+lambda'/2/pi*J';

Cp=1-(V/Vinf).^2

angles=min(beta):0.01:max(beta);

Cp_exact=1-4*sin(angles).^2;

figure(1);

plot(R*cos(0:0.01:2*pi),R*sin(0:0.01:2*pi),'b',X,Y,'r',x_mid,y_mid,'k^');

axis equal;legend('Exact Shape','Panel approximation','Control Points');

figure(2);

plot(angles,Cp_exact,'b',beta,Cp,'k^');grid;

legend('C_p (exact)', 'C_p (Source Panel Method)');

Answer 1

clc;clear all;close all;
Vinf=30;
R=1;%cylinder radius
n=100;%number of panels
dtheta=2*pi/n;
alfa=0;%angle of attack;
theta=pi+pi/n:-dtheta:-pi+pi/n;%central angle
X=R*cos(theta);
Y=R*sin(theta);
for index=1:n
%angle of flow with tangent of panel
phi(index)=-alfa+...
atan2((Y(index+1)-Y(index)),(X(index+1)-X(index)));
%angle of flow with normal of panel
beta(index)=phi(index)+pi/2;
midpoint_x(index)=(X(index+1)+X(index))/2;
midpoint_y(index)=(Y(index+1)+Y(index))/2;
S(index)=sqrt((Y(index+1)-Y(index))^2+...
(X(index+1)-X(index))^2);%length of panel
end
%The Source Panel Method
for p=1:n
neighbors(:,p)=[1:p-1 p+1:n];
xi=midpoint_x(p);
yi=midpoint_y(p);
for index=1:n-1
m=neighbors(index,p);
Xj=X(m);
Yj=Y(m);
Xj1=X(m+1);
Yj1=Y(m+1);
A=-(xi-Xj)*cos(phi(m))-(yi-Yj)*sin(phi(m));
B=(xi-Xj)^2+(yi-Yj)^2;
C=sin(phi(p)-phi(m));
D=(yi-Yj)*cos(phi(p))-(xi-Xj)*sin(phi(p));
E=sqrt(B-A^2);
Sj=S(m);
I(p,m)=C/2*log((Sj^2+2*A*Sj+B)/B)+...
(D-A*C)/E*(atan2((Sj+A),E)-atan2(A,E));
J(p,m)=(D-A*C)/2/E*log((Sj^2+2*A*Sj+B)/B)...
-C*(atan2((Sj+A),E)-atan2(A,E));
end
F(p,1)=Vinf*cos(beta(p));
end

M=I/2/pi+eye(n)/2;

lambda=-inv(M)*F;
fprintf('The sum of all sources is %f',lambda'*S');%check sum

%Recoving velocity at the nodes
V=Vinf*sin(beta)+lambda'/2/pi*J';

Cp=1-(V/Vinf).^2;
angles=min(beta):0.01:max(beta);
Cp_exact=1-4*sin(angles).^2;
subplot(1,2,1);plot(R*cos(0:0.01:2*pi),R*sin(0:0.01:2*pi),'b',...
X,Y,'r',midpoint_x,midpoint_y,'g^');axis equal;
legend('Exact Shape','Panel approximation','Control Points')
subplot(1,2,2);plot(angles,Cp_exact,'b',beta,Cp,'r^');grid;
legend('C_p (exact', 'C_p (Source Panel Method)');