TIPS for those working on SIMPLE
 

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!

i can see lot of swearing and banging your head against a sharp corner behind this post... :D but kudos for sharing.

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.

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

First of all,u should check if u apply the wall boundary appropriately.Because u are using staggered grid,B.C for u and v should be DIFFERENT.
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 under-relaxation 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.

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.

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.

Yeah I realized that later when you mentioned overall mass balance. Unfortunately the overall mass balance for my code is not resulting in zero.

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.

sad ...
 

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 Backward-facing step flow ;

The reason why people like lid-driven 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

hey houkensjtu,

I have kind of narrowed down the cause for divergence. When a moving lid is horizontal then the divergence occurs in the x-velocity but no divergence in the y-velocity. The vice-versa occurs when the moving lid is vertical. Any ideas why this is happening.

Thanks

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

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.

You should check again every coefficient in your discrete equation, especially for the boundary cell, they are extremely important! If you are following versteeg text book, you should definitely read the SIMPLE algorithm example and calculate every step by yourself on paper.

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!

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:nYp-1)=u(nX-1,2:nYp-1); % First extrapolate values from the domain
u(nX,2:nYp-1)=u(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1)); % This multiplier is the mass in divided by the mass out, which ensures continuity.
v(nXp,:)=v(nXp-1,:); % See Versteeg page 273
p(nXp,2:nYp-1)=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 u-momentum 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 u-momentum 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

I read through your description briefly and I think you should apply the boundary condition of Wall3(Which as you said, ensure the continuity condition) EVERY TIME before you solve the pressure correction equation. Coefficient of outlet velocity is cut as you said but they will enter the source term of pressure correction equation.
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!

this is source of error for you, you do not apply fixed pressure boundary at the walls. Walls are always neumann conditions.

emmm... I think it's not the problem.

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.

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

Thank-you so much for the reply! The pressure boundary condition here was indeed a mistake, however since my domain is I=2 to I=NI-1, and J=2 to NJ-1, those cells never actually came into the solution.

Thanks again for the help! so you are saying that after solving the momentum equations yielding u* and v*, that i should apply the outlet boundary condition to these guessed velocity fields u* and v*? That makes sense, as that velocity cv (at i=NI), is the location of the east face of the pressure control volume at I=NI-1, i.e. the last pressure cv in the domain.

On one quick note about the initial guesses, which initialise the solution... at the moment, I have a guessed pressure field p* which is a linear distribution from west to east, i.e. 10-8-6-4-2-0, and I have also set the intial velocity field u* equal to the inlet vlaue, and v*=0. Anderson suggests putting in a velocity spike in a v control volume to create a 2D flow also.
If i set the u velocity field to zero (i=3 onward to i=NI-1), my outlet b.c returns NaN, and the solution fails due to this piece of code:

u(nX,2:nYp-1)=u(nX-1,2:nYp-1); % First extrapolate values from the domain

u(nX,2:nYp-1)=u(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1)); % This multiplier is the mass in divided by the mass out, which ensures continuity

So I guess the question is, do i assign velocity values that correspond with p*, or is it just sufficient that the u-velocities be non-zero?

This is all much appreciated and I am making more progress than I have in a long time! Best Regards, Michael

Yeah you should ensure continuity on outlet every time before solving pressure correction.
As for NaN problem, I guess, because you applied 0 to all u velocity as initial guess, so when you apply
actually you are dividing by zero,or sth very near zero.
As far as I know...pressure initial guess is not so crucial problem.
good luck

I assume then that this is what you mean:

%% Apply Boundary Conditions
boundary_conditions();
%% Solve for u*
[ustar, diJ]=u_momentum();

%% Check continuity again
% u(nX,2:nYp-1)=ustar(nX-1,2:nYp-1); % First extrapolate values from the domain
% u(nX,2:nYp-1)=ustar(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1));
%% Solve for v*
[vstar, dIj]=v_momentum();

%% Solve for p'
[pdash] = pressure_correction(diJ, dIj);


Where now I have used the values of u* along i=NI-1, whereas in the boundary conditions, I had used the values of the initial guess.

Thank you very much for your notes!
My cavity code wasn't working till I used your suggested relaxation factors! And now it is working like a charm!!