How rho_ and psi_ are calculated in compressible solvers of OpenFOAM?

 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 10, 2012, 11:29 How rho_ and psi_ are calculated in compressible solvers of OpenFOAM? #1 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 6 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();` I tracked down this rho() function to see how density is calculated. For the compressible thermophysical models, in the base class that is basicPsiThermo.H, one can find this: Code: ``` //- Density [kg/m^3] - uses current value of pressure virtual tmp rho() const { return p_*psi(); }``` Is my understanding of how rho_ is updated correct?! 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]); psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);``` This is true when sensible enthalpy is used as the basis of thermophysical model. For scenarios where h or e are used similar procedure follows. Is my understanding correct?! 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 { return 1.0/(R()*T); }``` But I'm not sure, and this in itself raises additional questions. Is this where psi_ is calculated or I'm wrong?! 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... dongchao yang, hz283, Lisandro Maders and 1 others like this.

 July 10, 2012, 11:59 #2 Senior Member     Marco A. Turcios Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 727 Rep Power: 18 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 run-time 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...#x37-2050007.1 JBUNSW likes this.

 July 10, 2012, 21:24 #3 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 6 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. 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]);` Furthermore, considering the definition of psi function in perfectGasI.H, I can break my first question into the following three detailed questions: Code: ```inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const { return 1.0/(R()*T); }``` 1.1 This function has no dependence on pressure; then why psi(p, T), in hsPsiThermo.C, passes both "pressure" and temperature to the psi function?! 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?! Code: ```inline Foam::scalar Foam::perfectGas::rho(scalar p, scalar T) const { return p/(R()*T); } inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const { return 1.0/(R()*T); }``` Lisandro Maders likes this.

 July 11, 2012, 11:54 #4 Senior Member     Marco A. Turcios Join Date: Mar 2009 Location: Vancouver, BC, Canada Posts: 727 Rep Power: 18 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! PicklER likes this.

 July 11, 2012, 18:48 #5 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 6 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

July 24, 2012, 08:04
#6
New Member

Join Date: Oct 2011
Location: Sydney, Australia
Posts: 27
Rep Power: 6
Dear Foamers,

I found the answer to some of my questions:

1. rho update process:
Quote:
 Originally Posted by JBUNSW 1. In pEqn.H of all compressible solvers (dieselFOAM, sonicFOAM, reactingFOAM, etc.) there is a line that reads Code: `rho = thermo.rho();` I tracked down this rho() function to see how density is calculated. For the compressible thermophysical models, in the base class that is basicPsiThermo.H, one can find this: Code: ``` //- Density [kg/m^3] - uses current value of pressure virtual tmp rho() const { return p_*psi(); }``` Is my understanding of how rho_ is updated correct?!
Yes, this is how density is updated when the thermo.rho() is invoked.

2.1 the psi calculation procedure:
Quote:
 Originally Posted by JBUNSW 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]); psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);``` This is true when sensible enthalpy is used as the basis of thermophysical model. For scenarios where h or e are used similar procedure follows. Is my understanding correct?!
Yes. This is correct too, with one minor correction. The invoking function can be hsPsiMixtureThermo.C, depending on the thermophysical model chosen.

2.2 the psi function:
Quote:
 Originally Posted by JBUNSW 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 { return 1.0/(R()*T); }``` But I'm not sure, and this in itself raises additional questions. Is this where psi_ is calculated or I'm wrong?!
Yes. That is where psi_ is updated if perfectGas is taken for EquationOfState in building the thermophysical model.

Cheers,
Jalal

 December 11, 2012, 15:33 #7 New Member   Mhsn Join Date: Oct 2012 Posts: 24 Rep Power: 5 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?

January 11, 2013, 20:32
#8
Senior Member

Join Date: Nov 2012
Posts: 168
Rep Power: 5
Quote:
 Originally Posted by JBUNSW Dear Foamers, I found the answer to some of my questions: 1. rho update process: Yes, this is how density is updated when the thermo.rho() is invoked. 2.1 the psi calculation procedure: Yes. This is correct too, with one minor correction. The invoking function can be hsPsiMixtureThermo.C, depending on the thermophysical model chosen. 2.2 the psi function: Yes. That is where psi_ is updated if perfectGas is taken for EquationOfState in building the thermophysical model. Cheers, Jalal
Hi JBUNSW,

This post is very helpful for me to understand the psi.

As you know, the compressibility psi is calculated through:

inline Foam::scalar Foam:erfectGas:si(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.

best regards,
H

 January 13, 2013, 04:03 How gas constant R is calculated in OpenFOAM? #9 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 6 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 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; Info << "MixtureSoFar = " << mixture_ << endl;``` You have to put the above inside the following loop: Code: ```template const ThermoType& Foam::multiComponentMixture::cellMixture ( const label celli ) const { mixture_ = Y_[0][celli]/speciesData_[0].W()*speciesData_[0]; for (label n=1; n

 April 5, 2013, 08:14 #10 New Member   Cameron Join Date: Jul 2012 Posts: 22 Rep Power: 5 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?

 April 5, 2013, 08:18 #11 Senior Member     Armin Join Date: Feb 2011 Location: Helsinki, Finland Posts: 152 Rep Power: 10 all the thermo stuff like rho and psi is updated in Code: ` thermo.correct();` It usually is called in the hEqn.H and the function itself calls another one called 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();`

 April 5, 2013, 20:31 #12 New Member   Cameron Join Date: Jul 2012 Posts: 22 Rep Power: 5 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. Last edited by c_dowd; April 5, 2013 at 20:56.

 April 5, 2013, 21:11 #13 New Member   Cameron Join Date: Jul 2012 Posts: 22 Rep Power: 5 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?

April 5, 2013, 21:18
#14
New Member

Join Date: Oct 2011
Location: Sydney, Australia
Posts: 27
Rep Power: 6
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:
 Originally Posted by JBUNSW Code: ```inline Foam::scalar Foam::perfectGas::psi(scalar, scalar T) const { return 1.0/(R()*T); }```
As for the thermo.correct(), it updates the following:

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

 April 5, 2013, 23:55 #15 New Member   Cameron Join Date: Jul 2012 Posts: 22 Rep Power: 5 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.

 Tags compressibility, density, psi_, rho_, thermo.rho()

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules