请教有关流体计算的编程问题

krzn

运行后画出的图线明显不符合实际啊；
虚心求教。。。
format
v=10;
R=0.5;
m=1;
c1=-1;
for b=0:0.1:2*pi
    circlex(m)=R*sin(b);
    circley(m)=R*cos(b);      
    r1(m)=(((c1/(v*sin(b)))^2+4*R^2)^0.5-(c1/(v*sin(b)))^2)/2;
    x1(m)=(r1(m)/(1+(tan(b))^2))^0.5;
    y1(m)=tan(b)*(r1(m)/(1+(tan(b))^2))^0.5
    m=m+1;
end
c2=-2;
m=1;
for b=0:0.1:2*pi
    r2(m)=(((c2/(v*sin(b)))^2+4*R^2)^0.5-(c2/(v*sin(b)))^2)/2;
    x2(m)=(r2(m)/(1+(tan(b))^2))^0.5;
    y2(m)=tan(b)*(r2(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
c3=-3;
m=1;
for b=0:0.1:2*pi
    r3(m)=(((c3/(v*sin(b)))^2+4*R^2)^0.5-(c3/(v*sin(b)))^2)/2;
    x3(m)=(r3(m)/(1+(tan(b))^2))^0.5;
    y3(m)=tan(b)*(r3(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
c4=-4;
m=1;
for b=0:0.1:2*pi
    r4(m)=(((c4/(v*sin(b)))^2+4*R^2)^0.5-(c4/(v*sin(b)))^2)/2;
    x3(m)=(r4(m)/(1+(tan(b))^2))^0.5;
    y3(m)=tan(b)*(r4(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
title('理想流体绕过圆柱流动的流线分布图')
plot(circlex,circley,x1,y1,x2,y2,x3,y3)

krzn

format
v=5;
R=0.5;
c1=1;
m=1;
for b=0:0.1:(pi/2)
    r1(m)=(((c1/(v*sin(b)))^2+4*(R^2))^0.5+c1/(v*sin(b)))/2;
    x1(m)=(r1(m)/(1+(tan(b))^2))^0.5;
    y1(m)=tan(b)*(r1(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
c2=3;
m=1;
for b=0:0.1:(pi/2)
    r2(m)=(((c2/(v*sin(b)))^2+4*R^2)^0.5+c2/(v*sin(b)))/2;
    x2(m)=(r2(m)/(1+(tan(b))^2))^0.5;
    y2(m)=tan(b)*(r2(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
c3=5;
m=1;
for b=0:0.1:(pi/2)
    r3(m)=(((c3/(v*sin(b)))^2+4*R^2)^0.5+c3/(v*sin(b)))/2;
    x3(m)=(r3(m)/(1+(tan(b))^2))^0.5;
    y3(m)=tan(b)*(r3(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
c4=7;
m=1;
for b=0:0.1:(pi/2)
    r4(m)=(((c4/(v*sin(b)))^2+4*R^2)^0.5+c4/(v*sin(b)))/2;
    x4(m)=(r4(m)/(1+(tan(b))^2))^0.5;
    y4(m)=tan(b)*(r4(m)/(1+(tan(b))^2))^0.5;
    m=m+1;
end
m=m+1;
for b=0:0.1:2*pi
    circlex(m)=R*sin(b);
    circley(m)=R*cos(b);   
    m=m+1;
end
plot(x1,y1,-x1,y1,x1,-y1,-x1,-y1,x2,y2,-x2,y2,x2,-y2,-x2,-y2,x3,y3,-x3,y3,x3,-y3,-x3,-y3,x4,y4,-x4,y4,x4,-y4,-x4,-y4,circlex,circley)
axis ('equal','auto')
title('理想流体绕过圆柱流动的流线分布图')
gtext('圆柱表面')
呵呵，这是改过的，那几个没有错啊
只是觉得for语句太多了，可以分段将数组中元素取出来画线吗？

xjzuo

for b=...的用法相当别扭．

建议去掉循环，直接用向量形式计算．
plot可用一个循环实现，每次加一条语句 hold on.

krzn

我想用二维数组代替前面的循环
但是画线时用plot 可以只取其中一列来画吗？

xjzuo

"其中一列" 是指什么?

krzn

上面循环中的x1,x2,x3我想用x(m,n)来保存
可是画线时我写的是
plot(x(:,1),y(:,1))
提示出错，改怎么样正确表示呢？
万分感谢～～

xjzuo

再次建议去掉循环:

%%%%%%%%%%%%%%%%%%%
...
c1=...;
c2=...;
...
b=0:0.1:2*pi;      %%去掉循环,改为向量
x=zeros(3,length(b));
y=zeros(3,length(b));

circlex=R*sin(b);
circley=R*cos(b);      
r1=(((c1./(v*sin(b))).^2+4*R^2).^0.5-(c1./(v*sin(b))).^2)/2;
x(1,:)=...
...

for k=1:3
plot(x(k,:),y(k,:),-x(k,:),y(k,:),x(k,:),-y(k,:),-x(k,:),-y(k,:))
hold on
end

hold on
plot(circlex,circley,'LineWidth',2)
...
%%%%%%%%%%%%%%%%%%%

另: 请勿重复发贴!!!(已删除重复贴)

[ 本帖最后由 xjzuo 于 2006-11-25 17:12 编辑 ]

krzn

谢谢了
不太会用向量来代替双循环啊
算了，还是用笨办法吧，
谢谢

happy

原帖由 krzn 于 2006-11-25 17:21 发表
谢谢了
不太会用向量来代替双循环啊
算了，还是用笨办法吧，
谢谢

流体计算如果不改用向量那是很考验人的耐性的哦