clear;
%參數(shù)設(shè)置
re=input('請(qǐng)輸入雷諾數(shù)的值\nRe=');                        %雷諾數(shù)取100，500，1000
nu=900e-6;                    %粘性系數(shù)
L=0.1;                          %空穴幾何尺寸
I=50;
J=50;
dx=L/I;
dy=L/J;
U=re*nu/L;                    %平板速度 
dt=1e-4;					    %時(shí)間步長
%給定初值
u(1:I+1,J+1)=U;					%給定u的初值，
u(1:I+1,1:J)=0;
v=zeros(I+1,J+1);				%給定v的初值
omega=zeros(I+1,J+1);           %給定具有梯度的渦量初場
omega(1:I+1,J+1)=3*U/(2*dy);
ps=zeros(I+1,J+1);            %給定流函數(shù)初場
omega1=ones(I+1,J+1);
ps1=ones(I+1,J+1);
tol1=ps-ps1;
tol2=omega-omega1;
t=0;
n=0;
while norm(tol1)>1e-4|norm(tol2)>1e-4
    n=n+1;
    omega1=omega;
    ps1=ps;
    t=t+dt;
    for i=2:I
        for j=2:J
            omega(i,j)=omega(i,j)-dt*(u(i,j)*(omega(i+1,j)-omega(i-1,j))/(2*dx)+v(i,j)*(omega(i,j+1)-omega(i,j-1))/(2*dy))+dt/re*((omega(i+1,j)-2*omega(i,j)+omega(i-1,j))/dx^2+(omega(i,j+1)-2*omega(i,j)+omega(i,j-1))/dy^2);
        end
    end
    for i=2:I
        for j=2:J
            ps(i,j)=ps(i,j)+0.25*(ps(i+1,j)+ps(i-1,j)+ps(i,j+1)+ps(i,j-1)-4*ps(i,j)+dx^2*omega(i,j));
        end
    end
    for i=2:I
         for j=2:J
             u(i,j)=(ps(i,j+1)-ps(i,j-1))/(2*dy);
             v(i,j)=-(ps(i+1,j)-ps(i-1,j))/(2*dx);
         end
     end
    for j=2:J
       omega(1,j)=1/(2*dx^2)*(ps(3,j)-8*ps(2,j));
       omega(I+1,j)=1/(2*dx^2)*(ps(I-1,j)-8*ps(I,j));
    end
    for i=2:I                                                                          %由流函數(shù)值和內(nèi)點(diǎn)上的渦量值決定渦量的邊界值
       omega(i,1)=1/(2*dy^2)*(ps(i,3)-8*ps(i,2));
       omega(i,J+1)=1/(2*dy^2)*(ps(i,J-1)-8*ps(i,J))-3*U/dy;
    end
tol1=ps-ps1;
tol2=omega-omega1;
end

for i=1:I+1						%求各網(wǎng)格點(diǎn)坐標(biāo)
    for j=1:J+1
        x(i,j)=dx*(i-1);
        y(i,j)=dy*(j-1);
    end
end
for i=1:I
    counter(i)=ps(i,i);
end
counter(i+1)=ps(j-1,2);
counter(i+2)=ps(j-2,3);
counter(i+3)=ps(j-3,3);
counter(i+4)=ps(j-4,3);
counter(i+5)=ps(j-5,3);
contour(x,y,ps,counter)