K-epsilon source term fvOptions

Kbshariff · May 18, 2021, 03:28

Hello foamers,

I would like to add a source term (developed by El-kasmi, 2008) to k-epsilon model in open foam. I found a similar source term defined using the atmTurbSource in OFv2006.
I copied the same file and modified the relation but I have an error in compilation.

Here is the original source term defined:

Code:

template<class AlphaFieldType, class RhoFieldType>
void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceEpsilon
(
    const AlphaFieldType& alpha,
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn,
    const label fieldi
) const
{

    const auto* turbPtr =
        mesh_.findObject<turbulenceModel>
        (
            turbulenceModel::propertiesName
        );

    const volScalarField& epsilon = turbPtr->epsilon();

    // (Heuristically derived from RS:Eq. 4, rhs-term:5)
    eqn +=
//        fvm::Sp(alpha()*rho()*C2_*sqr(epsilonAmb_)/(kAmb_*epsilon()), epsilon);
        fvm::Sp(alpha()*rho()*C2_*sqr(epsilonAmb_)/(kAmb_));
}


template<class AlphaFieldType, class RhoFieldType>
void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceOmega
(
    const AlphaFieldType& alpha,
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn,
    const label fieldi
) const
{
    const auto* turbPtr =
        mesh_.findObject<turbulenceModel>
        (
            turbulenceModel::propertiesName
        );
    const volScalarField& omega = turbPtr->omega();
    const volScalarField::Internal& beta =
        mesh_.lookupObjectRef<volScalarField::Internal>
        (
            word(turbPtr->type() + ":beta")
        );

    // (RS:Eq. 4, rhs-term:5)
//    eqn += fvm::Sp(alpha()*rho()*Cmu_*beta*sqr(omegaAmb_)/omega(), omega);
    eqn += fvm::Sp(alpha()*rho()*Cmu_*beta*sqr(omegaAmb_)/omega());
}


template<class AlphaFieldType, class RhoFieldType>
void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceK
(
    const AlphaFieldType& alpha,
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn,
    const label fieldi
) const
{
    const auto* turbPtr =
        mesh_.findObject<turbulenceModel>
        (
            turbulenceModel::propertiesName
        );
    const volScalarField& k = turbPtr->k();

    if (isEpsilon_)
    {
        // (Heuristically derived from RS:Eq. 3, rhs-term:4)
//        eqn += fvm::Sp(alpha()*rho()*epsilonAmb_/k(), k);
        eqn += fvm::Sp(alpha()*rho()*epsilonAmb_);
    }
    else
    {
        // (RS:Eq. 3, rhs-term:4)
//        eqn += fvm::Sp(alpha()*rho()*Cmu_*omegaAmb_*kAmb_/k(), k);
        eqn += fvm::Sp(alpha()*rho()*Cmu_*omegaAmb_*kAmb_);
    }
}

I have this compilation error

Code:

In file included from myAtmAmbientTurbSource.H:259:0,
                 from myAtmAmbientTurbSource.C:29:
myAtmAmbientTurbSourceTemplates.C: In instantiation of ‘void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceK(const AlphaFieldType&, const RhoFieldType&, Foam::fvMatrix<double>&, Foam::label) const [with AlphaFieldType = Foam::geometricOneField; RhoFieldType = Foam::geometricOneField; Foam::label = int]’:
myAtmAmbientTurbSource.C:170:9:   required from here
myAtmAmbientTurbSourceTemplates.C:106:49: error: no matching function for call to ‘Sp(const Foam::dimensioned<double>&)’
         eqn += fvm::Sp(alpha()*rho()*epsilonAmb_);
                                                 ^
