[Kbshariff]

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:

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

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]

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]

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

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]

What is line number 155 where the error occurred

[Kbshariff]

Originally Posted by mAlletto

What is line number 155 where the error occurred

This is the line 155

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)

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)

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]

Did you try

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

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

[Kbshariff]

Originally Posted by mAlletto

Did you try

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

fvc::Sp(A,B)

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

       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]

Hello,
Here is my topoSet Dict

    // 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]

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

[Kbshariff]

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

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]

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

[Kbshariff]

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]

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:

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:

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


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

topoSet:

    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]

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

        // 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]

Originally Posted by jherb

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

        // 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.