build your own turbulence model with buoyancy
 

Hallo all,

I try to make my own turbulence model with buoyancy effects. I began making my own kEpsilonModel with a dynamic library and included this one in the controlDict-data with the line libs ("libmyBCs.so");. I copied the kEpsilonModel, renamed it and compiled it. It works. Now I want to modify this model.

For my case it's neccessary to make a new class of RASModel with the Temperaturefield:

RASModel::RASModel
(
const word& type,
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& T,
transportModel& lamTransportModel
)

So I modified the RASModel.C & RasModel.H, the transportModel.C&H and the turbulenceModel.C&H datas. I used for all of them namespace incompressible.

Compiling all of them works without problems. But using this turbulenceModel with I case I get following error message:

/ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g
Reading thermophysical properties

Reading field T

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Creating turbulence model

Selecting RAS turbulence model kEpsilonTom



lookup of RASProperties from objectRegistry region0 successful
but it is not a RASModel, it is a kEpsilonTom#0 Foam::error::printStack(Foam::Ostream&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#3 Foam::incompressible::RASModel const& Foam::objectRegistry::lookupObject<Foam::incompres sible::RASModel>(Foam::word const&) const in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#4 Foam::incompressible::RASModels::nutWallFunctionFv PatchScalarField::calcNut() const in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#5 Foam::incompressible::RASModels::nutWallFunctionFv PatchScalarField::updateCoeffs() in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#6 Foam::fvPatchField<double>::evaluate(Foam::Pstream ::commsTypes) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#7 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::evaluate() in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#8 Foam::incompressible::RASModels::kEpsilonTom::kEps ilonTom(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::transportModel&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libmyincompressibleRASModels.so"
#9 Foam::incompressible::RASModel::adddictionaryConst ructorToTable<Foam::incompressible::RASModels::kEp silonTom>::New(Foam::GeometricField<Foam::Vector<d ouble>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::transportModel&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libmyincompressibleRASModels.so"
#10 Foam::incompressible::RASModel::New(Foam::Geometri cField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#11 main in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#12 __libc_start_main in "/lib/libc.so.6"
#13 _start at /build/buildd/glibc-2.9/csu/../sysdeps/x86_64/elf/start.S:116


From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/dm2/henry/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 121.

FOAM aborting

So my problem is, it is a RASModel, made with the new constructor with more arguments (temperature-field T is included).
So is it neccessary to create my own namespace? Has anybody an idea?

Regards Thomas

[QUOTE=Thomas Baumann;235374]Hallo all,

I try to make my own turbulence model with buoyancy effects. I began making my own kEpsilonModel with a dynamic library and included this one in the controlDict-data with the line libs ("libmyBCs.so");. I copied the kEpsilonModel, renamed it and compiled it. It works. Now I want to modify this model.

For my case it's neccessary to make a new class of RASModel with the Temperaturefield:

RASModel::RASModel
(
const word& type,
const volVectorField& U,
const surfaceScalarField& phi,
const volScalarField& T,
transportModel& lamTransportModel
)

So I modified the RASModel.C & RasModel.H, the transportModel.C&H and the turbulenceModel.C&H datas. I used for all of them namespace incompressible.

Compiling all of them works without problems. But using this turbulenceModel with I case I get following error message:

/ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g
Reading thermophysical properties

Reading field T

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Creating turbulence model

Selecting RAS turbulence model kEpsilonTom



lookup of RASProperties from objectRegistry region0 successful
but it is not a RASModel, it is a kEpsilonTom#0 Foam::error::printStack(Foam::Ostream&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#3 Foam::incompressible::RASModel const& Foam::objectRegistry::lookupObject<Foam::incompres sible::RASModel>(Foam::word const&) const in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#4 Foam::incompressible::RASModels::nutWallFunctionFv PatchScalarField::calcNut() const in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#5 Foam::incompressible::RASModels::nutWallFunctionFv PatchScalarField::updateCoeffs() in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#6 Foam::fvPatchField<double>::evaluate(Foam::Pstream ::commsTypes) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#7 Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::evaluate() in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#8 Foam::incompressible::RASModels::kEpsilonTom::kEps ilonTom(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::transportModel&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libmyincompressibleRASModels.so"
#9 Foam::incompressible::RASModel::adddictionaryConst ructorToTable<Foam::incompressible::RASModels::kEp silonTom>::New(Foam::GeometricField<Foam::Vector<d ouble>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::transportModel&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libmyincompressibleRASModels.so"
#10 Foam::incompressible::RASModel::New(Foam::Geometri cField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&) in "/user/hi204/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libincompressibleRASModels.so"
#11 main in "/user/hi204/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/buoyantSimpleOne"
#12 __libc_start_main in "/lib/libc.so.6"
#13 _start at /build/buildd/glibc-2.9/csu/../sysdeps/x86_64/elf/start.S:116


From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/dm2/henry/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 121.

FOAM aborting

So my problem is, it is a RASModel, made with the new constructor with more arguments (temperature-field T is included). So is it neccessary to create my own namespace? Has anybody an idea?

Regards Thomas

