Why we need rhoEqn.H in rhoPimpleFoam
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 |
Anyone? Still waiting for some explanations.
|
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 |
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? |
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. |
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
Calculating transient without using simpleRho:
Things get clearer now. However my personal bottleneck. I hope I am formulating it in a way you all understand what I mean:
|
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 :). |
Tobi
a few additions. Quote:
Quote:
Quote:
Code:
volScalarField& p = thermo.p(); Quote:
Quote:
Quote:
Quote:
Quote:
|
Quote:
|
Quote:
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? :confused: Quote:
Quote:
Quote:
Code:
rho = thermo.rho(); Quote:
Code:
rho = thermo.rho(); Code:
thermo.correct(); |
Quote:
I haven't considered that, thanks for correcting me.:) |
1 Attachment(s)
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. |
Tobi, is this the commit you mention in your video but couldn't find on the fly? https://github.com/OpenFOAM/OpenFOAM...6023bec2b66b50
|
1 Attachment(s)
Quote:
Code:
void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho) Code:
void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho) |
Quote:
|
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
|
1 Attachment(s)
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:
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). |
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). |
Quote:
|
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 |
All times are GMT -4. The time now is 08:43. |