[HaZe]

Hi Foamers,

I'm trying to create new fields from the single terms of the k and Epsilon equations of the turbulence model for every time step of writeInterval.

Unfortunately this is my first attempt to do some OF progamming and C++ programming at all. I read a lot in the forum, user and programmer's guide and tried different things recommanded in [1]-[3].
But I did't make any progress until now.

To get an idea of what to do generally I'm trying to create a new field of the diffusion term of k

- fvm::laplacian(DkEff(), k_)

hoping this will work for all other terms like convection, time derivate and so on accordingly.

My first problem is where to implement the code. In the solvers createFields.H I tried to add

    volSymmTensorField kDiffusion
    (
        IOobject
        (
            "kDiffusion",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ
            IOobject::AUTO_WRITE
        ),
        - fvm::laplacian(DkEff(), k_)
    );

but there are DkEff and k_ not declared. So I tried to do it in kEpsilon.C in a way like

tmp<???> kEpsilon::kDiffusion() const
{
    return tmp<???>
    (
        new ???
        (
            IOobject
            (
                "kDiffusion",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
           - fvm::laplacian(DkEff(), k_)
        )
    );
}

Here are two problems: at first I don't understand of which class ??? is and second would this be the right way to do it anyway?

I already did read about swak4foam which could maybe help me but I'm bound to OF1.5-dev so I can't use it.

And just for my understanding: is it necessary to compile the solver after I made a change in a adjeced lib (libcompressibleRASModels in my special case)? Or do I just have to compile the lib and everything will be ok?

[1] http://openfoamwiki.net/index.php/In...IOobject_class
[2] http://openfoamwiki.net/index.php/Op...objectRegistry
[3] http://www.cfd-online.com/Forums/ope...ta-output.html

[mlaurita]

I think that since both k_ and DfEff() are volScalarFields, the laplacian should also return a volScalarField (i.e. this is what the ??? should be)

Second, using fvm::laplacian is an implicit operation that actually returns an fvScalarMatrix, not what you want here. If you use fvc::laplacian you perform an explicit calculation on the two existing volScalarFields, returning your kDiffusion also as a volScalarField.

-MJL

[HaZe]

Yeah thanks! Now it's at least compiling without errors and I did understand the difference of fvc/fvm a little bit more even if I had read the Programmer's guide several times :-)

Next problem: I expected a new file in my time directories called "kDiffusion". But there is none.

My code so far (green for my additions/changes):

kEpsilon.H

 //- Return the turbulence kinetic energy dissipation rate
        tmp<volScalarField> epsilon() const
        {
            return epsilon_;
        }
//k Epsilon Details by HaZe 25.9.2013

        //- Return the kinetic energy diffusion
       tmp<volScalarField> kDiffusion() const;
            
//HaZe end
        //- Return the Reynolds stress tensor
        tmp<volSymmTensorField> R() const;

kEpsilon.C

tmp<volSymmTensorField> kEpsilon::devRhoReff() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "devRhoReff",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
           -muEff()*dev(twoSymm(fvc::grad(U_)))
        )
    );
}
// HaZe 25.9.2013
tmp<volScalarField> kEpsilon::kDiffusion() const
{
    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "kDiffusion",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            -fvc::laplacian(DkEff(), k_)
        )
    );
}
//haze ende

tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(muEff(), U)
      - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
    );
}

After changing the code I compiled the lib and the solver. The solver runs but I'm missing the output of "kDiffusion". Is there anything else to change?

[mlaurita]

Well, you're part of the way there. You've successfully added a function to your turbulence model that calculates kDiffusion, but you haven't actually called it from anywhere.

For example, somewhere in your kEpsilon.C file (probably inside of the update() function), you could declare a variable, e.g.:

