www.gusucode.com > MATLAB2008应用程序接口编程技术源码程序 > MATLAB2008应用程序接口编程技术源码程序/code/第4章/4.5/示例6/Areoplot.m
function Areoplot %%%定义常数数值 close all V_i = 25; G = 10; a = 1 ; %%%半径 c =-a*2; b =a*2; n =a*50; % %%迭代次数 [x,y]=meshgrid([c:(b-c)/n:b],[c:(b-c)/n:b]'); warning off %准备数据 for i=1:length(x); for k=1:length(x); if sqrt(x(i,k).^2+y(i,k).^2)<a; x(i,k)=0; y(i,k)=0; end end end %%%定义极坐标的数值 rho=sqrt(x.^2+y.^2); theta=atan2(y,x); %%%创建流线型的函数 z=V_i.*sin(theta).*rho.*(1-(a^2./(rho.^2)))-G*log(rho)/(2*pi); % 创建图形 n=100; r=ones(1,n+1)*a; t=[0:2*pi/n:2*pi]; %%%流线型图形contour(x,y,z,25) contour(x,y,z,25) hold on polar(t,r,'-k') axis square title('Stream Lines') grid off figure(2) contour(x,y,z,15) %建立环绕在圆柱体的矢量场 x=[-a*2:a/3:a*2]; [x]=meshgrid(x); y=x'; for i=1:length(x); for k=1:length(x); if sqrt(x(i,k).^2+y(i,k).^2)<a; x(i,k)=0; y(i,k)=0; end end end r=sqrt(x.^2+y.^2); theta=atan2(y,x); ur=V_i*cos(theta).*(1-a^2./(r.^2)); ut=-V_i*sin(theta).*(1+a^2./(r.^2))+G./(2*pi*r); u=ur.*cos(theta)-ut.*sin(theta); v=ur.*sin(theta)+ut.*cos(theta); %%%创建填充后的图形 t_r = 0:.1:2*pi; xxx = a*cos(t_r); yyy = a*sin(t_r); %%%填充后的图形和矢量场图形 figure(2) hold on quiver(x,y,u,v) fill(xxx,yyy,'y') axis square title('Speed Vectors') grid off warning on t=0:.1:2*pi; cp = 1 - 4*sin(t).^2 + 2* G / (pi*a*V_i) *sin(t) - (2* G/ (pi*a*V_i) )^2 ; cp_sim = 1 - 4*sin(t).^2 ; L = - 1.225*V_i*G; L = num2str(L); L = strcat('Kutta Joukowski Lift: ',L,' [N]'); figure(3) plot(t,cp,t,cp_sim,'--r') axis([0 2*pi min(cp) max(cp_sim)]) title('Pressure coefficient around the surface (standard air density)') xlabel('Theta (angle with orizonthal)') ylabel('C_p') legend('Lifting solution','Symmetrical solution') grid on