myAtmAmbientTurbSourceTemplates.C:106:49: note: candidates are:
In file included from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/finiteVolume/lnInclude/fvmSup.H:165:0,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/finiteVolume/lnInclude/fvm.H:49,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude/linearViscousStress.C:30,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude/linearViscousStress.H:116,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude/Stokes.H:45,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude/laminarModel.C:30,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/TurbulenceModels/turbulenceModels/lnInclude/laminarModel.H:201,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/TurbulenceModels/incompressible/lnInclude/turbulentTransportModel.H:49,
                 from myAtmAmbientTurbSource.H:118,
                 from myAtmAmbientTurbSource.C:29:
/soft/OpenFOAM-v2006/OpenFOAM-v2006/src/finiteVolume/lnInclude/fvmSup.C:102:1: note: template<class Type> Foam::tmp<Foam::fvMatrix<Type> > Foam::fvm::Sp(const Internal&, const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&)
 Foam::fvm::Sp
 ^
/soft/OpenFOAM-v2006/OpenFOAM-v2006/src/finiteVolume/lnInclude/fvmSup.C:102:1: note:   template argument deduction/substitution failed:
In file included from myAtmAmbientTurbSource.H:259:0,
                 from myAtmAmbientTurbSource.C:29:
myAtmAmbientTurbSourceTemplates.C:106:49: note:   cannot convert ‘Foam::operator*<Foam::dimensioned<double> >((*(const Foam::one*)(& Foam::operator*((*(const Foam::oneField*)(&(& alpha)->Foam::geometricOneField::operator()())), (*(const Foam::oneField*)(&(& rho)->Foam::geometricOneField::operator()()))))), (* &((const Foam::fv::myAtmAmbientTurbSource*)this)->Foam::fv::myAtmAmbientTurbSource::epsilonAmb_))’ (type ‘const Foam::dimensioned<double>’) to type ‘const Internal& {aka const Foam::DimensionedField<double, Foam::volMesh>&}’
         eqn += fvm::Sp(alpha()*rho()*epsilonAmb_);

Kindly assist me in debugging please.

Kabir

mAlletto · May 18, 2021, 05:42

Code:

fvm::Sp(alpha()*rho()*C2_*sqr(epsilonAmb_)/(kAmb_))

Check out how implicit sources have to be defined. I didn't read your equation but the epsilon is missing in your term

Kbshariff · May 18, 2021, 06:28

I have now modified the equations. All variables like epsilonAmb_, kAmb_ were declared as constant in boundary conditions. All parameters are explicit now, I have suppressed all depending on implicit terms in.C file because it's not needed in the new equation.

Here is the error I am unable to debug now

Code:

