[stefan.gracik]

Has anyone implemented any zero equation turbulence models to OpenFOAM? That is, a turbulence model which uses algebraic equations to calculate turbulent viscosity rather than PDE's (k, ε etc.). The classic example of this is the Prandtl mixing length model.

I'm trying to implement a specialized model and am having trouble modifying the 2 and 1 eq models to remove the PDE's, so I figured having a zero equation model already implemented would help me understand what I need to do.

Thanks

[stefan.gracik]

I've made some progress on this, but am stuck again. It is an equation designed for modeling external airflow over buildings. The turbulence model I'm trying to implement is depended on this distance to wall, where H is the average height of the building (implemented as a dimensionedScalar)

for walldist <= 1.3*H, nut = const*exp(walldist/H)*U

and

for walldist >1.3*H, nut = const*walldist ... etc.

I thought the best way to write this would be with an if statement. Currently I have written it like this (where d_ is the walldist)

 
 
    if
    (d_ <= (1.3*H_))
    {
    nut_ = a_*zh_*exp(-b_*d_/H_)*U_*pow((d_/H_),2);
    nut_.correctBoundaryConditions();
    }
    else
    {
    nut_ = 0.16*Uh_*(d_ + z0_)/Foam::log((d_ + z0_)/z0_);
    nut_.correctBoundaryConditions();
    }

but receive a pretty long error