volScalarField kDiffusion
(
IOobject
(
"kDiffusion",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
kDiffusion()
);

Here, the value returned by your kDiffusion() function is assigned to a variable also called kDiffusion. If you do this after the k equation has been solved you'll get kDiffusion based on the latest values. Then you can write the field with:

kDiffusion.write();

There may also be a way to do it by adding kDiffusion_ as a member of the kEpsilon class and assigning it in the constructor. This could lead to it auto writing with the rest of the variables, but not having tried this myself I'm not sure of the details.

MJL

[HaZe]

Originally Posted by mlaurita

Well, you're part of the way there. You've successfully added a function to your turbulence model that calculates kDiffusion, but you haven't actually called it from anywhere.

This sounds logical to me.

I followed your instructions but it's still not working. There's no update() function as you mentioned so I'm not sure where to place it. I tried several locations inside kEpsilon.C but it was always wrong. This is my actual code for kEpsilon.C:

   
#include "kEpsilon.H"
#include "addToRunTimeSelectionTable.H"
#include "wallFvPatch.H"
            
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
            
namespace Foam
{
namespace compressible
{
namespace RASModels
{
            
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
               
defineTypeNameAndDebug(kEpsilon, 0);
addToRunTimeSelectionTable(RASModel, kEpsilon, dictionary);
                
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
     
kEpsilon::kEpsilon
(
    const volScalarField& rho,
    const volVectorField& U,
    const surfaceScalarField& phi,
    basicThermo& thermophysicalModel
)
:
    RASModel(typeName, rho, U, phi, thermophysicalModel),
      
    Cmu_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cmu",
            coeffDict_,
            0.09
        )
    ),
    C1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C1",
            coeffDict_,
            1.44
        )
    ),   
    C2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C2",
            coeffDict_,
            1.92
        )
    ),          
    C3_ 
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C3",
            coeffDict_,
            -0.33 
        )
    ),
    alphak_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphak",
            coeffDict_,
            1.0
        )
    ),
    alphaEps_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphaEps",
            coeffDict_,
            0.76923
        )
    ),
    alphah_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphah",
            coeffDict_,
            1.0
        )
    ),
                
    k_
    (
        IOobject
        (
            "k",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),
     
    epsilon_
    (
        IOobject
        (
            "epsilon",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),   
      
    mut_
    (
        IOobject
        (
            "mut",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
    )                

{
#   include "wallViscosityI.H"   

    volScalarField kDiffusion
    (
        IOobject
        (
            "kDiffusion",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        kDiffusion()
    );
    
    kDiffusion.write();

    printCoeffs();
}
    
     
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
      
tmp<volSymmTensorField> kEpsilon::R() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "R",
                runTime_.timeName(),  
                mesh_,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
            k_.boundaryField().types()
        )
    );  
}
        
tmp<volSymmTensorField> kEpsilon::devRhoReff() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "devRhoReff",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
           -muEff()*dev(twoSymm(fvc::grad(U_)))
        )
    );
}

//haze
tmp<volScalarField> kEpsilon::kDiffusion() const
{
    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "kDiffusion",
                runTime_.timeName(),
                mesh_,        
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            -fvc::laplacian(DkEff(), k_)
        )
    );
}
//haze ende

tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(muEff(), U)
      - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
    );
}
        
bool kEpsilon::read()
{
    if (RASModel::read())
    {
        Cmu_.readIfPresent(coeffDict_);   
        C1_.readIfPresent(coeffDict_);
        C2_.readIfPresent(coeffDict_);
        C3_.readIfPresent(coeffDict_);
        alphak_.readIfPresent(coeffDict_);
        alphaEps_.readIfPresent(coeffDict_);
        alphah_.readIfPresent(coeffDict_);

        return true;
    }
    else
    {
        return false;
    }
}


void kEpsilon::correct()
{
    // Bound in case of topological change.  Signalling currently not available
    // HJ, 22/Aug/2007
    bound(k_, k0_);
    bound(epsilon_, epsilon0_);

    if (!turbulence_)
    {
        // Re-calculate viscosity
        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
#       include "wallViscosityI.H"
        return;
    }

    RASModel::correct();

    volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));

    if (mesh_.moving())
    {
        divU += fvc::div(mesh_.phi());
    }

    tmp<volTensorField> tgradU = fvc::grad(U_);
    volScalarField G = mut_*(tgradU() && dev(twoSymm(tgradU())));
    tgradU.clear();

#   include "wallFunctionsI.H"

    // Dissipation equation
    tmp<fvScalarMatrix> epsEqn
    (
        fvm::ddt(rho_, epsilon_)
      + fvm::div(phi_, epsilon_)
      - fvm::laplacian(DepsilonEff(), epsilon_)
     ==
        C1_*G*epsilon_/k_
      - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
      - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
    );

#   include "wallDissipationI.H"

    epsEqn().relax();

    solve(epsEqn);
    bound(epsilon_, epsilon0_);


    // Turbulent kinetic energy equation

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(rho_, k_)
      + fvm::div(phi_, k_)
      - fvm::laplacian(DkEff(), k_)
     ==
        G
      - fvm::SuSp((2.0/3.0)*rho_*divU, k_)
      - fvm::Sp(rho_*epsilon_/k_, k_)
    );

    kEqn().relax();
    solve(kEqn);
    bound(k_, k0_);


    // Re-calculate viscosity
    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;