If you work through how the dispatch works, you'll understand why your code must segfault/abort on you.
For example, in kEpsilon.C you'll see it is registered via a dictionary constructor:
Which you can follow through to RASModel.C so see where the dispatch actually occurs in RASModel::New()

You are thus currently stuck with an identical number of parameters if you wish to hook into the current RASModel selection mechanism.

A workaround would be to retain the normal parameter list (like the other RASModels) and assign the temperature reference within the constructor using the objectRegistry lookup mechanism. You can use "T" as the default name for the temperature field, but also provision for an alternative name by using a 'readIfPresent' on the corresponding model coefficients sub-dictionary.

Hallo Mark,
I had a closer look at the different solvers (for example buoyantBoussinesqPisoFoam). There is the RASModel.H included. So I think making a dynamic bibliothek with new MyRASModel.H&C with a new constructor isn't the best solution.

So I'm trying your proposal with the old RASModel.H&C and include the reference for the temperaturefield in the new myTurbulenceModel.C&H like you discribed.

Regards Thomas

Finally I modified the kEpsilon model without modifying the class RASModel.

My modifications being able to include the temperature field in the turbulence modelling transport equations are:

kEpsilonBuoyancy.H (former kEpsilon.H)

volScalarField Tnew_;

virtual tmp<volScalarField> T() const
{
return Tnew_;
}

kEpsilonBuoyancy.C (former kEpsilon.C)

Tnew_
(
IOobject
(
"T",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh_
),

I modified the transport equations, implemented the turbulent Prandtl number, thermal expansion coefficient, some constants, too.
The first simulations show nice results.

Regards Thomas

Are you really sure that you want a tmp field for Temperature? Also, how will your Tnew_ have the current value of "T" at each time during the calculation?
My first impression is that it may somehow work (ie, not crash), but it isn't quite right either.

Hallo Mark,

thank you for the information.
I'm trying to implement the temperature field using the lookup-method.

But my experience with C++ and programming in OpenFOAM is not very high. I have looked in some codes and tried to implement it in my model, but nothing has matched until now.

I think the code must look like:
const volScalarField& T = ???.lookupObject<volScalarField>("T");

So I'm further searching for "how to include T" without writing a new constructor of the RASModel-class.

Regards Thomas

This looks about right. The ??? just needs to be replaced by an objectRegistry. Why not try using the one attached to another of your fields?
For example,
const volScalarField& T = U.db().lookupObject<volScalarField>("T");

/mark

Implementing

const volScalarField& T = U.db().lookupObject<volScalarField>("T");

and compiling the new turbulence-model I get following error message:

clude -I. -I/user/hi204/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude -I/user/hi204/OpenFOAM/OpenFOAM-1.6/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/kEpsilonBuoyancy.o
kEpsilonBuoyancy/kEpsilonBuoyancy.C:157: error: ‘U’ was not declared in this scope
kEpsilonBuoyancy/kEpsilonBuoyancy.C:157: error: expected primary-expression before ‘>’ token
make: *** [Make/linux64GccDPOpt/kEpsilonBuoyancy.o] Fehler 1


Regards Thomas

Okay, but I did say "For example" - I don't know which fields you have available in your class nor what names you've given them. Maybe you've got them named U_ or tkeValues or whatever. I don't even know what your class looks like or anything, how can I know what you should type to get at your fields?!?

Hallo Mark,

I think it works now. Your tips helped me very much.
By the way, you were right: the velocity field was U_, not U. It took me some time to understand your tips and I had to understand more about C++.
For the sake of completeness: I want to use the new turbulence model in simulations using the BuoyantSimpleFoam solver.


I implemented a pointer to get the temperature field without modifying the constructor of RASModel.


myTurbulenceModel.H
const volScalarField* Tnew_;

myTurbulenceModel.C within the constructor:
Tnew_ = &(U_.db().lookupObject<volScalarField>("T"));

Because of the overloaded Operator * in the turbulent kinetic energy and dissipation equation I created a function Tnew() to use the temperature field T for modelling the turbulent viscoity. So I implemented following member function, too:

myTurbulenceModel.H
virtual const volScalarField Tnew();

myTurbulenceModel.C
const volScalarField myTurbulenceModel::Tnew()
{
return *Tnew_;
}

It's not easy to test, if the right temperature field is used.
So I did the same method to get the velocity field U, modified the kEpsilon turbulence model by solving the k and epsilon equation using this new method to get the velocity field (the same as the T field). I simulated a basic problem with the basic kEpsilon and with the modified kEpsilon model.
Comparing the residuals showed no differences. So I think the temperature-field should be included in a right way.

So thank you very much for your help.

Regards Thomas

Servus Thomas,

It looks mostly okay, except that in OpenFOAM they mostly use references instead of pointers and it look more consistent if you used references too, but that is, I suppose, a matter of taste for your case.

Assuming that you stick with pointers, you should likely be writing it like this:

The return type is 'const volScalarField&' since you are returning a const (read-only) reference to the field. If you leave out the '&' on the return type, you are copying back the field, which isn't so great for performance.

You should also have a trailing 'const' on the Tnew() method, since calling it does not change anything inside of your model - ie, you can call the Tnew() method on a read-only object without any problem.

/mark