Joukowski Airfoil

One of the more important potential flow results obtained using conformal mapping are the solutions of the potential flows past a family of airfoil shapes known as Joukowski foils. Like some of the other solutions presented here, we begin with a known solution, namely the flow past a circular cylinder with circulationand use conformal mapping to transform the cylinder into an airfoil shape.

Download Ebook 
close all;
clear all;
clc;
%% inputs ::
row=1.225;
p_inf=1.023e5;
chord=1;
v_inf=100;
max_thickness=.12;
max_camber=.06;
alpha=[-4,-2,0,2,4,6,8,10,12,14,16,18,20,10]*pi/180;
b=chord/4;
e=max_thickness/1.3;
beta=2*max_camber;
a=b*(1+e)/cos(beta);
%% program ::
%% airfoil coordinates ::
theta_=0:pi/100:2*pi;
r_=a;
x_=r_*cos(theta_);
y_=r_*sin(theta_);
x=x_-b*e;
y=y_+a*beta;
x1_ex=x.*(1+b^2./(x.^2+y.^2));
y1_ex=y.*(1-b^2./(x.^2+y.^2));
x1_app=2*b*cos(theta_);
y1_app=2*b*e*(1-cos(theta_)).*sin(theta_)+2*b*beta*sin(theta_).^2;
figure;
whitebg;
plot(x1_ex,y1_ex,'g-','LineWidth',2);
hold on;
plot(x1_app,y1_app,'r-','LineWidth',2);
xlabel('$X_{B}$','interpreter','latex','FontSize',14);
ylabel('$Y_{B}$','interpreter','latex','FontSize',14);
title('$Airfoil$ $shape$ $approximate$ $,$ $exact$','interpreter','latex','FontSize',14);
legend('Exact','Approximate');
grid;
axis equal
hold off;
%% calculations of v and cp ::
clear theta_ r_ r
theta_=0:pi/100:2*pi;
r_=a;
theta=atan2(r_*sin(theta_)+a*beta,r_*cos(theta_)-b*e);
r=b*(1+e*(1-cos(theta))+beta*sin(theta));
i=1;     %index of alpha=-4
vt_=-v_inf*(2*sin(theta_-alpha(i))+2*sin(alpha(i)+beta));
V11=sqrt(vt_.^2./(1-2*(b^2./r.^2).*cos(2*theta)+(b^4./r.^4)));
figure;
plot(x1_app,y1_app,'g-','LineWidth',2);
hold on;
plot(x1_app,V11/v_inf,'r-','LineWidth',2);
%------------------------------%
i=6;     %index of alpha=6
vt_=-v_inf*(2*sin(theta_-alpha(i))+2*sin(alpha(i)+beta));
V16=sqrt(vt_.^2./(1-2*(b^2./r.^2).*cos(2*theta)+(b^4./r.^4)));
plot(x1_app,V16/v_inf,'b-','LineWidth',2);
%------------------------------%
i=8;     %index of alpha=10
vt_=-v_inf*(2*sin(theta_-alpha(i))+2*sin(alpha(i)+beta));
V18=sqrt(vt_.^2./(1-2*(b^2./r.^2).*cos(2*theta)+(b^4./r.^4)));
plot(x1_app,V18/v_inf,'y-','LineWidth',2);
%-------------------------------%
xlabel('$X_{B}$','interpreter','latex','FontSize',14);
ylabel('$V/V_{inf}$','interpreter','latex','FontSize',14);
title('$Velocity$ $over$ $airfoil$','interpreter','latex','FontSize',14);
legend('Airfoil','alpha=-4','alpha=6','alpha=10');
grid;
axis square
hold off;
%----------------------------------%
cp1=1-(V11/v_inf).^2;
cp6=1-(V16/v_inf).^2;
cp8=1-(V18/v_inf).^2;
figure;
plot(x1_app,y1_app,'g-','LineWidth',2);
hold on;
plot(x1_app,cp1,'r-','LineWidth',2);
plot(x1_app,cp6,'b-','LineWidth',2);
plot(x1_app,cp8,'y-','LineWidth',2);
xlabel('$X_{B}$','interpreter','latex','FontSize',14);
ylabel('$C_{p}$','interpreter','latex','FontSize',14);
title('$Pressure$ $over$ $airfoil$','interpreter','latex','FontSize',14);
legend('Airfoil','alpha=-4','alpha=6','alpha=10');
grid;
axis square
hold off;
%% calculation of CL CM Xc.p ::
% we will calculate cm numerically because of it hard to get analytical
% form for it.
alpha_=[-10:.1:10]*pi/180;
P=zeros(1,size(theta_,2));
Cm=zeros(1,size(alpha_,2));
Cl=zeros(1,size(alpha_,2));
Cd=zeros(1,size(alpha_,2));
Xcp=zeros(1,size(alpha_,2));
Xac=zeros(1,size(alpha_,2));
Xcp_alpha=zeros(1,size(alpha_,2));
cl_alpha=zeros(1,size(alpha_,2));
clear theta_ r_
theta_=0:pi/100:2*pi;
r_=a;
theta=atan2(r_*sin(theta_)+a*beta,r_*cos(theta_)-b*e);
r=b*(1+e*(1-cos(theta))+beta*sin(theta));
for i=1:size(alpha_,2)
    vt_=-v_inf*(2*sin(theta_-alpha_(i))+2*sin(alpha_(i)+beta));
    V=sqrt(vt_.^2./(1-2*(b^2./r.^2).*cos(2*theta)+(b^4./r.^4)));
    P=p_inf+.5*row*(v_inf.^2-V.^2);
    M=0;
    Ly=0;
    Lx=0;
   for j=1:size(theta_,2)-1 
    M=M-(P(j)+P(j+1))/2*((x1_app(j)+x1_app(j+1))/2*(x1_app(j+1)-x1_app(j))-...
        (y1_app(j)+y1_app(j+1))/2*(y1_app(j+1)-y1_app(j))*0);
    Ly=Ly+(P(j)+P(j+1))/2*(x1_app(j+1)-x1_app(j));
    Lx=Ly+(P(j)+P(j+1))/2*(y1_app(j+1)-y1_app(j));
   end
   D=Ly*sin(alpha_(i))+Lx*cos(alpha_(i));
   L=Ly*cos(alpha_(i))-Lx*sin(alpha_(i));
   Cm(i)=M/(.5*row*v_inf^2*chord^2);
   Cl(i)=L/(.5*row*v_inf^2*chord);
   Cd(i)=D/(.5*row*v_inf^2*chord);
   Xcp(i)=-M/Ly;