#   include "wallViscosityI.H"

}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam

I get the compiler error:

kEpsilon/kEpsilon.C: In constructor ‘Foam::compressible::RASModels::kEpsilon::kEpsilon(const Foam::volScalarField&, const Foam::volVectorField&, const Foam::surfaceScalarField&, Foam::basicThermo&)’:
kEpsilon/kEpsilon.C:172: error: no match for call to ‘(Foam::volScalarField) ()’

Among a lot of other variants where to place the code I also tried according to this post
and the following post.

  
#include "kEpsilon.H"
#include "addToRunTimeSelectionTable.H"
#include "wallFvPatch.H"
            
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
            
namespace Foam
{
namespace compressible
{
namespace RASModels
{
            
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
               
defineTypeNameAndDebug(kEpsilon, 0);
addToRunTimeSelectionTable(RASModel, kEpsilon, dictionary);
                
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
     
kEpsilon::kEpsilon
(
    const volScalarField& rho,
    const volVectorField& U,
    const surfaceScalarField& phi,
    basicThermo& thermophysicalModel
)
:
    RASModel(typeName, rho, U, phi, thermophysicalModel),
      
    Cmu_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cmu",
            coeffDict_,
            0.09
        )
    ),
    C1_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C1",
            coeffDict_,
            1.44
        )
    ),   
    C2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C2",
            coeffDict_,
            1.92
        )
    ),          
    C3_ 
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C3",
            coeffDict_,
            -0.33 
        )
    ),
    alphak_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphak",
            coeffDict_,
            1.0
        )
    ),
    alphaEps_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphaEps",
            coeffDict_,
            0.76923
        )
    ),
    alphah_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphah",
            coeffDict_,
            1.0
        )
    ),
                
    k_
    (
        IOobject
        (
            "k",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),
     
    epsilon_
    (
        IOobject
        (
            "epsilon",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),   
      
    mut_
    (
        IOobject
        (
            "mut",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
    ),
                
    kDiffusion_ 
    (
        IOobject
        (
            "kDiffusion",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        kDiffusion()
    )

{
#   include "wallViscosityI.H"   

if(runTime_.outputTime())
    {
        kDiffusion_.write();
    }
    printCoeffs();
}
    
     
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
      
tmp<volSymmTensorField> kEpsilon::R() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "R",
                runTime_.timeName(),  
                mesh_,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
            k_.boundaryField().types()
        )
    );  
}
        
tmp<volSymmTensorField> kEpsilon::devRhoReff() const
{
    return tmp<volSymmTensorField>
    (
        new volSymmTensorField
        (
            IOobject
            (
                "devRhoReff",
                runTime_.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
           -muEff()*dev(twoSymm(fvc::grad(U_)))
        )
    );
}

//haze
tmp<volScalarField> kEpsilon::kDiffusion() const
{
    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "kDiffusion",
                runTime_.timeName(),
                mesh_,        
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            -fvc::laplacian(DkEff(), k_)
        )
    );
}
//haze ende

tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(muEff(), U)
      - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
    );
}
        
bool kEpsilon::read()
{
    if (RASModel::read())
    {
        Cmu_.readIfPresent(coeffDict_);   
        C1_.readIfPresent(coeffDict_);
        C2_.readIfPresent(coeffDict_);
        C3_.readIfPresent(coeffDict_);
        alphak_.readIfPresent(coeffDict_);
        alphaEps_.readIfPresent(coeffDict_);
        alphah_.readIfPresent(coeffDict_);

        return true;
    }
    else
    {
        return false;
    }
}


void kEpsilon::correct()
{
    // Bound in case of topological change.  Signalling currently not available
    // HJ, 22/Aug/2007
    bound(k_, k0_);
    bound(epsilon_, epsilon0_);

    if (!turbulence_)
    {
        // Re-calculate viscosity
        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
#       include "wallViscosityI.H"
        return;
    }

    RASModel::correct();

    volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));

    if (mesh_.moving())
    {
        divU += fvc::div(mesh_.phi());
    }

    tmp<volTensorField> tgradU = fvc::grad(U_);
    volScalarField G = mut_*(tgradU() && dev(twoSymm(tgradU())));
    tgradU.clear();

