[lachuktr]

Currently, I am simulating the atmospheric boundary layer in OpenFOAM. I want to add the source term to the epsilon equation of the k-epsilon turbulence model.

Since z is not defined in kEpsilon.C, I am having trouble initiating the "z" value

"z" is the distance from the ground.
Kindly help.

[crubio.abujas]

Hi Lakshman,

I think you can try scalarCodedSource attached to epsilon equation. Depending on the OF version you're using the epsilon equation may or may not ask fvOption for additional sources. Ive tried this on OF7 and seems to work fine.

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

epsilonSource
{
    type    scalarCodedSource;
    name    internalField;
    selectionMode    all;
    fields    (epsilon);

    codeAddSup
    #{
        const vectorField& C = mesh().C();
        const scalarField& V = mesh_.V();

        const volScalarField& k = mesh_.lookupObject<volScalarField>("k");

        scalar C_mu(0.09);      // Coefficient of turbulent viscosity
        scalar z0(0.03);        // Aerodynamics roughness

        scalarField& source = eqn.source();

        forAll(C, i)
        {
           source[i] -= C_mu * pow(k[i], 2) / pow( C[i].z() + z0 , 2 ) * V[i];
        }

        Info<< "epsilonSource min/max: " 
            << gMin(source) << "/" << gMax(source)
            << endl;
  #};
 
}

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

I'm not totally sure if you have to multiply it by the volume or not, you should check if it work as expected.

I hope that helps.

[lachuktr]

Originally Posted by crubio.abujas

Hi Lakshman,

I think you can try scalarCodedSource attached to epsilon equation. Depending on the OF version you're using the epsilon equation may or may not ask fvOption for additional sources. Ive tried this on OF7 and seems to work fine.

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

epsilonSource
{
    type    scalarCodedSource;
    name    internalField;
    selectionMode    all;
    fields    (epsilon);

    codeAddSup
    #{
        const vectorField& C = mesh().C();
        const scalarField& V = mesh_.V();

        const volScalarField& k = mesh_.lookupObject<volScalarField>("k");

        scalar C_mu(0.09);      // Coefficient of turbulent viscosity
        scalar z0(0.03);        // Aerodynamics roughness

        scalarField& source = eqn.source();

        forAll(C, i)
        {
           source[i] -= C_mu * pow(k[i], 2) / pow( C[i].z() + z0 , 2 ) * V[i];
        }

        Info<< "epsilonSource min/max: " 
            << gMin(source) << "/" << gMax(source)
            << endl;
  #};
 
}

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

I'm not totally sure if you have to multiply it by the volume or not, you should check if it work as expected.

I hope that helps.

Thank you so much Carlos... it seems to be working fine.

Since I have not much experience in C++ coding and FvOption, can you please explain the codes you wrote here for better understanding. it will be very helpful for me in the future.

[crubio.abujas]

hi Lakshman,


I've added a commented version of the code. The object mesh_ has a lot of information about the simulation, as it acts as register you can recall different fields to be used in your code (for example the k field). Check out the fvMesh class for more information.


Once you get the information required (center of the cells, volumes, other fields and constants) you loop across all the cells in your domain to apply the given source term.


The system of equation is automatically read by the codedSource. Be aware that this eqn is a fvMatrix so you can see further details in the documentation. In short, it stores the coefficients of the matrix A and the terms of the source b of the system Ax=b. I think that it stores the source as Ax - b = 0, so you need to switch the sign on the source term (thats the by of the -= in the source).

        const vectorField& C = mesh().C();   // Read from the mesh the center of each cell.

        const scalarField& V = mesh_.V();  // Read from the mesh the volume of each cell.

        const volScalarField& k = mesh_.lookupObject<volScalarField>("k"); // lookup for the kinetic turbulence energy field.

        // Set the constant refereed inthe formula.

        scalar C_mu(0.09);      // Coefficient of turbulent viscosity
        scalar z0(0.03);        // Aerodynamics roughness

        // the codedOption will reach the system of ecuations associated 
        // with the field (Ax=b), the idea is to add an explicit source to it 
        // so the source is extracted (b).
        scalarField& source = eqn.source();

        // loop for each of the cells.
        forAll(C, i)
        {
           // For each cell add the required source term.
           // pow (x, 2) will raise x^2
           source[i] -= C_mu * pow(k[i], 2) / pow( C[i].z() + z0 , 2 ) * V[i];
        }

        // This just prints the max/min source terms for debugging purposes.
        Info<< "epsilonSource min/max: " 
            << gMin(source) << "/" << gMax(source)
            << endl;

[lachuktr]

Thanks a lot for your help...

Can you help me to add the same source term in KEpsilon.C file rather than using fvOption.

