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