#   include "wallFunctionsI.H"

    // Dissipation equation
    tmp<fvScalarMatrix> epsEqn
    (
        fvm::ddt(rho_, epsilon_)
      + fvm::div(phi_, epsilon_)
      - fvm::laplacian(DepsilonEff(), epsilon_)
     ==
        C1_*G*epsilon_/k_
      - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
      - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
    );

#   include "wallDissipationI.H"

    epsEqn().relax();

    solve(epsEqn);
    bound(epsilon_, epsilon0_);


    // Turbulent kinetic energy equation

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(rho_, k_)
      + fvm::div(phi_, k_)
      - fvm::laplacian(DkEff(), k_)
     ==
        G
      - fvm::SuSp((2.0/3.0)*rho_*divU, k_)
      - fvm::Sp(rho_*epsilon_/k_, k_)
    );

    kEqn().relax();
    solve(kEqn);
    bound(k_, k0_);


    // Re-calculate viscosity
    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;

#   include "wallViscosityI.H"

}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam

// ************************************************************************* //

with kEpsilon.H:

\*---------------------------------------------------------------------------*/
         
#ifndef compressiblekEpsilon_H
#define compressiblekEpsilon_H

#include "RASModel.H"
        
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
            
namespace Foam
{
namespace compressible
{
namespace RASModels
{
        
/*---------------------------------------------------------------------------*\
                           Class kEpsilon Declaration
\*---------------------------------------------------------------------------*/
                
class kEpsilon
:
    public RASModel
{
    // Private data
         
        // Model coefficients
             
//            dimensionedScalar Cmu;
            dimensionedScalar Cmu_;
            dimensionedScalar C1_;
            dimensionedScalar C2_;
            dimensionedScalar C3_;
            dimensionedScalar alphak_;
            dimensionedScalar alphaEps_;
            dimensionedScalar alphah_;

        // Fields

            volScalarField k_;
            volScalarField epsilon_;
            volScalarField mut_;
//Haze add
            volScalarField kDiffusion_;
//End Haze

public:

    //- Runtime type information
    TypeName("kEpsilon");

    // Constructors

        //- from components
        kEpsilon
        (
            const volScalarField& rho,
            const volVectorField& U,
            const surfaceScalarField& phi,
            basicThermo& thermophysicalModel
        );


    // Destructor

        ~kEpsilon()
        {}


    // Member Functions

      //- Return the turbulence viscosity
        tmp<volScalarField> mut() const
        {
            return mut_;
        }

        //- Return the effective diffusivity for k
        tmp<volScalarField> DkEff() const
        {
            return tmp<volScalarField>
            (
                new volScalarField("DkEff", alphak_*mut_ + mu())
            );
        }   

        //- Return the effective diffusivity for epsilon
        tmp<volScalarField> DepsilonEff() const
        {
            return tmp<volScalarField>
            (
                new volScalarField("DepsilonEff", alphaEps_*mut_ + mu())
            );
        }

        //- Return the effective turbulent thermal diffusivity
        tmp<volScalarField> alphaEff() const
        {
            return tmp<volScalarField>
            (
                new volScalarField("alphaEff", alphah_*mut_ + alpha())
            );
        }
             
        //- Return the turbulence kinetic energy
        tmp<volScalarField> k() const
        {
            return k_;
        }
            
        //- Return the turbulence kinetic energy dissipation rate
        tmp<volScalarField> epsilon() const
        {
            return epsilon_;
        }
//k Epsilon Details by HaZE 23.9.2013
            
        //- Return the kinetic energy diffusion
        tmp<volScalarField> kDiffusion() const;
            
//Details end
        //- Return the Reynolds stress tensor
        tmp<volSymmTensorField> R() const;

        //- Return the effective stress tensor including the laminar stress
        tmp<volSymmTensorField> devRhoReff() const;

        //- Return the source term for the momentum equation
        tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
        
        //- Solve the turbulence equations and correct the turbulence viscosity
        void correct();
            
        //- Read turbulenceProperties dictionary
        bool read();
};
          

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam

Produces:

kEpsilon/kEpsilon.C: In constructor ‘Foam::compressible::RASModels::kEpsilon::kEpsilon(const Foam::volScalarField&, const Foam::volVectorField&, const Foam::surfaceScalarField&, Foam::basicThermo&)’:
kEpsilon/kEpsilon.C:176: error: ‘((Foam::compressible::RASModels::kEpsilon*)this)->Foam::compressible::RASModels::kEpsilon::kDiffusion’ does not have class type

Now I'm really stuck and don't know what else to try.

Could you please help me another time. I also don't understand why I have to create a second IOobject for the variable "kDiffusion". Isn't this already done in my member function?

[mlaurita]

My mistake, I meant correct() (not update()).

But after thinking about this a little more I think you're almost there...
in the assignment section of the constructor try:

mut_ ( IOobject ( "mut", runTime_.timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_) ), kDiffusion_ ( IOobject ( "kDiffusion", runTime_.timeName(), mesh_, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh_ )
{
kDiffusion_ = kDiffusion();
kDiffusion_.updateBoundaryConditions();
}

Then, your kDiffusion() function can simply be:
tmp<volScalarField> kEpsilon::kDiffusion() const { return (
-fvc::laplacian(DkEff(), k_) );
}

And finally, at the end of the correct() function you can add:
mut_ = rho_*Cmu_*sqr(k_)/epsilon_; # include "wallViscosityI.H"

kDiffusion_ = kDiffusion();
kDiffusion_.updateBoundaryConditions();
}


This creates kDiffusion_ as a variable that is updated with the latest value of k_ and then automatically written to the disk whenever the solution is written. I'm not 100% sure you need the call to updateBoundaryConditions(), but I was basing this off of the way that nut is handled in the incompressible version. One catch is that you need to have a "kDiffusion" file in the 0/ directory, but you can specify an internal field of 0 and boundary fields as "calculated".

I tried this with the incompressible version of kEpsilon and it worked fine, though I don't know if the resulting "kDiffusion" values that are output are accurate, but that'll be up to you.

One last thing, it's always good practice to copy the turbulence model to a new folder and give the files a new name ("myKEpsilon.H", "myKEpsilon.C", for example). Also be sure to change the typename in both the .H and .C files. That way you still have the original one available and you can select your model by specifying "RASModel myKEpsilon;" in the RASProperties dictionary.

Good luck!

MJL

[mlaurita]

Sorry, let me try again with the code tags this time!