myAtmAmbientTurbSource.C: In constructor ‘Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSource(const Foam::word&, const Foam::word&, const Foam::dictionary&, const Foam::fvMesh&)’:
myAtmAmbientTurbSource.C:155:1: error: expected ‘{’ before ‘void’
 void Foam::fv::myAtmAmbientTurbSource::addSup // I have checked but there is no unclosed  {}
 ^
In file included from myAtmAmbientTurbSource.H:259:0,
                 from myAtmAmbientTurbSource.C:29:
myAtmAmbientTurbSourceTemplates.C: In instantiation of ‘void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceK(const AlphaFieldType&, const RhoFieldType&, Foam::fvMatrix<double>&, Foam::label) const [with AlphaFieldType = Foam::geometricOneField; RhoFieldType = Foam::geometricOneField; Foam::label = int]’:
myAtmAmbientTurbSource.C:170:9:   required from here
myAtmAmbientTurbSourceTemplates.C:108:49: error: no matching function for call to ‘Sp(const Foam::dimensioned<double>&)’
         eqn += fvm::Sp(alpha()*rho()*epsilonAmb_);

THank you,

Kabir

mAlletto · May 18, 2021, 07:50

What is line number 155 where the error occurred

Kbshariff · May 18, 2021, 09:35

Quote:

Originally Posted by mAlletto

What is line number 155 where the error occurred

This is the line 155

Code:

void Foam::fv::myAtmAmbientTurbSource::addSup //line 155
(
    fvMatrix<scalar>& eqn,
    const label fieldi
)

{
    if (fieldi == 1)
    {
        myAtmAmbientTurbSourceK
        (
            geometricOneField(),
            geometricOneField(),
            eqn,
            fieldi
        );
        return;
    }

    if (isEpsilon_)
    {
        myAtmAmbientTurbSourceEpsilon
        (
            geometricOneField(),
            geometricOneField(),
            eqn,
            fieldi
        );
    }
    else
    {
        myAtmAmbientTurbSourceOmega
        (
            geometricOneField(),
            geometricOneField(),
            eqn,
            fieldi
        );
    }
}

I bypass the error by using empty braket {}.

Now I have these error ( a & b) and are repeated in several lines in the code
(a)

Code:

myAtmAmbientTurbSourceTemplates.C:114:58: error: no matching function for call to ‘Sp(Foam::tmp<Foam::DimensionedField<double, Foam::volMesh> >)’
         eqn += fvm::Sp(alpha()*rho()*Cmu_*omegaAmb_*kAmb_);
// corresponing line 114 is shown in blue below

void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceK
(
    const AlphaFieldType& alpha,
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn,
    const label fieldi
) const
{
/*
    const auto* turbPtr =
        mesh_.findObject<turbulenceModel>
        (
            turbulenceModel::propertiesName
        );
    //const volScalarField& k = turbPtr->k();
*/
    if (isEpsilon_)
    {
        // (Heuristically derived from RS:Eq. 3, rhs-term:4)
 //       eqn += fvm::Sp(alpha()*rho()*epsilonAmb_/k(), k);
        eqn += fvm::Sp(alpha()*rho()*epsilonAmb_);
    }
    else
    {
        // (RS:Eq. 3, rhs-term:4)
//        eqn += fvm::Sp(alpha()*rho()*Cmu_*omegaAmb_*kAmb_/k(), k);
        eqn += fvm::Sp(alpha()*rho()*Cmu_*omegaAmb_*kAmb_);
    }
}

(b)

Code:

myAtmAmbientTurbSourceTemplates.C:83:53: note:   candidate expects 2 arguments, 1 provided
     eqn += fvm::Sp(alpha()*rho()*Cmu_*sqr(omegaAmb_));
// corresponding line 83 is shown in blue below

template<class AlphaFieldType, class RhoFieldType>
void Foam::fv::myAtmAmbientTurbSource::myAtmAmbientTurbSourceOmega
(
    const AlphaFieldType& alpha,
    const RhoFieldType& rho,
    fvMatrix<scalar>& eqn,
    const label fieldi
) const
{
/*
    const auto* turbPtr =
        mesh_.findObject<turbulenceModel>
        (
            turbulenceModel::propertiesName
        );
    const volScalarField& omega = turbPtr->omega();
    const volScalarField::Internal& beta =
        mesh_.lookupObjectRef<volScalarField::Internal>
        (
            word(turbPtr->type() + ":beta")
        );
*/
    // (RS:Eq. 4, rhs-term:5)
//    eqn += fvm::Sp(alpha()*rho()*Cmu_*beta*sqr(omegaAmb_)/omega(), omega);
    eqn += fvm::Sp(alpha()*rho()*Cmu_*sqr(omegaAmb_));
}

Merci

mAlletto · May 18, 2021, 10:30

Did you try

Code:

fvc::Sp(alpha()*rho()*C2_*sqr(epsilonAmb_)/(kAmb_))

Please have a look how explicit and implicit sources are implemented in OF

Kbshariff · May 26, 2021, 09:15

Quote:

Originally Posted by mAlletto

Did you try

Code:

fvc::Sp(alpha()*rho()*C2_*sqr(epsilonAmb_)/(kAmb_))

Please have a look at how explicit and implicit sources are implemented in OF

Sorry for my late reply.

I have check how implicit and explicit terms are applied.
the implicit is in the form of

Code:

fvc::Sp(A,B)

In my case since all the values are constants, I applied the source term as:

Code:

       eqn += alpha()*rho()*epsilonAmb_;;

And it compiled successfully.
The problem now is that the source term is applied in the entire domain not only the region I defined as topoSet.
Here are the results without source and with a source at disc region.
How do I limit the source term to the cellSet region only, please?

Thank you

Kbshariff · June 2, 2021, 06:41

Hello,
Here is my topoSet Dict

Code:

    // sourceDisk1
    {
        name    sourceDisk1CellSet;
        type    cellSet;
        action  new;
        source  cylinderToCell;
        sourceInfo
        {
            radius 0.135;
		p1 (5.9 2.5 0.225);
		p2 (6.1 2.5 0.225);
        }
    }
    {
        name    sourceDisk1;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set sourceDisk1CellSet;
        }
    }

