CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

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

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

Like Tree9Likes
  • 3 Post By JBUNSW
  • 1 Post By mturcios777
  • 1 Post By mturcios777
  • 3 Post By JBUNSW
  • 1 Post By JBUNSW

Reply
 
LinkBack Thread Tools Display Modes
Old   July 10, 2012, 12:29
Default How rho_ and psi_ are calculated in compressible solvers of OpenFOAM?
  #1
New Member
 
Jalal Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 23
Rep Power: 5
JBUNSW is on a distinguished road
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<volScalarField> 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 and PicklER like this.
JBUNSW is offline   Reply With Quote

Old   July 10, 2012, 12:59
Default
  #2
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 686
Rep Power: 17
mturcios777 will become famous soon enough
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.
mturcios777 is offline   Reply With Quote

Old   July 10, 2012, 22:24
Default
  #3
New Member
 
Jalal Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 23
Rep Power: 5
JBUNSW is on a distinguished road
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);
 }
JBUNSW is offline   Reply With Quote

Old   July 11, 2012, 12:54
Default
  #4
Senior Member
 
mturcios777's Avatar
 
Marco A. Turcios
Join Date: Mar 2009
Location: Vancouver, BC, Canada
Posts: 686
Rep Power: 17
mturcios777 will become famous soon enough
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.
mturcios777 is offline   Reply With Quote

Old   July 11, 2012, 19:48
Default
  #5
New Member
 
Jalal Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 23
Rep Power: 5
JBUNSW is on a distinguished road
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
JBUNSW is offline   Reply With Quote

Old   July 24, 2012, 09:04
Default
  #6
New Member
 
Jalal Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 23
Rep Power: 5
JBUNSW is on a distinguished road
Dear Foamers,

I found the answer to some of my questions:

1. rho update process:
Quote:
Originally Posted by JBUNSW View Post
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<volScalarField> 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 View Post
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 View Post
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
fabian_roesler, hz283 and PicklER like this.
JBUNSW is offline   Reply With Quote

Old   December 11, 2012, 15:33
Default
  #7
New Member
 
Mhsn
Join Date: Oct 2012
Posts: 24
Rep Power: 4
mhsn is on a distinguished road
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?
mhsn is offline   Reply With Quote

Old   January 11, 2013, 20:32
Default
  #8
Senior Member
 
Join Date: Nov 2012
Posts: 168
Rep Power: 4
hz283 is on a distinguished road
Quote:
Originally Posted by JBUNSW View Post
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.

I really appreciate it if you can give me some hints about this problem.

best regards,
H
hz283 is offline   Reply With Quote

Old   January 13, 2013, 04:03
Default How gas constant R is calculated in OpenFOAM?
  #9
New Member
 
Jalal Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 23
Rep Power: 5
JBUNSW is on a distinguished road
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;
Info << "MixtureSoFar = " << mixture_ << endl;
You have to put the above inside the following loop:

Code:
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellMixture
 (
     const label celli
 ) const
 {
     mixture_ = Y_[0][celli]/speciesData_[0].W()*speciesData_[0];
 
     for (label n=1; n<Y_.size(); n++)
     {
         mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];
     }
 
     return mixture_;
 }
Then wmake the reactionThermo directory for your changes to be applied. Then try running any solver that uses multiComponentMixture class. For example, all the solvers that use ODEChemistryModel can reveal this structure.

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
             const scalar cTot = sum(c);
             ThermoType mixture(0.0*specieThermo_[0]);
             for (label i=0; i<nSpecie_; i++)
             {
                 mixture += (c[i]/cTot)*specieThermo_[i];
             }
             Ti = mixture.TH(hi, Ti);
Now that you know how the molar weight of the mixture is calculated(i.e., mixture_.W()), you can see that gas constant for the mixture can be easily obtained by dividing the universal gas constant by molar weight in specieI.H:

Code:
inline scalar specie::R() const
 {
     return RR/molWeight_;
 }
Hope this helps.

Jalal
JBUNSW is offline   Reply With Quote

Old   April 5, 2013, 09:14
Default
  #10
New Member
 
Cameron
Join Date: Jul 2012
Posts: 22
Rep Power: 4
c_dowd is on a distinguished road
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?
c_dowd is offline   Reply With Quote

Old   April 5, 2013, 09:18
Default
  #11
Senior Member
 
dkxls's Avatar
 
Armin
Join Date: Feb 2011
Location: Helsinki, Finland
Posts: 118
Rep Power: 9
dkxls will become famous soon enough
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();
dkxls is offline   Reply With Quote

Old   April 5, 2013, 21:31
Default
  #12
New Member
 
Cameron
Join Date: Jul 2012
Posts: 22
Rep Power: 4
c_dowd is on a distinguished road
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 21:56.
c_dowd is offline   Reply With Quote

Old   April 5, 2013, 22:11
Default
  #13
New Member
 
Cameron
Join Date: Jul 2012
Posts: 22
Rep Power: 4
c_dowd is on a distinguished road
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?
c_dowd is offline   Reply With Quote

Old   April 5, 2013, 22:18
Default
  #14
New Member
 
Jalal Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 23
Rep Power: 5
JBUNSW is on a distinguished road
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 View Post
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
PicklER likes this.
JBUNSW is offline   Reply With Quote

Old   April 6, 2013, 00:55
Default
  #15
New Member
 
Cameron
Join Date: Jul 2012
Posts: 22
Rep Power: 4
c_dowd is on a distinguished road
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.
c_dowd is offline   Reply With Quote

Reply

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

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Need help with boundary conditions: open to atmosphere Wolle OpenFOAM 2 April 11, 2011 08:32


All times are GMT -4. The time now is 01:22.