[lorbekl]

Hello all,

Thanks to your help in this thread Advice regarding algorithm for solution of discretised equations in compressible flow I have managed to make a steady state solver for 1D flow based on a compressible SIMPLE algorithm as described by LuckyTran.

Now I am struggling to extend the steady state model to a transient model. Following this advice:

Lucky Tran

In transient solvers: one usually solves the continuity equation for the density at the new time (since it's wrong to assume that density is unchanged). Then one solves the momentum equation for the velocity. Then solve the energy equation. Then do the pressure correction. Finally there is also a density correction at the end based on the final T,P to get the correct density.

I would like your opinion on a couple of things.
- Why do we solve for density twice in one iteration (continuity eq. and eq. of state)? Isn't this a bit unusual?
- if continuity is supposed to be satisfied after the pressure correction equation, isn't it problematic to do a density correction afterwards as this changes the mass fluxes?
- Should the mass fluxes be recalculated after the pressure correction eq. (when cont. is satisfied) or at the end of the iteration?

I would in general like to hear any opinions on this matter. I noticed in the book by Ferziger and Peric that the EOS is used to calculate pressure but I couldn't figure out their algorithm. Any thoughts on that algorithm would be appreciated aswell!

[LuckyTran]

Originally Posted by lorbekl

I noticed in the book by Ferziger and Peric that the EOS is used to calculate pressure but I couldn't figure out their algorithm. Any thoughts on that algorithm would be appreciated aswell!

EOS is used to calculate pressure has the same meaning as when I said calculate density from T & P. It should be apparent that an EOS is used to calculate density at some point (not always strictly to calculate density, because sometimes you use the EOS to calculate pressure from density & temperature).

I would just add that SIMPLE is a predictor-corrector method. Saying that you "solve" anything has a different meaning. What you are doing is calculating a prediction and then applying correctors. It also means you can't fully satisfy all the governing equations in 1 sweep.

Remember you are implicit in time. You need a value for density before you can calculate anything else. Density appears in the continuity equation and momentum equation. You cannot solve the momentum equation for velocity unless you already have an initial guess for density. How you get this initial guess can be done a few ways.

Originally Posted by lorbekl

- if continuity is supposed to be satisfied after the pressure correction equation, isn't it problematic to do a density correction afterwards as this changes the mass fluxes?
- Should the mass fluxes be recalculated after the pressure correction eq. (when cont. is satisfied) or at the end of the iteration?

You satisfy continuity after the pressure correction but you satisfy the wrong continuity equation because you have the incorrect density.

If you don't update the density, then you already know that you do not satisfy continuity. Updating density means you no longer satisfy continuity that is true.

All fluxes need to be recalculated, and that's why you normally do more than 1 iteration of SIMPLE. The inconsistencies you have between density, pressure, & velocity (hopefully) get corrected with many iterations (just like in a steady algorithm). You have a choice. You can correct them immediately or you can take them into another iteration.

There are some exceptions, i.e. when you take really small timesteps so that variables do not change drastically and coupling between equations is limited.

In conclusion I agree it's still debatable whether you pull density out of the continuity equation at the start.

[lorbekl]

Updating density means you no longer satisfy continuity that is true.

But isn't this the condition you have to keep at the end of each iteration in a SIMPLE algorithm?

You cannot solve the momentum equation for velocity unless you already have an initial guess for density. How you get this initial guess can be done a few ways.

Would it be possible to just use the value from the previous time step/iteration and then use only the EOS to calculate density? In essence this would mean the algorithm would be the same as in the steady state case (momentum -> energy -> EOS solved for density -> press. cor).

In conclusion I agree it's still debatable whether you pull density out of the continuity equation at the start.

What would be the alternatives?

[FMDenaro]

Maybe the issue can be seen by looking at the continuity equation integrated in time:


rho(t1,x) =rho(t0,x) - Int [t0,t1] div (rho v) dt


This exact equation highlights that the condition at t0 is not sufficient to determine the density at t1 as you should know as both density and velocity change in the time interval (t0,t1). This change is fully coupled to the other variables.

The numerical integration can have some tricks. For example, using explicit methods, the RHS is computed without having information on the state at t1 and there is no apparent reason to perform a correction on the computed density. However, small numerical errors in the density can be strongly amplified (for example at low Mach) and a correction is sometimes applied.

[arjun]

Originally Posted by lorbekl

But isn't this the condition you have to keep at the end of each iteration in a SIMPLE algorithm?

Yes and no.

When your pressure solver finishes you get new pressure and thus you get new density.

But as soon as you called energy solver, now you shall have new density due to temperature change.

Now it is very important to update density at this point but as soon as you do, your continuity no longer satisfies.

Originally Posted by lorbekl

Would it be possible to just use the value from the previous time step/iteration and then use only the EOS to calculate density? In essence this would mean the algorithm would be the same as in the steady state case (momentum -> energy -> EOS solved for density -> press. cor).

Solver is very very sensitive to density, pressure and temperature trio here. You might solve some cases by the approach you said but anything difficult and expect to run into trouble.

Originally Posted by lorbekl

What would be the alternatives?

For the reasons mentioned, in my solver Wildkatze there is a special version of segregated solver available for combustion problems.

In this we solve

1. Momentum
2. Continuity (mainly due to velocity change)

3. Energy
4. Again continuity (from the system set up in step 3), I just update the source term and solve the system again, which got disturbed due to energy.

Now it is consistent for pressure, density and temperature.


The solver shows very good robustness so far.

[lorbekl]

Arjun, if I may ask, when is the equation of state used in your solver?

[lorbekl]

Unfortunately I still haven't been able to get a transient solution going. I would be grateful if I could get any opinions on it because it seems that I am not considering something correctly. This is the current algorithm that is causing me problems.

-> Load all variables from a steady-state solution (p,h0,h,rho,u) - let's call them PHI.

start time loop

-> State that PHI_0=PHI (variables at old time level are those currently in the memory)
-> State that PHI*=PHI (first guess for inner iteration is also equal to last value in memory)

start inner loop

-> Again state that PHI*=PHI (this line is skipped in first iteration)

-> Obtain density from the continuity equation using old values of velocity u* (continuity equation is discretized with central differencing scheme on regular mesh (u is stored on cell faces - staggered) and implicit in time)

-> Obtain velocity from momentum equation based on new density (for the transient term), old pressure, old mass fluxes (for convective terms) and old friction factor. Momentum eq. is discretized with 1st. order upwind on staggered mesh and implicit in time.

-> Obtain total enthalpy from energy equation based on new density (for transient term), old mass fluxes (for convective term). It's discretized with 1st. order upwind on regular mesh and is implicit in time.
-> new enthalpies are determined based on new total enthalpy and new velocity (h0=h+0.5*u^2)

-> Obtain pressure and velocity corrections by inserting momentum eq. into continuity eq. (SIMPLE). I only do a pressure correction p' and not a density correction (rho') because Ma is always < 0,3. It has done the trick in steady case but I will probably implement rho' aswell.