My problem now is the source is applied throughout the domain. How can I apply it to soureDisk only?

Thank you.

mAlletto · June 2, 2021, 08:31

You can specify for a fvoption the cellset where it is applied

Kbshariff · June 2, 2021, 09:04

Here is my fvOptions, I have specified the location defined in topoSet for myAtmTurbSource as follows:

Code:

disk1
{
    type            actuationDiskSource;
    variant         Froude;
    selectionMode   cellSet;
    cellSet         actuationDisk1;
    diskArea        0.50265;
    diskDir         (1 0 0);
    writeToFile     true;
    sink            true;
    Cp		    0.42;
    Ct		    0.83;


    
    monitorMethod   cellSet;
    cellSet         actuationDisk1;

}

atmAmbientTurbSource1
{
    type                atmAmbientTurbSource;
    atmAmbientTurbSourceCoeffs
    {
   	selectionMode   cellSet;
    	cellSet         actuationDisk1;
        kAmb            0.0216;   	// +15%
        epsilonAmb      3.725e-3;
    }
}

The code disk1 applies the actuationDisk term at the defined location but the ambTurbSource applies it in the whole domain as shown in the previous post.

mAlletto · June 5, 2021, 04:28

How do you know the fvoption applies in the whole domain and not only in the cellset you specified

Kbshariff · June 8, 2021, 09:00

Quote:

Originally Posted by mAlletto

How do you know the fvoption applies in the whole domain and not only in the cellset you specified

The figure below shows the values of TI, k, and epsilon for two simulations (a) without source and (b) with source applied to the disc

(a) The profile without source is decaying as expected
(b) the profile with a source at the disc maintains the kAmb and epsilonAmb throughout the domain.
Moreover, when I used the 'select all' instead of 'cellset' (actuatorDisc) the result is also the same as the source in the disc region only.

Kbshariff · June 11, 2021, 10:59

Hi,

I think I found the reason why the source term is applied in the whole domain instead of the region defined.

I run a simulation of the channel with the source term but without the actuator disc to see the effect of the source term.

From my simpleFoam log, I found this:

Code:

Creating finite volume options from "constant/fvOptions"

Selecting finite volume options type atmAmbientTurbSource
    Source: atmAmbientTurbSource1
    - selecting cells using cellZone zone1
    - selected 12325 cell(s) with volume 0.0929195759224
    Applying atmAmbientTurbSource to: epsilon and k

Starting time loop

turbulenceFields turbulenceFields: storing fields:
    turbulenceProperties:I

From the log file, the source is applied to epsilon and k,

How can I change it to

Applying atmAmbientTurbSource to: zone1 ?

Here is my fvOptions and topoSet

fvOptions:

Code:

    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


atmAmbientTurbSource1
{
    type                atmAmbientTurbSource;
    atmAmbientTurbSourceCoeffs
    {
   	selectionMode   cellZone;
    	cellZone       	zone1;
        kAmb            0.0096;   	// +10%
        epsilonAmb      1.1e-3;
    }
}

topoSet:

Code:

    object      topoSetDict;
}

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

