CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Convergence problem with pressure correction equation in SIMPLE

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   October 13, 2012, 08:18
Default Convergence problem with pressure correction equation in SIMPLE
  #1
Senior Member
 
Join Date: Dec 2011
Location: Madrid, Spain
Posts: 133
Rep Power: 5
michujo is on a distinguished road
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:

p'=(1-\omega)\cdot p'+\omega \cdot(a_E\cdot p'_E-a_W\cdot p'_W+a_N\cdot p'_N-a_S\cdot p'_S+b)/a_p

where the suscripts E,W,N and S stand for east, west, north and south sides, b is the mass imbalance term and a_P=a_E+a_W+a_S+a_N.

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
michujo is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Having troubles in understanding the correction term in pressure equation kaifu OpenFOAM 1 October 8, 2012 05:19
Pressure correction BC using SIMPLE André Gaathaug Main CFD Forum 2 August 3, 2007 01:38
pressure correction in SIMPLE manoj Main CFD Forum 7 July 19, 2007 08:21
Preconditioning for Pressure Correction Equation cfd101 Main CFD Forum 1 February 23, 2006 14:34
Pressure correction question - Please help! Barry Main CFD Forum 1 February 5, 2003 00:27


All times are GMT -4. The time now is 19:02.