[rks171]

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(i-1) + 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 1-D 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'(i-1) - 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'(i-1) + 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(i-1), 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_new-rho_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 non-physical, but my methodology seems to follow perfectly with the description of the SIMPLE algorithm given by Patankar.

[arjun]

sorry I did not read your full post, its too large for me.

Originally Posted by rks171

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 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.

[rks171]

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, 1E-5 seconds to 1E-8 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.

[umassche]

...........

[arjun]

Originally Posted by rks171

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, 1E-5 seconds to 1E-8 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.

It is only in the beginning of the calculation. As you approach convergence things settle down.

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.

[madhukar_m_rao]

It is also a good idea to have mass conserving initial and boundary conditions for unsteady simulations of low Mach number or incompressible flows.

[rks171]

Upon convergence your mass imbalance will be approaching 0 and hence no reasons for any pressure correction.

Yes, I understand this. However, my problem isn't beginning of simulation vs. end of simulation. It's the fact that if I run the simulation with a timestep size of 1E-5 seconds, then my first pressure correction is 3000 Pa (for example). If I run the simulation again, but with a timestep size of 1E-8 seconds, which is three orders of magnitude smaller, my pressure correction is now 3,000,000 Pa. I'm sure that, in either case, the pressure will get smaller on the next step and subsequent steps (if the extreme pressure correction doesn't move outside of my steam tables). It's just that it doesn't seem that running the simulation with a smaller timestep should lead to a bigger pressure correction. If anything, I would think that dropping the timestep size to 1E-8 seconds would cause the pressure correction to go from 3000 Pa to 3 Pa, or somewhere in that order.

Thank you for your assistance.

[arjun]

Originally Posted by rks171

Yes, I understand this. However, my problem isn't beginning of simulation vs. end of simulation. It's the fact that if I run the simulation with a timestep size of 1E-5 seconds, then my first pressure correction is 3000 Pa (for example). If I run the simulation again, but with a timestep size of 1E-8 seconds, which is three orders of magnitude smaller, my pressure correction is now 3,000,000 Pa. I'm sure that, in either case, the pressure will get smaller on the next step and subsequent steps (if the extreme pressure correction doesn't move outside of my steam tables). It's just that it doesn't seem that running the simulation with a smaller timestep should lead to a bigger pressure correction. If anything, I would think that dropping the timestep size to 1E-8 seconds would cause the pressure correction to go from 3000 Pa to 3 Pa, or somewhere in that order.

Thank you for your assistance.

I thought we already settled up that part that decreasing the time step by 3 orders would increase the pressure correction by 3 orders. Which is what you observe. As I explained this is related to Ap.

Originally Posted by rks171

If anything, I would think that dropping the timestep size to 1E-8 seconds would cause the pressure correction to go from 3000 Pa to 3 Pa, or somewhere in that order.

Sorry but I do not agree that it should happen.
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.

[rks171]

I thought we already settled up that part that decreasing the time step by 3 orders would increase the pressure correction by 3 orders. Which is what you observe. As I explained this is related to Ap.

Yes, I'm sorry, we have. It really perplexed me, though, because so often in my work, if the solution was changing too much in a timestep, the remedy was always to reduce the timestep size in order to reduce the changes and prevent the solution from going unstable and diverging. In this case, reducing the timestep made the problems bigger and I didn't expect that at all. It thought it had to be a mistake in my implementation of the algorithm. But now that I can see that it is actually the behavior of the algorithm, how do I correct this problem? Have my code reduce the under-relaxation factor when the pressure correction exceeds a certain threshold?

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.

[arjun]

Originally Posted by rks171

Yes, I'm sorry, we have. It really perplexed me, though, because so often in my work, if the solution was changing too much in a timestep, the remedy was always to reduce the timestep size in order to reduce the changes and prevent the solution from going unstable and diverging. In this case, reducing the timestep made the problems bigger and I didn't expect that at all. It thought it had to be a mistake in my implementation of the algorithm. But now that I can see that it is actually the behavior of the algorithm, how do I correct this problem? Have my code reduce the under-relaxation factor when the pressure correction exceeds a certain threshold?

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.

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.

[rks171]

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.

[arjun]

Originally Posted by rks171

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.

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.

[rks171]

If it is your own code then how about remove the time contribution from Ap that is needed for pressure correction.

That's an interesting idea... does the time contribution need to be in the pressure correction equation? I wonder if it will affect the converged pressure field solution. I'll test that out.

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 under-relax the solution.

Since the pressure correction equation is source-dominated, 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 pressure-velocity 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 inter-phase 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.

[arjun]

Originally Posted by rks171

That's an interesting idea... does the time contribution need to be in the pressure correction equation?

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.

[madhukar_m_rao]

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.

[rks171]

It is a two-phase 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 re-calculation 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 re-calculation and get the simulation to converge with those models on.

[madhukar_m_rao]

You may want to try implementing the compressibility modifications to the pressure-correction as described in chapter 10 of the book by Ferziger and Peric.

[rks171]

You may want to try implementing the compressibility modifications to the pressure-correction as described in chapter 10 of the book by Ferziger and Peric.

I don't have that book right now, actually, but I do have a question about this temporal term. When we discretize the temporal term in, say, the continuity equation, we get:

(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?

[madhukar_m_rao]

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

[rks171]

I expect that you have the corresponding terms for the
momentum and energy equations also.

Yes, I have the terms setup in the same way for the momentum equation also... I didn't put my energy equation in yet because I'm still having trouble getting the "momentum -> pressure correction -> velocity/density correction -> continuity" iteration loop to work.

In addition, your mass fluxes will depend on the density correction as well as the
pressure correction.

My code has steam tables. After the pressure correction gives the new absolute pressures, I use that as the new saturation pressure and use the saturated liquid enthalpy to get liquid density via correlation and the saturated vapor enthalpy to get the vapor density via correlation.

If you cannot get access to Ferziger
and Peric to work out the compressible pressure correction equation (with associated
boundary conditions) let me know.

I checked and I will have access to this book. I will take a look at this book tomorrow. Right now, I'm stuck with getting pressure to converge on a value. If I force phase void fraction consistency (i.e. sum of all voids = 1) at the end of my continuity equation step, the pressure and density enter into a vicious cycle where a small change in pressure makes density change in the opposite direction and then that density change causes an even larger pressure change which causes an even larger density change, etc., until the pressure diverges to unrealistic values. If I let the void fractions go where they please (and the sum of voids isn't forced to equal 1 each iteration), the pressure slowly goes up, then down and drops really fast as soon as the relative pressure becomes negative in my system. It seems to me that the void fraction change damps out the the pressure/density relationship, but my code still isn't converging even when I don't manipulate the voids.

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.