actions
(
    {
        name    box1b;
        type    cellSet;
        action  new;
        source  boxToCell;
        box     (3.8 2.4 0.4) (4.2 2.6 1.6);
    }

    {
        name    zone1;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        set     box1b;
    }



);

Merci

- Kabir

jherb · June 12, 2021, 16:16

The documentation https://www.openfoam.com/documentati...urbSource.html says:

Code:

        // Mandatory (inherited) entries (unmodifiable)
        selectionMode    all;

If you look at the source code of other source terms, e.g. https://github.com/OpenFOAM/OpenFOAM...sSource.C#L128 there are loops over the cells, the source should be applied to.
So you can expect, that your source is applied in the whole flow domain.

Are you using OpenFOAM.com or OpenFOAM.org? For the first one, here is another example: https://develop.openfoam.com/Develop...gSource.C#L193

Kbshariff · June 15, 2021, 06:42

Quote:

Originally Posted by jherb

The documentation https://www.openfoam.com/documentati...urbSource.html says:

Code:

        // Mandatory (inherited) entries (unmodifiable)
        selectionMode    all;

If you look at the source code of other source terms, e.g. https://github.com/OpenFOAM/OpenFOAM...sSource.C#L128 there are loops over the cells, the source should be applied to.
So you can expect, that your source is applied in the whole flow domain.

Are you using OpenFOAM.com or OpenFOAM.org? For the first one, here is another example: https://develop.openfoam.com/Develop...gSource.C#L193

Hi,

I am using OpenFOAM.com, specifically v2006. Thanks, I will check the new example.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Adding source term with fvOptions to k-e model	Tibo99	OpenFOAM Programming & Development	20	March 15, 2023 07:12
[swak4Foam] swak4foam for OpenFOAM 4.0	mnikku	OpenFOAM Community Contributions	80	May 17, 2022 08:06
Source term for EVAPORATION in Energy Equ. - technical difficulty ?	Kummi	OpenFOAM	1	September 9, 2019 09:32
[foam-extend.org] Problems installing foam-extend-4.0 on openSUSE 42.2 and Ubuntu 16.04	ordinary	OpenFOAM Installation	19	September 3, 2019 18:13
DxFoam reader update	hjasak	OpenFOAM Post-Processing	69	April 24, 2008 01:24

May 18, 2021, 05:42		#2
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 615 Rep Power: 15	Code: fvm::Sp(alpha()rho()C2_*sqr(epsilonAmb_)/(kAmb_)) Check out how implicit sources have to be defined. I didn't read your equation but the epsilon is missing in your term

May 18, 2021, 07:50		#4
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 615 Rep Power: 15	What is line number 155 where the error occurred

May 18, 2021, 10:30		#6
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 615 Rep Power: 15	Did you try Code: fvc::Sp(alpha()rho()C2_*sqr(epsilonAmb_)/(kAmb_)) Please have a look how explicit and implicit sources are implemented in OF

June 2, 2021, 08:31		#9
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 615 Rep Power: 15	You can specify for a fvoption the cellset where it is applied

June 5, 2021, 04:28		#11
mAlletto Senior Member Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 615 Rep Power: 15	How do you know the fvoption applies in the whole domain and not only in the cellset you specified

June 12, 2021, 16:16		#14
jherb Senior Member Joachim Herb Join Date: Sep 2010 Posts: 650 Rep Power: 21	The documentation https://www.openfoam.com/documentati...urbSource.html says: Code: // Mandatory (inherited) entries (unmodifiable) selectionMode all; If you look at the source code of other source terms, e.g. https://github.com/OpenFOAM/OpenFOAM...sSource.C#L128 there are loops over the cells, the source should be applied to. So you can expect, that your source is applied in the whole flow domain. Are you using OpenFOAM.com or OpenFOAM.org? For the first one, here is another example: https://develop.openfoam.com/Develop...gSource.C#L193