
[Sponsors] 
August 20, 2011, 22:26 
TIPS for those working on SIMPLE

#1 
Member
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 8 
Hey guys!
I am a Master student and i'm programming SIMPLE method on myself.(duct flow) I have been confused for a very very long time,and now I eventually get familiar with SIMPLE and I'd like to share some of experience to all of you. 1.Use Poiseuille flow to check your momentum equation.Apply the exact solution of velocity and pressure drop as boundary and initial condition.If your momentum equation is correct, it should never change the velocity field.Especially pay attention to the wall boundary condition. 2.Check mass balance over all flow field.This is VERY important and have NOT been stressed in most books,like Patankar's and HK.Versteeg's text book. For example,u assumed a pressure field to be even as an initial condition.Then u solve the momentum equation,the result is that velocity field will definitely NOT fulfill the overall mass balance,i.e.,the mass flux on outlet will not equal flux you set on the inlet. Terrible thing is,If u don't correct the velocity on outlet to make overall mass balance,the pressure correction equation will BLOW UP!You can never get a pressure correction solution, and the pressure correction will grow up and never stop. By the way,you may wonder, that why all the sample codes are written for a cavity flow, one important reason I guess, cavity flow is INTERNAL flow and u never worry about overall mass balance,so there is no necessary to correct. However, if u modify them to a duct flow you will definitely FAIL,because of the blow up of pressure correction equation. Hope this can help someone,good luck! 

August 21, 2011, 05:10 

#2  
Senior Member
Vieri Abolaffio
Join Date: Jul 2010
Location: Always on the move.
Posts: 308
Rep Power: 10 
Quote:


August 21, 2011, 21:41 

#3 
Senior Member
Cean
Join Date: Feb 2010
Posts: 128
Rep Power: 9 
I think that is just the difference between duct flow & cavity flow.
for duct flow, you need to add outlet type of BC which check the global mass balance. 

August 23, 2011, 16:41 

#4 
New Member
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 8 
hey houkensjtu
I am coding the driven cavity using SIMPLE in C++ using the staggered grid and ghost cell approach, however my pressure keeps growing and the residuals start diverging after the 15th iteration or so. Any ideas why this might be happening. Thanks 

August 23, 2011, 21:56 

#5  
Member
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 8 
Quote:
Then I still suggest u to check the overall mass balance,i.e.,(u(i,j+1)u(i,j))+(v(i+1,j+1)v(i,j+1)),is the sum of this term being zero?If so,u should get a solution for pressure correction,and CHECK after u apply the correction to velocity,the velocity field should be divergence free. Also,u need to pay attention to underrelaxation factor,usually,factor for pressure should be less than 0.2,and 0.6~0.8 for momentum equation.This also might be a reason cause diverging. 

August 23, 2011, 22:10 

#6 
New Member
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 8 
Thanks houkensjtu
I had taken care of the under relaxation. However the sum of (u(i,j+1)u(i,j))+(v(i+1,j+1)v(i,j+1)) for each cell is not zero. I will try fixing it. Thanks again. 

August 23, 2011, 22:42 

#7 
Member
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 8 
Be careful, do not misunderstand my point.The "sum" doesn't mean this term should be zero for each cell, but the overall summation of all cell should be zero.


August 23, 2011, 22:46 

#8 
New Member
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 8 
Yeah I realized that later when you mentioned overall mass balance. Unfortunately the overall mass balance for my code is not resulting in zero.


August 23, 2011, 23:35 

#9 
New Member
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 8 
The other problem that I am facing is when i update either u or v then I get the u or v to look like it should. How ever when I update both u and v then i get junk values.


August 24, 2011, 09:30 
sad ...

#10 
Member

I guess you come from shanghai jiao tong university ....China
well, as for your tips ... i would like to share some words with you .... the all stuff you mentioned is really too simple and basic ... the essence is how to treat open boundary conditions , that is all ... if you want to talk about something useful, please read carefully Tao 's book named 'recent development of numerical heat transfer' (2001, pp234  248) or you can also read J.H. Fergizer, M. Peric, Computational Methods for Fluid Dynamics, Springer, 2002. pp. and its code _ the implemented details for Backwardfacing step flow ; The reason why people like liddriven cavity flow can be summarized as, This problem is a nice one for testing for several reasons. First, as mentioned above, there is a great deal of literature to compare with. Second, the (laminar) solution is steady. Third, the boundary conditions are simple and compatible with most numerical approaches. as for open boundary condition and its treatment for mass conservation is another topic ... best wishes ... and it is so nice to meet people from china who would like to write some codes by himself, if possible, we can share some more in CFD, xikuanghit@gmail.com 

August 26, 2011, 14:38 

#11 
New Member
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 8 
hey houkensjtu,
I have kind of narrowed down the cause for divergence. When a moving lid is horizontal then the divergence occurs in the xvelocity but no divergence in the yvelocity. The viceversa occurs when the moving lid is vertical. Any ideas why this is happening. Thanks 

July 7, 2012, 12:53 

