[JBUNSW]

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

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:

             //- 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.

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

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...

[mturcios777]

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]

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:

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:

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?!

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);
 }

[mturcios777]

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!

[JBUNSW]

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]

Dear Foamers,

I found the answer to some of my questions:

1. rho update process:

Originally Posted by JBUNSW

1. In pEqn.H of all compressible solvers (dieselFOAM, sonicFOAM, reactingFOAM, etc.) there is a line that reads

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:

             //- 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:

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.

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:

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

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

[mhsn]

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?

[hz283]

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.

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

best regards,
H

[JBUNSW]

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.

Info << "speciesData_[" << n << "] = " << speciesData_[n] << endl;
Info << "MixtureSoFar = " << mixture_ << endl;

You have to put the above inside the following loop:

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:

// 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:

inline scalar specie::R() const
 {
     return RR/molWeight_;
 }

Hope this helps.

Jalal

[c_dowd]

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?

[dkxls]

all the thermo stuff like rho and psi is updated in

    thermo.correct();

It usually is called in the hEqn.H and the function itself calls another one called

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:

rho = thermo.rho();

[c_dowd]

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.

[c_dowd]

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?

[JBUNSW]

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:

Originally Posted by JBUNSW

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

[c_dowd]

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.

[KeiJun]

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

[bbita]

Dear Foamers,

I add Y to the main program in compressibleMultiphaseInterFoam through following post:
Define simple multicomponent composition
I changed Y and I see in the output results, Y of different components in different phases are changing but density and viscosity are not updating. Any suggestion?

Many thanks,

[Kedar G. Bhide]

I am using OpenFoam 6, and therefore unable to locate the place of calculation of rho as
given in this thread. The structure has changes I think between the versions. I am trying to locate the calculation of rho in tutorial case 'rhoPimpleFoam/RAS/cavity' in OpenFoam 6.
Anybody has idea how it is done?
I have looked inside 'equationOfState/perfectGas' but it is not the one.

Thank you,
Kedar