计算流体力学RAI(2)(2)

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])

流

函数图像对比：

速度图像对比：