CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   Modify turbulence kEpsilon for incompressible multiphase (

franciscofelis March 27, 2014 11:50

Modify turbulence kEpsilon for incompressible multiphase
Hello to everyone.

I'm trying to build a modify version of the incompressible kEpsilon turbulence introducing a variable density formulation for the transport equations to describe a two phase mean flow (air+water).

I'am using the twoLiquidMixingFoam as a base and I have created my personalized k-epsilon model from the incompressible part. My idea is to change the ddt(k) with ddt(rho, k) and add other source terms also as a function of "rho".

I'm rather new to this, but until now every problem that I have encountered I have managed to solve it by reading, reading and reading... until now:

I couldn't find a way to inherit the volScalarField rho and surfaceScalarField rhoPhi to the kEpsilon class.

Is there any efficient way to do that? Here is what I want:

To go from this:


      + fvm::div(phi_, k_)
      - fvm::laplacian(DkEff(), k_)
      - fvm::Sp(epsilon_/k_, k_)

To this:


        fvm::ddt(rho_, k_)
      + fvm::div(rhoPhi_, k_)
      - fvm::laplacian(rho_*DkEff(), k_)
      - fvm::Sp(rho*epsilon_/k_, k_)

And for that, to my understanding, rho and rhoPhi must be inherited from turbulencemodel->RASModel->mykEpsilon.

I would really appreciate your help!

- Francisco

kmooney March 31, 2014 12:02

Hi Francisco,

A few things:
1. When you mention inheriting a field, what you probably want to do is use a object registry lookup. Its an extremely powerful and easy thing that makes OpenFOAM super cool. If you wanted to pull in the density field, as long as its declared as a regIOobject, you can do this:


volScalarField& rho = mesh().lookupObject<volScalarField>("rho")
Basically, as long as you can get a reference to the mesh class, you can load up pretty much any field which is registered.

2. Compressible turbulence models already assume a non uniform density so that might be a better place to start then the incompressible turbulence classes.

3. This whole conversation might be moot because foam 2.3 has a multiphase kEp model :).


franciscofelis April 1, 2014 06:10

Hello Kyle,

Thanks for your answer. Effectively I looked also into the incompressible side of the Turbulence model constructor, however, It seemed more simple to add "rho" than modify the thermophysical part to work along the Multiphase transport one.

I've started also to look into the 2.3 version and the new universal turbulence models for multiphase, however I don't feel myself so "strong" in OpenFOAM for using it yet...

For now, I've tried the Mesh lookup method that you mentioned and it didn't work. Here is the output:


rhokEpsilon.C:144:33: error: ‘mesh’ was not declared in this scope
    volScalarField& rho_ = mesh().lookupObject<volScalarField>("rho");

To my understanding, the way that the turbulence model is constructed, he doesn't know that there is any mesh, a part from the constructors, any other variable that works beneath the "turbulence model" has to be inherited or read directly from an IO file, like k for instance:


        autoCreateK("k", mesh_)

From the constructor inherited from RASModel and Incompressible Turbulence Model these are the only available:


    const volVectorField& U, // hoping to add rho somewhere here
    const surfaceScalarField& phi,
    transportModel& transport,
    const word& turbulenceModelName,
    const word& modelName
    RASModel(modelName, U, phi, transport, turbulenceModelName),

So, to my understanding, it seems a little more complicated... Do you have any other inside on a way to call rho into the KEpsilon class ?

Many thanks in advance.


alexeym April 1, 2014 07:07


as kEpsilon is a child of RASModel class, which in turn is a child of trubulenceModel class it has mesh_ property. So your call should be


const volScalarField& rho = mesh_.lookupObject<volScalarField>("rho");
or, if we forget about inheritance, you can use


const volScalarField& rho = U_.mesh().lookupObject<volScalarField>("rho");
maybe in final implementation you'll move this call into the constructor of derived kEpsilon model.

franciscofelis April 1, 2014 10:36

Hello Alexey,

Thanks, it works like a charm!

So, finally, what I've done was to change the transport model for k and epsilon to a variable density formulation, which I think is more real when the difference between rho1 and rho2 is important. Here is the final code:

1) change in the production term of k: Do you think that this form of the variable density boussinesq is good? (rho is added at the end into the equations)


volScalarField G
        (((2.0/3.0)*I)*k_ - 2.0*nut_*(symm(fvc::grad(U_))-(1.0/3.0)*I*(fvc::div(U_))))
        && fvc::grad(U_)

2) Read rho and phiRho from the mesh:


    const volScalarField& rho_ = mesh_.lookupObject<volScalarField>("rho");
    const surfaceScalarField& rhoPhi_ = mesh_.lookupObject<surfaceScalarField>("rhoPhi");

3) Modify k and epsilon equations:


// Dissipation equation
    tmp<fvScalarMatrix> epsEqn
        fvm::ddt(rho_, epsilon_)
      + fvm::div(rhoPhi_, epsilon_)
      - fvm::laplacian(rho_*DepsilonEff(), epsilon_)
      - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)

// Turbulent kinetic energy equation
    tmp<fvScalarMatrix> kEqn
        fvm::ddt(rho_, k_)
      + fvm::div(rhoPhi_, k_)
      - fvm::laplacian(rho_*DkEff(), k_)
      - fvm::Sp(rho_*epsilon_/k_, k_)

All of this works inside a slight variation of the twoLiquidMixingFoam, which describes the dispersed phase (water droplets into air) into a single equation (alpha1). I hope that'll do the work...

Many thanks to everyone.

-- Francisco

All times are GMT -4. The time now is 19:32.