CFD Online Discussion Forums - HELP!Problems on SIMPLE method

I'm a master student and trying to do a laminar simulation by SIMPLE method.
I use matlab and the result seems not good.
The most serious problem I think is,It seems very difficult to get convergence solution on Pressure correction.(I have tried Jacobi and ADI)NOT ONLY the absolute value of pressure grows with iteration，but ALSO the pressure drop grows,while the contour shape of pressure drop seems NOT strange.

The code is writing by matlab,based on staggered grid by Patankar.(I didn't do residual check in the program just for easy to read.The fact is even I check for residual,it can't get a convergence result)

clear all;
%initialization
nx=30;
ny=10;
%velcity field
u=zeros(ny,nx+1,3);
v=zeros(ny+1,nx,3);
%pressure
p=zeros(ny,nx,3);
%co. for momentum equation
aeu=zeros(ny,nx);
awu=zeros(ny,nx);
asu=zeros(ny,nx);
anu=zeros(ny,nx);
apu=zeros(ny,nx);
aev=zeros(ny,nx);
awv=zeros(ny,nx);
asv=zeros(ny,nx);
anv=zeros(ny,nx);
apv=zeros(ny,nx);
%co. for pressure correction equation
pce=zeros(ny,nx);
pcw=zeros(ny,nx);
pcs=zeros(ny,nx);
pcn=zeros(ny,nx);
aep=zeros(ny,nx);
awp=zeros(ny,nx);
asp=zeros(ny,nx);
anp=zeros(ny,nx);
app=zeros(ny,nx);
pc=zeros(ny,nx,3);
%mass unbalance
mp=zeros(ny,nx);
errp=zeros(ny,nx);
raw=1;
dx=0.1;
dy=0.1;
d=1;
u(2:ny-1,:,:)=30;
for q=1:5000
for ll=1:30
for i=2:ny-1
for j=2:nx
aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j-1,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i,j-1,1)+0.5*v(i,j,1)));
anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j-1,1)+0.5*v(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j-1,1)+0.5*v(i+1,j,1)));
apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i,j-1,1)-p(i,j,1))*dy)/apu(i,j);
end
end
for i=2:ny
for j=2:nx-1
aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i-1,j+1,1)+0.5*u(i,j+1,1))/d))^5)+max(0,-raw*dy*(0.5*u(i-1,j+1,1)+0.5*u(i,j+1,1)));
awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i-1,j,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i-1,j,1)+0.5*u(i,j,1)));
asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j,1)-p(i,j,1))*dy)/apv(i,j);
end
end
end

for ll=1:100
for i=2:ny-1
for j=2:nx-1
cep(i,j)=raw*dy*(0.375*(u(i,j,1)+u(i,j+1,1))+(u(i-1,j,1)+u(i+1,j,1)+u(i-1,j+1,1)+u(i+1,j+1,1))/16);
cwp(i,j)=raw*dy*(0.375*(u(i,j,1)+u(i,j-1,1))+(u(i-1,j,1)+u(i+1,j,1)+u(i-1,j-1,1)+u(i+1,j-1,1))/16);
cnp(i,j)=raw*dx*(0.375*(v(i+1,j,1)+v(i,j,1))+(v(i, j-1,1)+v(i,j+1,1)+v(i+1,j-1,1)+v(i+1,j+1,1))/16);
csp(i,j)=raw*dx*(0.375*(v(i-1,j,1)+v(i,j,1))+(v(i,j-1,1)+v(i,j+1,1)+v(i-1,j-1,1)+v(i-1,j+1,1))/16);
pcep(i,j)=cep(i,j)/d;
pcwp(i,j)=cwp(i,j)/d;
pcnp(i,j)=cnp(i,j)/d;
pcsp(i,j)=csp(i,j)/d;
aep(i,j)=d*(max(0,(1-0.1*abs(pcep(i,j)))^5)+max(-pcep(i,j),0));
awp(i,j)=d*(max(0,(1-0.1*abs(pcwp(i,j)))^5)+max(pcwp(i,j),0));
anp(i,j)=d*(max(0,(1-0.1*abs(pcnp(i,j)))^5)+max(-pcnp(i,j),0));
asp(i,j)=d*(max(0,(1-0.1*abs(pcsp(i,j)))^5)+max(pcsp(i,j),0));
app(i,j)=aep(i,j)+awp(i,j)+asp(i,j)+anp(i,j);
end
end
for i=2:ny-1
for j=2:nx-1
aepc(i,j)=2*raw*dy^2/(app(i,j)+app(i,j+1));
awpc(i,j)=2*raw*dy^2/(app(i,j)+app(i,j-1));
anpc(i,j)=2*raw*dy^2/(app(i,j)+app(i+1,j));
aspc(i,j)=2*raw*dy^2/(app(i,j)+app(i-1,j));
end
end
for i=2:ny-1
for j=2:nx-1
mp(i,j)=raw*dx*((u(i,j+1,1)-u(i,j-1,1))/2+(v(i+1,j,1)-v(i-1,j,1))/2);
end
end

for i=2:ny-1
aepc(i,nx-1)=raw*dy^2/app(i,nx-1);
end
for i=2:ny-1
awpc(i,2)=raw*dy^2/app(i,2);
end
for j=2:nx-1
anpc(ny-1,j)=raw*dy^2/app(ny-1,j);
end
for j=2:nx-1
aspc(2,j)=raw*dy^2/app(2,j);
end

for i=2:ny-1
for j=2:nx-1
appc(i,j)=aepc(i,j)+awpc(i,j)+anpc(i,j)+aspc(i,j);
end
end

for i=2:ny-1
for j=2:nx-1
pc(i,j,1)=(pc(i,j+1,1)*aepc(i,j)+pc(i,j-1,1)*awpc(i,j)+pc(i+1,j,1)*anpc(i,j)+pc(i-1,j,1)*aspc(i,j)-mp(i,j))/appc(i,j);
end
end
end

for i=2:ny-1
pc(i,1,1)=pc(i,2,1);
pc(i,nx,1)=pc(i,nx-1,1);
end
for j=2:nx-1
pc(1,j,1)=pc(2,j,1);
pc(ny,j,1)=pc(ny-1,j,1);
end
pc(1,1,1)=(pc(1,2,1)+pc(2,1,1))/2;
pc(1,nx,1)=(pc(1,nx-1,1)+pc(2,nx,1))/2;
pc(ny,1,1)=0.5*(pc(ny,2,1)+pc(ny-1,1,1));
pc(ny,nx,1)=0.5*(pc(ny-1,nx,1)+pc(ny,nx-1,1));

p(:,:,1)=p(:,:,1)+0.01*pc(:,:,1);
for i=2:ny-1
for j=2:nx
u(i,j,1)=u(i,j,1)+dx/apu(i,j)*(pc(i,j-1,1)-pc(i,j,1));
end
end
for i=2:ny
for j=2:nx-1
v(i,j,1)=v(i,j,1)+dy/apv(i,j)*(pc(i-1,j,1)-pc(i,j,1));
end
end
err=norm(mp(:,:))
if err<10^-10
break
end

end