    mut_
    (
        IOobject
        (
            "mut",              
            runTime_.timeName(),             
            mesh_,             
            IOobject::NO_READ,              
            IOobject::NO_WRITE         
        ), 
        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)     
    ),                     
    kDiffusion_     
    (         
        IOobject
        (              
            "kDiffusion", 
            runTime_.timeName(),
            mesh_,              
            IOobject::MUST_READ,             
            IOobject::AUTO_WRITE         
        ), 
        mesh_     
    )
{
         kDiffusion_ = kDiffusion();
         kDiffusion_.updateBoundaryConditions();
}

tmp<volScalarField> kEpsilon::kDiffusion() const 
{     
    return     
    (
                -fvc::laplacian(DkEff(), k_)     
    ); 
}

mut_ = rho_*Cmu_*sqr(k_)/epsilon_;  

#   include "wallViscosityI.H"

    kDiffusion_ = kDiffusion();
    kDiffusion_.updateBoundaryConditions();

[hiuluom]

Hi all.

I have just used openfoam. I has LES model and it runs good. And now, I would like to calculate some terms in LES as first term1 = nu{p*[d2(ui/xi)]} and second term2 = p*[d((uiuj)/xj)]. After that, get average term1 - average term2.
d2 is second derivative and d is first derivative.

Anybody could help me how to setup all math in openfoam?

Thanks.

Thanh.

[HaZe]

@hiuluom: Would you please start a new thread for your question. Otherwise it's getting complicated with more than one discussion per thread.

@mlaurita: Thank you so much! Now it's working!

At first the compiler complained about "volScalarField" has no member kDiffusion_.updateBoundaryConditions(). But I changed kDiffusion_.updateBoundaryConditions() to kDiffusion_.write and now it's working.

I will do some testing and then post the complete working code.

Thanks for everything and keep FOAMing.

[HaZe]

Thank you mlaurita! Here's the extended kEpsilon Model so far. It's working with OF-1.5-dev.

[abbott.hn]

Originally Posted by mlaurita

kDiffusion_
(
IOobject
(
"kDiffusion",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
)
{
kDiffusion_ = kDiffusion();
kDiffusion_.updateBoundaryConditions();
}

when define the member kDiffusion_ like this ,how should I specify the BC for it? thanks for your reply

[HaZe]

Sorry for not mentioning you have to make a new file in your time directory for every new k** and epsilon** field.
The boundary condition is

type calculated;

for every patch. The dimension for k** is

dimensions      [1 -1 -3 0 0 0 0];

and for epsilon**

dimensions      [1 -1 -4 0 0 0 0];

And at last I have to say I didn't take the "updateCoefficients" line into my code. It didn't compile with it, so I don't know yet if it has any effects on the results.

It would be great if anyone using this code would share his/her experiences with it.

[aki_yafuji]

Hi Haze, and everyone!

Thank you for sharing the extended kEpsilon code! It's working with OF-2.2 by modifying a little, and I can output new fields such as kProduction!
I would like to use these new fields in the interFoam solver, how should I do??

In interFoam.C, I add
/////////
volScalarField A1 = turbulence->nut(); // OK
volScalarField A2 = turbulence->k(); // OK
volSymmTensorField A3 = turbulence->R(); // OK
volScalarField A4 = turbulence->kProduction(); // did not work ...
/////////

but, I recieved the following error message. ..

interFoam.C:106:40: error: ‘class Foam::incompressible::turbulenceModel’ has no member named ‘kProduction’

If you wouldn't mind, would you please give me some advice?

aki

[HaZe]

That should be easy. My code works for the compressible Turbulence Models. If you want to use it with incompressible solvers you have to adjust it for the incompressible k and epsilon function and add it to the incompressible Turbulence Models, too.

Meanwhile I figured out the dimensions of the k** and Epsilon** terms are scaled by density (look at the dimensions). I don't know why and where this happens. I already completed my studies and it was possible for me to fix this bug post processing.

If you want to use the fields directly in the solver, don't call the new functions from the solver. I tried this and didn't make it this way. For me it was more elegant to get the values from the new calculated fields.
This works something like:

volScalarField& kProductionPtr = const_cast<volScalarField&>(mesh.lookupObject<volScalarField>("kProduction"));

[Nitin Minocha]

Hello
First of all, i would like to thank you for this post. I tried to use the modified kEpsilonDetailed model, in order to solve additional terms of turbulent kinetic energy and dissipation rate. I am using OpenFOAM 2.2.0. I am using buoyantBoussinesqPisoFoam solver for incompressible fluid. I got the following error message while compiling the code.
KEmake: Warning: File `options' has modification time 1.8e+03 s in the future
make: warning: Clock skew detected. Your build may be incomplete.
make: Warning: File `linux64Gcc46DPOpt/options' has modification time 15 s in the future
make: warning: Clock skew detected. Your build may be incomplete.
make: Warning: File `Make/linux64Gcc46DPOpt/dontIncludeDeps' has modification time 15 s in the future
make: *** No rule to make target `kEpsilonDetailed/kEpsilonDetailed.dep', needed by `Make/linux64Gcc46DPOpt/dependencies'. Stop.
I am sending you my .c and .h files alongwith make directory. Please have a look and help me solve this problem.

[HaZe]

Originally Posted by Nitin Minocha

Hello
First of all, i would like to thank you for this post. I tried to use the modified kEpsilonDetailed model, in order to solve additional terms of turbulent kinetic energy and dissipation rate. I am using OpenFOAM 2.2.0. I am using buoyantBoussinesqPisoFoam solver for incompressible fluid. I got the following error message while compiling the code.
KEmake: Warning: File `options' has modification time 1.8e+03 s in the future
make: warning: Clock skew detected. Your build may be incomplete.
make: Warning: File `linux64Gcc46DPOpt/options' has modification time 15 s in the future
make: warning: Clock skew detected. Your build may be incomplete.
make: Warning: File `Make/linux64Gcc46DPOpt/dontIncludeDeps' has modification time 15 s in the future
make: *** No rule to make target `kEpsilonDetailed/kEpsilonDetailed.dep', needed by `Make/linux64Gcc46DPOpt/dependencies'. Stop.
I am sending you my .c and .h files alongwith make directory. Please have a look and help me solve this problem.

Hi there,
it's great to hear somebody is using my work. Since I didn't work with OF for the last months I'm not sure if I can provide any helpful hints.

But there are some points you can check before further problems occur:
- Did you recognise my model is for compressible fluids?
- Did you recognise my model is made for OF-1.5-dev (I don't know what changes are made between OF-1.5-dev and your OF-version but I'm sure it's a lot!)
- Were you able to compile my model without problems and use it?
- Did you read my last post? In my special case I was able to correct my calculations post processing - but that's no nice way to do so.
- At last I didn't understand why you are trying to extend my model and not extending the existing one from OF-2.2.0. What I did is only to write some new fields to visualize all terms of the k-Epsilon-Model - nothing more.

And last of all would you please start a new thread with your questions.

Regards,
HaZe