
[Sponsors] 
April 24, 2012, 19:44 
Timestep and Pressure Correction Relationship in SIMPLE

#1 
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
I'm noticing some behavior that seems strange to me in a program I wrote that utilizes the SIMPLE algorithm. During the first iteration of the first timestep in my simulation, and due to a jump in void fraction between two mesh cells, there was a rather large net outflow of mass from a cell. As one would expect, this led to a predicted drop in pressure. But the pressure drop that was predicted was much too large, so I figured that, because there was a lot happening right in the beginning, I should decrease my timestep size. I did this, and the predicted pressure drop in the cell INCREASED. This seems completely illogical to me because less time in one step means less time for mass to go in and out of the cell and, thus, should mean less change in pressure.
But then I was looking at the equations for SIMPLE and I can't tell if this behavior is actually a feature of the SIMPLE algorithm. Our form of the momentum equations is: a*u(i) = b*u(i1) + c*u(i+1) + s a, b, and c are coefficients and s is the source due to drag and pressure drop and so on. The u terms are the velocities at three adjacent grid points (this is a 1D problem). The a and s terms include the time step, dt. In a, we have rho*dx/dt. This multiplies by u(i), which is the new time velocity. s contains the old time density and velocity and also has a dt in it. So we see that by decreasing dt, the value of a will increase with respect to the b and c coefficients. This will result in u changing by smaller amounts over each timestep. This is logical. But now we form the pressure correction equation. A velocity correction is obtained from the momentum equation: a*u'(i) = p'(i1)  p'(i+1) In other words, the change in the pressure around the node will change the velocity. Now our velocity is defined u = u* + u', where u* was our guessed velocity and u' is the correction to make it the true velocity, u. We can substitute this definition of the velocity into the continuity equation and now we have an equation in terms of the pressure corrections throughout the mesh. x*p'(i) = y*p'(i1) + z*p'(i+1) + S The x, y, and z coefficients contain the previous a coefficients from the momentum equation due to our definition of the velocity correction term in terms of the pressure corrrections. So x is rho/a(i), y is rho/a(i1), and z is rho/a(i+1). All these a coefficients have dt in the denominator. So what we end up having is that the behavior of the x, y, z coefficients is exactly the opposite of the behavior of the a, b, and c coefficients from before. a, b, and c got big when dt went small. x, y, and z will go small when dt goes small. Now, there is also a dt showing up in the S term... we have the change in density over the timestep i.e. dx*(rho_new  rho_old)/dt. So the S term should get smaller too, but this is my first iteration and first timestep, so my (rho_newrho_old) is zero and the dt in the source term has no effect at all. S also contains the convective terms, i.e. rho*u_in  rho*u_out. And because my void fraction increases a lot in one of the cells next door, the u_out is multiplying by a much larger density than the u_in. So what I see is that decreasing dt is decreasing my x, y, z coefficients, but it has no effect on my source term, S. Thus, my pressure correction is increasing by 2 extra orders of magnitude when I decrease my dt by two orders of magnitude. This just seems nonphysical, but my methodology seems to follow perfectly with the description of the SIMPLE algorithm given by Patankar. 

April 24, 2012, 20:03 

#2  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 390
Rep Power: 10 
sorry I did not read your full post, its too large for me.
Quote:
This is actually not illogical. It sounds illogical but based on my experience with coding SIMPLE it is not. This is assuming that you are talking about beginning of the calculation. Upon convergence, it will settle to correct values. Why this happens. As you decrease your time step. The momentum coefficient increases. Ap = Ap_time + Ap_due_to_convection_diffusion As deltA decrease Ap_time increases and thus Ap increase. This leads to higher pressure corrections in the start , when solution not converged. So if you check after starting few iterations pressure drop shall be higher than you expected. If you are talking about fully converged field then I agree , it is strange to have pressure drop dependent on time step. 

April 24, 2012, 20:35 

#3 
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Thank you for replying. I'm sorry for the long last post, but I was giving my reasoning for thinking this behavior was correct. But the short question would be, in the SIMPLE algorithm, if you drop the timestep size, will you see a corresponding increase in the pressure corrections?
So if I read your response right, you're saying that dropping the time step from, say, 1E5 seconds to 1E8 seconds will cause the temporal portion of the momentum coefficient to increase. Then this, in turn, causes the pressure correction to increase, too. I actually did this, and, decreasing the timestep by 3 orders of magnitude caused an increase in my pressure corrections by about 3 orders of magnitude. If this is how the algorithm works, though, then whether it is beginning of the timestep or fully converged, the pressure correction will always react in exactly the opposite direction as the change in timestep size, right? This seems like it would cause instability to me. 

