fluid 2 mechanical engineering homework
ymalharbi
airfoil_codes.zipairfoil codes/CalculateMAC.m
function [MAC,yMAC]=CalculateMAC(y,c) %-------------------------------------------------------------------------- %CalculateMAC %Version 1.00 %Created by Stepen %Created 25 March 2011 %-------------------------------------------------------------------------- %CalculateMAC calculates the mean aerodynamic chord length of a wing based %on its spanwise chord length distribution. %-------------------------------------------------------------------------- %Syntax: %MAC=CalculateMAC(y,c) %Input argument: %- y (m x 1 num) specifies the spanwise location of all wing stations. %- c (m x 1 num) specifies the chord length of all wing stations. %Output argument: %- MAC (1 x 1 num) specifies the mean aerodynamic chord length of wing with % given chord length distribution. %- yMAC (p x 1 num) specifies the spanwise location of mean aerodynamic % chord. %-------------------------------------------------------------------------- %Example: %[MAC,yMAC]=CalculateMAC([0,10],[3,1]) % calculates MAC and MAC's location of a straight tapered wing. %[MAC,yMAC]=CalculateMAC([0,4],[6,0]) % calculates MAC and MAC's location of a delta wing. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input camber line vertexes size_y=size(y); size_c=size(c); if (size_y(1)~=size_c(1))||(size_y(1)~=size_c(1))||... (numel(y)~=numel(c)) error('Input y and c must have the same size!') end planformcount=numel(y)-1; %Calculating wing area wing_area=0; for i=1:planformcount wing_area=wing_area+(0.5*(y(i+1)-y(i))*(c(i+1)+c(i))); end %Calculating MAC length MAC=0; for i=1:planformcount m=(c(i+1)-c(i))/(y(i+1)-y(i)); c0=c(i); dy=y(i+1)-y(i); MAC=MAC+((1/wing_area)*(((m^2)*(dy^3)/3)+... (m*c0*(dy^2))+... ((c0^2)*dy))); clear m c0 dy end %Finding MAC location(s) yMACcount=0; for i=1:planformcount if c(i)~=c(i+1) if (MAC<max([c(i),c(i+1)]))&&(MAC>min([c(i),c(i+1)])) yMACcount=yMACcount+1; end elseif c(i)==c(i+1) if c(i)==MAC yMACcount=yMACcount+1; end end end yMAC=zeros(yMACcount,1); counter=0; for i=1:planformcount if c(i)~=c(i+1) if (MAC<max([c(i),c(i+1)]))&&(MAC>min([c(i),c(i+1)])) counter=counter+1; yMAC(counter)=interp1([c(i),c(i+1)],[y(i),y(i+1)],MAC); end elseif c(i)==c(i+1) if c(i)==MAC counter=counter+1; yMAC(counter)=0.5*(y(i)+y(i+1)); end end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromDoublet2D (1).m
function [u,v,un,vn]=CalculateVelocityFromDoublet2D(s,xc,yc,x,y,dtheta) %-------------------------------------------------------------------------- %CalculateVelocityFromDoublet2D %Version 1.00 %Created by Stepen %Created 8 January 2011 %Last modified 2 August 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromDoublet2D calculates the flow velocity influenced by %several given doublet flow at several given locations in 2D domain. %CalculateVelocityFromDoublet2D assumes that the doublets that influence %the flow are point doublets. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromDoublet2D(s,xc,yc,x,y) %Input argument: %- s (m x 1 num) specifies the strength of all point doublets. %- xc (m x 1 num) specifies the x axis location of all point doublets. %- yc (m x 1 num) specifies the y axis location of all point doublets. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %- dtheta (m x 1) specifies the rotation of point doublets. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input point doublets s_size=size(s); xc_size=size(xc); yc_size=size(yc); dtheta_size=size(dtheta); if (s_size(2)~=1)||... (xc_size(2)~=1)||(yc_size(2)~=1)||... (dtheta_size(2)~=1) error('Input point doublets must be a i x 1 array!') end if (numel(s)~=numel(xc))||(numel(s)~=numel(yc))||... (numel(s)~=numel(dtheta)) error('Input s, xc, yc, and dtheta do not have the same size!') end doubletcount=numel(s); %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Preallocating array for speed u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount for k=1:1:doubletcount if (x(i,j)==xc(k))&&(y(i,j)==yc(k)) u(i,j)=u(i,j); v(i,j)=v(i,j); else r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]); theta=atan2(y(i,j)-yc(k),x(i,j)-xc(k))+deg2rad(dtheta); sintheta=sin(theta); costheta=cos(theta); vr=(1/(2*pi()*r*r))*s(k)*costheta; vtheta=-(1/(2*pi()*r*r))*s(k)*sintheta; u(i,j)=u(i,j)+((vr*((x(i,j)-xc(k))/r))+... (vtheta*((y(i,j)-yc(k))/r))); v(i,j)=v(i,j)+((vr*((y(i,j)-yc(k))/r))+... (vtheta*((xc(k)-x(i,j))/r))); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromDoublet2D.m
function [u,v,un,vn]=CalculateVelocityFromDoublet2D(s,xc,yc,x,y,dtheta) %-------------------------------------------------------------------------- %CalculateVelocityFromDoublet2D %Version 1.00 %Created by Stepen %Created 8 January 2011 %Last modified 2 August 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromDoublet2D calculates the flow velocity influenced by %several given doublet flow at several given locations in 2D domain. %CalculateVelocityFromDoublet2D assumes that the doublets that influence %the flow are point doublets. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromDoublet2D(s,xc,yc,x,y) %Input argument: %- s (m x 1 num) specifies the strength of all point doublets. %- xc (m x 1 num) specifies the x axis location of all point doublets. %- yc (m x 1 num) specifies the y axis location of all point doublets. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %- dtheta (m x 1) specifies the rotation of point doublets. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input point doublets s_size=size(s); xc_size=size(xc); yc_size=size(yc); dtheta_size=size(dtheta); if (s_size(2)~=1)||... (xc_size(2)~=1)||(yc_size(2)~=1)||... (dtheta_size(2)~=1) error('Input point doublets must be a i x 1 array!') end if (numel(s)~=numel(xc))||(numel(s)~=numel(yc))||... (numel(s)~=numel(dtheta)) error('Input s, xc, yc, and dtheta do not have the same size!') end doubletcount=numel(s); %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Preallocating array for speed u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount for k=1:1:doubletcount if (x(i,j)==xc(k))&&(y(i,j)==yc(k)) u(i,j)=u(i,j); v(i,j)=v(i,j); else r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]); theta=atan2(y(i,j)-yc(k),x(i,j)-xc(k))+(pi/180).*dtheta; sintheta=sin(theta); costheta=cos(theta); vr=(1/(2*pi()*r*r))*s(k)*costheta; vtheta=-(1/(2*pi()*r*r))*s(k)*sintheta; u(i,j)=u(i,j)+((vr*((x(i,j)-xc(k))/r))+... (vtheta*((y(i,j)-yc(k))/r))); v(i,j)=v(i,j)+((vr*((y(i,j)-yc(k))/r))+... (vtheta*((xc(k)-x(i,j))/r))); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromSource2D (1).m
function [u,v,un,vn]=CalculateVelocityFromSource2D(s,xc,yc,x,y) %-------------------------------------------------------------------------- %CalculateVelocityFromSource2D %Version 1.00 %Created by Stepen %Created 8 January 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromSource2D calculates the flow velocity influenced by %several given source flow at several given locations in 2D domain. %CalculateVelocityFromSource2D assumes that the sources that influence the %flow are point sources. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromSource2D(s,xc,yc,x,y) %Input argument: %- s (m x 1 num) specifies the strength of all point sources. %- xc (m x 1 num) specifies the x axis location of all point sources. %- yc (m x 1 num) specifies the y axis location of all point sources. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input point sources s_size=size(s); xc_size=size(xc); yc_size=size(yc); if (s_size(2)~=1)||(xc_size(2)~=1)||(yc_size(2)~=1) error('Input point sources must be a i x 1 array!') end if (numel(s)~=numel(xc))||(numel(s)~=numel(yc)) error('Input s, xc, and yc do not have the same size!') end sourcecount=numel(s); %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Preallocating array for speed u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount for k=1:1:sourcecount if (x(i,j)==xc(k))&&(y(i,j)==yc(k)) u(i,j)=u(i,j); v(i,j)=v(i,j); else r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]); vr=s(k)/(2*pi()*r); u(i,j)=u(i,j)+(vr*((x(i,j)-xc(k))/r)); v(i,j)=v(i,j)+(vr*((y(i,j)-yc(k))/r)); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromSource2D (2).m
function [u,v,un,vn]=CalculateVelocityFromSource2D(s,xc,yc,x,y) %-------------------------------------------------------------------------- %CalculateVelocityFromSource2D %Version 1.00 %Created by Stepen %Created 8 January 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromSource2D calculates the flow velocity influenced by %several given source flow at several given locations in 2D domain. %CalculateVelocityFromSource2D assumes that the sources that influence the %flow are point sources. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromSource2D(s,xc,yc,x,y) %Input argument: %- s (m x 1 num) specifies the strength of all point sources. %- xc (m x 1 num) specifies the x axis location of all point sources. %- yc (m x 1 num) specifies the y axis location of all point sources. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input point sources s_size=size(s); xc_size=size(xc); yc_size=size(yc); if (s_size(2)~=1)||(xc_size(2)~=1)||(yc_size(2)~=1) error('Input point sources must be a i x 1 array!') end if (numel(s)~=numel(xc))||(numel(s)~=numel(yc)) error('Input s, xc, and yc do not have the same size!') end sourcecount=numel(s); %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Preallocating array for speed u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount for k=1:1:sourcecount if (x(i,j)==xc(k))&&(y(i,j)==yc(k)) u(i,j)=u(i,j); v(i,j)=v(i,j); else r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]); vr=s(k)/(2*pi()*r); u(i,j)=u(i,j)+(vr*((x(i,j)-xc(k))/r)); v(i,j)=v(i,j)+(vr*((y(i,j)-yc(k))/r)); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromSource2D.m
function [u,v,un,vn]=CalculateVelocityFromSource2D(s,xc,yc,x,y) %-------------------------------------------------------------------------- %CalculateVelocityFromSource2D %Version 1.00 %Created by Stepen %Created 8 January 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromSource2D calculates the flow velocity influenced by %several given source flow at several given locations in 2D domain. %CalculateVelocityFromSource2D assumes that the sources that influence the %flow are point sources. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromSource2D(s,xc,yc,x,y) %Input argument: %- s (m x 1 num) specifies the strength of all point sources. %- xc (m x 1 num) specifies the x axis location of all point sources. %- yc (m x 1 num) specifies the y axis location of all point sources. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input point sources s_size=size(s); xc_size=size(xc); yc_size=size(yc); if (s_size(2)~=1)||(xc_size(2)~=1)||(yc_size(2)~=1) error('Input point sources must be a i x 1 array!') end if (numel(s)~=numel(xc))||(numel(s)~=numel(yc)) error('Input s, xc, and yc do not have the same size!') end sourcecount=numel(s); %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Preallocating array for speed u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount for k=1:1:sourcecount if (x(i,j)==xc(k))&&(y(i,j)==yc(k)) u(i,j)=u(i,j); v(i,j)=v(i,j); else r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]); vr=s(k)/(2*pi()*r); u(i,j)=u(i,j)+(vr*((x(i,j)-xc(k))/r)); v(i,j)=v(i,j)+(vr*((y(i,j)-yc(k))/r)); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromSourcePanel2D.m
function [u,v,un,vn]=CalculateVelocityFromSourcePanel2D(s,xp,yp,x,y,dotcount) %-------------------------------------------------------------------------- %CalculateVelocityFromSourcePanel2D %Version 1.00 %Created by Stepen %Created 9 January 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromSourcePanel2D calculates the flow velocity influenced %by a panel of uniformly distributed source at several given locations in %2D domain. CalculateVelocityFromSourcePanel2D requires %CalculateVelocityFromSource2D m-file function in the same directories to %run. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromSourcePanel2D(s,xp,yp,x,y) %Input argument: %- s (1 x 1 num) specifies the strength of source of the panel per unit % length. %- xp (2 x 1 num) specifies the x axis location of the panel's edge. %- yp (2 x 1 num) specifies the y axis location of the panel's edge. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %- dotcount (1 x 1 num) specifies the number of vertexes to be generated on % each the panel. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input panel xp_size=size(xp); yp_size=size(yp); if numel(s)~=1 error('Input source strength must be a scalar!') end if (xp_size(1)~=2)||(yp_size(1)~=2)||(xp_size(2)~=1)||(yp_size(2)~=1) error('Input panel edges location must be a 2 x 1 array!') end %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Checking input dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Preallocating array for speed xg=zeros(dotcount-1,1); yg=zeros(dotcount-1,1); xcp=zeros(dotcount-1,1); ycp=zeros(dotcount-1,1); l=zeros(dotcount-1,1); u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating panel properties length=norm([xp(2)-xp(1),yp(2)-yp(1)]); sine=(yp(2)-yp(1))/length; cosine=(xp(2)-xp(1))/length; nsine=cosine; ncosine=-sine; %Defining panel grids and its control point dx=(xp(2)-xp(1))/(dotcount-1); dy=(yp(2)-yp(1))/(dotcount-1); for i=1:1:dotcount xg(i)=xp(1)+((i-1)*dx); yg(i)=yp(1)+((i-1)*dy); end for i=1:1:dotcount-1 xcp(i)=0.5*(xg(i+1)+xg(i)); ycp(i)=0.5*(yg(i+1)+yg(i)); l(i)=norm([(xg(i+1)-xg(i)),(yg(i+1)-yg(i))]); end %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount if (x(i,j)==(0.5*(xp(2)+xp(1))))&&(y(i,j)==(0.5*(yp(2)+yp(1)))) vn=0.5*s; u(i,j)=vn*ncosine; v(i,j)=vn*nsine; else for k=1:1:dotcount-1 [du,dv]=CalculateVelocityFromSource2D(s,xcp(k),ycp(k),x(i,j),y(i,j)); u(i,j)=u(i,j)+(du*l(k)); v(i,j)=v(i,j)+(dv*l(k)); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/CalculateVelocityFromVortex2D.m
function [u,v,un,vn]=CalculateVelocityFromVortex2D(s,xc,yc,x,y) %-------------------------------------------------------------------------- %CalculateVelocityFromVortex2D %Version 1.00 %Created by Stepen %Created 8 January 2011 %-------------------------------------------------------------------------- %CalculateVelocityFromVortex2D calculates the flow velocity influenced by %several given vortex flow at several given locations in 2D domain. %CalculateVelocityFromVortex2D assumes that the vortexs that influence the %flow are point vortexes. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn]=CalculateVelocityFromvortex2D(s,xc,yc,x,y) %Input argument: %- s (m x 1 num) specifies the strength of all point vortexes. %- xc (m x 1 num) specifies the x axis location of all point vortexes. %- yc (m x 1 num) specifies the y axis location of all point vortexes. %- x (p x q num) specifies the x axis location of the corresponding % velocity vector. %- y (p x q num) specifies the y axis location of the corresponding % velocity vector. %Output argument: %- u (m x n) specifies the horizontal velocity component at the % corresponding location. %- u (m x n) specifies the vertical velocity component at the corresponding % location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input point vortexes s_size=size(s); xc_size=size(xc); yc_size=size(yc); if (s_size(2)~=1)||(xc_size(2)~=1)||(yc_size(2)~=1) error('Input point vortexs must be a i x 1 array!') end if (numel(s)~=numel(xc))||(numel(s)~=numel(yc)) error('Input s, xc, and yc do not have the same size!') end vortexcount=numel(s); %Checking input location x_size=size(x); y_size=size(y); if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2)) error('Input x and y do not have the same size!') end rowcount=x_size(1); colcount=x_size(2); %Preallocating array for speed u=zeros(rowcount,colcount); v=zeros(rowcount,colcount); un=zeros(rowcount,colcount); vn=zeros(rowcount,colcount); %Calculating influenced velocity for i=1:1:rowcount for j=1:1:colcount for k=1:1:vortexcount if (x(i,j)==xc(k))&&(y(i,j)==yc(k)) u(i,j)=u(i,j); v(i,j)=v(i,j); else r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]); vtheta=-s(k)/(2*pi()*r); u(i,j)=u(i,j)+(vtheta*((y(i,j)-yc(k))/r)); v(i,j)=v(i,j)+(vtheta*((xc(k)-x(i,j))/r)); end end end end %Calculating normalized vector for i=1:1:rowcount for j=1:1:colcount V=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V; vn(i,j)=v(i,j)/V; end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/ClassicThinAirfoilAnalysis.m
function [Cl,Cm_ac,x_cp]=ClassicThinAirfoilAnalysis(xc,yc,AoA,IntMod) %-------------------------------------------------------------------------- %ClassicThinAirfoilAnalysis %Version 1.00 %Created by Stepen %Created 6 February 2011 %-------------------------------------------------------------------------- %ClassicThinAirfoilAnalysis calculates aerodynamic characteristic of an %airfoil with given camber line by distributing vortex sheet at airfoil's %camber line and and adjusting its strength so that airfoil's camber line %is a streamline of the flow. ClassicThinAirfoilAnalysis results only valid %for potential flow cases (incompressible, inviscid, and irrotational %flow). %-------------------------------------------------------------------------- %Syntax: %[Cl,Cm_LE,x_cp]=ClassicThinAirfoilAnalysis(xc,yc) %Input argument: %- xc (m x 1 num) specifies the x axis location of the airfoil's camber % line vertexes in fraction of chord. %- yc (m x 1 num) specifies the y axis location of the airfoil's camber % line vertexes in fraction of chord. %- AoA (n x 1 num) specifies the evaluation angle of attack of the airfoil % in degree. %- IntMod (str) specifies the integral evaluation method. Enter 'theta' to % use constant theta interval for integration or 'xc' to use constant xc % interval for integration. %Output argument: %- Cl (1 x 1 num) specifies the airfoil's lift coefficient at given % evaluation angles of attack. %- Cm_ac (1 x 1 num) specifies the airfoils' pitching moment coeficient % around airfoil's aerodynamic center (theoritically at 0.25 chord). %- x_cp (1 x 1 num) specifies the airfoil's center of pressure location % from leading-edge in fraction of chord at given evaluation angles of % attack. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input camber line vertexes size_x=size(xc); size_y=size(yc); if (size_x(1)~=size_y(1))||(size_x(1)~=size_y(1))||... (numel(xc)~=numel(yc)) error('Input xc and yc must have the same size!') end %Checking input IntMod if nargin==3 IntMod='theta'; end if (~strcmpi(IntMod,'theta'))&&(~strcmpi(IntMod,'xc')) error('Integration method input must be theta or xc!') end %Defining integration size and interval if strcmpi(IntMod,'theta') theta=0:pi/1000:pi; dotcount=numel(theta); elseif strcmpi(IntMod,'xc') theta=acos(1-(2*xc)); dotcount=numel(xc); end panelcount=dotcount-1; %Preallocating array for speed theta_panel=zeros(1,panelcount); dydx=zeros(1,panelcount); Cl=zeros(1,numel(AoA)); x_cp=zeros(1,numel(AoA)); %Calculating theta_panel and dy/dx at each panel if strcmpi(IntMod,'theta') xtemp=0.5*(1-cos(theta)); ytemp=interp1(xc,yc,xtemp); xc=xtemp; yc=ytemp; end for i=1:1:panelcount theta_panel(i)=acos(1-(2*(mean([xc(i) xc(i+1)])))); dydx(i)=(yc(i+1)-yc(i))/(xc(i+1)-xc(i)); end %Calculating A0, A1, A2, and Cl_0 from integration A0=0; A1=0; A2=0; Cl_0=0; for i=1:1:panelcount dtheta=theta(i+1)-theta(i); A0=A0-((1/pi)*dydx(i)*dtheta); A1=A1+((2/pi)*dydx(i)*cos(theta_panel(i))*dtheta); A2=A2+((2/pi)*dydx(i)*cos(2*theta_panel(i))*dtheta); Cl_0=Cl_0+(2*dydx(i)*(cos(theta_panel(i))-1)*dtheta); end %Calculating lift coefficient and center of pressure location for i=1:1:numel(AoA) Cl(i)=Cl_0+((pi/180).*(AoA(i))*2*pi); x_cp(i)=0.25*(1+((pi/Cl(i))*(A1-A2))); end %Calculating pitching moment coefficient Cm_ac=(pi/4)*(A2-A1); %CodeEnd------------------------------------------------------------------- end
airfoil codes/GenerateAirfoil.m
function [xu,yu,xl,yl]=GenerateAirfoil(xc,yc,yt) %-------------------------------------------------------------------------- %GenerateAirfoil %Version 1.10 %Created by Stepen %Created 3 December 2010 %Last modified 17 December 2010 %-------------------------------------------------------------------------- %GenerateAirfoil generates the airfoil coordinates from the given camber %line and thickness distribution. GenerateAirfoil generates airfoil's %vertexes by assuming that airfoil thickness at each chordwise loation is %measured perpendicular to the camber line. %-------------------------------------------------------------------------- %Syntax: %[xu,yu,xl,yl]=GenerateAirfoil(xc,yc,yt) %Input argument: %- xc (m x 1 num) specifies the x axis location of the airfoil's camber % line vertexes. %- yc (m x 1 num) specifies the y axis location of the airfoil's camber % line vertexes. %- yt (m x 1 num) specifies the airfoil's thickness at the corresponding % camber line vertexes. %Output argument: %- xu (i x 1 num) specifies the x axis location of airfoil's upper surface % vertexes in fraction of chord. The airfoil's upper surface vertex are % arranged from leading edge (the first element of xu) to the trailing % edge (the last element of xu). %- yu (i x 1 num) specifies the y axis location of airfoil's upper surface % vertexes in fraction of chord. The airfoil's upper surface vertex are % arranged from leading edge (the first element of yu) to the trailing % edge (the last element of yu). %- xl (i x 1 num) specifies the x axis location of airfoil's lower surface % vertexes in fraction of chord. The airfoil's lower surface vertex are % arranged from leading edge (the first element of xl) to the trailing % edge (the last element of xl). %- yl (i x 1 num) specifies the y axis location of airfoil's lower surface % vertexes in fraction of chord. The airfoil's lower surface vertex are % arranged from leading edge (the first element of yl) to the trailing % edge (the last element of yl). %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input xc size_xc=size(xc); if ~iscolumn(xc) error('Input xc must be one column array!') end %Checking input yc if ~iscolumn(yc) error('Input yc must be one column array!') end %Checking input yt if ~iscolumn(yt) error('Input yt must be one column array!') end if yt(1)~=0 error('Airfoil thickness at leading edge must be zero!') end %Checking input consistency if (numel(xc)~=numel(yc))||(numel(xc)~=numel(yt)) error('Input xc, yc, and yt must have the same size!') end %Preallocating array for speed gc=zeros(size_xc(1),1); xu=zeros(size_xc(1),1); xl=zeros(size_xc(1),1); yu=zeros(size_xc(1),1); yl=zeros(size_xc(1),1); %Calculating camber line's gradients for i=2:1:size_xc(1)-1 gc(i)=(yc(i+1)-yc(i-1))/(xc(i+1)-xc(i-1)); end sc=atand(gc); %Generating airfoil vertexes for i=1:1:size_xc(1) xu(i)=xc(i)-yt(i)*sind(sc(i)); yu(i)=yc(i)+yt(i)*cosd(sc(i)); xl(i)=xc(i)+yt(i)*sind(sc(i)); yl(i)=yc(i)-yt(i)*cosd(sc(i)); end %CodeEnd------------------------------------------------------------------- end
airfoil codes/GenerateBezierCurve.m
function [x,y,z]=GenerateBezierCurve(xcp,ycp,zcp,dotcount) %-------------------------------------------------------------------------- %GenerateBezierCurve %Version 1.00 %Created by Stepen %Created 19 January 2011 %Last modified 31 January 2011 %-------------------------------------------------------------------------- %GenerateBezierCurve generates Bezier curve from given control points and %returns the vertexes of the Bezier curve. %-------------------------------------------------------------------------- %Syntax: %[x,y,z]=GenerateBezierCurve(xcp,ycp,zcp,dotcount) %Input argument: %- xcp (m x 1 num) specifies the x axis location of Bezier curve's control % points. %- ycp (m x 1 num) specifies the y axis location of Bezier curve's control % points. %- zcp (m x 1 num) specifies the z axis location of Bezier curve's control % points. %- dotcount (1 x 1 int) specifies the number of vertexes to be generated on % the Bezier curve. %Output argument: %- x (n x 1 num) specifies the x axis location of the Bezier curve's % vertexes. %- y (n x 1 num) specifies the y axis location of the Bezier curve's % vertexes. %- z (n x 1 num) specifies the z axis location of the Bezier curve's % vertexes. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input control points numel_x=numel(xcp); numel_y=numel(ycp); numel_z=numel(zcp); if (numel_x~=numel_y)||(numel_x~=numel_z) error('Input xcp, ycp, and zcp do not have the same size!') end bc=numel_x-1; %Checking input dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Calculating segment line length and position parameter length=1/(dotcount-1); s=0:length:1; %Preallocating array for speed x=zeros(dotcount,1); y=zeros(dotcount,1); z=zeros(dotcount,1); %Generating vertexes on Bezier Curve for i=1:1:numel(s) for j=1:1:numel_x x(i)=x(i)+((factorial(bc)/(factorial(j-1)*factorial(bc-j+1)))*((1-s(i))^(bc+1-j))*(s(i)^(j-1))*xcp(j)); y(i)=y(i)+((factorial(bc)/(factorial(j-1)*factorial(bc-j+1)))*((1-s(i))^(bc+1-j))*(s(i)^(j-1))*ycp(j)); z(i)=z(i)+((factorial(bc)/(factorial(j-1)*factorial(bc-j+1)))*((1-s(i))^(bc+1-j))*(s(i)^(j-1))*zcp(j)); end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/GenerateCSpline.m
function [x,y,z]=GenerateCSpline(xcp,ycp,zcp,dx,dy,dz,dotcount) %-------------------------------------------------------------------------- %GenerateCSpline %Version 1.00 %Created by Stepen %Created 23 January 2011 %Last modified 31 January 2011 %-------------------------------------------------------------------------- %GenerateCSpline generates cubic Hermite curve from given control points %and the tangent vector on each control points and returns the vertexes of %the curve. %-------------------------------------------------------------------------- %Syntax: %[x,y,z]=GenerateCSpline(xcp,ycp,zcp,dx,dy,dz,dotcount) %Input argument: %- xcp (m x 1 num) specifies the x axis location of cubic Hermite curve's % control points. %- ycp (m x 1 num) specifies the y axis location of cubic Hermite curve's % control points. %- zcp (m x 1 num) specifies the z axis location of cubic Hermite curve's % control points. %- dx (m x 1 num) specifes the unit vector component in x axis direction of % the corresponding tangent vector at each control point. %- dy (m x 1 num) specifes the unit vector component in y axis direction of % the corresponding tangent vector at each control point. %- dz (m x 1 num) specifes the unit vector component in z axis direction of % the corresponding tangent vector at each control point. %- dotcount (1 x 1 int) specifies the number of vertexes to be generated on % the polynomial curve. %Output argument: %- x (n x 1 num) specifies the x axis location of the cubic Hermite curve's % vertexes. %- y (n x 1 num) specifies the y axis location of the cubic Hermite curve's % vertexes. %- z (n x 1 num) specifies the z axis location of the cubic Hermite curve's % vertexes. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input control points and its tangent numel_x=numel(xcp); numel_y=numel(ycp); numel_z=numel(zcp); numel_dx=numel(dx); numel_dy=numel(dy); numel_dz=numel(dz); if (numel_x~=numel_y)||(numel_x~=numel_z)||(numel_x~=numel_dx)||... (numel_x~=numel_dy)||(numel_x~=numel_dz) error('Input xcp, ycp, and tangent vector do not have the same size!') end segmentcount=numel_x-1; %Checking input dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Preallocating array for speed x=zeros(segmentcount,dotcount); y=zeros(segmentcount,dotcount); z=zeros(segmentcount,dotcount); veclength=zeros(size(xcp)); %Calculating segment line length and position parameter length=1/(dotcount-1); s=0:length:1; %Normalizing tangent vector for i=1:1:numel_x veclength(i)=norm([dx(i),dy(i),dz(i)]); if veclength(i)==0 error('Input tangential vector must not be 0!') end end ndx=dx./veclength; ndy=dy./veclength; ndz=dz./veclength; %Generating vertexes on CSpline curve for i=1:1:segmentcount for j=1:1:dotcount x(i,j)=((2*s(j)^3)-(3*s(j)^2)+1)*xcp(i)+... ((-2*s(j)^3)+(3*s(j)^2))*xcp(i+1)+... ((s(j)^3)-(2*s(j)^2)+s(j))*ndx(i)*abs(xcp(i+1)-xcp(i))+... ((s(j)^3)-(s(j)^2))*ndx(i+1)*abs(xcp(i+1)-xcp(i)); y(i,j)=((2*s(j)^3)-(3*s(j)^2)+1)*ycp(i)+... ((-2*s(j)^3)+(3*s(j)^2))*ycp(i+1)+... ((s(j)^3)-(2*s(j)^2)+s(j))*ndy(i)*abs(ycp(i+1)-ycp(i))+... ((s(j)^3)-(s(j)^2))*ndy(i+1)*abs(ycp(i+1)-ycp(i)); z(i,j)=((2*s(j)^3)-(3*s(j)^2)+1)*zcp(i)+... ((-2*s(j)^3)+(3*s(j)^2))*zcp(i+1)+... ((s(j)^3)-(2*s(j)^2)+s(j))*ndz(i)*abs(zcp(i+1)-zcp(i))+... ((s(j)^3)-(s(j)^2))*ndz(i+1)*abs(zcp(i+1)-zcp(i)); end end %Rearranging vertexes x=reshape(x',numel(x),1); y=reshape(y',numel(y),1); z=reshape(z',numel(y),1); for j=segmentcount-1:-1:1 x(j*dotcount)=[]; y(j*dotcount)=[]; z(j*dotcount)=[]; end %CodeEnd------------------------------------------------------------------- end
airfoil codes/GenerateNACACamberLine.m
function [xc,yc]=GenerateNACACamberLine(cl,a,dotcount) %-------------------------------------------------------------------------- %GenerateNACACamberLine %Version 1.00 %Created by Stepen ([email protected] %Created 3 December 2010 %-------------------------------------------------------------------------- %GenerateNACACamberLine generates the camber line coordinate of a NACA %airfoil with given design lift coefficient and uniform load distribution %characteristic. The math equation used to generates the airfoil %coordinates is based on Theory of Wing Section Chapter 4 by Abbott and %Doenhoff. %-------------------------------------------------------------------------- %Syntax: [xc,yc]=GenerateNACACamberLine(cl,a,dotcount) %Input argument: %- cl (1 x 1 num) specifies the airfoil's design lift coefficient. %- a (1 x 1 num) specifies the chordwise location in fraction of chord % where load distribution is uniform from leading-edge to this point and % linearly decreasing to zero behind this point. %- dotcount (1 x 1 int) specifies the number of vertexes to be generated on % the airfoil's camber/mean line. %Output argument: %- xc (m x 1 num) specifies the x axis location of the camber line's % vertexes. %- yc (m x 1 num) specifies the y axis location of the camber line's % vertexes. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input cl if numel(cl)~=1 error('Lift Coefficient must be scalar!') end %Checking input a if numel(a)~=1 error('Uniform load location must be a scalar!') end if (a>1)||(a<0) error('Uniform load location values must be between 0 and 1!') end %Checking input dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Calculating x-axis location of camber line vertexes panelcount=dotcount-1; panellength=1/panelcount; xc=(0:panellength:1)'; %Preallocating array for speed yc=zeros(dotcount,1); %Calculating y-axis location of camber line vertexes for i=2:1:dotcount-1 tempval1=(-1/(1-a))*(((a^2)*((0.5*log(a))-0.25))+0.25); tempval2=((1/(1-a))*((0.5*((1-a)^2)*log(1-a))-(0.25*((1-a)^2))))+tempval1; if xc(i)==a xc(i)=xc(i)+0.0001*xc(i); end if a==1.0 yc(i)=(-cl/(4*pi()))*(((1-xc(i))*log(1-xc(i)))+(xc(i)*log(xc(i)))); else yc(i)=(cl/(2*pi()*(a+1)))*(((1/(1-a))*((0.5*((a-xc(i))^2)*log(abs(a-xc(i))))-... (0.5*((1-xc(i))^2)*log(abs(1-xc(i))))-... (0.25*((a-xc(i))^2))+... (0.25*((1-xc(i))^2))))-... (xc(i)*log(xc(i)))+... (tempval1)-... (tempval2*xc(i))); end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/GenerateNACASeries4Airfoil.m
function [xu,yu,xl,yl,xc,yc]=GenerateNACASeries4Airfoil(afid,dotcount,xmod) %-------------------------------------------------------------------------- %GenerateNACASeries4Airfoil %Version 1.20 %Created by Stepen ([email protected]) %Created 20 November 2010 %Last modified 30 November 2011 %-------------------------------------------------------------------------- %GenerateNACASeries4Airfoil generates the airfoil vertexes' coordinate of %the given NACA Series 4 Airfoil. The math equation used to generates the %airfoil coordinates is based on Theory of Wing Section Chapter 6 by Abbott %and Doenhoff. %-------------------------------------------------------------------------- %Syntax: %[xu,yu,xl,yl,xc,yc]=GenerateNACASeries4Airfoil(afid,dotcount,xmod) %Input argument: %- afid (1 x 4 str) specifies NACA Series 4 Airfoil identifier. %- dotcount (1 x 1 int) specifies the number of vertexes to be generated on % the airfoil's camber/mean line. %- xmod (str) specifies the mode of airfoil vertex distribution. Enter % 'Uniform' to create uniform vertex distribution or 'Cosine' to create % vertex distribution based on cosine function (More vertex on the leading % edge region). %Output argument: %- xu (i x 1 num) specifies the x axis location of airfoil's upper surface % vertexes in fraction of chord. The airfoil's upper surface vertex are % arranged from leading edge (the first element of xu) to the trailing % edge (the last element of xu). %- yu (i x 1 num) specifies the y axis location of airfoil's upper surface % vertexes in fraction of chord. The airfoil's upper surface vertex are % arranged from leading edge (the first element of yu) to the trailing % edge (the last element of yu). %- xl (i x 1 num) specifies the x axis location of airfoil's lower surface % vertexes in fraction of chord. The airfoil's lower surface vertex are % arranged from leading edge (the first element of xl) to the trailing % edge (the last element of xl). %- yl (i x 1 num) specifies the y axis location of airfoil's lower surface % vertexes in fraction of chord. The airfoil's lower surface vertex are % arranged from leading edge (the first element of yl) to the trailing % edge (the last element of yl). %- xc (i x 1 num) specifies the x axis location of airfoil's camber line % vertexes in fraction of chord. The airfoil's camber line vertex are % arranged from leading edge (the first element of xc) to the trailing % edge (the last element of xc). %- yc (i x 1 num) specifies the y axis location of airfoil's camber line % vertexes in fraction of chord. The airfoil's camber line vertex are % arranged from leading edge (the first element of yc) to the trailing % edge (the last element of yc). %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input afid if ~ischar(afid) error('Airfoil identifier must be a string!') end if numel(afid)~=4 error('Airfoil identifier must be a 4 digit number!') end if isempty(str2double(afid)) error('Airfoil identifier must be a 4 digit number!') end %Checking input dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Checking input xmod if nargin<3 xmod='Cosine'; end if (~strcmpi(xmod,'Uniform'))&&(~strcmpi(xmod,'Cosine')) error('Vertex distribution input must be Uniform or Cosine!') end %Assigning identifier to equation coefficient id1=str2double(afid(1)); id2=str2double(afid(2)); id3=str2double(afid([3,4])); m=id1*(1/100); %Maximum camber p=id2*(10/100); %Maximum camber location t=id3*(1/100); %Maximum thickness if t==0 warning(['Zero thickness airfoil!',... ' Airfoil will be just a camber line!']) end %Calculating x-axis location of camber line vertexes panelcount=dotcount-1; if strcmpi(xmod,'Uniform') panellength=1/panelcount; xc=(0:panellength:1)'; elseif strcmpi(xmod,'Cosine') deltadeg=90/panelcount; xc=1-cosd(0:deltadeg:90)'; end %Preallocating array for speed yc=zeros(dotcount,1); gc=zeros(dotcount,1); yt=zeros(dotcount,1); xu=zeros(dotcount,1); xl=zeros(dotcount,1); yu=zeros(dotcount,1); yl=zeros(dotcount,1); %Calculating y-axis location of camber line vertexes and camber gradient if m~=0 for i=1:1:dotcount if xc(i)<=p yc(i)=(m/(p^2))*((2*p*xc(i))-(xc(i)^2)); elseif xc(i)>p yc(i)=(m/((1-p)^2))*(1-(2*p)+(2*p*xc(i))-(xc(i)^2)); end gc(i)=(m/p^2)*(2*p+2*xc(i)); end end %Converting camber gradient to camber slope and normal sc=atand(gc); %Calculating thickness distribution for i=1:1:dotcount yt(i)=5*t*((0.29690*(xc(i)^0.5))-... (0.12600*xc(i))-... (0.35160*(xc(i)^2))+... (0.28430*(xc(i)^3))-... (0.10150*(xc(i)^4))); end %Generating airfoil vertexes for i=1:1:dotcount xu(i)=xc(i)-yt(i)*sind(sc(i)); yu(i)=yc(i)+yt(i)*cosd(sc(i)); xl(i)=xc(i)+yt(i)*sind(sc(i)); yl(i)=yc(i)-yt(i)*cosd(sc(i)); end %CodeEnd------------------------------------------------------------------- end
airfoil codes/GenerateNACASeries5Airfoil.m
function [xu,yu,xl,yl,xc,yc]=GenerateNACASeries5Airfoil(afid,dotcount,xmod) %-------------------------------------------------------------------------- %GenerateNACASeries5Airfoil %Version 1.20 %Created by Stepen ([email protected]) %Created 20 November 2010 %Last modified 30 November 2011 %-------------------------------------------------------------------------- %GenerateNACASeries5Airfoil generates the airfoil vertexes' coordinate of %the given NACA Series 5 Airfoil. The math equation used to generates the %airfoil coordinates is based on Theory of Wing Section Chapter 6 by Abbott %and Doenhoff. %-------------------------------------------------------------------------- %Syntax: %[xu,yu,xl,yl,xc,yc]=GenerateNACASeries5Airfoil(afid,dotcount,xmod) %Input argument: %- afid (1 x 5 str) specifies NACA Series 5 Airfoil identifier. %- dotcount (1 x 1 int) specifies the number of vertexes to be generated on % the airfoil's camber/mean line. %- xmod (str) specifies the mode of airfoil vertex distribution. Enter % 'Uniform' to create uniform vertex distribution or 'Cosine' to create % vertex distribution based on cosine function (More vertex on the leading % edge region). %Output argument: %- xu (i x 1 num) specifies the x axis location of airfoil's upper surface % vertexes in fraction of chord. The airfoil's upper surface vertex are % arranged from leading edge (the first element of xu) to the trailing % edge (the last element of xu). %- yu (i x 1 num) specifies the y axis location of airfoil's upper surface % vertexes in fraction of chord. The airfoil's upper surface vertex are % arranged from leading edge (the first element of yu) to the trailing % edge (the last element of yu). %- xl (i x 1 num) specifies the x axis location of airfoil's lower surface % vertexes in fraction of chord. The airfoil's lower surface vertex are % arranged from leading edge (the first element of xl) to the trailing % edge (the last element of xl). %- yl (i x 1 num) specifies the y axis location of airfoil's lower surface % vertexes in fraction of chord. The airfoil's lower surface vertex are % arranged from leading edge (the first element of yl) to the trailing % edge (the last element of yl). %- xc (i x 1 num) specifies the x axis location of airfoil's camber line % vertexes in fraction of chord. The airfoil's camber line vertex are % arranged from leading edge (the first element of xc) to the trailing % edge (the last element of xc). %- yc (i x 1 num) specifies the y axis location of airfoil's camber line % vertexes in fraction of chord. The airfoil's camber line vertex are % arranged from leading edge (the first element of yc) to the trailing % edge (the last element of yc). %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input afid if ~ischar(afid) error('Airfoil identifier must be a string!') end if numel(afid)~=5 error('Airfoil identifier must be a 5 digit number!') end if isempty(str2double(afid)) error('Airfoil identifier must be a 5 digit number!') end id=str2double(afid([2,3])); if (id~=10)&&(id~=20)&&(id~=30)&&(id~=40)&&(id~=50) error(['Airfoil identifier must start with X10,',... ' X20, X30, X40, or X50!']) end %Checking input dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Checking input xmod if nargin<3 xmod='Cosine'; end if (~strcmpi(xmod,'Uniform'))&&(~strcmpi(xmod,'Cosine')) error('Vertex distribution input must be Uniform or Cosine!') end %Declaring look-up table mtable=[0.0580,0.1260,0.2025,0.2900,0.3910]; ktable=[361.4000,51.6400,15.9570,6.6430,3.2300]; %Assigning identifier to equation coefficient id1=str2double(afid(1)); id2=str2double(afid([2,3])); id3=str2double(afid([4,5])); cl=id1*(3/2)*(10/100); %Design lift coefficient p=id2*(1/2)*(1/100); %Maximum camber location t=id3*(1/100); %Maximum thickness if t==0 warning(['Zero thickness airfoil!',... ' Airfoil will be just a camber line!']) end m=mtable(str2double(afid(2))); k=ktable(str2double(afid(2))); %Calculating x-axis location of camber line vertexes panelcount=dotcount-1; if strcmpi(xmod,'Uniform') panellength=1/panelcount; xc=(0:panellength:1)'; elseif strcmpi(xmod,'Cosine') deltadeg=90/panelcount; xc=1-cosd(0:deltadeg:90)'; end %Preallocating array for speed yc=zeros(dotcount,1); gc=zeros(dotcount,1); yt=zeros(dotcount,1); xu=zeros(dotcount,1); xl=zeros(dotcount,1); yu=zeros(dotcount,1); yl=zeros(dotcount,1); %Calculating y-axis location of camber line vertexes and camber gradient if m~=0 for i=1:1:dotcount if xc(i)<=p yc(i)=(k/6)*((xc(i)^3)-... (3*m*(xc(i)^2))+... ((m^2)*(3-m)*xc(i))); gc(i)=(k/6)*((2*(xc(i)^2))-... (6*m*xc(i))+... (((3*(m^2))-(m^3)))); elseif xc(i)>p yc(i)=(k/6)*(m^3)*(1-xc(i)); gc(i)=-1*(k/6)*(m^3); end end end %Correcting camber line vertexes and camber gradient for extended series yc=yc*(id1/2); gc=gc*(id1/2); %Converting camber gradient to camber slope and normal sc=atand(gc); %Calculating thickness distribution for i=1:1:dotcount yt(i)=5*t*((0.29690*(xc(i)^0.5))-... (0.12600*xc(i))-... (0.35160*(xc(i)^2))+... (0.28430*(xc(i)^3))-... (0.10150*(xc(i)^4))); end %Generating airfoil vertexes for i=1:1:dotcount xu(i)=xc(i)-yt(i)*sind(sc(i)); yu(i)=yc(i)+yt(i)*cosd(sc(i)); xl(i)=xc(i)+yt(i)*sind(sc(i)); yl(i)=yc(i)-yt(i)*cosd(sc(i)); end %CodeEnd------------------------------------------------------------------- end
airfoil codes/lecture.m
clear all close all figure(1) %%create points on airfoil using GenerateNACASeries4Airfoil [xu,yu,xl,yl,xc,yc]=GenerateNACASeries4Airfoil('0012',5,'Cosine') plot(xc,yc) axis('equal') %% run n cases for n=1:1 %%set angle of attack alpha=(n-3)*5 vx=cosd(alpha); vy=sind(alpha); %%calculate results for each case using sourcepanelmethod2d %% note that I had to rearrange the points given by generate airfoil to go from the LE on top to the TE and back along the bottom [lambda,cp,cl,cd]=SourcePanelMethod2D([xu; flipud(xl)],[yu; flipud(yl)],vx,vy,10); %%save array of results cla(n)=cl; cda(n)=cd; aoa(n)=alpha; end figure(2) plot(aoa,cla)
airfoil codes/PotentialFlow2DSim.m
function [u,v,un,vn,V]=PotentialFlow2DSim(elmnflow,x,y) %-------------------------------------------------------------------------- %PotentialFlow2DSim %Version 1.10 %Created by Stepen ([email protected]) %Created 24 November 2010 %Last modofied 3 August 2011 %-------------------------------------------------------------------------- %PotentialFlow2DSim generates the velocity field of a potential flow %(irrotational, inviscid, and incompressible) from the superposition of %the given potential elementary flow such as uniform flow, source, sink, %doublets, and vortex. %-------------------------------------------------------------------------- %Syntax: %[u,v,un,vn,V]=PotentialFlow2DSim(elmnflow,x,y) %Input argument: %- elmnflow (i x 5 num) specifies the elementary potential flow that % contribute to the flow. Each row of elmnflow array specifies one % elementary flow. The first column specifies elementary flow's type, % where 1 for uniform flow, 2 for source/sinks, 3 for doublets, and 4 for % vortex. % For uniform flow, the second and third column specifies the uniform % velocity on x and y axis direction. The fourth column is not used but % must be assigned to an arbitrary number for array completion. % For other flow, the second and third column specifies the x and y % axis location of elementary flow center and the fourth column % specifies its strength. % For doublet flow, the fifth column specifies the rotation of doublet. %- x (m x n num) specifies the x axis location for velocity vector % evaluation. %- y (m x n num) specifies the y axis location for velocity vector % evaluation. %Output argument %- u (m x n num) specifies the horizontal velocity component at the % corresponding location. %- v (m x n num) specifies the vertical velocity component at the % corresponding location. %- un (m x n num) specifies the normalized horizontal velocity component of % the corresponding location. un output was generated for plotting % purpose. Use this output instead of u for velocity quiver plotting. %- vn (m x n num) specifies the normalized vertical velocity component of % the corresponding location. vn output was generated for plotting % purpose. Use this output instead of v for velocity quiver plotting. %- V (m x n num) specifies the total velocity magnitude of the % corresponding loation. %-------------------------------------------------------------------------- %Example: %[x,y]=meshgrid(-1:0.1:1,-1:0.1:1); %[u,v,un,vn,V]=PotentialFlow2DSim([1,1,0,0;3,0,0,-10],x,y); % generates a combination of uniform flow with doublet centered at 0,0 % which forms a 2D flow around cylinder. The resulting flow field can be % displayed by plotting the normalized vector using quiver(x,y,un,vn) % command. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input elmnflow size_elmnflow=size(elmnflow); if numel(size_elmnflow)~=2 error('Elementary flow array must have size of i x 5!') end if size_elmnflow(2)~=5 error('Elementary flow array must have size of i x 5!') end %Checking input x and y size_x=size(x); dim_x=numel(size_x); size_y=size(y); dim_y=numel(size_y); if dim_x~=2 error('Input x must be a 2D array!') end if dim_y~=2 error('Input y must be a 2D array!') end if (size_x(1)~=size_y(1))||(size_x(2)~=size_y(2)) error('Input x and y do not have the same size!') end %Preallocating array for speed u=zeros(size_x); v=zeros(size_x); un=zeros(size_x); vn=zeros(size_x); V=zeros(size_x); %Calculating velocity vector from all elementary flow for i=1:1:size_x(1) for j=1:1:size_x(2) for k=1:1:size_elmnflow(1) if elmnflow(k,1)==1 utemp=elmnflow(k,2); vtemp=elmnflow(k,3); else xc=elmnflow(k,2); yc=elmnflow(k,3); sc=elmnflow(k,4); r=elmnflow(k,5); if elmnflow(k,1)==2 [utemp,vtemp]=CalculateVelocityFromSource2D... (sc,xc,yc,x(i,j),y(i,j)); elseif elmnflow(k,1)==3 [utemp,vtemp]=CalculateVelocityFromDoublet2D... (sc,xc,yc,x(i,j),y(i,j),r); elseif elmnflow(k,1)==4 [utemp,vtemp]=CalculateVelocityFromVortex2D... (sc,xc,yc,x(i,j),y(i,j)); end end u(i,j)=u(i,j)+utemp; v(i,j)=v(i,j)+vtemp; clear utemp clear vtemp clear xc yc sc r end end end %Calculating normalized vector for i=1:1:size_x(1) for j=1:1:size_x(2) V(i,j)=norm([u(i,j),v(i,j)]); un(i,j)=u(i,j)/V(i,j); vn(i,j)=v(i,j)/V(i,j); end end %CodeEnd------------------------------------------------------------------- end
airfoil codes/PotentialFlow2DSim_GUI.m
function PotentialFlow2DSim_GUI %-------------------------------------------------------------------------- %PotentialFlowSim_GUI %Version 1.00 %Created by Stepen %Created 31 March 2011 %-------------------------------------------------------------------------- %PotentialFlowSim_GUI starts a GUI program for PotentialFlow2DSim.m %-------------------------------------------------------------------------- %How to use PotentialFlowSim_GUI: %Input your elementary potential flow in the 'Flow Control' panel to %generate the quiver plot of your potential flow's velocity field. Use %'PlotControl' panel to pan or zoom your plot. If your quiver plot looks %poorly scaled, use 'Switch Plot' button to switch to quiver plot of %velocity direction and contour plot of velocity magnitude to get better %plot for your flow. If you want to save your velocity field data, use %'Save' button. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Preparing MATLAB environment close all clear all %Declaring global variable ScreenSize=get(0,'ScreenSize'); global x y u v un vn V elmnflow=zeros(0,5); elmnflowcount=0; indexview=1; xpos=0; ypos=0; zoomstate=1; resolution=10; plottype=1; %Preparing GUI mainwindow=figure('Name','PotentialFlow2DSim',... 'NumberTitle','Off',... 'Menubar','none',... 'Resize','off',... 'Units','pixels',... 'Position',[0.5*(ScreenSize(3)-1000),... 0.5*(ScreenSize(4)-600),... 1000,600]); axes('Parent',mainwindow,... 'Units','pixels',... 'Position',[25,25,550,550]); uicontrol('Parent',mainwindow,... 'Style','pushbutton',... 'String','Switch Plot',... 'Units','pixels',... 'Position',[0.775 0.05 0.1 0.05],... 'Callback',@switchplotfcn); uicontrol('Parent',mainwindow,... 'Style','pushbutton',... 'String','Save Data',... 'Units','pixels',... 'Position',[0.775 0.125 0.1 0.05],... 'Callback',@savedatafcn); flowcontrolpanel=uipanel('Parent',mainwindow,... 'Title','FlowControl',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[590,250,380,325]); uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','U:',... 'HorizontalAlignment','left',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[30,275,25,20]); uedit=uicontrol('Parent',flowcontrolpanel,... 'Style','edit',... 'String','0',... 'Units','pixels',... 'Position',[55,275,40,25],... 'Callback',@updateplotfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'Units','pixels',... 'Position',[95,288,20,12],... 'Callback',@uplusfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'Units','pixels',... 'Position',[95,275,20,12],... 'Callback',@uminfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','V:',... 'HorizontalAlignment','left',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[150,275,25,20]); vedit=uicontrol('Parent',flowcontrolpanel,... 'Style','edit',... 'String','0',... 'Units','pixels',... 'Position',[175,275,40,25],... 'Callback',@updateplotfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'Units','pixels',... 'Position',[215,288,20,12],... 'Callback',@vplusfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'Units','pixels',... 'Position',[215,275,20,12],... 'Callback',@vminfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','Elementary Flow',... 'HorizontalAlignment','center',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[30,235,60,30]); uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','Strength',... 'HorizontalAlignment','center',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[100,242.5,50,15]); uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','Center Location',... 'HorizontalAlignment','center',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[160,242.5,120,15]); uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','Rotation',... 'HorizontalAlignment','center',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[290,242.5,50,15]); moveelmnflowupbutton=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'HorizontalAlignment','center',... 'Enable','off',... 'Units','pixels',... 'Position',[2,235,21,10],... 'Callback',@moveelmnflowupfcn); moveelmnflowdownbutton=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'HorizontalAlignment','center',... 'Enable','off',... 'Units','pixels',... 'Position',[2,25,21,10],... 'Callback',@moveelmnflowdownfcn); uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','Add Elementary Flow',... 'HorizontalAlignment','center',... 'Units','pixels',... 'Position',[250,5,120,20],... 'Callback',@addelmnflowfcn); elmnindextext=zeros(8,1); elmnflowpopup=zeros(8,1); strengthedit=zeros(8,1); strengthplusbutton=zeros(8,1); strengthminbutton=zeros(8,1); xcedit=zeros(8,1); xcplusbutton=zeros(8,1); xcminbutton=zeros(8,1); ycedit=zeros(8,1); ycplusbutton=zeros(8,1); ycminbutton=zeros(8,1); selectcenterbutton=zeros(8,1); thetaedit=zeros(8,1); thetaplusbutton=zeros(8,1); thetaminbutton=zeros(8,1); rmvelmnflowbutton=zeros(8,1); for i=1:8 elmnindextext(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','text',... 'String','0',... 'HorizontalAlignment','center',... 'BackgroundColor',[0.8,0.8,0.8],... 'Visible','off',... 'Enable','on',... 'Units','pixels',... 'Position',[2,235-(i*25),21,20],... 'Callback',@updateplotfcn); elmnflowpopup(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','popup',... 'String',{'Source',... 'Doublet',... 'Vortex'},... 'Enable','off',... 'Units','pixels',... 'Position',[30,235-(i*25),60,25],... 'Callback',@updateplotfcn); strengthedit(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','edit',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[100,239-(i*25),30,20],... 'Callback',@updateplotfcn); strengthplusbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[130,249-(i*25),20,10]); strengthminbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[130,239-(i*25),20,10]); xcedit(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','edit',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[160,239-(i*25),30,20],... 'Callback',@updateplotfcn); xcplusbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[190,249-(i*25),20,10]); xcminbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[190,239-(i*25),20,10]); ycedit(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','edit',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[215,239-(i*25),30,20],... 'Callback',@updateplotfcn); ycplusbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[245,249-(i*25),20,10]); ycminbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[245,239-(i*25),20,10]); selectcenterbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[270,239-(i*25),10,20]); thetaedit(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','edit',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[290,239-(i*25),30,20],... 'Callback',@updateplotfcn); thetaplusbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[320,249-(i*25),20,10]); thetaminbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','',... 'Enable','off',... 'Units','pixels',... 'Position',[320,239-(i*25),20,10]); rmvelmnflowbutton(i)=uicontrol('Parent',flowcontrolpanel,... 'Style','pushbutton',... 'String','X',... 'HorizontalAlignment','center',... 'Enable','off',... 'Units','pixels',... 'Position',[350,239-(i*25),... 25,20]); end set(strengthplusbutton(1),'Callback',@strength1plusfcn) set(strengthplusbutton(2),'Callback',@strength2plusfcn) set(strengthplusbutton(3),'Callback',@strength3plusfcn) set(strengthplusbutton(4),'Callback',@strength4plusfcn) set(strengthplusbutton(5),'Callback',@strength5plusfcn) set(strengthplusbutton(6),'Callback',@strength6plusfcn) set(strengthplusbutton(7),'Callback',@strength7plusfcn) set(strengthplusbutton(8),'Callback',@strength8plusfcn) set(strengthminbutton(1),'Callback',@strength1minfcn) set(strengthminbutton(2),'Callback',@strength2minfcn) set(strengthminbutton(3),'Callback',@strength3minfcn) set(strengthminbutton(4),'Callback',@strength4minfcn) set(strengthminbutton(5),'Callback',@strength5minfcn) set(strengthminbutton(6),'Callback',@strength6minfcn) set(strengthminbutton(7),'Callback',@strength7minfcn) set(strengthminbutton(8),'Callback',@strength8minfcn) set(xcplusbutton(1),'Callback',@xc1plusfcn) set(xcplusbutton(2),'Callback',@xc2plusfcn) set(xcplusbutton(3),'Callback',@xc3plusfcn) set(xcplusbutton(4),'Callback',@xc4plusfcn) set(xcplusbutton(5),'Callback',@xc5plusfcn) set(xcplusbutton(6),'Callback',@xc6plusfcn) set(xcplusbutton(7),'Callback',@xc7plusfcn) set(xcplusbutton(8),'Callback',@xc8plusfcn) set(xcminbutton(1),'Callback',@xc1minfcn) set(xcminbutton(2),'Callback',@xc2minfcn) set(xcminbutton(3),'Callback',@xc3minfcn) set(xcminbutton(4),'Callback',@xc4minfcn) set(xcminbutton(5),'Callback',@xc5minfcn) set(xcminbutton(6),'Callback',@xc6minfcn) set(xcminbutton(7),'Callback',@xc7minfcn) set(xcminbutton(8),'Callback',@xc8minfcn) set(ycplusbutton(1),'Callback',@yc1plusfcn) set(ycplusbutton(2),'Callback',@yc2plusfcn) set(ycplusbutton(3),'Callback',@yc3plusfcn) set(ycplusbutton(4),'Callback',@yc4plusfcn) set(ycplusbutton(5),'Callback',@yc5plusfcn) set(ycplusbutton(6),'Callback',@yc6plusfcn) set(ycplusbutton(7),'Callback',@yc7plusfcn) set(ycplusbutton(8),'Callback',@yc8plusfcn) set(ycminbutton(1),'Callback',@yc1minfcn) set(ycminbutton(2),'Callback',@yc2minfcn) set(ycminbutton(3),'Callback',@yc3minfcn) set(ycminbutton(4),'Callback',@yc4minfcn) set(ycminbutton(5),'Callback',@yc5minfcn) set(ycminbutton(6),'Callback',@yc6minfcn) set(ycminbutton(7),'Callback',@yc7minfcn) set(ycminbutton(8),'Callback',@yc8minfcn) set(selectcenterbutton(1),'Callback',@selectcenter1fcn) set(selectcenterbutton(2),'Callback',@selectcenter2fcn) set(selectcenterbutton(3),'Callback',@selectcenter3fcn) set(selectcenterbutton(4),'Callback',@selectcenter4fcn) set(selectcenterbutton(5),'Callback',@selectcenter5fcn) set(selectcenterbutton(6),'Callback',@selectcenter6fcn) set(selectcenterbutton(7),'Callback',@selectcenter7fcn) set(selectcenterbutton(8),'Callback',@selectcenter8fcn) set(thetaplusbutton(1),'Callback',@theta1plusfcn) set(thetaplusbutton(2),'Callback',@theta2plusfcn) set(thetaplusbutton(3),'Callback',@theta3plusfcn) set(thetaplusbutton(4),'Callback',@theta4plusfcn) set(thetaplusbutton(5),'Callback',@theta5plusfcn) set(thetaplusbutton(6),'Callback',@theta6plusfcn) set(thetaplusbutton(7),'Callback',@theta7plusfcn) set(thetaplusbutton(8),'Callback',@theta8plusfcn) set(thetaminbutton(1),'Callback',@theta1minfcn) set(thetaminbutton(2),'Callback',@theta2minfcn) set(thetaminbutton(3),'Callback',@theta3minfcn) set(thetaminbutton(4),'Callback',@theta4minfcn) set(thetaminbutton(5),'Callback',@theta5minfcn) set(thetaminbutton(6),'Callback',@theta6minfcn) set(thetaminbutton(7),'Callback',@theta7minfcn) set(thetaminbutton(8),'Callback',@theta8minfcn) set(rmvelmnflowbutton(1),'Callback',@rmvelmnflow1fcn) set(rmvelmnflowbutton(2),'Callback',@rmvelmnflow2fcn) set(rmvelmnflowbutton(3),'Callback',@rmvelmnflow3fcn) set(rmvelmnflowbutton(4),'Callback',@rmvelmnflow4fcn) set(rmvelmnflowbutton(5),'Callback',@rmvelmnflow5fcn) set(rmvelmnflowbutton(6),'Callback',@rmvelmnflow6fcn) set(rmvelmnflowbutton(7),'Callback',@rmvelmnflow7fcn) set(rmvelmnflowbutton(8),'Callback',@rmvelmnflow8fcn) flowvisualizationpanel=uipanel('Parent',mainwindow,... 'Title','FlowVisualization',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[590,140,185,100]); uicontrol('Parent',flowvisualizationpanel,... 'Style','pushbutton',... 'String','Distribute particles',... 'Units','pixels',... 'Position',[42.5,10,100,20],... 'Callback',@flowvisfcn) uicontrol('Parent',flowvisualizationpanel,... 'Style','pushbutton',... 'String','Place particles',... 'Units','pixels',... 'Position',[42.5,60,100,20],... 'Callback',@flowviscustomfcn) plotcontrolpanel=uipanel('Parent',mainwindow,... 'Title','PlotControl',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[785,140,185,100]); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','+',... 'Units','pixels',... 'Position',[45 45 20 20],... 'Callback',@increaseresolutionfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','-',... 'Units','pixels',... 'Position',[65 45 20 20],... 'Callback',@decreaseresolutionfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Left',... 'Units','pixels',... 'Position',[5 45 40 20],... 'Callback',@panleftfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Up',... 'Units','pixels',... 'Position',[45 65 40 20],... 'Callback',@panupfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Down',... 'Units','pixels',... 'Position',[45 25 40 20],... 'Callback',@pandownfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Right',... 'Units','pixels',... 'Position',[85 45 40 20],... 'Callback',@panrightfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Zoom in',... 'Units','pixels',... 'Position',[125 55 50 20],... 'Callback',@zoominfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Zoom out',... 'Units','pixels',... 'Position',[125 35 50 20],... 'Callback',@zoomoutfcn); uicontrol('Parent',plotcontrolpanel,... 'Style','pushbutton',... 'String','Switch View',... 'Units','pixels',... 'Position',[15 5 155 20],... 'Callback',@switchviewfcn); uicontrol('Parent',mainwindow',... 'Style','pushbutton',... 'String','Save Data',... 'Units','pixels',... 'Position',[590 100 380 25],... 'Callback',@savedatafcn) uicontrol('Parent',mainwindow',... 'Style','pushbutton',... 'String','Save Plot',... 'Units','pixels',... 'Position',[590 75 380 25],... 'Callback',@saveplotfcn) msgpanel=uipanel('Parent',mainwindow,... 'Title','Message Box',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[590 25 380 40]); msgtext=uicontrol('Parent',msgpanel,... 'Style','text',... 'String',['Start creating your potential flow',... ' by modifying FlowControl panel...'],... 'HorizontalAlignment','left',... 'BackgroundColor',[0.8,0.8,0.8],... 'Units','pixels',... 'Position',[5 5 370 15]); %Listing local function %Start of updateplot function updateplot %Creating plot grid xmin=xpos-0.5*zoomstate; xmax=xpos+0.5*zoomstate; ymin=ypos-0.5*zoomstate; ymax=ypos+0.5*zoomstate; delta=zoomstate/resolution; [x,y]=meshgrid(xmin:delta:xmax,ymin:delta:ymax); x=round(x/delta)*delta; y=round(y/delta)*delta; %Reading flow properties readpanel; %Simulating flow [u,v,un,vn,V]=PotentialFlow2DSim(elmnflow,x,y); %Updating plot cla hold on if plottype==2 contourf(x,y,V,20,'LineStyle','none') end if plottype==1 quiver(x,y,u,v,'k') elseif plottype==2 quiver(x,y,un,vn,'k') end axis equal axis ([xmin xmax ymin ymax]) %Displaying message set(msgtext,'String','Plot updated!') end %End of updateplot %Start of readpanel function readpanel %Reading uniform flow panel elmnflow(1,1)=1; temp=str2double(get(uedit,'String')); if ~isnan(temp) elmnflow(1,2)=temp; end temp=str2double(get(vedit,'String')); if ~isnan(temp) elmnflow(1,3)=temp; end %Reading elementary flow data maxview=min(8,elmnflowcount); for panel=1:maxview for stat=1:5 if stat==1 temp=get(elmnflowpopup(panel),'Value')+1; elseif stat==4 temp=str2double(get(strengthedit(panel),'String')); elseif stat==2 temp=str2double(get(xcedit(panel),'String')); elseif stat==3 temp=str2double(get(ycedit(panel),'String')); elseif stat==5 temp=str2double(get(thetaedit(panel),'String')); end if ~isnan(temp) elmnflow(panel+indexview,stat)=temp; else elmnflow(panel+indexview,stat)=0; end end end end %End of readpanel %Start of updatepanel function updatepanel %Updating uniform flow edit boxes set(uedit,'String',elmnflow(1,2)) set(vedit,'String',elmnflow(1,3)) %Updating view changer button if elmnflowcount<=8 set(moveelmnflowupbutton,'Enable','off') set(moveelmnflowdownbutton,'Enable','off') else set(moveelmnflowupbutton,'Enable','on') set(moveelmnflowdownbutton,'Enable','on') end %Resetting elementary flow edit boxes set(elmnindextext,'Visible','off') set(elmnflowpopup,'Enable','off') set(strengthedit,'String','') set(strengthedit,'Enable','off') set(strengthplusbutton,'Enable','off') set(strengthminbutton,'Enable','off') set(xcedit,'String','') set(xcedit,'Enable','off') set(xcplusbutton,'Enable','off') set(xcminbutton,'Enable','off') set(ycedit,'String','') set(ycedit,'Enable','off') set(ycplusbutton,'Enable','off') set(ycminbutton,'Enable','off') set(selectcenterbutton,'Enable','off') set(thetaedit,'String','') set(thetaedit,'Enable','off') set(thetaplusbutton,'Enable','off') set(thetaminbutton,'Enable','off') set(rmvelmnflowbutton,'Enable','off') %Updating elementary flow edit boxes maxview=min(8,elmnflowcount); for panel=1:maxview set(elmnindextext(panel),'String',num2str(indexview+panel-1)) set(elmnindextext(panel),'Visible','on') set(elmnflowpopup(panel),'Value',elmnflow(indexview+panel,1)-1) set(elmnflowpopup(panel),'Enable','on') set(strengthedit(panel),'String',elmnflow(indexview+panel,4)) set(strengthedit(panel),'Enable','on') set(strengthplusbutton(panel),'Enable','on') set(strengthminbutton(panel),'Enable','on') set(xcedit(panel),'String',elmnflow(indexview+panel,2)) set(xcedit(panel),'Enable','on') set(xcplusbutton(panel),'Enable','on') set(xcminbutton(panel),'Enable','on') set(ycedit(panel),'String',elmnflow(indexview+panel,3)) set(ycedit(panel),'Enable','on') set(ycplusbutton(panel),'Enable','on') set(ycminbutton(panel),'Enable','on') set(selectcenterbutton(panel),'Enable','on') set(thetaedit(panel),'String',elmnflow(indexview+panel,5)) set(thetaedit(panel),'Enable','on') set(thetaplusbutton(panel),'Enable','on') set(thetaminbutton(panel),'Enable','on') set(rmvelmnflowbutton(panel),'Enable','on') end end %End of updatepanel %Listing callback function %Start of panleftfcn function panleftfcn(~,~) xpos=xpos-0.1*zoomstate; updateplot end %End of panleftfcn %Start of panupfcn function panupfcn(~,~) ypos=ypos+0.1*zoomstate; updateplot end %End of panupfcn %Start of pandownfcn function pandownfcn(~,~) ypos=ypos-0.1*zoomstate; updateplot end %End of pandownfcn %Start of panrightfcn function panrightfcn(~,~) xpos=xpos+0.1*zoomstate; updateplot end %End of panrightfcn %Start of zoominfcn function zoominfcn(~,~) zoomstate=zoomstate/2; updateplot end %End of zoominfcn %Start of zoomoutfcn function zoomoutfcn(~,~) zoomstate=zoomstate*2; updateplot end %End of zoomoutfcn %Start of increaseresolutionfcn function increaseresolutionfcn(~,~) resolution=resolution*2; updateplot end %End of increaseresolutionfcn %Start of decreaseresolutionfcn function decreaseresolutionfcn(~,~) resolution=resolution/2; updateplot end %End of decreaseresolutionfcn %Start of updateplotfcn function updateplotfcn(~,~) readpanel updateplot end %End of updateplotfcn %Start of moveelmnflowupfcn function moveelmnflowupfcn(~,~) if indexview>1 indexview=indexview-1; updatepanel end end %End of moveelmnflowupfcn %Start of moveelmnflowdownfcn function moveelmnflowdownfcn(~,~) if indexview+7<elmnflowcount indexview=indexview+1; updatepanel end end %End of moveelmnflowdownfcn %Start of uplusfcn function uplusfcn(~,~) set(uedit,'String',num2str(str2double(get(uedit,'String'))+0.1)) readpanel updateplot end %End of uplusfcn %Start of vplusfcn function vplusfcn(~,~) set(vedit,'String',num2str(str2double(get(vedit,'String'))+0.1)) readpanel updateplot end %End of vplusfcn %Start of strengthplus1button function strength1plusfcn(~,~) set(strengthedit(1),... 'String',... num2str(str2double(get(strengthedit(1),'String'))+0.1)) readpanel updateplot end %End of strengthplus1button %Start of strengthplus2button function strength2plusfcn(~,~) set(strengthedit(2),... 'String',... num2str(str2double(get(strengthedit(2),'String'))+0.1)) readpanel updateplot end %End of strengthplus2button %Start of strengthplus3button function strength3plusfcn(~,~) set(strengthedit(3),... 'String',... num2str(str2double(get(strengthedit(3),'String'))+0.1)) readpanel updateplot end %End of strengthplus3button %Start of strengthplus4button function strength4plusfcn(~,~) set(strengthedit(4),... 'String',... num2str(str2double(get(strengthedit(4),'String'))+0.1)) readpanel updateplot end %End of strengthplus4button %Start of strengthplus5button function strength5plusfcn(~,~) set(strengthedit(5),... 'String',... num2str(str2double(get(strengthedit(5),'String'))+0.1)) readpanel updateplot end %End of strengthplus5button %Start of strengthplus6button function strength6plusfcn(~,~) set(strengthedit(6),... 'String',... num2str(str2double(get(strengthedit(6),'String'))+0.1)) readpanel updateplot end %End of strengthplus6button %Start of strengthplus7button function strength7plusfcn(~,~) set(strengthedit(7),... 'String',... num2str(str2double(get(strengthedit(7),'String'))+0.1)) readpanel updateplot end %End of strengthplus7button %Start of strengthplus8button function strength8plusfcn(~,~) set(strengthedit(8),... 'String',... num2str(str2double(get(strengthedit(8),'String'))+0.1)) readpanel updateplot end %End of strengthplus8button %Start of xc1plusfcn function xc1plusfcn(~,~) set(xcedit(1),... 'String',num2str(str2double(get(xcedit(1),'String'))+0.1)) readpanel updateplot end %End of xc1plusfcn %Start of xc2plusfcn function xc2plusfcn(~,~) set(xcedit(2),... 'String',num2str(str2double(get(xcedit(2),'String'))+0.1)) readpanel updateplot end %End of xc2plusfcn %Start of xc3plusfcn function xc3plusfcn(~,~) set(xcedit(3),... 'String',num2str(str2double(get(xcedit(3),'String'))+0.1)) readpanel updateplot end %End of xc3plusfcn %Start of xc4plusfcn function xc4plusfcn(~,~) set(xcedit(4),... 'String',num2str(str2double(get(xcedit(4),'String'))+0.1)) readpanel updateplot end %End of xc4plusfcn %Start of xc5plusfcn function xc5plusfcn(~,~) set(xcedit(5),... 'String',num2str(str2double(get(xcedit(5),'String'))+0.1)) readpanel updateplot end %End of xc5plusfcn %Start of xc6plusfcn function xc6plusfcn(~,~) set(xcedit(6),... 'String',num2str(str2double(get(xcedit(6),'String'))+0.1)) readpanel updateplot end %End of xc6plusfcn %Start of xc7plusfcn function xc7plusfcn(~,~) set(xcedit(7),... 'String',num2str(str2double(get(xcedit(7),'String'))+0.1)) readpanel updateplot end %End of xc7plusfcn %Start of xc8plusfcn function xc8plusfcn(~,~) set(xcedit(8),... 'String',num2str(str2double(get(xcedit(8),'String'))+0.1)) readpanel updateplot end %End of xc8plusfcn %Start of yc1plusfcn function yc1plusfcn(~,~) set(ycedit(1),... 'String',num2str(str2double(get(ycedit(1),'String'))+0.1)) readpanel updateplot end %End of yc1plusfcn %Start of yc2plusfcn function yc2plusfcn(~,~) set(ycedit(2),... 'String',num2str(str2double(get(ycedit(2),'String'))+0.1)) readpanel updateplot end %End of yc2plusfcn %Start of yc3plusfcn function yc3plusfcn(~,~) set(ycedit(3),... 'String',num2str(str2double(get(ycedit(3),'String'))+0.1)) readpanel updateplot end %End of yc3plusfcn %Start of yc4plusfcn function yc4plusfcn(~,~) set(ycedit(4),... 'String',num2str(str2double(get(ycedit(4),'String'))+0.1)) readpanel updateplot end %End of yc4plusfcn %Start of yc5plusfcn function yc5plusfcn(~,~) set(ycedit(5),... 'String',num2str(str2double(get(ycedit(5),'String'))+0.1)) readpanel updateplot end %End of yc5plusfcn %Start of yc6plusfcn function yc6plusfcn(~,~) set(ycedit(6),... 'String',num2str(str2double(get(ycedit(6),'String'))+0.1)) readpanel updateplot end %End of yc6plusfcn %Start of yc7plusfcn function yc7plusfcn(~,~) set(ycedit(7),... 'String',num2str(str2double(get(ycedit(7),'String'))+0.1)) readpanel updateplot end %End of yc7plusfcn %Start of yc8plusfcn function yc8plusfcn(~,~) set(ycedit(8),... 'String',num2str(str2double(get(ycedit(8),'String'))+0.1)) readpanel updateplot end %End of yc8plusfcn %Start of theta1plusfcn function theta1plusfcn(~,~) set(thetaedit(1),... 'String',num2str(str2double(get(thetaedit(1),'String'))+0.1)) readpanel updateplot end %End of theta1plusfcn %Start of theta2plusfcn function theta2plusfcn(~,~) set(thetaedit(2),... 'String',num2str(str2double(get(thetaedit(2),'String'))+0.1)) readpanel updateplot end %End of theta2plusfcn %Start of theta3plusfcn function theta3plusfcn(~,~) set(thetaedit(3),... 'String',num2str(str2double(get(thetaedit(3),'String'))+0.1)) readpanel updateplot end %End of theta3plusfcn %Start of theta4plusfcn function theta4plusfcn(~,~) set(thetaedit(4),... 'String',num2str(str2double(get(thetaedit(4),'String'))+0.1)) readpanel updateplot end %End of theta4plusfcn %Start of theta5plusfcn function theta5plusfcn(~,~) set(thetaedit(5),... 'String',num2str(str2double(get(thetaedit(5),'String'))+0.1)) readpanel updateplot end %End of theta5plusfcn %Start of theta6plusfcn function theta6plusfcn(~,~) set(thetaedit(6),... 'String',num2str(str2double(get(thetaedit(6),'String'))+0.1)) readpanel updateplot end %End of theta6plusfcn %Start of theta7plusfcn function theta7plusfcn(~,~) set(thetaedit(7),... 'String',num2str(str2double(get(thetaedit(7),'String'))+0.1)) readpanel updateplot end %End of theta7plusfcn %Start of theta8plusfcn function theta8plusfcn(~,~) set(thetaedit(8),... 'String',num2str(str2double(get(thetaedit(8),'String'))+0.1)) readpanel updateplot end %End of theta8plusfcn %Start of uminfcn function uminfcn(~,~) set(uedit,'String',num2str(str2double(get(uedit,'String'))-0.1)) readpanel updateplot end %End of uminfcn %Start of vminfcn function vminfcn(~,~) set(vedit,'String',num2str(str2double(get(vedit,'String'))-0.1)) readpanel updateplot end %End of vminfcn %Start of strengthmin1button function strength1minfcn(~,~) set(strengthedit(1),... 'String',... num2str(str2double(get(strengthedit(1),'String'))-0.1)) readpanel updateplot end %End of strengthmin1button %Start of strengthmin2button function strength2minfcn(~,~) set(strengthedit(2),... 'String',... num2str(str2double(get(strengthedit(2),'String'))-0.1)) readpanel updateplot end %End of strengthmin2button %Start of strengthmin3button function strength3minfcn(~,~) set(strengthedit(3),... 'String',... num2str(str2double(get(strengthedit(3),'String'))-0.1)) readpanel updateplot end %End of strengthmin3button %Start of strengthmin4button function strength4minfcn(~,~) set(strengthedit(4),... 'String',... num2str(str2double(get(strengthedit(4),'String'))-0.1)) readpanel updateplot end %End of strengthmin4button %Start of strengthmin5button function strength5minfcn(~,~) set(strengthedit(5),... 'String',... num2str(str2double(get(strengthedit(5),'String'))-0.1)) readpanel updateplot end %End of strengthmin5button %Start of strengthmin6button function strength6minfcn(~,~) set(strengthedit(6),... 'String',... num2str(str2double(get(strengthedit(6),'String'))-0.1)) readpanel updateplot end %End of strengthmin6button %Start of strengthmin7button function strength7minfcn(~,~) set(strengthedit(7),... 'String',... num2str(str2double(get(strengthedit(7),'String'))-0.1)) readpanel updateplot end %End of strengthmin7button %Start of strengthmin8button function strength8minfcn(~,~) set(strengthedit(8),... 'String',... num2str(str2double(get(strengthedit(8),'String'))-0.1)) readpanel updateplot end %End of strengthmin8button %Start of xc1minfcn function xc1minfcn(~,~) set(xcedit(1),... 'String',num2str(str2double(get(xcedit(1),'String'))-0.1)) readpanel updateplot end %End of xc1minfcn %Start of xc2minfcn function xc2minfcn(~,~) set(xcedit(2),... 'String',num2str(str2double(get(xcedit(2),'String'))-0.1)) readpanel updateplot end %End of xc2minfcn %Start of xc3minfcn function xc3minfcn(~,~) set(xcedit(3),... 'String',num2str(str2double(get(xcedit(3),'String'))-0.1)) readpanel updateplot end %End of xc3minfcn %Start of xc4minfcn function xc4minfcn(~,~) set(xcedit(4),... 'String',num2str(str2double(get(xcedit(4),'String'))-0.1)) readpanel updateplot end %End of xc4minfcn %Start of xc5minfcn function xc5minfcn(~,~) set(xcedit(5),... 'String',num2str(str2double(get(xcedit(5),'String'))-0.1)) readpanel updateplot end %End of xc5minfcn %Start of xc6minfcn function xc6minfcn(~,~) set(xcedit(6),... 'String',num2str(str2double(get(xcedit(6),'String'))-0.1)) readpanel updateplot end %End of xc6minfcn %Start of xc7minfcn function xc7minfcn(~,~) set(xcedit(7),... 'String',num2str(str2double(get(xcedit(7),'String'))-0.1)) readpanel updateplot end %End of xc7minfcn %Start of xc8minfcn function xc8minfcn(~,~) set(xcedit(8),... 'String',num2str(str2double(get(xcedit(8),'String'))-0.1)) readpanel updateplot end %End of xc8minfcn %Start of yc1minfcn function yc1minfcn(~,~) set(ycedit(1),... 'String',num2str(str2double(get(ycedit(1),'String'))-0.1)) readpanel updateplot end %End of yc1minfcn %Start of yc2minfcn function yc2minfcn(~,~) set(ycedit(2),... 'String',num2str(str2double(get(ycedit(2),'String'))-0.1)) readpanel updateplot end %End of yc2minfcn %Start of yc3minfcn function yc3minfcn(~,~) set(ycedit(3),... 'String',num2str(str2double(get(ycedit(3),'String'))-0.1)) readpanel updateplot end %End of yc3minfcn %Start of yc4minfcn function yc4minfcn(~,~) set(ycedit(4),... 'String',num2str(str2double(get(ycedit(4),'String'))-0.1)) readpanel updateplot end %End of yc4minfcn %Start of yc5minfcn function yc5minfcn(~,~) set(ycedit(5),... 'String',num2str(str2double(get(ycedit(5),'String'))-0.1)) readpanel updateplot end %End of yc5minfcn %Start of yc6minfcn function yc6minfcn(~,~) set(ycedit(6),... 'String',num2str(str2double(get(ycedit(6),'String'))-0.1)) readpanel updateplot end %End of yc6minfcn %Start of yc7minfcn function yc7minfcn(~,~) set(ycedit(7),... 'String',num2str(str2double(get(ycedit(7),'String'))-0.1)) readpanel updateplot end %End of yc7minfcn %Start of yc8minfcn function yc8minfcn(~,~) set(ycedit(8),... 'String',num2str(str2double(get(ycedit(8),'String'))-0.1)) readpanel updateplot end %End of yc8minfcn %Start of theta1minfcn function theta1minfcn(~,~) set(thetaedit(1),... 'String',num2str(str2double(get(thetaedit(1),'String'))-0.1)) readpanel updateplot end %End of theta1minfcn %Start of theta2minfcn function theta2minfcn(~,~) set(thetaedit(2),... 'String',num2str(str2double(get(thetaedit(2),'String'))-0.1)) readpanel updateplot end %End of theta2minfcn %Start of theta3minfcn function theta3minfcn(~,~) set(thetaedit(3),... 'String',num2str(str2double(get(thetaedit(3),'String'))-0.1)) readpanel updateplot end %End of theta3minfcn %Start of theta4minfcn function theta4minfcn(~,~) set(thetaedit(4),... 'String',num2str(str2double(get(thetaedit(4),'String'))-0.1)) readpanel updateplot end %End of theta4minfcn %Start of theta5minfcn function theta5minfcn(~,~) set(thetaedit(5),... 'String',num2str(str2double(get(thetaedit(5),'String'))-0.1)) readpanel updateplot end %End of theta5minfcn %Start of theta6minfcn function theta6minfcn(~,~) set(thetaedit(6),... 'String',num2str(str2double(get(thetaedit(6),'String'))-0.1)) readpanel updateplot end %End of theta6minfcn %Start of theta7minfcn function theta7minfcn(~,~) set(thetaedit(7),... 'String',num2str(str2double(get(thetaedit(7),'String'))-0.1)) readpanel updateplot end %End of theta7minfcn %Start of theta8minfcn function theta8minfcn(~,~) set(thetaedit(8),... 'String',num2str(str2double(get(thetaedit(8),'String'))-0.1)) readpanel updateplot end %End of theta8minfcn %Start of selectcenter1fcn function selectcenter1fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 1 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(1),'String',num2str(xp(1))) set(ycedit(1),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter1fcn %Start of selectcenter2fcn function selectcenter2fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 2 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(2),'String',num2str(xp(1))) set(ycedit(2),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter2fcn %Start of selectcenter3fcn function selectcenter3fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 3 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(3),'String',num2str(xp(1))) set(ycedit(3),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter3fcn %Start of selectcenter4fcn function selectcenter4fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 4 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(4),'String',num2str(xp(1))) set(ycedit(4),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter4fcn %Start of selectcenter5fcn function selectcenter5fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 5 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(5),'String',num2str(xp(1))) set(ycedit(5),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter5fcn %Start of selectcenter6fcn function selectcenter6fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 6 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(6),'String',num2str(xp(1))) set(ycedit(6),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter6fcn %Start of selectcenter7fcn function selectcenter7fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 7 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(7),'String',num2str(xp(1))) set(ycedit(7),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter7fcn %Start of selectcenter8fcn function selectcenter8fcn(~,~) %Displaying message set(msgtext,... 'String',['Select a point on plot to move',... ' the center of elementary flow no 8 to...']) %Prompting user input [xp,yp]=getpts; %Applying changes set(xcedit(8),'String',num2str(xp(1))) set(ycedit(8),'String',num2str(yp(1))) %Updating plot readpanel updateplot end %End of selectcenter8fcn %Start of addelmnflowfcn function addelmnflowfcn(~,~) %Adding new elementary flow if elmnflowcount>=8 indexview=elmnflowcount-6; end elmnflowcount=elmnflowcount+1; elmnflow(elmnflowcount+1,:)=[2,0,0,1,0]; updatepanel updateplot %Displaying message set(msgtext,'String','New element flow has been added!') end %End of addelmndflowfcn %Start of rmvelmnflow1fcn function rmvelmnflow1fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+1,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview),... ' deleted!']) end %End of rmvelmnflow1fcn %Start of rmvelmnflow2fcn function rmvelmnflow2fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+2,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+1),... ' deleted!']) end %End of rmvelmnflow2fcn %Start of rmvelmnflow3fcn function rmvelmnflow3fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+3,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+2),... ' deleted!']) end %End of rmvelmnflow3fcn %Start of rmvelmnflow4fcn function rmvelmnflow4fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+4,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+3),... ' deleted!']) end %End of rmvelmnflow4fcn %Start of rmvelmnflow5fcn function rmvelmnflow5fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+5,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+4),... ' deleted!']) end %End of rmvelmnflow5fcn %Start of rmvelmnflow6fcn function rmvelmnflow6fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+6,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+5),... ' deleted!']) end %End of rmvelmnflow6fcn %Start of rmvelmnflow7fcn function rmvelmnflow7fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+7,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+6),... ' deleted!']) end %End of rmvelmnflow7fcn %Start of rmvelmnflow8fcn function rmvelmnflow8fcn(~,~) %Deleting elementary flow elmnflowcount=elmnflowcount-1; elmnflow(indexview+8,:)=[]; if (indexview+7>elmnflowcount)&&(indexview>1) indexview=indexview-1; end %Updating panel updatepanel updateplot %Displaying message set(msgtext,'String',['Elementary flow no ',... num2str(indexview+7),... ' deleted!']) end %End of rmvelmnflow8fcn %Start of flowvisfcn function flowvisfcn(~,~) xp=[x(:,1);x(:,end);x(1,:)';x(end,:)']; yp=[y(:,1);y(:,end);y(1,:)';y(end,:)']; %Initiating particle scatter pscatter=scatter(xp,yp,'filled','k'); %Displaying message set(msgtext,'String','Performing flow visualization...') %Visualizing particle movement for timestep=1:1000 %Finding particle speed up=interp2(x,y,u,xp,yp); vp=interp2(x,y,v,xp,yp); xp=xp+up*0.05*zoomstate; yp=yp+vp*0.05*zoomstate; %Deleting out-of-bounds particles xnan=find(isnan(xp)); xp(xnan)=[]; yp(xnan)=[]; ynan=find(isnan(yp)); xp(ynan)=[]; yp(ynan)=[]; %Deleting stagnation point particles stag=(up==0)&(vp==0); xp(stag)=[]; yp(stag)=[]; %Plotting particle movement set(pscatter,'XData',xp) set(pscatter,'YData',yp) pause(0.001) %Anticipating lost of all particles if isempty(xp) break end end clear pscatter %Displaying message set(msgtext,'String','Flow visualization finished!') end %End of flowvisfcn %Start of flowviscustomfcn function flowviscustomfcn(~,~) %Displaying message set(msgtext,... 'String','Select multiple points on plot to place particle...') %Prompting user input [xp,yp]=getpts; %Initiating particle scatter pscatter=scatter(xp,yp,'filled','k'); %Displaying message set(msgtext,'String','Performing flow visualization') %Visualizing particle movement for timestep=1:1000 %Finding particle speed up=interp2(x,y,u,xp,yp); vp=interp2(x,y,v,xp,yp); xp=xp+up*0.01; yp=yp+vp*0.01; %Deleting out-of-bounds particles xnan=find(isnan(xp)); xp(xnan)=[]; yp(xnan)=[]; ynan=find(isnan(yp)); xp(ynan)=[]; yp(ynan)=[]; %Deleting stagnation point particles stag=(up==0)&(vp==0); xp(stag)=[]; yp(stag)=[]; %Plotting particle movement set(pscatter,'XData',xp) set(pscatter,'YData',yp) pause(0.001) %Anticipating lost of all particles if isempty(xp) break end end set(pscatter,'XData',[]) set(pscatter,'YData',[]) clear pscatter %Displaying message set(msgtext,'String','Flow visualization finished!') end %End of flowviscustomfcn %Start of switchviewfcn function switchviewfcn(~,~) %Switching view plottype=3-plottype; updateplot %Displaying message if plottype==1 set(msgtext,'String','Switching to absolute plot!') else set(msgtext,'String','Switching to normalized plot!') end end %End of switchviewfcn %Start of savedatafcn function savedatafcn(~,~) %Saving plot data save('PotentialFlowSimData.mat','x','y','u','v','V'); %Displaying message set(msgtext,'String',['Data successfully saved to',... ' PotentialFlowSimData.mat!']) end %End of savedatafcn %Start of saveplotfcn function saveplotfcn(~,~) %Generating dummy axes tempfig=figure(); tempaxes=axes('Parent',tempfig); readpanel updateplot %Saving plot figure saveas(tempaxes,'PotentialFlowSim.fig') %Deleting dummy axes close(tempfig) %Displaying message set(msgtext,'String',['Data successfully saved to',... ' PotentialFlowSim.fig!']) end %End of saveplotfcn %CodeEnd------------------------------------------------------------------- end
airfoil codes/SourcePanelMethod2D.m
function [lambda,cp,cl,cd]=SourcePanelMethod2D(xp,yp,uf,vf,dotcount) %-------------------------------------------------------------------------- %SourcePanelMethod2D %Version 1.10 %Created by Stepen %Created 13 August 2011 %-------------------------------------------------------------------------- %SourcePanelMethod2D performs 2D source panel method the given geometry and %finds the source strength at all panel that satisfies boundary condition %at all panels. SourcePanelMethod2D requires %CalculateVelocityFromSourcePanel2D and CalculateVelocityFromSource2D %m-file function in the same directories to run. %-------------------------------------------------------------------------- %Syntax: %lambda=SourcePanelMethod2D(x,y,uf,vf,dotcount) %Input argument: %- xp (m x 1 num) specifies the x axis location of object's panels. To % define closed loop panels, xp must be arranged in clockwise direction % (for normal vector determination) and the first value of xp must be the % same with the last value of xp. %- yp (m x 1 num) specifies the y axis location of object's panels. To % define closed loop panels, yp must be arranged in clockwise direction % (for normal vector determination) and the first value of yp must be the % same with the last value of yp. %- uf (1 x 1 num) specifies the horizontal velocity component of % freestream. %- vf (1 x 1 num) specifies the vertital velocity component of freestream. %- dotcount (1 x 1 int) specifies the number of vertexes to be generated on % each the panel. %Output argument: %- lambda (p x 1 num) specifies the source strength per unit length at the % corresponding panel. %- cp (p x 1 num) specifies the pressure coefficient on the corresponding % panel. %- cl (1 x 1 num) specifies the lift coefficient of the object. %- cd (1 x 1 num) specifies the drag coefficient of the object. %-------------------------------------------------------------------------- %Example: %[lambda,cp,cl,cd]=SourcePanelMethod2D(10*sind(0:15:360)',... % 10*cosd(0:15:360)',... % 10,0,11) % calculates source strength distribution and cp over a circle and % displays its resulting flow field. %[lambda,cp,cl,cd]=SourcePanelMethod2D([0;0.25;1;0.25;0],... % [0;0.15;0;-0.1;0],... % 10,0,11) % calculates source strength distribution and cp over a simple diamond % wedge airfoil and displays its resulting flow field. %-------------------------------------------------------------------------- %CodeStart----------------------------------------------------------------- %Checking input object's panels x_size=size(xp); y_size=size(yp); if (x_size(2)~=1)||(y_size(2)~=1) error('Input panel location must be a i x 1 array!') end if (numel(xp)~=numel(yp)) error('Input xp and yp do not have the same size!') end panelcount=numel(xp)-1; %Checking input u and v if numel(uf)~=1 error('Freestream horizontal velocity component must be a scalar!') end if numel(vf)~=1 error('Freestream vertital velocity component must be a scalar!') end %Checking dotcount if numel(dotcount)~=1 error('Number of vertex must be scalar!') end if (mod(dotcount,1~=0))||(dotcount<0) error('Number of vertex must be positive integer!') end %Preallocating array for speed xcp=zeros(panelcount,1); ycp=zeros(panelcount,1); ucp=zeros(panelcount,1); vcp=zeros(panelcount,1); Vcp=zeros(panelcount,1); cp=zeros(panelcount,1); cl=0; cd=0; u=zeros(panelcount); v=zeros(panelcount); l=zeros(panelcount,1); sine=zeros(panelcount,1); cosine=zeros(panelcount,1); nsine=zeros(panelcount,1); ncosine=zeros(panelcount,1); InfluenceMatrix=zeros(panelcount); FreeStreamMatrix=zeros(panelcount,1); %Calculating panel's properties for i=1:1:panelcount %Calculating panel's control point location xcp(i,1)=0.5*(xp(i)+xp(i+1)); ycp(i,1)=0.5*(yp(i)+yp(i+1)); %Calculating panel's length l(i,1)=norm([xp(i+1)-xp(i),yp(i+1)-yp(i)]); %Calculating panel's orientation and its normal sine(i,1)=(yp(i+1)-yp(i))/l(i); cosine(i,1)=(xp(i+1)-xp(i))/l(i); nsine(i,1)=cosine(i,1); ncosine(i,1)=-sine(i,1); end %Calculating source influence on all panel and generating InfluenceMatrix for i=1:1:panelcount for j=1:1:panelcount [u(i,j),v(i,j)]=CalculateVelocityFromSourcePanel2D(1,[xp(j);xp(j+1)],[yp(j);yp(j+1)], xcp(i),ycp(i),dotcount); InfluenceMatrix(i,j)=(u(i,j)*(ncosine(i)))+(v(i,j)*(nsine(i))); end end %Generating FreeStreamMatrix for i=1:1:panelcount FreeStreamMatrix(i,1)=-(uf*(ncosine(i)))-((vf*(nsine(i)))); end %Calculating LambdaMatrix LambdaMatrix=(InfluenceMatrix^-1)*FreeStreamMatrix; lambda=LambdaMatrix; %Calculating velocity and pressure coefficient on all panel Vf=norm([uf,vf]); for i=1:1:panelcount for j=1:1:panelcount [utemp,vtemp]=CalculateVelocityFromSourcePanel2D... (LambdaMatrix(j),[xp(j);xp(j+1)],... [yp(j);yp(j+1)],... xcp(i),ycp(i),dotcount); ucp(i)=ucp(i)+utemp; vcp(i)=vcp(i)+vtemp; end ucp(i)=ucp(i)+uf; vcp(i)=vcp(i)+vf; Vcp(i)=norm([ucp(i),vcp(i)]); cp(i)=1-((Vcp(i)/Vf)^2); end %Calculating lift and coefficient for i=1:1:panelcount cl=cl+(cp(i)*nsine(i)*l(i)); cd=cd+(cp(i)*ncosine(i)*l(i)); end %Plotting result %Generating meshgrid xmin=min(xp);xmax=max(xp); dx=xmax-xmin; ymin=min(yp);ymax=max(yp); dy=xmax-xmin; [xg,yg]=meshgrid(xmin-dx:dx/10:xmax+dx,ymin-dy:dy/10:ymax+dy); [inpanel,onpanel]=inpolygon(xg,yg,xp,yp); grid_size=size(xg); %Calculating velocity at all grid points ug=zeros(grid_size); vg=zeros(grid_size); for i=1:1:grid_size(1) for j=1:1:grid_size(2) if (inpanel(i,j)==1)||(onpanel(i,j)==1) ug(i,j)=NaN; yg(i,j)=NaN; else for k=1:1:panelcount [utemp,vtemp]=CalculateVelocityFromSourcePanel2D... (LambdaMatrix(k),[xp(k);xp(k+1)],... [yp(k);yp(k+1)],... xg(i,j),yg(i,j),... dotcount); ug(i,j)=ug(i,j)+utemp; vg(i,j)=vg(i,j)+vtemp; end ug(i,j)=ug(i,j)+uf; vg(i,j)=vg(i,j)+vf; end end end %Generating plot hold on axis equal quiver(xg,yg,ug,vg) patch(xp,yp,'w') %CodeEnd------------------------------------------------------------------- end
