Commit ab0dd1f1 authored by Kutalmis Bercin's avatar Kutalmis Bercin Committed by Andrew Heather
Browse files

DEFEATURE: deprecate v2f model in favour of kEpsilonPhitF

  - kEpsilonPhitF is a kEpsilon-based model which originated
    from (Durbin, 1995)’s v2-f methodology. However, the majority of
    v2-f model variants proved to be numerically stiff for segregated
    solution algorithms due to the coupled formulations of v2 and f fields,
    particularly on wall boundaries.

    The v2-f variant (i.e. OpenFOAM’s v2f model) due to
    (Lien and Kalitzin, 2001) reformulated the original v2-f model to enable
    segregated computations; however, a number of shortcomings regarding
    the model fidelity were reported in the literature.

    To overcome the shortcomings of the v2-f methodology, the v2-f approach
    was re-evaluated by (Laurence et al., 2005) by transforming v2 scale into
    its equivalent non-dimensional form, i.e. phit, to reduce the numerical
    stiffness.
    This variant, i.e. kEpsilonPhitF, is believed to provide numerical
    robustness, and insensitivity to grid anomalies while retaining the
    theoretical model fidelity of the original v2-f model.

    Accordingly the v2f RANS model is deprecated in favour of the variant
    kEpsilonPhitF model.
