%--------------------------------------------------------%
%???????SIMPLE?Algorithm?-?staggered?mesh
%???????Incompressible?Lid-driven?Flow?in?Cavity
%--------------------------------------------------------%
close?all
clear
clc
disp（‘SIMPLE?Algorithm‘）;

%%?Setting
n=51;
m=n;
dx=1/（n-1）;
dy=1/（m-1）;
Re=100;
dt=0.001;
tmax=10000;
k1max=300;
k2max=300;
alpha=1;
alpha_p=1;
relax=1
r_p_critial=1e-2;
r_u_critial=1e-2;
r_v_critial=1e-2;
r_pc_critial=1e-2;

%%?Initialization
u=zeros（n-1m）;
u（:m）=1;??????%?lid
v=zeros（nm-1）;
p=zeros（nm）;
r_p=zeros（tmax1）;
r_u=zeros（tmax1）;
r_v=zeros（tmax1）;

%%?Loops
tic
fprintf（‘Iter\tP-residual\t\tU-residual\t\tV-residual\n‘）;
A1=1/dx^2;
A2=1/dy^2;
A3=1/（2*A1+2*A2）;
A=1/dt+2/Re*（A1+A2）;
for?k=1:tmax
????%%?Compute?tentative?velocities?uv
????u1=zeros（n-1m）;
????u1（:m）=1;
????v1=zeros（nm-1）;
????for?k1=1:k1max
????????for?i=2:n-2
????????????for?j=2:m-1
????????????????vn=（v（ij）+v（i+1j））/2;
????????????????vs=（v（ij-1）+v（i+1j-1））/2;
????????????????a1=-u（i+1j）/（2*dx）+A1/Re;
????????????????a2=u（i-1j）/（2*dx）+A1/Re;
????????????????a3=-vn/（2*dy）+A2/Re;
????????????????a4=vs/（2*dy）+A2/Re;
????????????????u1（ij）=（1-relax）*u1（ij）+relax*（u（ij）/dt+a1*u1（i+1j）+a2*u1（i-1j）+a3*u1（ij+1）+a4*u1（ij-1）-（p（i+1j）-p（ij））/dx）/A;
????????????end
????????end
????????r_u1=norm（u-u1）;
????????if?r_u1????end
????for?k1=1:k1max
????????for?i=2:n-1
????????????for?j=2:m-2
????????????????ue=（u（ij）+u（ij+1））/2;
????????????????uw=（u（i-1j）+u（i-1j+1））/2;
????????????????b1=-v（ij+1）/（2*dy）+A2/Re;
????????????????b2=v（ij-1）/（2*dy）+A2/Re;
????????????????b3=-ue/（2*dx）+A1/Re;
????????????????b4=uw/（2*dx）+A1/Re;
????????????????v1（ij）=（1-relax）*v1（ij）+relax*（v（ij）/dt+b1*v1（ij+1）+b2*v1（ij-1）+b3*v1（i+1j）+b4*v1（i-1j）-（p（ij+1）-p（ij））/dy）/A;
????????????end
????????end
????????r_v1=norm（v-v1