[cryabroad]

Hello Foamers,

This has been a question that bothers me, why do we need to solve for density in the rhoPimpleFoam solver? I searched a lot on this forum but the threads I found are not that clear.

For the rhoPimpleFoam solver, here's what I understand. The solver first solves the momentum equation (UEqn.H), then energy equation (EEqn.H), then the pressure equation (pEqn.H) and additional equations like turbulent ones. At the end of the EEqn.H, there's thermo.correct(), which means properties like viscosity, thermal diffusivity and temperature are updated. At the beginning of the pEqn.H, there's rho = thermo.rho(), so density is updated before the pressure is solved.

However, after the pressure is solved, there's #include "rhoEqn.H" (line 84 of pEqn.H, OpenFOAM-4.x), then there's another rho = thermo.rho() line (line 91 of pEqn.H, OpenFOAM-4.x). Why do we need to solve an equation for density when density has already been (and should be) determined by the equation of state? Won't this give you a density field that is not consistent with the equation of state? Also, why do we need to update the density twice in the pressure equation? In rhoSimpleFoam the rho = thermo.rho() line only appears once at the end of the pEqn.H.

I feel like this can be a problem related to compressible flows. I've mainly dealt with incompressible flows so it's hard for me to go through these compressible solvers, although they use almost the same pressure-velocity coupling strategy.

Thanks in advance,

Ruiyan

[cryabroad]

Anyone? Still waiting for some explanations.

[cryabroad]

I think I can answer my own question, and maybe someone find this useful.