end
Mslope=(Cm(135)-Cm(120))/(alpha_(135)-alpha_(120))
for i=1:size(alpha_,2)-1
    cl_alpha(i)=(Cl(i+1)-Cl(i))/(alpha_(i+1)-alpha_(i));
    Xcp_alpha(i)=(Xcp(i+1)-Xcp(i))/(alpha_(i+1)-alpha_(i));
    Xac(i)=Xcp(i)+Cl(i)/cl_alpha(i)*Xcp_alpha(i);
end
figure;
plot(alpha_,Cd,'g-','LineWidth',2);
xlabel('$alpha$','interpreter','latex','FontSize',14);
ylabel('$C_{D}$','interpreter','latex','FontSize',14);
title('$C_{D}$','interpreter','latex','FontSize',14);
grid on;
%---------------%
figure;
plot(alpha_,Cl,'r-','LineWidth',2);
hold on;
Cl=2*pi*(1+e)*sin(alpha_+beta);
plot(alpha_,Cl,'g-','LineWidth',2);
xlabel('$alpha$','interpreter','latex','FontSize',14);
ylabel('$C_{L}$ $Exact$ $,$ $approximate$','interpreter','latex','FontSize',14);
title('$C_{L}$','interpreter','latex','FontSize',14);
legend('numerical','analytical');
grid on;
%---------------%
figure;
plot(alpha_,Cm,'r-','LineWidth',2);
xlabel('$alpha$','interpreter','latex','FontSize',14);
ylabel('$C_{m}$','interpreter','latex','FontSize',14);
title('$C_{m}$ $over$ $airfoil$ $using$ $numerical$ $method$','interpreter','latex','FontSize',14);
grid on;
figure;
plot(alpha_,Xcp,'r-','LineWidth',2);
xlabel('$alpha$','interpreter','latex','FontSize',14);
ylabel('$X_{c.p}$','interpreter','latex','FontSize',14);
title('$X_{c.p}$ $variation$ $with$ $alpha$','interpreter','latex','FontSize',14);
grid on;
figure;
plot(alpha_,Xac,'r-','LineWidth',2);
xlabel('$alpha$','interpreter','latex','FontSize',14);
ylabel('$X_{a.c}$','interpreter','latex','FontSize',14);
title('$X_{a.c}$ $variation$ $with$ $alpha$','interpreter','latex','FontSize',14);
grid on;
%%  draw of streeamlines ::
%%  draw of streeamlines ::
clear theta_ r_ theta vr_ vt_ x_ y_ x y x1 y1 r r_ alpha_
i=14;
theta_=0:pi/100:2*pi;
r_=a:a/80:6*a; 
[theta_,r_]=meshgrid(theta_,r_);
alpha__=[0:1:360]*pi/180;
figure('Renderer','zbuffer');
axis tight manual
set(gca,'NextPlot','replaceChildren');
ved= VideoWriter('airfoil');
ved.Quality=100;
ved.FrameRate=10;
open(ved);
for i=1:size(alpha__,2)
theta=atan2(r_.*sin(theta_)+a*beta,r_.*cos(theta_-alpha__(i))-b*e);
r=b*(1+e*(1-cos(theta))+beta*sin(theta));
x_=r_.*cos(theta_);
y_=r_.*sin(theta_);
x=x_-b*e;
y=y_+a*beta;
x1=x.*(1+b^2./(x.^2+y.^2));
y1=y.*(1-b^2./(x.^2+y.^2));
epsi=v_inf*(r_.*sin(theta_-alpha__(i))+a^2./r_.*sin(alpha__(i)-theta_)+2*a*sin(alpha__(i)+beta).*log(r_));
%--------------transform--------------%
xx1=x1*cos(alpha__(i))+y1*sin(alpha__(i));
yy1=-x1*sin(alpha__(i))+y1*cos(alpha__(i));
colordef black
contour(xx1,yy1,epsi,50,'w','LineWidth',1);
hold on;
x1_exx=x1_ex*cos(alpha__(i))+y1_ex*sin(alpha__(i));
y1_exx=-x1_ex*sin(alpha__(i))+y1_ex*cos(alpha__(i));
plot(x1_exx,y1_exx,'c-','LineWidth',2);
axis([-1 1 -1 1]);
hold off;
%---------------%
F(i) = getframe;
writeVideo(ved, F(i));
end
close(ved);
%__________________________________%
% vr_=v_inf*(1-a^2./r_.^2).*cos(theta_-alpha(i));
% vt_=-v_inf*(sin(theta_-alpha(i)).*(1+a^2./r_.^2)+2*(a./r_).*sin(alpha(i)+beta));
% A=vr_.*cos(theta_)-vt_.*sin(theta_);
% B=-(vr_.*sin(theta_)+vt_.*cos(theta_));
% C=1-(b^2./r_.^2).*cos(2*theta);
% D=(b^2./r_.^2).*sin(2*theta);
% u1=(A.*C+B.*D)./(C.^2+D.^2);
% v1=(D.*A-B.*C)./(C.^2+D.^2);
% ystart=linspace(-1,1,50);
% xstart=-ones(1,50);
% figure;
% x11=x1*cos(alpha(i))+y1*sin(alpha(i));
% y11=-x1*sin(alpha(i))+y1*cos(alpha(i));
% quiver(x11,y11,u1,v1,2,'-');
% hold on;
% axis([-1 1 -1 1]);
% x1_exx=x1_ex*cos(alpha(i))+y1_ex*sin(alpha(i));
% y1_exx=-x1_ex*sin(alpha(i))+y1_ex*cos(alpha(i));
% plot(x1_exx,y1_exx,'b-','LineWidth',2);
% %streamslice(x1,y1,u1,v1,xstart,ystart,'w');

Leave a Reply

Your email address will not be published. Required fields are marked *