parent 0ff95341
......@@ -85,9 +85,6 @@ makeRASModel(kOmegaSSTSAS);
#include "kOmegaSSTLM.H"
makeRASModel(kOmegaSSTLM);
#include "v2f.H"
makeRASModel(v2f);
#include "LRR.H"
makeRASModel(LRR);
......
......@@ -81,9 +81,6 @@ makeRASModel(kOmegaSSTSAS);
#include "kOmegaSSTLM.H"
makeRASModel(kOmegaSSTLM);
#include "v2f.H"
makeRASModel(v2f);
#include "LRR.H"
makeRASModel(LRR);
......
turbulenceModel.C
RAS/v2f/v2fBase.C
RAS/kEpsilonPhitF/kEpsilonPhitFBase.C
LESdelta = LES/LESdeltas
......@@ -57,12 +55,6 @@ kqRWallFunctions = $(wallFunctions)/kqRWallFunctions
$(kqRWallFunctions)/kqRWallFunction/kqRWallFunctionFvPatchFields.C
$(kqRWallFunctions)/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C
v2WallFunctions = $(wallFunctions)/v2WallFunctions
$(v2WallFunctions)/v2WallFunction/v2WallFunctionFvPatchScalarField.C
fWallFunctions = $(wallFunctions)/fWallFunctions
$(fWallFunctions)/fWallFunction/fWallFunctionFvPatchScalarField.C
RASBCs = RAS/derivedFvPatchFields
......
......@@ -119,7 +119,6 @@ kEpsilonPhitF<BasicTurbulenceModel>::kEpsilonPhitF
transport,
propertiesName
),
kEpsilonPhitFBase(),
Cmu_
(
......
......@@ -78,14 +78,13 @@ SourceFiles
kEpsilonPhitF.C
SeeAlso
v2f.H
kEpsilon.C
\*---------------------------------------------------------------------------*/
#ifndef kEpsilonPhitF_H
#define kEpsilonPhitF_H
#include "kEpsilonPhitFBase.H"
#include "RASModel.H"
#include "eddyViscosity.H"
......@@ -103,8 +102,7 @@ namespace RASModels
template<class BasicTurbulenceModel>
class kEpsilonPhitF
:
public eddyViscosity<RASModel<BasicTurbulenceModel>>,
public kEpsilonPhitFBase
public eddyViscosity<RASModel<BasicTurbulenceModel>>
{
// Private Member Functions
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "kEpsilonPhitFBase.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
defineTypeNameAndDebug(kEpsilonPhitFBase, 0);
} // End namespace RASModels
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::RASModels::kEpsilonPhitFBase
Group
grpRASTurbulence
Description
Abstract base-class for the \c k-epsilon-phit-f model to provide boundary
condition access to the \c phit and \c f fields.
See also
Foam::RASModels::kEpsilonPhitF
SourceFiles
kEpsilonPhitFBase.C
\*---------------------------------------------------------------------------*/
#ifndef kEpsilonPhitFBase_H
#define kEpsilonPhitFBase_H
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class kEpsilonPhitFBase Declaration
\*---------------------------------------------------------------------------*/
class kEpsilonPhitFBase
{
public:
//- Runtime type information
TypeName("kEpsilonPhitFBase");
// Constructors
kEpsilonPhitFBase()
{}
//- Destructor
virtual ~kEpsilonPhitFBase()
{}
// Member Functions
//- Return the normalised wall-normal fluctuating velocity scale field
virtual tmp<volScalarField> phit() const = 0;
//- Return the elliptic relaxation factor field
virtual tmp<volScalarField> f() const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2012-2017 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "v2f.H"
#include "fvOptions.H"
#include "bound.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class BasicTurbulenceModel>
tmp<volScalarField> v2f<BasicTurbulenceModel>::Ts() const
{
// SAF: limiting thermo->nu(). If psiThermo is used rho might be < 0
// temporarily when p < 0 then nu < 0 which needs limiting
return
max
(
k_/epsilon_,
6.0*sqrt
(
max
(
this->nu(),
dimensionedScalar(this->nu()().dimensions(), Zero)
)
/ epsilon_
)
);
}
template<class BasicTurbulenceModel>
tmp<volScalarField> v2f<BasicTurbulenceModel>::Ls() const
{
// SAF: limiting thermo->nu(). If psiThermo is used rho might be < 0
// temporarily when p < 0 then nu < 0 which needs limiting
return
CL_
* max
(
pow(k_, 1.5)/epsilon_,
Ceta_*pow025
(
pow3
(
max
(
this->nu(),
dimensionedScalar(this->nu()().dimensions(), Zero)
)
)/epsilon_
)
);
}
template<class BasicTurbulenceModel>
void v2f<BasicTurbulenceModel>::correctNut()
{
this->nut_ = min(CmuKEps_*sqr(k_)/epsilon_, this->Cmu_*v2_*Ts());
this->nut_.correctBoundaryConditions();
fv::options::New(this->mesh_).correct(this->nut_);
BasicTurbulenceModel::correctNut();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
v2f<BasicTurbulenceModel>::v2f
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaRhoPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName,
const word& type
)
:
eddyViscosity<RASModel<BasicTurbulenceModel>>
(
type,
alpha,
rho,
U,
alphaRhoPhi,
phi,
transport,
propertiesName
),
v2fBase(),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
this->coeffDict_,
0.22
)
),
CmuKEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CmuKEps",
this->coeffDict_,
0.09
)
),
C1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C1",
this->coeffDict_,
1.4
)
),
C2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C2",
this->coeffDict_,
0.3
)
),
CL_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CL",
this->coeffDict_,
0.23
)
),
Ceta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ceta",
this->coeffDict_,
70.0
)
),
Ceps2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ceps2",
this->coeffDict_,
1.9
)
),
Ceps3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ceps3",
this->coeffDict_,
-0.33
)
),
sigmaK_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaK",
this->coeffDict_,
1.0
)
),
sigmaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaEps",
this->coeffDict_,
1.3
)
),
k_
(
IOobject
(
IOobject::groupName("k", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
epsilon_
(
IOobject
(
IOobject::groupName("epsilon", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
v2_
(
IOobject
(
IOobject::groupName("v2", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
f_
(
IOobject
(
IOobject::groupName("f", alphaRhoPhi.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
v2Min_(dimensionedScalar("v2Min", v2_.dimensions(), SMALL)),
fMin_(dimensionedScalar("fMin", f_.dimensions(), Zero))
{
bound(k_, this->kMin_);
bound(epsilon_, this->epsilonMin_);
bound(v2_, v2Min_);
bound(f_, fMin_);
if (type == typeName)
{
this->printCoeffs(type);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
bool v2f<BasicTurbulenceModel>::read()
{
if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
{
Cmu_.readIfPresent(this->coeffDict());
CmuKEps_.readIfPresent(this->coeffDict());
C1_.readIfPresent(this->coeffDict());
C2_.readIfPresent(this->coeffDict());
CL_.readIfPresent(this->coeffDict());
Ceta_.readIfPresent(this->coeffDict());
Ceps2_.readIfPresent(this->coeffDict());
Ceps3_.readIfPresent(this->coeffDict());
sigmaK_.readIfPresent(this->coeffDict());
sigmaEps_.readIfPresent(this->coeffDict());
return true;
}
return false;
}
template<class BasicTurbulenceModel>
void v2f<BasicTurbulenceModel>::correct()
{
if (!this->turbulence_)
{
return;
}
// Local references
const alphaField& alpha = this->alpha_;
const rhoField& rho = this->rho_;
const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
const volVectorField& U = this->U_;
volScalarField& nut = this->nut_;
fv::options& fvOptions(fv::options::New(this->mesh_));
eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
// Use N=6 so that f=0 at walls
const dimensionedScalar N("N", dimless, 6.0);
const volTensorField gradU(fvc::grad(U));
const volScalarField S2(2*magSqr(dev(symm(gradU))));
const volScalarField G(this->GName(), nut*S2);
const volScalarField Ts(this->Ts());
const volScalarField L2(type() + ":L2", sqr(Ls()));
const volScalarField v2fAlpha
(
type() + ":alpha",
1.0/Ts*((C1_ - N)*v2_ - 2.0/3.0*k_*(C1_ - 1.0))
);
const volScalarField Ceps1
(
"Ceps1",
1.4*(1.0 + 0.05*min(sqrt(k_/v2_), scalar(100)))
);
// Update epsilon (and possibly G) at the wall
epsilon_.boundaryFieldRef().updateCoeffs();
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(alpha, rho, epsilon_)
+ fvm::div(alphaRhoPhi, epsilon_)
- fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
==
Ceps1*alpha*rho*G/Ts
- fvm::SuSp(((2.0/3.0)*Ceps1 + Ceps3_)*alpha*rho*divU, epsilon_)
- fvm::Sp(Ceps2_*alpha*rho/Ts, epsilon_)
+ fvOptions(alpha, rho, epsilon_)
);
epsEqn.ref().relax();
fvOptions.constrain(epsEqn.ref());
epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
solve(epsEqn);
fvOptions.correct(epsilon_);
bound(epsilon_, this->epsilonMin_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k_)