计算流体力学RAI(2)

计算流体力学第二次实验

小雷诺数：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;

for j=1:99;

ap(i,j)=max(u(i,j),-u(i,j))/h+max(v(i,j),-v(i,j))/h+4/(h*h*re); aw(i,j)=max(0,u(i,j)/h)+1/(h*h*re); ae(i,j)=max(0,-u(i,j)/h)+1/(h*h*re); an(i,j)=max(0,-v(i,j)/h)+1/(h*h*re); as(i,j)=max(0,v(i,j)/h)+1/(h*h*re); end; end;

for i=2:100; for j=2:100;

w2(i,j)=(ae(i-1,j-1)*w1(i,j+1)+aw(i-1,j-1)*w1(i,j-1)+an(i-1,j-1)*w1(i-1,j)+as(i-1,j-1)*w1(i+1,j))/ap(i-1,j-1); 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.00001) 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;

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; for j=1:99;

ap(i,j)=max(u(i,j),-u(i,j))/h+max(v(i,j),-v(i,j))/h+4/(h*h*re); aw(i,j)=max(0,u(i,j)/h)+1/(h*h*re); ae(i,j)=max(0,-u(i,j)/h)+1/(h*h*re); an(i,j)=max(0,-v(i,j)/h)+1/(h*h*re);

as(i,j)=max(0,v(i,j)/h)+1/(h*h*re); end; end;

for i=2:100; for j=2:100;

w2(i,j)=(ae(i-1,j-1)*w1(i,j+1)+aw(i-1,j-1)*w1(i,j-1)+an(i-1,j-1)*w1(i-1,j)+as(i-1,j-1)*w1(i+1,j))/ap(i-1,j-1); 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])

流

函数图像对比：

速度图像对比：

本文来源：https://www.bwwdw.com/article/mwb6.html