[cryabroad]

Hi Guys,

First of all, thank you for taking the time to read this thread.

I've been learning OF for the past several months and for the sake of practicing programming, I am trying to understand how OF calculates a volScalarField.

For example, for compressible flows using the psiThermo class and the ideal gas law, according to OpenFOAM, the density rho (which is a volScalarField) is calculated by rho = psi*p, where psi = 1/(R*T). This is documented in the psiThermo.C file, line 93 to 102. This works fine because all the variables are declared as volScalarField. So basically OpenFOAM is calculating a volScalarField (rho) based on other volScalarField (psi and p).

My questions is, how does OpenFOAM actually process the calculation? It still need to do something like rho[celli] = psi[celli]*p[celli] right?

To test this idea, I decide to replace the original rho() function (line 160) in the rhoThermo.C file.The reasons why I chose to modify the rhoThermo class instead of the psiThermo class are as following. First, mesh is passed into the rho() function as its argument in the rhoThermo class, this way I can use information associated with mesh; second, rho_ is a protected data in the rhoThermo class , while rho_ is not decleared in the psiThermo class. I replaced the rho() function with the following code

forAll(mesh.V(), celli)
{
    rho_[celli] = psi_[celli]*p_[celli]
}

When I recompiled the class I get the following error:

assignment of read-only location ‘((const Foam::testrhoThermo*)this)->Foam::testrhoThermo::rho_.Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::<anonymous>.Foam::DimensionedField<double, Foam::volMesh>::<anonymous>.Foam::Field<double>::<anonymous>.Foam::List<double>::<anonymous>.Foam::UList<T>::operator[]<double>(celli)’
         rho_[celli] = psi_[celli]*p_[celli];

and the error is pointing right at the assignment operator =. It seems like rho_[celli] can only be read but not be assigned?

This question can be generalized to how can one assign values to a volScalarField cell-by-cell instead of doing calculations on the volScalarField level? e.g., not p = x*y (where p, x, y are any volScalarField) but p[celli] = x[celli]*y[celli]?

I'm ruuning OF-4.x. on Ubuntu 16.04. And thank you again for reading.

Thanks

[Blumenkind]

The syntax seems to be right.

rho seems to be declared as a const. (e.g. "read-only")




For such projects and experiments the use of Openfoams test cases is recommendable.

[cryabroad]

Originally Posted by Blumenkind

The syntax seems to be right.

rho seems to be declared as a const. (e.g. "read-only")




For such projects and experiments the use of Openfoams test cases is recommendable.

Hi Blumenkind,

Thank you for the response! I think you are right about rho being declared as a const. My confusion still exists though: on the very bottom level of the code OF still has to do calculations cell-by-cell (this is my guess), which means OF itself has the ability to assign all the cell values. If this is the case, we as users should also be able to do similar operations right?

What "Openfoams test cases" are you referring to? Is it something I'm missing? I was planning to use the modified class to those tutorial cases of compressible flows.

Thanks.

[Blumenkind]

Test applications can be found in application/test in your openfoam installation.

To your problem: Your code does not lead to any further conclusions. rhoThermo has two functions rho(). One of those is constant. Be sure to not use this one. Furthermore I would suggest you write a new function for your calculation and base it one the function correctRho(volScalarField).

[cryabroad]

Originally Posted by Blumenkind

Test applications can be found in application/test in your openfoam installation.

To your problem: Your code does not lead to any further conclusions. rhoThermo has two functions rho(). One of those is constant. Be sure to not use this one. Furthermore I would suggest you write a new function for your calculation and base it one the function correctRho(volScalarField).

Thank you for the quick response, really appreciate it!

I think you are right, I'll create a new function instead of modifying the existing ones, that should overcome the restriction of rho being "read-only". Will work on that and see how it goes.

Thanks.

[cryabroad]

Hey Bluemenkind,

I hope you are doing well. Just want to share with you what I've been working on recently, it seems like I've found something interesting.

The code I'm working on now is not in file psiThermo.C, instead, it is in file hePsiThermo.C. In the description of its header file, this hePsiThermo class is used for calculating "Energy for a mixture based on compressibility".

In the hePsiThermo.C file, line 41 to 57, the code reads

forAll(TCells, celli)
{
const typename MixtureType::thermoType& mixture_ =
this->cellMixture(celli);

TCells[celli] = mixture_.THE
(
hCells[celli],
pCells[celli],
TCells[celli]
);

psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);

muCells[celli] = mixture_.mu(pCells[celli], TCells[celli]);
alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
}

What I did is to change the line associated with psiCells to something like

psiCells[celli] = C/pCells[celli];

where C is a constant. Ideally it should a function but for my case I used a constant for testing, and I randomly give C a value of 0.8. After this minor modification, I recompiled the psiThermo and hePsiThermo class, and used the rhoSimpleFoam solver (which uses these two classes) to solve the squareBend case in the tutorial. After several iterations the density converges to 0.8, which is the value I set it to.

Note that I didn't change the psiThermo class at all, in that class the density is still calculated as psi*p. However, because now psi is calculated as 0.8/p (from the modified hePsiThermo.C file), rho is (0.8/p)*p, which is 0.8.

I want to thank you for your previous discussions, and hope this little test is helpful to other interested readers (if there are any) as well.