[carlDahmen]

Hello Foamer!

I am working with CFDEM, but my question is OpenFOAM related.
My aim is to modify the Smagorinsky turbulence model and implement a viscosity that depends on the local voidfraction according to Batchelor and Green (1972) as: nu_t = nu*(a*(1-voidfraction)+b*(1-voidfraction)^2)

Since my OpenFOAM knowledge is lacking, have I tried to follow the Smagorinsky model as closely as I can and followed guides to implement a new turbelence model. But my problem comes down to accessing the voidfraction, I assume this can be done by using IOobject but I do not know how.

Here is my code (Protected Member Functions):

template<class BasicTurbulenceModel>
tmp<volScalarField> BatchelorGreen<BasicTurbulenceModel>::k
(
    const tmp<volScalarField>& voidfrac
) const
{
    volScalarField solidfrac(scalar(1)-voidfrac);

    volScalarField c1(nu_*a_*solidfrac);
    volScalarField c2(nu_*b_*sqr(solidfrac));

    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("k", this->U_.group()),
                this->runTime_.timeName(),
                this->mesh_
            ),
            c1+c2
        )
    );
}

template<class BasicTurbulenceModel>
void BatchelorGreen<BasicTurbulenceModel>::correctNut()
{
    new volScalarField
    (
	IOobject
	(
	    IOobject::groupName("voidfraction", this->U_.group()),
	    this->runTime_.timeName(),
	    this->mesh_,
	    IOobject::MUST_READ,
	    IOobject::NO_WRITE
	)
    )
	
    volScalarField k(this->k(this->voidfraction_)); 

    this->nut_ = k;
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

[carlDahmen]

Hey!


My code runs now, need to verify the result. It is attached for anyone interested.

\*---------------------------------------------------------------------------*/

#include "BatchelorGreen.H"
#include "fvOptions.H"

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

namespace Foam
{
namespace LESModels
{

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

template<class BasicTurbulenceModel>
tmp<volScalarField> BatchelorGreen<BasicTurbulenceModel>::k
(
    const tmp<volScalarField>& voidfrac
) const
{
    volScalarField solidfrac(scalar(1)-voidfrac);

    volScalarField c1(a_*solidfrac+b_*sqr(solidfrac));

    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("k", this->U_.group()),
                this->runTime_.timeName(),
                this->mesh_
            ),
            c1
        )
    );
}

template<class BasicTurbulenceModel>
void BatchelorGreen<BasicTurbulenceModel>::correctNut()
{
    
    volScalarField k(voidfraction_); 

    this->nut_ = this->nu()*k;
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
BatchelorGreen<BasicTurbulenceModel>::BatchelorGreen
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    LESeddyViscosity<BasicTurbulenceModel>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    a_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "a",
            this->coeffDict_,
            2.5
        )
    ),

    b_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "b",
            this->coeffDict_,
            7.6
        )
    ),

    voidfraction_
    (
    IOobject
    (
        "voidfraction",
        this->runTime_.timeName(),
            this->mesh_,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    ),
    this->mesh_
    )

{
    if (type == typeName)
    {
        this->printCoeffs(type);
    }
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
bool BatchelorGreen<BasicTurbulenceModel>::read()
{
    if (LESeddyViscosity<BasicTurbulenceModel>::read())
    {
        a_.readIfPresent(this->coeffDict());
    b_.readIfPresent(this->coeffDict());

        return true;
    }
    else
    {
        return false;
    }
}



template<class BasicTurbulenceModel>
void BatchelorGreen<BasicTurbulenceModel>::correct()
{
    LESeddyViscosity<BasicTurbulenceModel>::correct();
    correctNut();
}


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

} // End namespace LESModels
} // End namespace Foam

[carlDahmen]

Seems like the voidfraction is not updated, how do I declare it in the Header file?

[joshwilliams]

Originally Posted by carlDahmen

Seems like the voidfraction is not updated, how do I declare it in the Header file?

In your header file you can put

volVectorField UsTilde_;

Or in your *.C file you can put somewhere

    const volScalarField& voidFraction = mesh.lookupObject<volScalarField>("voidFraction");

[carlDahmen]

Thank you for your answer!

I have put in my correctNut() in the *.C file:

const volScalarField& voidfraction_ = this->mesh().lookupObject<volScalarField>("voidfraction");

But I get the error message:

error: expected primary-expression before '>' token

Do you now what's wrong?

[carlDahmen]

It's a template class and I needed:

const volScalarField& voidfraction_ = this->mesh().template lookupObject<volScalarField>("voidfraction");