April 24, 2012, 20:47 
test

#4 
New Member
John Damnson
Join Date: Feb 2011
Posts: 15
Rep Power: 7 
...........


April 24, 2012, 21:23 

#5  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 390
Rep Power: 10 
Quote:
You should remember that now since time step is decreased velocity changes slowly even if pressure is changing very much. Upon convergence your mass imbalance will be approaching 0 and hence no reasons for any pressure correction. Ideal condition mass imbalance = 0 , would produce 0 pressure corrections. Last edited by arjun; April 24, 2012 at 22:29. 

April 25, 2012, 05:06 

#6 
New Member
Dr. Madhukar M Rao
Join Date: Nov 2010
Location: Bangalore, India
Posts: 20
Rep Power: 7 
It is also a good idea to have mass conserving initial and boundary conditions for unsteady simulations of low Mach number or incompressible flows.


April 25, 2012, 08:35 

#7  
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Quote:
Thank you for your assistance. 

April 25, 2012, 09:32 

#8  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 390
Rep Power: 10 
Quote:
Quote:
In pressure correction Ap appears as 1/Ap So say you had A * p = src p_old = A^1 . src Now if you increase Ap by say factor of 1000 you will decrease A by factor of 1000 , so (A/1000) * p = src p = (A^1 . src) * 1000 which is what you observe. PS: (this is to give rough idea). Good luck. 

April 25, 2012, 09:45 

#9  
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Quote:
Also, another reason I got hung up on this behavior is that I can't see a physical reasoning why making the window of time that you're doing the analysis for would make the pressure change in the cell larger. If you are looking at a smaller window of time, then less mass can cross over into your cell. If less mass is coming in or going out, then the amount by which the pressure changes due to that mass transfer should also be less. Numerically, by looking at the SIMPLE equations, yes, I can see why the pressure change gets larger with a smaller timestep. But it just doesn't make sense to me physically, and it seems to me that, above all else, our solution should produce physically realistic results. 

April 25, 2012, 20:01 

#10  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 390
Rep Power: 10 
Quote:
I can't be much help but I can tell you this: by chosing smaller time step even if pressure corrections are larger the velocity field changes slowly. So it is more stable in transient case. If you observe opposite of this then there are only two possiblity 1. You code has some bug. 2. Your mesh is too bad at some place. Velocity field resposes to pressure correction like this u' = (volume/Ap) dp'dx So you see p' is now 1000 times but so does Ap that is u' = (volume/Ap) dp'dx = u' = (volume/ (Ap*1000) ) d(p'*1000)dx So your velocity correction should be of similar order and not the larger than first case. 

April 25, 2012, 20:33 

#11 
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Indeed, the velocity corrections are small, but the pressure corrections change the absolute pressure in my cells. The absolute pressure is used to calculate the state properties, like density, specific heat, surface tension, etc. It only takes a handful of large pressure corrections during successive iterations to put the absolute pressure in the cell right outside of a realistic value. Because the pressure change has only a small impact on the velocity field (due to the large value of Ap), the code keeps applying a large pressure change each iteration to try and correct the mass imbalance caused by the incorrect velocity field.


April 25, 2012, 20:45 

#12  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 390
Rep Power: 10 
Quote:
If it is your own code then how about remove the time contribution from Ap that is needed for pressure correction. I am not sure it would work or not but my feeling is that it might just work. I could try but i am in transition from version 2 to version 3 of my code and lots of things are mixed up. Edited to add: Another option is that first run the simulation with really small time steps and switch off all the other models. Once you get proper pressure field switch the models on . Is it possible. 

April 25, 2012, 21:15 

