Commit ffa363b3 authored by Henry Weller's avatar Henry Weller
Browse files

epsilonWallFunction: Updated to work with both low- and high-Reynolds number turbulence models

This boundary condition provides a turbulence dissipation wall constraint
for low- and high-Reynolds number turbulence models.

The condition can be applied to wall boundaries for which it
- calculates \c epsilon and \c G
- specifies the near-wall epsilon value

where

\vartable
    epsilon | turblence dissipation field
    G       | turblence generation field
\endvartable

The model switches between laminar and turbulent functions based on the
laminar-to-turbulent y+ value derived from kappa and E.

Recent tests have shown that this formulation is more accurate than
the standard high-Reynolds number form for 10 < y+ < 30 with both
standard and continuous wall-functions.

Replaces epsilonLowReWallFunction and should be used for all
low-Reynolds number models for which the epsilonLowReWallFunction BC was
recommended.
parent 1405a310
......@@ -43,7 +43,6 @@ $(nutWallFunctions)/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarFiel
epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
$(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
$(epsilonWallFunctions)/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
omegaWallFunctions = $(wallFunctions)/omegaWallFunctions
$(omegaWallFunctions)/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
......
......@@ -41,7 +41,7 @@ Description
Wall boundary conditions are:
k = kLowReWallFunction
epsilon = epsilonLowReWallFunction
epsilon = epsilonWallFunction
v2 = v2WallFunction
f = fWallFunction
......@@ -90,7 +90,7 @@ See also
Foam::RASModels::v2fBase
Foam::RASModels::kEpsilon
Foam::kLowReWallFunctionFvPatchScalarField
Foam::epsilonLowReWallFunctionFvPatchScalarField
Foam::epsilonWallFunctionFvPatchScalarField
Foam::v2WallFunctionFvPatchScalarField
Foam::fWallFunctionFvPatchScalarField
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "epsilonLowReWallFunctionFvPatchScalarField.H"
#include "turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::scalar Foam::epsilonLowReWallFunctionFvPatchScalarField::yPlusLam
(
const scalar kappa,
const scalar E
)
{
scalar ypl = 11.0;
for (int i=0; i<10; i++)
{
ypl = log(max(E*ypl, 1))/kappa;
}
return ypl;
}
void Foam::epsilonLowReWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbModel,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G0,
scalarField& epsilon0
)
{
const label patchi = patch.index();
const scalarField& y = turbModel.y()[patchi];
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
const tmp<volScalarField> tk = turbModel.k();
const volScalarField& k = tk();
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
const scalarField magGradUw(mag(Uw.snGrad()));
const DimensionedField<scalar, volMesh>& G =
db().lookupObject<DimensionedField<scalar, volMesh>>
(
turbModel.GName()
);
// Set epsilon and G
forAll(nutw, facei)
{
label celli = patch.faceCells()[facei];
scalar yPlus = Cmu25*sqrt(k[celli])*y[facei]/nuw[facei];
scalar w = cornerWeights[facei];
if (yPlus > yPlusLam_)
{
epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
G0[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(kappa_*y[facei]);
}
else
{
epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
G0[celli] += G[celli];
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::epsilonLowReWallFunctionFvPatchScalarField::
epsilonLowReWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
epsilonWallFunctionFvPatchScalarField(p, iF),
yPlusLam_(yPlusLam(kappa_, E_))
{}
Foam::epsilonLowReWallFunctionFvPatchScalarField::
epsilonLowReWallFunctionFvPatchScalarField
(
const epsilonLowReWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
epsilonWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
yPlusLam_(ptf.yPlusLam_)
{}
Foam::epsilonLowReWallFunctionFvPatchScalarField::
epsilonLowReWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
epsilonWallFunctionFvPatchScalarField(p, iF, dict),
yPlusLam_(yPlusLam(kappa_, E_))
{}
Foam::epsilonLowReWallFunctionFvPatchScalarField::
epsilonLowReWallFunctionFvPatchScalarField
(
const epsilonLowReWallFunctionFvPatchScalarField& ewfpsf
)
:
epsilonWallFunctionFvPatchScalarField(ewfpsf),
yPlusLam_(ewfpsf.yPlusLam_)
{}
Foam::epsilonLowReWallFunctionFvPatchScalarField::
epsilonLowReWallFunctionFvPatchScalarField
(
const epsilonLowReWallFunctionFvPatchScalarField& ewfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
epsilonWallFunctionFvPatchScalarField(ewfpsf, iF),
yPlusLam_(ewfpsf.yPlusLam_)
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
epsilonLowReWallFunctionFvPatchScalarField
);
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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::epsilonLowReWallFunctionFvPatchScalarField
Group
grpWallFunctions
Description
This boundary condition provides a turbulence dissipation wall function
condition for low- and high-Reynolds number turbulent flow cases.
The condition can be applied to wall boundaries, whereby it inserts near
wall epsilon values directly into the epsilon equation to act as a
constraint.
The model operates in two modes, based on the computed laminar-to-turbulent
switch-over y+ value derived from kappa and E.
Usage
\table
Property | Description | Required | Default value
Cmu | model coefficient | no | 0.09
kappa | Von Karman constant | no | 0.41
E | model coefficient | no | 9.8
\endtable
Example of the boundary condition specification:
\verbatim
<patchName>
{
type epsilonLowReWallFunction;
}
\endverbatim
See also
Foam::epsilonWallFunctionFvPatchScalarField
SourceFiles
epsilonLowReWallFunctionFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef epsilonLowReWallFunctionFvPatchScalarField_H
#define epsilonLowReWallFunctionFvPatchScalarField_H
#include "epsilonWallFunctionFvPatchScalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class epsilonLowReWallFunctionFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class epsilonLowReWallFunctionFvPatchScalarField
:
public epsilonWallFunctionFvPatchScalarField
{
protected:
// Protected data
//- Y+ at the edge of the laminar sublayer
scalar yPlusLam_;
// Protected Member Functions
//- Calculate the Y+ at the edge of the laminar sublayer
scalar yPlusLam(const scalar kappa, const scalar E);
//- Calculate the epsilon and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
);
public:
//- Runtime type information
TypeName("epsilonLowReWallFunction");
// Constructors
//- Construct from patch and internal field
epsilonLowReWallFunctionFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
epsilonLowReWallFunctionFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// epsilonLowReWallFunctionFvPatchScalarField
// onto a new patch
epsilonLowReWallFunctionFvPatchScalarField
(
const epsilonLowReWallFunctionFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
epsilonLowReWallFunctionFvPatchScalarField
(
const epsilonLowReWallFunctionFvPatchScalarField&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField>
(
new epsilonLowReWallFunctionFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
epsilonLowReWallFunctionFvPatchScalarField
(
const epsilonLowReWallFunctionFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new epsilonLowReWallFunctionFvPatchScalarField(*this, iF)
);
}
//- Destructor
virtual ~epsilonLowReWallFunctionFvPatchScalarField()
{}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "epsilonWallFunctionFvPatchScalarField.H"
#include "nutWallFunctionFvPatchScalarField.H"
#include "turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
......@@ -173,7 +174,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
scalarField& epsilon0
)
{
// accumulate all of the G and epsilon contributions
// Accumulate all of the G and epsilon contributions
forAll(cornerWeights_, patchi)
{
if (!cornerWeights_[patchi].empty())
......@@ -186,7 +187,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
}
}
// apply zero-gradient condition for epsilon
// Apply zero-gradient condition for epsilon
forAll(cornerWeights_, patchi)
{
if (!cornerWeights_[patchi].empty())
......@@ -201,48 +202,62 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
void Foam::epsilonWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const turbulenceModel& turbModel,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
scalarField& G0,
scalarField& epsilon0
)
{
const label patchi = patch.index();
const scalarField& y = turbulence.y()[patchi];
const scalarField& y = turbModel.y()[patchi];
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
const tmp<volScalarField> tk = turbulence.k();
const tmp<volScalarField> tk = turbModel.k();
const volScalarField& k = tk();
const tmp<scalarField> tnuw = turbulence.nu(patchi);
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const tmp<scalarField> tnutw = turbulence.nut(patchi);
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchi];
const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
const scalarField magGradUw(mag(Uw.snGrad()));
typedef DimensionedField<scalar, volMesh> FieldType;
const FieldType& G = db().lookupObject<FieldType>(turbModel.GName());
// Set epsilon and G
forAll(nutw, facei)
{
label celli = patch.faceCells()[facei];
const label celli = patch.faceCells()[facei];
scalar w = cornerWeights[facei];
const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
epsilon[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
const scalar w = cornerWeights[facei];
G[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(kappa_*y[facei]);
if (yPlus > yPlusLam_)
{
epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
G0[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(kappa_*y[facei]);
}
else
{
epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
G0[celli] += w*G[celli];
}
}
}
......@@ -260,6 +275,7 @@ epsilonWallFunctionFvPatchScalarField
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
epsilon_(),
initialised_(false),
......@@ -283,6 +299,7 @@ epsilonWallFunctionFvPatchScalarField
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
yPlusLam_(ptf.yPlusLam_),
G_(),
epsilon_(),
initialised_(false),
......@@ -305,6 +322,7 @@ epsilonWallFunctionFvPatchScalarField
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
epsilon_(),
initialised_(false),
......@@ -313,7 +331,7 @@ epsilonWallFunctionFvPatchScalarField
{