Hello everyone,
let's try to start the discussions in this group by posting about a topic that might be of general interest for people using thermodynamic libraries in their solvers. You will see that some things are still unclear so please comment and help finding the answers.
Since the release of OF 1.6 I was asking myself the above question and the documentation in the header files "Basic thermodynamic properties based on compressibility/density" did not really help to clarify things for me.
So I took a look into the source code and started from the end:
1) The solvers reactingFoam and rhoReactingFoam use psiChemistryModel and rhoChemistryModel which then use
hsCombustionThermo and hsReactionThermo, respectively. Both pairs of classes are different only by the name not by the way the functions are implemented.
2) They use hsPsiMixtureThermo and hsRhoMixtureThermo to calculate the values for T, psi, mu and alpha as a mixture of the components' cell values. How the cell mixture itself is calculated depends on multiComponentMixture, dieselMixture, etc. which is used by both hsPsiMixtureThermo and hsRhoMixtureThermo. The two classes are actually different only in one line: hsRhoMixtureThermo calculates also the density as a cell mixture of the component densities.
3) The actual difference can be found in basicPsiThermo and basicRhoThermo:
When you call thermo->rho(), basicPsiThermo will return rho=p*psi. Psi is the ideal gas compressibility psi=1/RT if you use "perfectGas" or zero if you use icoPolynomial. As a result, icoPolynomials should always return rho=0 in psiThermo!!
On the other hand, basicRhoThermo will return a rho field that is stored as a private variable. This means that when using rhoThermo, you might have two different density fields, one in the solver and one in the rhoThermo class. To avoid confusing both fields in the object registry, the rhoThermo class names its rho field "rhoThermo". The rho field is updated when you call thermo->correct() and rho is calculated as a cell mixture of the components' rho values. If you select perfect gas model, then the component's rho will be rho= p/RT.
So in the case of ideal (perfect) gases, the difference between psiThermo and rhoThermo is: nothing!
(Well, there is a small difference, since rhoThermo is updated once every time step and psiThermo will always use the up-to-date pressure field)
The difference arises when using icoPolynimial for simulating liquids: Here, psiThermo will return rho=0! (Why is psiChemistryModel also defined for icoPoly8ThermoPhysics??? It should not work!). rhoThermo will just return the value for rho from the defined polynomial. As I understand this is the only difference between psiThermo and rhoThermo, since other fluid properties are calculated as cell mixtures e.g. by multiComponentMixture.
Another difference is in the PISO loop of the solver application: In pEqn.H of rhoReactingFoam, p*psi is first subtracted from rho, then the p equation is solved with an additional correction term and then p*psi is again added to rho. The purpose of this in unclear to me. Is it for numerical stability purposes??
As a conclusion, I will prefer rhoThermo for my solvers, since only this one seems to be valid for both liquids and gases...
Which one are you using?
Best Regards,
Dominik |
