计算流体力学第二次实验
小雷诺数：Re=100
程序：
f1=zeros(101,101); f2=zeros(101,101); w1=zeros(101,101); w2=zeros(101,101); re=100; h=0.01; u0=0.1; for i=2:100; w1(1,i)=-2*u0/h; end;
for i=2:100; for j=2:100;
f2(i,j)=(f1(i-1,j)+f1(i+1,j)+f1(i,j-1)+f1(i,j+1)+h*h*w1(i,j))/4; end end
a=zeros(99,99); b=zeros(99,99); for i=1:99; for j=1:99;
a(i,j)=(f2(i,j+1)-f2(i+2,j+1))/4; b(i,j)=(f2(i+1,j+2)-f2(i+1,j))/4; end end
for i=2:100; for j=2:100;
w2(i,j)=((1-a(i-1,j-1)*re)*w1(i,j+1)+(1+a(i-1,j-1)*re)*w1(i,j-1)+(1+b(i-1,j-1)*re)*w1(i-1,j)+(1-b(i-1,j-1)*re)*w1(i+1,j))/4; end; end;
for i=2:100;
w2(i,101)=-2*f2(i,100)/(h*h); w2(i,1)=-2*f2(i,2)/(h*h); w2(1,i)=-2*f2(2,i)/(h*h)-2*u0/h; w2(101,i)=-2*f2(100,j)/(h*h);
end; z=max(w1-w2,w2-w1); c=max(z(:,:));
d=max(c);
while(d>=0.0000001) f1=f2;w1=w2;
for i=2:100; for j=2:100;
f2(i,j)=(f1(i-1,j)+f1(i+1,j)+f1(i,j-1)+f1(i,j+1)+h*h*w1(i,j))/4; end;
end; for i=1:99; for j=1:99;
a(i,j)=(f2(i,j+1)-f2(i+2,j+1))/4; b(i,j)=(f2(i+1,j+2)-f2(i+1,j))/4; end;
end;
for i=2:100; for j=2:100;
w2(i,j)=((1-a(i-1,j-1)*re)*w1(i,j+1)+(1+a(i-1,j-1)*re)*w1(i,j-1)+(1+b(i-1,j-1)*re)*w1(i-1,j)+(1-b(i-1,j-1)*re)*w1(i+1,j))/4; end; end;
for i=2:100;
w2(i,101)=-2*f2(i,100)/(h*h); w2(i,1)=-2*f2(i,2)/(h*h); w2(1,i)=-2*f2(2,i)/(h*h)-2*u0/h; w2(101,i)=-2*f2(100,j)/(h*h);
end; z=max(w1-w2,w2-w1); c=max(z(:,:)); d=max(c); end;
u=zeros(101,101); v=zeros(101,101); for i=2:100; for j=2:100;
u(i,j)=(f2(i-1,j)-f2(i+1,j))/(2*h); v(i,j)=(f2(i,j-1)-f2(i,j+1))/(2*h); end; end;
for i=2:100; u(1,i)=0.1; end;
y=1:-0.01:0;x=0:0.01:1;
figure(1) contour(x,y,f2) figure(2) quiver(x,y,u,v) axis([0 1 0 1])
流函数图像对比：
速度图像对比：
大雷诺数：Re=400
程序：
f1=zeros(101,101); f2=zeros(101,101); w1=zeros(101,101); w2=zeros(101,101); re=400; h=0.01; u0=0.1; for i=2:100; w1(1,i)=-2*u0/h; end;
for i=2:100; for j=2:100;
f2(i,j)=(f1(i-1,j)+f1(i+1,j)+f1(i,j-1)+f1(i,j+1)+h*h*w1(i,j))/4; end end
ap=zeros(99,99);ae=zeros(99,99);aw=zeros(99,99);an=zeros(99,99);as=zeros(99,99);u=zeros(99,99);v=zeros(99,99); for i=1:99; for j=1:99;
u(i,j)=(f2(i,j+1)-f2(i+2,j+1))/(2*h); v(i,j)=-(f2(i+1,j+2)-f2(i+1,j))/(2*h); end end
for i=1:99;