#12 
Member
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 7 
Hi everyone!
I know that I am being an absolute pest on this forum, but desperate times call for desperate measures! I am about 80% complete on my implementation of the SIMPLE algorithm in MATLAB for steady flow between two flat plates, and similar to this thread, I am having issues with the pressure correction equation... also, I am struggling to understand how to implement residuals and what is involved... as with most CFD textbooks, everything always seems a little vague, well, to me anyway! I'd greatly appreciate any assistance, the purpose of my project is for educatinal reasons. the code follows versteeg, and is easily readable to both beginners to cfd AND to Matlab, and I fully intend on making it a "go to code" for anyone with no experience, as I (and everyone else on this thread) have experienced the blood sweat and tears behinf a cfd code!! (mine continues! ) best regards, Michael 0565113@gmail.com 

July 7, 2012, 14:16 

#13  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 
Quote:
Read peric's book. There is nothing vague in it. There are all the information in it to write a perfectly working code and all these information are provided systematically and in a language that is easy to understand. As far as writing codes is concerned, here is small but the most important aspect that i have learned in last 8 years. If you do everything correctly, it will all work correctly. Doing everything correctly is hardest part , get this and you will see that things do work for you too like they work for others. 

July 7, 2012, 22:09 

#14  
Member
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 8 
Quote:
After that, my first question is, what kind of outlet are you applying? NOTING in Versteeg text book, all SIMPLE algorithm examples are done with a pressure outlet, NOT velocity outlet. Velocity outlet need extra modification to insure pressure correction equation being solvable. Good luck! 

July 12, 2012, 05:49 

#15 
Member
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 7 
Hello everyone, thanks for the advice! I went through it again and found a mistake, but things are still not right.. I am trying to solve for steady Poiseulle flow on a backward staggered grid using FVM and SIMPLE, these are the initial boundary conditions that I set...
%% Wall 1,The inlet u(2,=uin; % U values are stored at i=2 v(1,=vin; % V values are stored at I=1 p(1:2,=pstaticin; % Pressure is stored at I=1, and using dp/dx=0 enters the domain at I=2 %% Wall 2, The top plate %For Poiseulle flow; No slip boundary condition: u(:,nYp)=0; v(:,nY)=0; p(:,nYp)=0; %% Wall 3, Outlet u(nX,2:nYp1)=u(nX1,2:nYp1); % First extrapolate values from the domain u(nX,2:nYp1)=u(nX1,2:nYp1)*sum(u(2,2:nYp1))/sum(u(nX,2:nYp1)); % This multiplier is the mass in divided by the mass out, which ensures continuity. v(nXp,=v(nXp1,; % See Versteeg page 273 p(nXp,2:nYp1)=0; % Constant static pressure value of '0' set at the outlet (I = NI) %% Wall 4, Plate %For Poiseulle flow; No slip boundary condition: u(:,1)=0; v(:,2)=0; p(:,1)=0; From here, I solve for the momentum equations, noting that the source term is the pressure gradient, as well as the source terms in appendix C of Versteeg. The wall shear stress also enters the umomentum equations in the central coefficient, Versteeg: Sp = (mu/delta yp)*Acell ... eqn 9.11 pg 275. Since I am using a backward staggered grid, there is no need to incorporate the inlet and outlet boundary conditions as source terms (right hand side of equation), since the hybrid differencing scheme ustilizes the already stored inlet value of uin at i=2, and the umomentum equations are calculated from i=3. Once u* and v* are calculated, solve for the pressure correction equaiton, being sure to cut the coefficients to the inlet (aW), outlet (aE), and the top and bottom plates (aN and aS), noting also that the coefficients are simply the central coeficients of the u & v momentum equations, multiplied by the density and the cellface area (rho*d*a)i,J etc... and that the source term is the continuity imbalance accross a pressure volume, i.e West face  East face + south face north face. I am using the MATLAB function BICGSTAB to solve the sets of equations... can anyone point out something wrong? REgards, 0565113@gmail.com 

July 14, 2012, 02:27 

#16  
Member
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 8 
Quote:
Also, I think BICGSTAB is ok but not necessary, because your grid number is certainly not so huge(maybe several hundreds or thousands). Usually backslash is far more quick for small problem. Matlab's backslash is very capable. good luck! 

July 14, 2012, 03:01 

#17  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 
Quote:
this is source of error for you, you do not apply fixed pressure boundary at the walls. Walls are always neumann conditions. 

July 14, 2012, 03:13 

#18  
Member
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 8 
Quote:
For pressure correction equation, as you said it is neumann boundary condition, so pressure correction on wall side is cut. For the source term because velocity on wall is 0 so it also doesn't matter. 

July 14, 2012, 03:47 

#19 
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 723
Rep Power: 19 
Thats what you believe.
in reality it matters a lot. But in case someone is fond of debugging programs, it is good thing to believe in. :D 

July 14, 2012, 11:48 

#20 
Member
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 7 
Thankyou so much for the reply! The pressure boundary condition here was indeed a mistake, however since my domain is I=2 to I=NI1, and J=2 to NJ1, those cells never actually came into the solution.


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
LidDriven cavity flow with SIMPLE method  luckyxu  Main CFD Forum  6  November 9, 2011 08:17 
Some tips for Snappyhexmesh required  basilwatson  OpenFOAM  3  December 9, 2010 03:53 
The correction on pressure equation of SIMPLE algorithm in MRFSimpleFOAM solver  renyun0511  OpenFOAM Running, Solving & CFD  0  November 10, 2010 02:47 
Help me on SIMPLE  L. Aouanouk  Main CFD Forum  6  April 17, 2003 05:08 
flow over a 2D cyl using SIMPLE  T  Main CFD Forum  1  January 27, 2001 08:32 