Convergence problem with pressure correction equation in SIMPLE
 

Hi all, I'm having some problems with my 2D, incompressible solver using the SIMPLE algorithm (I'm following Patankar and Versteeg&Malalasekera).
The pressure correction equation never converges!
If I plot the error of the pressure correction equation after every iteration, calculated as the difference between the old and the new p' solutions, it gets to a point where it does not decrease any more ¿?.

I'm using succesive over-relaxation, so the computed p' is:



where the suscripts E,W,N and S stand for east, west, north and south sides, b is the mass imbalance term and .

Does anyone have any hint as to what I am doing wrong? your help will be highly appreciated.

Thanks a lot in advance.
Cheers,
Michujo.



P.S: I include here below the piece of code that is not working (I'm using MATLAB).








uprime=zeros(Nx+1,Ny);
vprime=zeros(Nx,Ny+1);
pprime=zeros(Nx,Ny);
pprime_old=zeros(Nx,Ny);
iter=0;itermax=200;

%Calculate pressure correction pprime

b=zeros(Nx,Ny);
a=zeros(Nx,Ny);
ax=zeros(Nx,Ny);
ay=zeros(Nx,Ny);
du=zeros(Nx,Ny);
dv=zeros(Nx,Ny);

for I=1:Nx
for J=1:Ny
b(I,J)=rho*deltay(J)*ustar(I,J)-rho*deltay(J)*ustar(I+1,J)+rho*deltax(I)*vstar(I,J )-rho*deltax(I)*vstar(I,J+1); %Mass imbalance term
du(I,J)=deltay(J)*omegauv/ap_u(I,J);
dv(I,J)=deltax(I)*omegauv/ap_v(I,J);
ax(I,J)=rho*du(I,J)*deltay(J);
ay(I,J)=rho*dv(I,J)*deltax(I);
end
end

for I=1:Nx
for J=1:Ny
if I>=2 && I<=(Nx-1) && J>=2 && J<=(Ny-1) %Interior points
a(I,J)=ax(I+1,J)+ax(I,J)+ay(I,J+1)+ay(I,J);
elseif I==1 && J>1 && J<Ny %West face
a(I,J)=ax(I+1,J)+ay(I,J+1)+ay(I,J);
elseif I==Nx && J>1 && J<Ny %East face
a(I,J)=ax(I,J)+ay(I,J+1)+ay(I,J);
elseif I>1 && I<Nx && J==1 %South face
a(I,J)=ax(I+1,J)+ax(I,J)+ay(I,J+1);
elseif I>1 && I<Nx && J==Ny %North face
a(I,J)=ax(I+1,J)+ax(I,J)+ay(I,J);
elseif I==1 && J==1 %Southwest corner
a(I,J)=ax(I+1,J)+ay(I,J+1);
elseif I==Nx && J==1 %Southeast corner
a(I,J)=ax(I,J)+ay(I,J+1);
elseif I==1 && J==Ny %Northwest corner
a(I,J)=ax(I+1,J)+ay(I,J);
elseif I==Nx && J==Ny %Northeast corner
a(I,J)=ax(I,J)+ay(I,J);
end %End cell type
end
end

while error>tol %&& iter<itermax
for I=1:Nx
for J=1:Ny
if I>=2 && I<=(Nx-1) && J>=2 && J<=(Ny-1) %Interior points
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I+1,J)*pprime(I+1,J )+ax(I,J)*pprime(I-1,J)+ay(I,J+1)*pprime(I,J+1)+ay(I,J)*pprime(I,J-1)+b(I,J))/a(I,J);
elseif I==1 && J>1 && J<Ny %West face
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I+1,J)*pprime(I+1,J )+ay(I,J+1)*pprime(I,J+1)+ay(I,J)*pprime(I,J-1)+b(I,J))/a(I,J);
elseif I==Nx && J>1 && J<Ny %East face
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I,J)*pprime(I-1,J)+ay(I,J+1)*pprime(I,J+1)+ay(I,J)*pprime(I,J-1)+b(I,J))/a(I,J);
elseif I>1 && I<Nx && J==1 %South face
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I+1,J)*pprime(I+1,J )+ax(I,J)*pprime(I-1,J)+ay(I,J+1)*pprime(I,J+1)+b(I,J))/a(I,J);
elseif I>1 && I<Nx && J==Ny %North face
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I+1,J)*pprime(I+1,J )+ax(I,J)*pprime(I-1,J)+ay(I,J)*pprime(I,J-1)+b(I,J))/a(I,J);
elseif I==1 && J==1 %Southwest corner
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I+1,J)*pprime(I+1,J )+ay(I,J+1)*pprime(I,J+1)+b(I,J))/a(I,J);
elseif I==Nx && J==1 %Southeast corner
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I,J)*pprime(I-1,J)+ay(I,J+1)*pprime(I,J+1)+b(I,J))/a(I,J);
elseif I==1 && J==Ny %Northwest corner
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I+1,J)*pprime(I+1,J )+ay(I,J)*pprime(I,J-1)+b(I,J))/a(I,J);
elseif I==Nx && J==Ny %Northeast corner
pprime(I,J)=(1-omegap)*pprime(I,J)+omegap*(ax(I,J)*pprime(I-1,J)+ay(I,J)*pprime(I,J-1)+b(I,J))/a(I,J);
end %End cell type if condition

end
end %End cells loop

error=max(max(abs(pprime-pprime_old)));
pprime_old=pprime;
iter=iter+1
end %End iteration loop