clear;
clc;
syms V c alpha
alpha=2*pi/45; % the attack angle is 8 degree and change to radia;
c=1; % the nondimensional chord length;
n=200; % the half number of numerical panel quantity;
t=0.12*c; % the max thickness of the airfoil
X(1)=c;
for i=1:n/2+1;
X(i)=c*(n/2+1-i)/(n/2);
X(n+2-i)=X(i);
end
% the following is to form camber line for NACA 2412;
% m = maximum camber as fraction of chord,p = chord wise position of maximum
% camber. e.g. NACA 2412 means m = 0.02, p = 0.4 and tmax = 0.12.
m=0.02;
p=0.4;
for i=1:n+1
if X(i)<=p
YC(i)=(m/p^2)*(2*p*X(i)-X(i)^2);
else
YC(i)=(m/(1-p)^2)*((1-2*p)+2*p*X(i)-X(i)^2);
end
end
for i=1:n+1
if i<=n/2+1
Y(i)=-5*t*(0.29690*X(i)^0.5-0.12600*X(i)-0.35160*X(i)^2+0.28430*X(i)^3-0.10150*X(i)^4)+YC(i);
else
Y(i)=+5*t*(0.29690*X(i)^0.5-0.12600*X(i)-0.35160*X(i)^2+0.28430*X(i)^3-0.10150*X(i)^4)+YC(i);
end
end
for i=1:n
x(i)=(X(i+1)+X(i))/2;
y(i)=(Y(i+1)+Y(i))/2;
end
for i=1:n
phi(i)=atan2((Y(i+1)-Y(i)),(X(i+1)-X(i)));
if phi(i)<0
phi(i)=phi(i)+2*pi;
end
end
for i=1:n
for j=1:n
A(i,j)=-(x(i)-X(j))*cos(phi(j))-(y(i)-Y(j))*sin(phi(j));
B(i,j)=(x(i)-X(j))^2+(y(i)-Y(j))^2;
C1(i,j)=sin(phi(i)-phi(j));
C2(i,j)=cos(phi(i)-phi(j));
D(i,j)=(x(i)-X(j))*cos(phi(i))+(y(i)-Y(j))*sin(phi(i));
S(j)=((X(j+1)-X(j))^2+(Y(j+1)-Y(j))^2)^0.5;
E(i,j)=(x(i)-X(j))*sin(phi(j))-(y(i)-Y(j))*cos(phi(j)); % (B(i,j)-A(i,j)^2)^0.5;
end
end
for i=1:n
for j=1:n
F(i,j)=log((S(j)^2+2*A(i,j)*S(j)+B(i,j))/B(i,j));
G(i,j)=atan((S(j)+A(i,j))/E(i,j))-atan(A(i,j)/E(i,j));
J(i,j)=-(C2(i,j)/2)*F(i,j)-C1(i,j)*G(i,j);
K(i,j)=-C2(i,j)+(1/S(j))*(F(i,j)*(A(i,j)*C2(i,j)+D(i,j)/2)+G(i,j)*(A(i,j)*C1(i,j)+E(i,j)*C2(i,j)));
end
end
for i=1:n
for j=1:n
if j==1
JK(i,j)=J(i,j)-K(i,j);
else
JK(i,j)=K(i,(j-1))+J(i,j)-K(i,j);
end
end
JK(i,n+1)=K(i,n);
JK(n+1,1)=1;
JK(n+1,n+1)=1;
end
for i=1:n
RHS(i,1)=2*pi*sin(alpha-phi(i)); % 2*pi*V*sin(alpha-phi(i));
RHS(n+1,1)=0;
end
gamma=(inv(JK)*RHS)';
for i=1:n+1
Cp(i)=1-(gamma(i))^2; % 1-(gamma(i)/V)^2;
vpa(Cp,6);
end
hold on
plot(X,Y);axis equal
xlabel('X');
ylabel('Y');
title('Figure of the NACA 2412 airfoil');
figure
hold off
hold on
plot(X,Cp,'k-+');
xlabel('X');
ylabel('Cp');
title('x~Cp');
hold off
% the following is to calculate the whole circulation Gamma;
for i=1:n
Gamma=sum(((gamma(i)+gamma(i+1))/2)*S(i));
end
% % Validation on an elliptic airfoil
t=10;
k=(t-1)/(t+1);
for i=1:n/2+1
thetal(i)=pi*(n/2+1-i)/(n/2);
thetal(n+2-i)=thetal(i);
v(i)=-2*sin(thetal(i))/((1-k)^2+4*k*sin(thetal(i))^2)^0.5;
Cp1(i)=1-v(i)^2;
end
for i=1:n+1
v(i)=-2*sin(thetal(i))/((1-k)^2+4*k*sin(thetal(i))^2)^0.5;
Cp1(i)=1-v(i)^2;
end
for i=1:n+1
thetau(i)=pi*(i-1)/n;
v(i)=-2*sin(thetau(i))/((1-k)^2+4*k*sin(thetau(i))^2)^0.5;
Cpu(i)=1-v(i)^2;
end
figure
hold on
plot(thetal,Cp1,'R-*');
plot(thetal,v,'k:.');
legend('thetal~Cpl','thetal~velocity');
% plot(thetau,Cpu);
xlabel('theta');
ylabel('Cp');
hold off作者: zhoulin96100193 时间: 2007-4-30 21:17
不错，谢谢分享了作者: 异域 时间: 2007-8-27 17:38
嗯，很好！支持一下！作者: citizenc16 时间: 2016-5-2 13:51
看不懂，但是支持一下楼主！ 作者: lical19368888 时间: 2016-6-19 00:51
支持一下