#13  
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Quote:
I did try out something different with the pressure correction equation today. Have you read Patankar's "Numerical Heat Transfer"? Well, Section 7.2 talks about source term linearization, so that it's form is S = Sc + Sp * phi, where phi is the dependent variable. In the case of the pressure correction equation, I take this to be P'. Of course, the source term has no P' in it, so I originally made Sp zero, but by choosing how you define Sc and Sp, you can underrelax the solution. Since the pressure correction equation is sourcedominated, it can be approximated by Sc + Sp * phi = 0. We can put a limit on what our new phi value is by setting the source coefficients: Sc = S * phi~ / (phi~  phi*) Sp =  S * phi~ / (phi~  phi*) phi~ is the limit that phi can change to and phi* is the current value of phi. If S is positive, it will increase phi so we use a max limit for phi~. If S is negative, it will decrease phi so we use a min limit for phi~. So I put my max and min allowable pressures in for phi~. I figured phi* was zero, since there is no "current pressure correction". Defining my source term in this way prevented the "steam table out of bounds" error I got before. That, and I used a larger timestep. This allowed the solution to converge to a pressurevelocity field. However, it only works when I shut off my continuity equations for the phases. So I'm only running momentum and pressure correction, but I'm still calling my steam tables each iteration in order to calculate the interphase terms. I think this is akin to your suggestion to run with other models off in order to get a correct pressure field. I think there's a problem with my continuity equation formulation, so I'll investigate that tomorrow. 

April 25, 2012, 23:13 

#14  
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 390
Rep Power: 10 
Quote:
As far as need is concerned, my opinion is that it is not (I could be wrong). but this is not checked. The reason it is not checked is that it might be unstable because Ap_time is very large compared to Ap_others. So it may not work out well. But it is no harm to test it once. 

April 26, 2012, 01:25 

#15 
New Member
Dr. Madhukar M Rao
Join Date: Nov 2010
Location: Bangalore, India
Posts: 20
Rep Power: 7 
Since you seem to have a compressible flow, are you including the terms due to the density correction in the pressure correction equation (See the description of the compressible SIMPLE algorithm in Patankar's book or in the book by Ferziger and Peric.


April 26, 2012, 08:41 

#16 
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
It is a twophase flow that I'm modeling, so there is vapor in there. The vapor speed doesn't get exceedingly high, though, so I can't say it's compressible. At any rate, I do have the dx/dt * ( (alpha*rho)_new  (alpha*rho)_old ) term present in my pressure correction equation, but for now, I have the density recalculation at the end of the iteration turned off. alpha is the void fraction of the phase. Today, I'm going to try turning on my continuity equation and density recalculation and get the simulation to converge with those models on.


April 27, 2012, 00:57 

#17 
New Member
Dr. Madhukar M Rao
Join Date: Nov 2010
Location: Bangalore, India
Posts: 20
Rep Power: 7 
You may want to try implementing the compressibility modifications to the pressurecorrection as described in chapter 10 of the book by Ferziger and Peric.


April 29, 2012, 12:26 

#18  
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Quote:
(rho_new * alpha_new  rho_old * alpha_old) / dt What I've been doing in my code is to set rho_old and alpha_old to the converged values from the last timestep. Then I enter the new timestep calculations, which will iterate until a level of convergence is reached. Each iteration, I update rho_new and alpha_new to their new iteration values, but I leave rho_old and alpha_old always set to their last timestep value. Is this the right behavior? Last edited by rks171; April 29, 2012 at 17:09. 

April 30, 2012, 02:08 

#19 
New Member
Dr. Madhukar M Rao
Join Date: Nov 2010
Location: Bangalore, India
Posts: 20
Rep Power: 7 
Yes, that is correct. I expect that you have the corresponding terms for the
momentum and energy equations also. In addition, your mass fluxes will depend on the density correction as well as the pressure correction. The density correction can be written in terms of the pressure correction by invoking the equation of state. If you cannot get access to Ferziger and Peric to work out the compressible pressure correction equation (with associated boundary conditions) let me know. Best regards, Madhu 

April 30, 2012, 16:06 

#20  
Member
Join Date: Dec 2011
Location: State College, PA
Posts: 87
Rep Power: 6 
Quote:
Quote:
Quote:
I'm pretty much out of ideas at this point, so I'll have to check out this book and see if there's something wrong with the way I formed my pressure correction equation. 

Tags 
pressure correction, simple algorithm, timestep size 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Simulation of a single bubble with a VOFmethod  Suzzn  CFX  18  October 2, 2009 04:18 
udf to export the pressure on a surface thread every timestep?  zeitistgeld  Fluent UDF and Scheme Programming  0  April 30, 2009 04:10 
pressure can be varied at different timestep ?  SKARP  CFX  5  July 20, 2007 15:39 
Neumann pressure BC and velocity field  Antech  Main CFD Forum  0  April 25, 2006 02:15 
hydrostatic and static pressure  pi  CFX  3  May 26, 2004 10:35 