in the pEqn.H from rhoPimpleFoam, the first rho = thermo.rho() is to update the density and use it to calculate pressure. After pressure is solved, a "test" density field is solved based on the velocity field (#include "rhoEqn.H"), which is used just to obtain the residuals for the continuity equation. In other words, this "test" density field is compared with thermo.rho(). Then, the density field is overwritten again with the second rho = thermo.rho() operation.

But there's still a not-so-clear question, why do we need to update the density before solving the pressure equation? If we just update the density after solving the pressure equation, does it cause problems?

Ruiyan

[Tobi]

Hey Ruiyan, very interesting question. I do have the same questions in mind. However, in the latest openfoam releases, the rhoEqn is solved first and I am not sure why we need it.

The other topic regarding the density update before the pressure is solved might be related to the fact, that for an EOS idealGas, the density is based on pressure and temperature and hence, we will get an updated density with respect to the energy inside the system. Even though, the density is not fine until the pressure is corrected, we still have a better assumption based on an updated energy field.

By the way, the Thermo.correct() should also update the density field (did not check that peace already).

After the pressure equation is solved, we again recalculate the density and take it with (now) an updated energy and pressure field. Additionally, using an incompressible perfect gas EOS, the update of the density field before solving the pressure field does make definitly more sense than updating it only after the pressure is updated (as the density only depends on the energy field).


However, the rhoEqn is used twice sometimes (before any transport equation is solved and after the pressure is corrected).


However, there is still some parts missing. You mentioned that rhoEqn will just calculate an error? Actually, what does this gain or provide?

[Tobermory]

Hi Tobi, I tried commenting on your YT video, but for some reason it would not show up ... maybe there is some delay for long posts? Anyway, let me drop my thoughts here instead/as well. I have also been thinking through the pressure corrector method in the last few months, so it has been useful for me to put some thoughts down on paper - please correct me if you think I am wrong in my discussion below.

First - just a reminder that the continuity equation is always included when using the pressure corrector method, compressible or incompressible ... it's just that the continuity equation is very simple in the incompressible case and it is absorbed into the pressure correction calculation (div phi = 0). For the compressible case, we now need another equation to solve for rho, and that is the equation of state (EOS).

The real crux of the matter is the density value in the momentum predictor stage. For compressible flow cases, you are now predicting the value of rho.U, and so you start by making an estimate for U ... but what about rho? Then you form the pressure correction equation ... again, with what rho? Choices for these rho values include: the old value from the previous time step (ie consistent with the old pressure and temperature); a value consistent with the old pressure and the new temperature; a value that reflects the new source term in the density equation ... etc. The answer is ... it depends on the scheme that you are running!

For example, if you are running a steady calculation or an unsteady calc with simpleRho active then you only update rho at the end of the pressure corrector loop, using the EOS & latest p, T. The continuity equation is again simple (div phi = 0), so you can use the same continuity error measure (is div phi~0?) as the incompressible case. Haven't worked out yet how it copes with a density source term here, though ... (thoughts?)

If you are running transient (not simpleRho), then you calculate rho at the start of the first pimple using EOS and p,T from last iter, then calculate the new T, then for each p-corrector loop you update rho using the EOS and the latest available p and the new T, solve the pressure field and update the mass fluxes, and then you solve the rho equation with the new mass fluxes and latest source term. In other words, the rho value is kept "up to date" at each stage. If the solution is well converged then the difference between this final density (from rhoEqn.H) and the EOS density (rho.thermo()) should be negligible ... and this is the basis of the compressibleContinuityErrs.H error measure (is rho.rhoEqn - rho.thermo ~ 0?).

At least, that is what I have divined so far. I have not found any references that demonstrate why this is the optimum approach, or the problems that you encounter if you don't follow this approach .. again pls comment anyone if you have found something useful.

[Tobi]

Hi, thank you for your reply. Probably you refer to the last video I published regarding the continuity equation and pressure-corrector within compressible solvers in OpenFOAM, right? So I did not receive any comment so no idea why things are not working out here. Maybe you are right and long comments are somehow re-checked - no idea.

However, to the topic. I agree on most things you said. I recall and reference to rhoPimpleFoam in the latest dev version that is on commit a3c8514:

Transient (using simpleRho() or steady-state

• Rho Simple should be an acronym that we simplify the calculation and assume that the density will not be influenced by the pressure much. E.g., buoyancy of a few degree (only T dependent) vs shock waves (definitely p dependent)
• The procedure inside each PIMPLE loop
  • Construct the momentum matrix
  • Solve it if you use the predictor
  • Construct the energy matrix and solve it
  • Update the thermodynamic properties based on the new T-field
    --> thermo.correct()
    However, the rho field is not updated in the solver (just in the thermodynamic lib)
  • Update the density field only if we use SIMPLEC
  • Constructing the pressure matrix and solve it
  • Correct the fluxes
  • Calculate the continuity errors (just using the fluxes)
  • relax the pressure field
  • Recalculate the velocity with the relaxed pressure and updated HbyA and rAAtU guy
  • Limit the pressure if activated
  • Update the density (this field should not include any updated pressure data as we never call thermo.correct() again or, right?)

This procedure is fine for me, if the density is not changing much related to the pressure (hence simpleRho activation) and of course, in steady-state analysis, we are not really worried about any intermediate state in which the density and pressure field should fit 100 % together. However, it seems that using SIMPLEC requires and intermediate update of the density before we construct the pressure matrix (this density will only include the change due to energy but maybe brings stability while using the SIMPLEC algorithm).

Calculating transient without using simpleRho:

• Summarizing the procedure here, I would say, that we incorporate both, temperature changes as well as pressure changes into the density calculation.
• The procedure (in my opinion) within each PIMPLE loop:
  • Solve the continuity equation at the first pimple iteration (probably to include any source terms right at the beginning, ??? ) to get a new density field
  • Construct the momentum matrix
    Solve the predictor step
  • Construct the energy equation and solve it
    Update all thermodynamic data (thermo.correct())
  • Update the density field (grab it from the thermodynamic)
  • Build the pressure matrix and solve it
  • Limit the pressure
  • Explicitly correct the density inside the thermodynamic library (thermo.correctRho())
    We incorporate any density change based on the pressure now which - in my opinion - is required for shock waves / transonic / high-speed flows in which compression and expansion effects happen; maybe traveling sound waves etc. However, the rho field in the solver is still not updated
  • Update the rho field only if we limit the pressure (not sure why we don't update it, if we are not limiting the pressure? Which is in contrary to your statement. You said, we updated it always but the pressureControl() function will only return true if we explicitly limit the pressure
    https://cpp.openfoam.org/v8/pressure...ce.html#l00248
  • Calculate the continuity equation and obtain a new density field
    Again not sure: If we don't limit the pressure, this is the old rho field without pressure influence, otherwise it is an up-to-date density field
  • Checking the continuity error
    I am sorry, I did not checked the header-file and yes of course, we calculate the difference between the actual rho field in the solver and the updated one in the thermo-library.
  • Relax the pressure
  • Recalculate the U field
  • Relax the pressure field if we are not in a transonic regime
  • After finishing the pimple loop, update the rho field with the latest density inside the thermodynamic library

Things get clearer now. However my personal bottleneck. I hope I am formulating it in a way you all understand what I mean:

1. Talking about transient cases (without the usage of simpleRho), I am wondering, why we only update the rho field after the pressure correction if we limit the pressure. The continuity equation will be different if we would have high influence of the pressure onto the density field. But maybe this is just a small error, not sure about that. However, I would expect to have an updated rho field in any case (including T and p) while solving the continuity error.
2. Talking about steady-state cases, I am wondering why we first relax the pressure field and after that we limit the pressure. I would expect to have it vice versa as we have it in the transient case.
3. The if condition to update the rho field for steady-state, simpleRho or adjustMass can be modified. The adjustMass could be removed in my opinion, as this option will only be true if we have steady-state, which is almost checked by using the first argument.

[T. J. P.]

Hi,

I came here through the YouTube Video https://www.youtube.com/watch?v=C9oXLtvwtwo as well.

I always thought, that the rho Equation in these solvers is there as safety measure to ensure mass conversation. So although in theory the pressure equation should conserve mass by definition the same is not necessarily true for the discretized equations especially if more advanced stuff like moving or overlayed meshes or like you mentioned source terms are present. So to ensure conservation the corrected rho (after the p equation) are solved again. If all worked in the coupling beforehand, the rho equation is already fullfilled, hence usually you read residuals of 0 in 0 iterations in the log for rho, so the equation is not really solved here. Thus I assume it's more of a safety measure for rare cases which usually does nothing but check if the residual of rho is 0 (hence the 0 iterations). But this is just my interpretation of this euqation, if somebody knows more about it and the cases it is required please correct me .

[Tobermory]

Tobi

a few additions.

Transient using simpleRho() or steady-state
...
[*]Update the thermodynamic properties based on the new T-field
--> thermo.correct()
However, the rho field is not updated in the solver (just in the thermodynamic lib)

Edit: agreed.

[*]Update the density field only if we use SIMPLEC

yes - absolutely - I had left out the discussion of SIMPLEC.

[*]Update the density (this field should not include any updated pressure data as we never call thermo.correct() again or, right?)

I don't think that I agree with this point. There is only one pressure field, thermo.p(); in createFields we define the pointer

volScalarField& p = thermo.p();

and so p and thermo.p() are the same thing. The call to thermo.correct() ends up actioning function correct() in the chosen thermoType (eg hePsiThermo), and just calls the calculate() function in the same class. As discussed above, this calculate() function updates T first, based on latest energy and p values, and then updates psi, mu and rho if using heRhoThermo, all using latest p, T values. So the latest p field IS used for density here, I think.Edit okay, as is usual with OpenFOAM, the correct answer is more complex than it first seems. If you run psiThermo, then thermo.rho_ is based on the latest pressure; if you run rhoThermo then it is unchanged

However, it seems that using SIMPLEC requires and intermediate update of the density before we construct the pressure matrix (this density will only include the change due to energy but maybe brings stability while using the SIMPLEC algorithm).

Agreed.

Calculating transient without using simpleRho:

...
[*]Construct the energy equation and solve it
Update all thermodynamic data (thermo.correct())
[*]Update the density field (grab it from the thermodynamic)
[*]Build the pressure matrix and solve it
[*]Limit the pressure
[*]Explicitly correct the density inside the thermodynamic library (thermo.correctRho())

We incorporate any density change based on the pressure now which - in my opinion - is required for shock waves / transonic / high-speed flows in which compression and expansion effects happen; maybe traveling sound waves etc. However, the rho field in the solver is still not updated

Again, I read this differently. Function thermo.correctRho() is declared in class fluidThermo.H (https://cpp.openfoam.org/v8/fluidThermo_8H_source.html) as a virtual function, and set to zero ... funnily enough I cannot find any other redefinition, so I assume that this line is vestigial! The comment line above it in fluidThermo says "Add the given density correction to the density field. Used to update the density field following pressure solution", but I don't see any code in v8 or dev that does this!

[*]Update the rho field only if we limit the pressure (not sure why we don't update it, if we are not limiting the pressure? Which is in contrary to your statement. You said, we updated it always but the pressureControl() function will only return true if we explicitly limit the pressure
https://cpp.openfoam.org/v8/pressure...ce.html#l00248

Strictly speaking I didn't say that the density is "updated always", but at each stage (eg each pressure corrector of each PIMPLE) - I was obviously misleading though. If you make a pressure adjustment by limiting the pressure, then you are probably nudging the p/rho/phi balance out of balance, hence maybe the need to update rho at this stage. It gets calculated properly, of course, in the continuity equation straight afterwards, so the only benefit is in the initial rho field and maybe in the source term - may enhance convergence and/or stability.

[*]Calculate the continuity equation and obtain a new density field
Again not sure: If we don't limit the pressure, this is the old rho field without pressure influence, otherwise it is an up-to-date density field

no - see above, the density field has been updated with the latest pressure. Edit depends on thermoType

[*]After finishing the pimple loop, update the rho field with the latest density inside the thermodynamic library

I am confused by the discussion of fields "inside the thermodynamic library" ... I do not think that there are duplicate rho and p fields, one in the thermo library and one outside ... the ONLY fields are thermo.p and thermo.rho.Edit This last part of my comment is junk. There is only one pressure field, but there ARE separate rho fields - one (rho) is generated in createFields, the other belongs to the thermo (psiThermo or rhoThermo) and is only updated when thermo.correct() is called. So this is even more complex!

[Tobermory]

Originally Posted by T. J. P.

I always thought, that the rho Equation in these solvers is there as safety measure to ensure mass conversation. So although in theory the pressure equation should conserve mass by definition the same is not necessarily true for the discretized equations especially if more advanced stuff like moving or overlayed meshes or like you mentioned source terms are present. So to ensure conservation the corrected rho (after the p equation) are solved again. If all worked in the coupling beforehand, the rho equation is already fullfilled, hence usually you read residuals of 0 in 0 iterations in the log for rho, so the equation is not really solved here. Thus I assume it's more of a safety measure for rare cases which usually does nothing but check if the residual of rho is 0 (hence the 0 iterations). But this is just my interpretation of this euqation, if somebody knows more about it and the cases it is required please correct me .

Edit On reflection, a zero residual does not mean that the solver is doing nothing - refer to this post rho residual always zero, why? sonicFoam, OF 5.x. Because it is a diagonal solver, it just doesn't need to iterate - it has the solution already. So the rhoEqn.H solver IS updating the rho field (which as Tobi points out is different to the thermo.rho field) ... and that is why you need check the difference between the two (i.e. rho based on fluxes, and rho based on thermo) in compressibleContinuityErrs.H, as a check on the continuity error.

[Tobi]

Originally Posted by T. J. P.

Hi,

I came here through the YouTube Video https://www.youtube.com/watch?v=C9oXLtvwtwo as well.

I always thought, that the rho Equation in these solvers is there as safety measure to ensure mass conversation. So although in theory the pressure equation should conserve mass by definition the same is not necessarily true for the discretized equations especially if more advanced stuff like moving or overlayed meshes or like you mentioned source terms are present. So to ensure conservation the corrected rho (after the p equation) are solved again. If all worked in the coupling beforehand, the rho equation is already fullfilled, hence usually you read residuals of 0 in 0 iterations in the log for rho, so the equation is not really solved here. Thus I assume it's more of a safety measure for rare cases which usually does nothing but check if the residual of rho is 0 (hence the 0 iterations). But this is just my interpretation of this euqation, if somebody knows more about it and the cases it is required please correct me .

Hmmm, I thought a diagonal solver will directly give the solution, without iterating over the linear system as we directly can calculating the solution using the inverse matrix (https://www.openfoam.com/documentati...e-solvers.html). Hence, we have 0 iterations always.

Just wanted to change the rho-solver but it is not possible?

I don't think that I agree with this point. There is only one pressure field, thermo.p(); in createFields we define the pointer

      volScalarField& p = thermo.p();

Agree. I missed the reference of the pressure field to the thermodynamic By the way, it is a reference and not a pointer even though references are handled similar to pointers inside the code. Nothing special here, I know that you know that

Again, I read this differently. Function thermo.correctRho() is declared in class fluidThermo.H (https://cpp.openfoam.org/v8/fluidThermo_8H_source.html) as a virtual function, and set to zero ... funnily enough I cannot find any other redefinition, so I assume that this line is vestigial! The comment line above it in fluidThermo says "Add the given density correction to the density field. Used to update the density field following pressure solution", but I don't see any code in v8 or dev that does this!

It is defined in psiThermo.H or rhoThermo.H

no - see above, the density field has been updated with the latest pressure.

Okay. The pressure is a reference but we are not calling the

rho = thermo.rho();

if we don't limit the pressure, right? At least before we calculate the continuity errors.

I am confused by the discussion of fields "inside the thermodynamic library" ... I do not think that there are duplicate rho and p fields, one in the thermo library and one outside ... the ONLY fields are thermo.p and thermo.rho. Edit This last bit is junk. There is only one pressure field, but there ARE separate rho fields - one (rho) is generated in createFields, the other belongs to the thermo (psiThermo or rhoThermo) and is only updated when thermo.correct() is called. So this is even more complex!

Exactly, I was missing the solver p field which is indeed the same as in the thermo library. However, as you pointed out (and what I mentioned in my first post), the rho field is a scalar field inside the solver which gets updated ONLY if we do:

rho = thermo.rho();

Otherwise, we are not grabbing the data. In addition, to ensure that the density is really up-to-date, we first need to call:

thermo.correct();

In both, heRhoThermo and hePsiThermo.

[T. J. P.]

Originally Posted by Tobermory

Edit On reflection, a zero residual does not mean that the solver is doing nothing - refer to this post rho residual always zero, why? sonicFoam, OF 5.x. Because it is a diagonal solver, it just doesn't need to iterate - it has the solution already. So the rhoEqn.H solver IS updating the rho field (which as Tobi points out is different to the thermo.rho field) ... and that is why you need check the difference between the two (i.e. rho based on fluxes, and rho based on thermo) in compressibleContinuityErrs.H, as a check on the continuity error.

I haven't considered that, thanks for correcting me.

[Tobermory]

I always find it easier to draw out a flow diagram ... this is what I think the main structure of rhoPimpleFoam is, with different "columns" for the different options of steady, unsteady, rhoSimple etc.

The added complication here is that whether thermo.rho() depends on the latest pressure or not depends on your choice of thermoType ... if you choose psiThermo then it is; if you choose rhoThermo, then it is not.

[Krapf]

Tobi, is this the commit you mention in your video but couldn't find on the fly? https://github.com/OpenFOAM/OpenFOAM...6023bec2b66b50

[Tobermory]

Originally Posted by Tobi

Function thermo.correctRho() is declared in class fluidThermo.H (https://cpp.openfoam.org/v8/fluidThermo_8H_source.html) as a virtual function, and set to zero ... funnily enough I cannot find any other redefinition, so I assume that this line is vestigial! The comment line above it in fluidThermo says "Add the given density correction to the density field. Used to update the density field following pressure solution", but I don't see any code in v8 or dev that does this!
It is defined in psiThermo.H or rhoThermo.H

Aaah - what confused me is that the definition (in the .C file) is empty for psiThermo:

 void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
 {}

and only has any action for rhoThermo, where it adds the deltaRho to thermo.rho_:

 void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
 {
     rho_ += deltaRho;
 }

I updated my schematic to correct this.

[Tobi]

Originally Posted by Tobermory

Aaah - what confused me is that the definition (in the .C file) is empty for psiThermo:

 void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
 {}

and only has any action for rhoThermo, where it adds the deltaRho to thermo.rho_:

 void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
 {
     rho_ += deltaRho;
 }

I updated my schematic to correct this.

Probably in psiThermo there is nothing to do as we already calculate rho based on psi whereas for the rho Thermo, we get it from the EOS. Well, I need to have a closer look at this - especially psi and others.

[Tobi]

By the way, your chart is not 100 % in terms of "limit p and update thermo:rho_". Here, the update of rho is only performed if someone sets any value into the fvSolution file to trigger the if-function

[Tobermory]

Correct - but this rho value gets overwritten almost immediately when the rhoEqn is solved, so all it does is provide a better rho value for the source term in the rhoEqn, if there is one. I've updated the schematic to note this.

Interestingly, if you unpick the logic, then we have the following:

1. steady, SIMPLE run: T, psi and thermo.rho_ are updated after the Energy equation, using the old value for p; rho is set to thermo.rho_ at the end of each pressure corrector stage; so for true SIMPLE operation (1 pressure corrector), the pressure solution and momentum correction is carried out with the old value of density, from the previous time step. This is consistent with the approach in rhoSimpleFoam. Importantly, if the thermoType is hePsiThermo, then the rho value at the end of the time step is based on the new pressure; but if the thermoType is heRhoThermo then it is based on the old pressure.
2. steady SIMPLEC run: if running SIMPLEC, then rho is updated with the value of thermo.rho_ at the start of each pressure corrector loop, so that each pressure corrector is carried out using a density value that is based on: the updated temperature, and either the updated pressure (hePsiThermo) or the old pressure (heRhoThermo).
3. unsteady simpleRho run: as per the steady SIMPLE/SIMPLEC run options above, but with one difference: once the new pressure field has been solved for, and the mass fluxes updated, the continuity equation is solved to update rho; and so the momentum correction is carried out with a rho value that is consistent with the latest fluxes; this rho value then gets overwritten with the thermo.rho_ value at the end of the pressure correction loop, as per the steady run cases.
4. unsteady PIMPLE run: same approach as the unsteady, SIMPLEC, simpleRho case above, i.e.: update T, psi, thermo.rho_ in the energy equation, i.e. in each PIMPLE iteration; at the start of each pressure corrector, set rho to rho.thermo_ (using latest p if hePsiThermo, or old p if heRhoThermo), solve for p and phi; solve continuity equation for rho and calculate momentum correction; set rho to rho.thermo_ (again using latest p if hePsiThermo, or old p if heRhoThermo) and move to next pressure corrector.

and so we can see that the unsteady PIMPLE approach uses a consistently more up-to-date version of the density as it loops through the pressure correctors and PIMPLEs, wherease the SIMPLE steady approach is happy to use the old value (unless SIMPLEC is activated).

[Tobi]

Very well summarized.


For the simpleRho option, I found this commit https://github.com/OpenFOAM/OpenFOAM...6023bec2b66b50


Which states, that simpleRho is for SIMPLE to improve convergence and stability for cases and large-courant numbers/time-steps (probably if we do PIMPLE with Co >> 1).

[Tobi]

Originally Posted by Krapf

Tobi, is this the commit you mention in your video but couldn't find on the fly? https://github.com/OpenFOAM/OpenFOAM...6023bec2b66b50

Yes. This was indeed the guy I was searching for - I missed your post, unfortunately.

[Jaco_97]

Hi everyone, I wish somebody could shed some light on the doubts I've been struggling with during the last few days. I am particularly interested in understanding the substantial difference between rhoThermo class and psiThermo class.

As far as I understand, in both cases thermo.correct() function calls a calculate() function which is responsible for updating temperature and other thermophysical properties. In particular, differently from psiThermo class, rhoThermo updates density as well within the aforementioned function according to a specific expression which derives from one of the available equations of state, e.g., rho=p/(RT) in case of perfect gas EOS. Reading previous comments under this post, I understood that the pressure appearing in the above expression is kept fixed and this is what makes the difference between "variable-density" and "compressible" solvers.

Provided that my reasoning is correct (please correct me if I am wrong), what is the specific part of the code where it can be actually read that the pressure does not change while updating the density by means of thermo.rho() in rhoThermo-based solvers?

Thanks in advance,
Jacopo