I have seen the YAP correction source term (attached file along with this reply )added to KEpsilon.C as follows.


template<class BasicTurbulenceModel>
tmp<volScalarField> YkEpsilon<BasicTurbulenceModel>::sYap() const
{

const volScalarField ell_eps = pow(Cmu_,-0.75)*kappa_*y_;
const volScalarField ell = pow(k_,1.5)/epsilon_;

//create a volScalrField=0 with the same dimensions of sYap
//this will be used to enforce sYap=max(sYap,0)
volScalarField zero_field
(
IOobject
(
"zero_field",
this->runTime_.timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
this->mesh_,
dimensionedScalar("0", dimensionSet(0,2,-4,0,0,0,0), 0.0)
);

return
max(Cyap_*pow(epsilon_,2)/k_*(ell/ell_eps-1)*pow((ell/ell_eps),2),zero_field);
}




Can you help me in adding the souce term which i attached in the first message of this forum simillar to the one given above.


I am also attaching the .C and .H files of yap correction implemented kEpsilon model for reference.

[crubio.abujas]

In the turbulence model the equation is solved when you call the correct() method. There you can see where the equations are implemented:

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

Here you can define an additional source term or modify a source already existing. I've tried to modify the epsilonSource as follows:

template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> YkEpsilon<BasicTurbulenceModel>::epsilonSource() const
{
    tmp<fvScalarMatrix> epsSource
    (
        new fvScalarMatrix
        (
            epsilon_,
            dimVolume*this->rho_.dimensions()*epsilon_.dimensions()
            /dimTime
        )
    );

    const vectorField& C = this->mesh().C();
    const scalarField& V = this->mesh().V();

    scalar C_mu(0.09);
    scalar z0(0.03);

    scalarField& source = epsSource.ref().source();

    forAll(C, i)
    {
        source[i] -=
            C_mu * pow(k_[i] , 2) / 
            pow( C[i].z() + z0 , 2 ) * V[i];
    }

    return epsSource;
}

which modified the term on the equation to include the correction factor you mention. You have there how to access the position and volume fields. In this case no lookup of the k field is necessary because is already refered in the class.

[lachuktr]

Originally Posted by crubio.abujas

In the turbulence model the equation is solved when you call the correct() method. There you can see where the equations are implemented:

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

Here you can define an additional source term or modify a source already existing. I've tried to modify the epsilonSource as follows:

template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> YkEpsilon<BasicTurbulenceModel>::epsilonSource() const
{
    tmp<fvScalarMatrix> epsSource
    (
        new fvScalarMatrix
        (
            epsilon_,
            dimVolume*this->rho_.dimensions()*epsilon_.dimensions()
            /dimTime
        )
    );

    const vectorField& C = this->mesh().C();
    const scalarField& V = this->mesh().V();

    scalar C_mu(0.09);
    scalar z0(0.03);

    scalarField& source = epsSource.ref().source();

    forAll(C, i)
    {
        source[i] -=
            C_mu * pow(k_[i] , 2) / 
            pow( C[i].z() + z0 , 2 ) * V[i];
    }

    return epsSource;
}

which modified the term on the equation to include the correction factor you mention. You have there how to access the position and volume fields. In this case no lookup of the k field is necessary because is already refered in the class.

thanks for the help again Carlos....

So can we do it like this to add a separate new source term (by giving a new name) rather than modifying the available one?...

if possible can you help me with that too

[crubio.abujas]

hi Lakshman,

I think you already have enough information to try it yourself. You just have to mix the solution you already have with the YkEpsilon source you shared with me.

If you find any particular problem you can figure out come with specific doubts.

Good luck!

[lachuktr]

Thanks a lot, carlos ... that helps a lot for me.

[lachuktr]

Originally Posted by crubio.abujas

hi Lakshman,

I think you already have enough information to try it yourself. You just have to mix the solution you already have with the YkEpsilon source you shared with me.

If you find any particular problem you can figure out come with specific doubts.

Good luck!

in that YkEpsilon source, Cmu is considered as a constant (0.09), but for my work, I want to change it to a variable as shown in the attachment.

I am not succeeding in the implementation. Kindly help

[crubio.abujas]

Hello Lakshman,

I think the most straightforward way to do it is creating a scalarField for C_mu. That should be quite simple to implement.

    scalar u_tau( 0.511);
    scalar z0(0.03);
    scalarField C_mu( pow(u_tau, 4) / pow(k_ , 2) );

    scalarField& source = epsSource.ref().source();

    forAll(C, i)
    {
        source[i] -=
            C_mu[i] * pow(k_[i] , 2) / 
            pow( C[i].z() + z0 , 2 ) * V[i];
     }

P.D: Ensure that you don't miss the C_mu[i] in the surce formula!