-> compute mass fluxes from new density and new velocity


-> update density based on last h and p with the EOS


end inner loop
t=t+dt
end time loop


BC's: Mass flux m'' and total enthalpy h0 at inlet face and pressure at the last node (not face) so that it coincides with the last face of the staggered grid. I extrapolate pressure to inlet/outlet. At the inlet p,h0 and m'' give sufficient info to determine all other variables, at the outlet I also extrapolate h0 and determine m'' from the continuity eq. The boundary values are updated everytime p,h0 or m'' change.

I notice that the problem begins after the pressure correction, because it is too large for a case where no change should be present anymore. The velocity correction then affects mass flux reconstruction and everything goes down from there. Beside not having a density correction in the pressure correction eq. does anyone see anything else I could be missing (and that could effect the magnitude of the pressure correction)?

[lorbekl]

Remember you are implicit in time. You need a value for density before you can calculate anything else. Density appears in the continuity equation and momentum equation. You cannot solve the momentum equation for velocity unless you already have an initial guess for density. How you get this initial guess can be done a few ways.

I've been following up on this topic and as I understand it now, the continuity equation is solved as a transport equation if you are using a (coupled) density based approach and not a (seggregated) pressure based approach. In the latter the continuity is worked into the pressure correction equation. Am I right? And also does that mean that the initial guess for density can be something as trivial as saying it's the same value as at the end of last step?

And in regard to that, I read an article by F. Moukalled where he states that the compressible SIMPLE algorithm is as follows:

-> solve momentum with old pressure and density
-> solve pressure correction eq. and update u,p,rho
-> solve energy equation
-> update density

Would you consider this a valid procedure?