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

TIPS for those working on SIMPLE

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

Like Tree1Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 20, 2011, 22:26
Default TIPS for those working on SIMPLE
  #1
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
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!
houkensjtu is offline   Reply With Quote

Old   August 21, 2011, 05:10
Default
  #2
Senior Member
 
sail's Avatar
 
Vieri Abolaffio
Join Date: Jul 2010
Location: Always on the move.
Posts: 308
Rep Power: 8
sail is on a distinguished road
Quote:
Originally Posted by houkensjtu View Post
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... but kudos for sharing.
sail is offline   Reply With Quote

Old   August 21, 2011, 21:41
Default
  #3
Senior Member
 
Cean
Join Date: Feb 2010
Posts: 128
Rep Power: 7
shirazbj is on a distinguished road
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.
shirazbj is offline   Reply With Quote

Old   August 23, 2011, 16:41
Default
  #4
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 6
illuminati5288 is on a distinguished road
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
illuminati5288 is offline   Reply With Quote

Old   August 23, 2011, 21:56
Default
  #5
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
Quote:
Originally Posted by illuminati5288 View Post
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.
houkensjtu is offline   Reply With Quote

Old   August 23, 2011, 22:10
Default
  #6
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 6
illuminati5288 is on a distinguished road
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.
illuminati5288 is offline   Reply With Quote

Old   August 23, 2011, 22:42
Default
  #7
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
Quote:
Originally Posted by illuminati5288 View Post
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.
houkensjtu is offline   Reply With Quote

Old   August 23, 2011, 22:46
Default
  #8
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 6
illuminati5288 is on a distinguished road
Yeah I realized that later when you mentioned overall mass balance. Unfortunately the overall mass balance for my code is not resulting in zero.
illuminati5288 is offline   Reply With Quote

Old   August 23, 2011, 23:35
Default
  #9
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 6
illuminati5288 is on a distinguished road
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.
illuminati5288 is offline   Reply With Quote

Old   August 24, 2011, 09:30
Default sad ...
  #10
Member
 
Wu Jian
Join Date: Jun 2009
Location: Poitiers
Posts: 33
Rep Power: 8
harbinyg is on a distinguished road
Send a message via MSN to harbinyg
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
harbinyg is offline   Reply With Quote

Old   August 26, 2011, 14:38
Default
  #11
New Member
 
Vignesh Suresh
Join Date: Jul 2011
Posts: 10
Rep Power: 6
illuminati5288 is on a distinguished road
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
illuminati5288 is offline   Reply With Quote

Old   July 7, 2012, 12:53
Default
  #12
Member
 
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 5
michaelmoor.aero is on a distinguished road
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
michaelmoor.aero is offline   Reply With Quote

Old   July 7, 2012, 14:16
Default
  #13
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 374
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by michaelmoor.aero View Post
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.
arjun is offline   Reply With Quote

Old   July 7, 2012, 22:09
Default
  #14
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
Quote:
Originally Posted by michaelmoor.aero View Post
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
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!
houkensjtu is offline   Reply With Quote

Old   July 12, 2012, 05:49
Default
  #15
Member
 
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 5
michaelmoor.aero is on a distinguished road
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
michaelmoor.aero is offline   Reply With Quote

Old   July 14, 2012, 02:27
Default
  #16
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
Quote:
Originally Posted by michaelmoor.aero View Post
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!
houkensjtu is offline   Reply With Quote

Old   July 14, 2012, 03:01
Default
  #17
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 374
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by michaelmoor.aero View Post
%% 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;

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

Old   July 14, 2012, 03:13
Default
  #18
Member
 
HouKen
Join Date: Jul 2011
Posts: 66
Rep Power: 6
houkensjtu is on a distinguished road
Quote:
Originally Posted by arjun View Post
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.
houkensjtu is offline   Reply With Quote

Old   July 14, 2012, 03:47
Default
  #19
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 374
Rep Power: 10
arjun is on a distinguished road
Quote:
Originally Posted by houkensjtu View Post
emmm... I think it's not the problem.
Thats what you believe.

Quote:
Originally Posted by houkensjtu View Post
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.
in reality it matters a lot. But in case someone is fond of debugging programs, it is good thing to believe in. :-D
arjun is offline   Reply With Quote

Old   July 14, 2012, 11:48
Default
  #20
Member
 
Michael Moor
Join Date: May 2012
Location: Ireland
Posts: 30
Rep Power: 5
michaelmoor.aero is on a distinguished road
Quote:
Originally Posted by arjun View Post
this is source of error for you, you do not apply fixed pressure boundary at the walls. Walls are always neumann conditions.
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.
michaelmoor.aero 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
Lid-Driven 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


All times are GMT -4. The time now is 15:04.