How rho_ and psi_ are calculated in compressible solvers of OpenFOAM?
Dear FOAMers,
I need to know how the density and compressibility are updated in OpenFOAM? The question has a general scope and applies to many compressible thermophysical models and solvers. I divide my question in two parts: 1. In pEqn.H of all compressible solvers (dieselFOAM, sonicFOAM, reactingFOAM, etc.) there is a line that reads Code:
rho = thermo.rho(); Code:
// Density [kg/m^3]  uses current value of pressure 2. Assuming that my impression of rho_ calculation is correct, we need to know psi() to find density, that is we need to know the compressibility. This raises two more questions: 2.1 My survey led to the speculation that psi_ is updated in hsPsiThermo.C. Code:
TCells[celli] = mixture_.THs(hsCells[celli], TCells[celli]); 2.2 If above is correct, then I need to know how the psi(p, T) works. I could find THs function in specieThermoI.H which basically updates temperature based on a given sensible ehthalpy and known composition mixture by an iterative method. However, I couldn't locate such a function for psi_ calculation. The only thing I could think of was the equation of state for perfect gas in perfectGasI.H Code:
inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const I appreciate it, if you could answer any of the above three questions. Any clue will be of great value to me. Thanks in advance... 
You are nearly 100% correct. This is where psi is calculated for a compressible thermophysical model. That is the reason its located in perfectGasI.H. If you look at the other equations of state, there are definitions of how psi should be calculated based on their assumptions.
The reason everything is scattered about is because the model is runtime selectable, and thus you need to be able to construct the thermo model from its various components. The user manual explains is fairly well on this page: http://www.openfoam.org/docs/user/th...#x372050007.1 
Hi,
Sincere thanks to Marco for quick reply. I'm afraid I have a few more questions. I would be grateful, if you or other OF experts could get involved in the discussion, as these issues are fundamental and advanced features of OF. :cool: The reason I am skeptical about my impression of psi_ calculation is due to the following: 1. As you can see in the code below, the function psi(p, T) in hsPsiThermo.C takes "pressure" and temperature of the current cell i and returns the psi value: Code:
psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]); Code:
inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const 1.2 Having said this, you can see from the above code that the syntax of this function is psi(scalar, scalar T). This means that the first argument is redundant or dummy (I don't know what they call this in C++). If we don't need the first argument, then why we mention it at all?! I guess this has something to do with maintaining a generic syntax for run time selection tables. Perhaps there should be some other implementations of psi(?, T) that do accept something as their first argument, but I don't know where are they and what are they?! It helps if you can name a few. 1.3 The same redundant first argument is preserved for the implementation of psi in the two incompressible equations of state models in IncompressibleI.H and icoPolynomialI.H. The psi implementation for both of these classes return zero, hence incompressible. For the latter, both arguments are dummy or redundant (or whatever they call such arguments). Why we have to keep these redundancies in syntax of our functions at all?! 2. The second major confusion I have about psi_ and rho_ update process lies in the presence of a rho_ calculation implementation in perfectGasI.H, as you can see below. If rho_ will be updated in basicPsiThermo.H, as I mentioned in my opening post in this thread, then why there should be another implementation of rho in perfectGasI.H?! :confused: Code:
inline Foam::scalar Foam::perfectGas::rho(scalar p, scalar T) const 
I've got another project I'm working on, so the answers will have to be a little short (just so you know I'm not being rude).
1. You have answered your own question (the FOAM is strong with this one). Maintaining the generic interface is what its all about, and depending on the types that can be created from the template class, the compiler will sort out what needs to be called. That is the only reason we keep the two parameters even if the first one (or none at all!) is discarded. That way we shouldn't have to recompile the code just because we want to use a different thermo model. 2. If you look at what started this discussion, you were wondering where psi comes from in the calculation rho = p*psi (in the particular solver you were looking at). My feeling is that because of the way the solver is structured, we need the boundary conditions on p (as well as the internal values) to properly calculate rho on the boundaries as well. Having both versions of the psi function gives us the flexibility to do the rho = p*psi. What I'm less sure about is how the compiler decided which one to call; that is a whole other level of code voodoo I haven't reached yet (but somehow I feel that hash tables and template classes are involved). Hope this helps! 
Thanks dear Marco for your comments. I appreciate it.
I will work on these questions, and will report here if any progress made. Will try to keep the explanatory format of the post, making it worth reading for future generations! Meanwhile, it will help a lot, if anybody likes to step forward and carve his signature on this stone. Jalal 
Dear Foamers,
I found the answer to some of my questions: 1. rho update process: Quote:
2.1 the psi calculation procedure: Quote:
2.2 the psi function: Quote:
Cheers, Jalal 
Hi guys,
I have a question regarding density calculation in reactingFoam as well. When I simulate a premixed combustion in a channel flow I should see a considerable change in density where the temperature has increased. Since pressure is constant, when I see the temperature has increased per say 5 times, then following the ideal gas law, the density should decrease 5 times, then following the continuity, the flow velocity should increase 5 times. But I don't see that difference in velocity. It looks like an incompressible flow!! I think reactingFoam is missing something and density doesn't seem being updated!!! Correct me if I'm wrong. ANY CLUE? 
Quote:
This post is very helpful for me to understand the psi. As you know, the compressibility psi is calculated through: inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const { return 1.0/(R()*T); } In combustion situation, the gas constant differs at each cell center and is a function of composition at each point. Do you how R() is calculated in combustion solver like reactingFoam? According to my understanding, the R=RR/W, here W is mixture average molecular weight and if we want to get it we need to know the mass fraction of all the (at least major) species. However, I failed to find this line of reasoning in Openfoam due to my not very proficient C++ skills. I really appreciate it if you can give me some hints about this problem. best regards, H 
How gas constant R is calculated in OpenFOAM?
Hi,
This is a very good and fundamental question. The answer, however, is not that easy and straightforward. I can only give some general tips, and you should work out the details yourself. OpenFoam uses a data structure named speciesData_ which is of type PtrList<ThermoType> declared in multiComponentMixture.H. This PtrList is filled in with janafThermo information in read function of multiComponentMixture.C. I forgot much of the details, but there is a chemkin lexer that reads in the species related data at run time. Depending on the species involved in the reaction mechanism, their data is read from therm.dat and chem.inp input files containing thermal and chemical information needed in a chemkin standard format. Atomic weight of each species are calculated by using the reference file "atomicWeights" under "species" directory. A good exercise would be Infoing the speciesData_ content. For example, you can insert this in cellMixture function of the multiComponentMixture.C. Code:
Info << "speciesData_[" << n << "] = " << speciesData_[n] << endl; Code:
template<class ThermoType> By scrutinising the contents of speciesData_ and MixtureSoFar as the loop progresses through the species in the chemical mechanism (I recommend a global mechanism for simplicity), you can understand that the second number in these arrays is the molecular weight. It's a bit tricky, but by looking closely you can unravel the pattern and underlying algorithm for calculating mixture weighted quantities. This is because OF is written in a highly efficient way to avoid reimplementing repeated pieces of code. The beauty of this simple += operator in the cellMixture function above is that, on the one hand, it preserves the mole based nature of the a[0] to a[5] coefficients of janafThermo functions in janafThermoI.H, by multiplying the final result by mixture W. On the other hand, it is used to calculate the molar weight of the mixture as well. Note that the function cellMixture calculates mixture properties for cell i. But the same can be applied to any arbitrary thermodynamic state independent of the mesh. An example of such use of += operator can be found in ODEChemistryModel: Code:
// update the temperature Code:
inline scalar specie::R() const Jalal 
Hi,
Thanks for you post about how rho is calculated, it cleared some things up for me. But where exactly in the solver does it call for psi to be updated? I can find the function for it in the thermophysical model, but I don't see where the model itself updates the psi value. In particular I'm looking at the fireFoam model, but I can't find it in any other compressible model either. Also, how would I force the model to update the psi value? 
all the thermo stuff like rho and psi is updated in
Code:
thermo.correct(); Code:
calculate() That's then finally the place where the magic happens. ;) Maybe to add: In some solvers rho is also updated in the pEqn.H, often something like: Code:
rho = thermo.rho(); 
Thanks for the answer, however thermo.correct() doesn't seem to be updating my psi. If I understand it correctly, the psi and T values are used to update the pressure in the thermo part? Basically, I'm trying to run fireFoam but hold the density constant, which so far has resulted in my pressure being constant as well, despite the 1800 degree temperatures. I've had some cases where it's obvious that T is being affected by the pressure (low pressures resulted in low temperatures), yet not the other way around.
Does thermo.correct only update psi? In this case, at which point does the psi update the pressure? I think there's something I'm missing here. 
Okay, answering my own question here. Looking further into it, psi is used in the typical compressible pressure equation in calculating pressure. thermo.correct() corrects psi based on the current temperature, which in turn affects pressure as the pressure equation is solved. Is there any straight forward way to use the incompressible pressure equation yet update pressure from the thermo properties based on T and psi?

Hi,
I have already answered this in post #6 above. I don't know about FireFoam, but as long as you're using "perfectGas" model the psi_ is updated by invoking the psi(p, T) function in perfectGasI.H. My understanding is that psi_, compressibility, actually doesn't depend on p, and is a function of mixture composition and temperature only: Quote:
1. T based on the value of h, e, or hs calculated in energy eq. 2. having updated T, it updates psi, mu, and alpha both for the internal and boundary cells and faces. Hope this helps, Jalal 
Thanks yes that cleared some things up. I now realise my constant pressure is resulting from the boundary conditions, nothing to do with the thermo. Good to know how it works though.

Dear Sir,
Above threads are very helpful for me to use reactingFoam. And then, I would like to ask you a question (it may be fundamental for you). I am using myReactingFoam solver, which I customize as I can calculate the additional source term of YEqn.H. When I checked the calculated rho at one cell, the value was different from the value which I calculated using the equation of state. For example, at one cell, when p = 40000, R = 288 (which is calculated from RR/W, W is calculated from Σ(Yi/Mi)), T = 293, rho calculated using perfectGas equation is 0.474. However, calculated rho by OF at the same cell is about 0.0657, which is about seven times smaller than 0.474. Why does this happen? Does this have relation to the equation of continuity; rhoEqn.H? Even if so, is it strange that rho calculated by OF is inconsistent a lot with that calculated by equation of state? I am looking forward to your comment. Sincerely, KeiJun 
All times are GMT 4. The time now is 08:38. 