wmakeLnInclude: linking include files to ./lnInclude
Making dependency list for source file mykEpsilon.C
SOURCE=mykEpsilon.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-100 -I/opt/openfoam211/src/turbulenceModels -I/opt/openfoam211/src/transportModels -I/opt/openfoam211/src/finiteVolume/lnInclude -I/opt/openfoam211/src/meshTools/lnInclude -I/opt/openfoam211/src/turbulenceModels/incompressible/RAS/lnInclude -IlnInclude -I. -I/opt/openfoam211/src/OpenFOAM/lnInclude -I/opt/openfoam211/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/mykEpsilon.o
mykEpsilon.H: In constructor ‘Foam::incompressible::RASModels::mykEpsilon::mykEpsilon(const volVectorField&, const surfaceScalarField&, Foam::transportModel&, const Foam::word&, const Foam::word&)’:
mykEpsilon.H:81:31: warning: ‘Foam::incompressible::RASModels::mykEpsilon::H_’ will be initialized after [-Wreorder]
mykEpsilon.H:77:31: warning:   ‘Foam::dimensionedScalar Foam::incompressible::RASModels::mykEpsilon::Uh_’ [-Wreorder]
mykEpsilon.C:47:1: warning:   when initialized here [-Wreorder]
mykEpsilon.C:160:19: error: no match for ‘operator<=’ in ‘((Foam::incompressible::RASModels::mykEpsilon*)this)->Foam::incompressible::RASModels::mykEpsilon::d_ <= Foam::operator*(const Foam::dimensioned<double>&, const Foam::dimensioned<Type>&) [with Type = double]((*(const Foam::dimensioned<double>*)(&((Foam::incompressible::RASModels::mykEpsilon*)this)->Foam::incompressible::RASModels::mykEpsilon::H_)))’
mykEpsilon.C:160:19: note: candidates are:
/opt/openfoam211/src/OpenFOAM/lnInclude/UList.C:224:6: note: bool Foam::UList<T>::operator<=(const Foam::UList<T>&) const [with T = double]
/opt/openfoam211/src/OpenFOAM/lnInclude/UList.C:224:6: note:   no known conversion for argument 1 from ‘Foam::dimensioned<double>’ to ‘const Foam::UList<double>&’
/opt/openfoam211/src/OpenFOAM/lnInclude/VectorSpaceI.H:693:13: note: template<class Form, class Cmpt, int nCmpt> bool Foam::operator<=(const Foam::VectorSpace<Form, Cmpt, nCmpt>&, const Foam::VectorSpace<Form, Cmpt, nCmpt>&)
mykEpsilon.C:162:50: error: no match for ‘operator=’ in ‘((Foam::incompressible::RASModels::mykEpsilon*)this)->Foam::incompressible::RASModels::mykEpsilon::nut_ = Foam::operator*(const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&, const Foam::tmp<Foam::GeometricField<double, PatchField, GeoMesh> >&) [with Type = Foam::Vector<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]((*(const Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >*)(& Foam::pow(const Foam::tmp<Foam::GeometricField<double, PatchField, GeoMesh> >&, const scalar&) [with PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh, Foam::scalar = double]((* &2.0e+0)))))’
mykEpsilon.C:162:50: note: candidates are:
/opt/openfoam211/src/OpenFOAM/lnInclude/GeometricField.C:1083:6: note: void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(const Foam::GeometricField<Type, PatchField, GeoMesh>&) [with Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/opt/openfoam211/src/OpenFOAM/lnInclude/GeometricField.C:1083:6: note:   no known conversion for argument 1 from ‘Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> >’ to ‘const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&’
/opt/openfoam211/src/OpenFOAM/lnInclude/GeometricField.C:1108:6: note: void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&) [with Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/opt/openfoam211/src/OpenFOAM/lnInclude/GeometricField.C:1108:6: note:   no known conversion for argument 1 from ‘Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> >’ to ‘const Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >&’
/opt/openfoam211/src/OpenFOAM/lnInclude/GeometricField.C:1144:6: note: void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(const Foam::dimensioned<Form>&) [with Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/opt/openfoam211/src/OpenFOAM/lnInclude/GeometricField.C:1144:6: note:   no known conversion for argument 1 from ‘Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> >’ to ‘const Foam::dimensioned<double>&’
make: *** [Make/linux64GccDPOpt/mykEpsilon.o] Error 1
stefan@stefan-OpenFOAM:~/OpenFOAM/stefan-2.1.1/src/turbulenceModels/incompressible/RAS/mykEpsilon$

Any idea what the issue is? I have a feeling it has something to do with comparing walldist to a dimensionedScalar, but I'm really not sure.

[cosimobianchini]

If I did not misunderstand you should use your if statement within a loop on the internal cells of your volScalarField nut_. What the error is saying is that it makes no sense a statement like <= for a volScalarField in fact for some cells it will be > for others <.

In order to make it work write something like

forAll(nut_,cellI)
{
if(d_[cellI] <= (1.3*H_))
{
nut_[cellI] = a_*zh_*exp(-b_*d_[cellI]/H_)*mag(U_)*pow((d_[cellI]/H_),2);
}
else
{
nut_[cellI] = 0.16*Uh_*(d_[cellI] + z0_)/Foam::log((d_[cellI] + z0_)/z0_);
}
}

assuming that all other variables except for d_ and U_ are of type scalar this should work.
Another thing you should pay attention to is in the use of vector variables (I guessed U_ it is velocity vector) in a statement with scalars on the LHS, I added the operator mag(U_) to evaluate the magnitude of velocity but I'm not sure this is what you want. Change the operator accordingly to your needs but check that it returns a scalar.

[stefan.gracik]

Thanks for the help!! but now I have a new error. d_ is generated from walldist.H, so it's the distance to the nearest wall. All of the other constants I've defined as dimensionedScalar

a_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "a",
            coeffDict_,
            0.09
        )
    ),
    b_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "b",
            coeffDict_,
            0.09
        )
    ),
    H_
    (
        dimensioned<scalar>//::lookupOrAddToDict
        (
            "H",
            dimensionSet(0, 1, 0, 0, 0, 0, 0),
            1.92
        )
    ),

    Uh_
    (
        dimensioned<scalar>//::lookupOrAddToDict
        (
            "Uh",
        dimensionSet(0, 1, -1, 0, 0, 0, 0),
            1.44
        )
    ),

    zh_
    (
        dimensioned<scalar>//::lookupOrAddToDict
        (
            "zh",
            dimensionSet(0, 1, 0, 0, 0, 0, 0),
            1.92
        )
    ),
    z0_
    (
        dimensioned<scalar>//::lookupOrAddToDict
        (
            "z0",
            dimensionSet(0,1,0,0,0,0,0),
            1.3
        )
    ),

but I still receive the error
"mykEpsilon.C:180:83: error: cannot convert ‘Foam::dimensioned<double>’ to ‘double’ in assignment"

Anyone have any idea why this could be?