6.2 粘性流体涡量-流函数法求解流场

clear;%参数设置Re=10;                        L=1;                          n=100;dh=L/n;psi=zeros(n+1,n+1);xi=zeros(n+1,n+1);rho=1.0;   for k=1:100000    err=0;    for i=2:n    xi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;    xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2;    end         for j=2:n    xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2;    xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2;    end        for i=2:n    for j=2:n        u(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh);        v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));        err1=(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+xi(i,j)*dh^2)/4-psi(i,j);        psi(i,j)=psi(i,j)+rho*err1;        err2=(xi(i+1,j)+xi(i-1,j)+xi(i,j+1)+xi(i,j-1))/4 ...            -Re*dh*(u(i,j)*(xi(i+1,j)-xi(i-1,j))+v(i,j)*(xi(i,j+1)-xi(i,j-1)))/8-xi(i,j);        xi(i,j)=xi(i,j)+rho*err2;        temp=max(abs(err1),abs(err2));        if err<temp        err=temp;        end    end    end    if err<1e-6    break;endend  kerrrho%psicontour(psi,100);