Files
openfoam/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
Andrew Heather 3b70a82bb7 BUG: Partial revert of commit fd87d0af1
omegaWallFunction
- re-instated behaviour when not using 'blended'
- turbulence generation always included when using 'blended'
- 'blended' now true by default

epsilonWallFunction
- re-instated low-Re switching
2017-12-12 08:48:38 +00:00

635 lines
15 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2017 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 "omegaWallFunctionFvPatchScalarField.H"
#include "nutWallFunctionFvPatchScalarField.H"
#include "turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "wallFvPatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
scalar omegaWallFunctionFvPatchScalarField::tolerance_ = 1e-5;
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void omegaWallFunctionFvPatchScalarField::checkType()
{
if (!isA<wallFvPatch>(patch()))
{
FatalErrorInFunction
<< "Invalid wall function specification" << nl
<< " Patch type for patch " << patch().name()
<< " must be wall" << nl
<< " Current patch type is " << patch().type() << nl << endl
<< abort(FatalError);
}
}
void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
{
os.writeEntry("Cmu", Cmu_);
os.writeEntry("kappa", kappa_);
os.writeEntry("E", E_);
os.writeEntry("beta1", beta1_);
os.writeEntry("blended", blended_);
}
void omegaWallFunctionFvPatchScalarField::setMaster()
{
if (master_ != -1)
{
return;
}
const volScalarField& omega =
static_cast<const volScalarField&>(this->internalField());
const volScalarField::Boundary& bf = omega.boundaryField();
label master = -1;
forAll(bf, patchi)
{
if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
{
omegaWallFunctionFvPatchScalarField& opf = omegaPatch(patchi);
if (master == -1)
{
master = patchi;
}
opf.master() = master;
}
}
}
void omegaWallFunctionFvPatchScalarField::createAveragingWeights()
{
const volScalarField& omega =
static_cast<const volScalarField&>(this->internalField());
const volScalarField::Boundary& bf = omega.boundaryField();
const fvMesh& mesh = omega.mesh();
if (initialised_ && !mesh.changing())
{
return;
}
volScalarField weights
(
IOobject
(
"weights",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
DynamicList<label> omegaPatches(bf.size());
forAll(bf, patchi)
{
if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchi]))
{
omegaPatches.append(patchi);
const labelUList& faceCells = bf[patchi].patch().faceCells();
forAll(faceCells, i)
{
label celli = faceCells[i];
weights[celli]++;
}
}
}
cornerWeights_.setSize(bf.size());
forAll(omegaPatches, i)
{
label patchi = omegaPatches[i];
const fvPatchScalarField& wf = weights.boundaryField()[patchi];
cornerWeights_[patchi] = 1.0/wf.patchInternalField();
}
G_.setSize(internalField().size(), 0.0);
omega_.setSize(internalField().size(), 0.0);
initialised_ = true;
}
omegaWallFunctionFvPatchScalarField&
omegaWallFunctionFvPatchScalarField::omegaPatch(const label patchi)
{
const volScalarField& omega =
static_cast<const volScalarField&>(this->internalField());
const volScalarField::Boundary& bf = omega.boundaryField();
const omegaWallFunctionFvPatchScalarField& opf =
refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchi]);
return const_cast<omegaWallFunctionFvPatchScalarField&>(opf);
}
void omegaWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
const turbulenceModel& turbModel,
scalarField& G0,
scalarField& omega0
)
{
// accumulate all of the G and omega contributions
forAll(cornerWeights_, patchi)
{
if (!cornerWeights_[patchi].empty())
{
omegaWallFunctionFvPatchScalarField& opf = omegaPatch(patchi);
const List<scalar>& w = cornerWeights_[patchi];
opf.calculate(turbModel, w, opf.patch(), G0, omega0);
}
}
// apply zero-gradient condition for omega
forAll(cornerWeights_, patchi)
{
if (!cornerWeights_[patchi].empty())
{
omegaWallFunctionFvPatchScalarField& opf = omegaPatch(patchi);
opf == scalarField(omega0, opf.patch().faceCells());
}
}
}
void omegaWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbModel,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G0,
scalarField& omega0
)
{
const label patchi = patch.index();
const scalarField& y = turbModel.y()[patchi];
const scalar Cmu25 = pow025(Cmu_);
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()));
// Set omega and G
forAll(nutw, facei)
{
const label celli = patch.faceCells()[facei];
const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
const scalar w = cornerWeights[facei];
const scalar omegaVis = 6*nuw[facei]/(beta1_*sqr(y[facei]));
const scalar omegaLog = sqrt(k[celli])/(Cmu25*kappa_*y[facei]);
// Switching between the laminar sub-layer and the log-region rather
// than blending has been found to provide more accurate results over a
// range of near-wall y+.
//
// For backward-compatibility the blending method is provided as an
// option
// Generation contribution is included using the blended option, or
// when using the switching option if operating in the laminar
// sub-layer
bool includeG = true;
if (blended_)
{
omega0[celli] += w*sqrt(sqr(omegaVis) + sqr(omegaLog));
}
else
{
if (yPlus > yPlusLam_)
{
omega0[celli] += w*omegaLog;
}
else
{
omega0[celli] += w*omegaVis;
includeG = false;
}
}
if (includeG)
{
G0[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(kappa_*y[facei]);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
beta1_(0.075),
blended_(true),
yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
beta1_(ptf.beta1_),
blended_(ptf.blended_),
yPlusLam_(ptf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
blended_(dict.lookupOrDefault<Switch>("blended", true)),
yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
// apply zero-gradient condition on start-up
this->operator==(patchInternalField());
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& owfpsf
)
:
fixedValueFvPatchField<scalar>(owfpsf),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
blended_(owfpsf.blended_),
yPlusLam_(owfpsf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& owfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(owfpsf, iF),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
blended_(owfpsf.blended_),
yPlusLam_(owfpsf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField& omegaWallFunctionFvPatchScalarField::G(bool init)
{
if (patch().index() == master_)
{
if (init)
{
G_ = 0.0;
}
return G_;
}
return omegaPatch(master_).G();
}
scalarField& omegaWallFunctionFvPatchScalarField::omega(bool init)
{
if (patch().index() == master_)
{
if (init)
{
omega_ = 0.0;
}
return omega_;
}
return omegaPatch(master_).omega(init);
}
void omegaWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
(
IOobject::groupName
(
turbulenceModel::propertiesName,
internalField().group()
)
);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbModel, G(true), omega(true));
}
const scalarField& G0 = this->G();
const scalarField& omega0 = this->omega();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G = db().lookupObjectRef<FieldType>(turbModel.GName());
FieldType& omega = const_cast<FieldType&>(internalField());
forAll(*this, facei)
{
label celli = patch().faceCells()[facei];
G[celli] = G0[celli];
omega[celli] = omega0[celli];
}
fvPatchField<scalar>::updateCoeffs();
}
void omegaWallFunctionFvPatchScalarField::updateWeightedCoeffs
(
const scalarField& weights
)
{
if (updated())
{
return;
}
const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
(
IOobject::groupName
(
turbulenceModel::propertiesName,
internalField().group()
)
);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbModel, G(true), omega(true));
}
const scalarField& G0 = this->G();
const scalarField& omega0 = this->omega();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G = db().lookupObjectRef<FieldType>(turbModel.GName());
FieldType& omega = const_cast<FieldType&>(internalField());
scalarField& omegaf = *this;
// only set the values if the weights are > tolerance
forAll(weights, facei)
{
scalar w = weights[facei];
if (w > tolerance_)
{
label celli = patch().faceCells()[facei];
G[celli] = (1.0 - w)*G[celli] + w*G0[celli];
omega[celli] = (1.0 - w)*omega[celli] + w*omega0[celli];
omegaf[facei] = omega[celli];
}
}
fvPatchField<scalar>::updateCoeffs();
}
void omegaWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix
)
{
if (manipulatedMatrix())
{
return;
}
matrix.setValues(patch().faceCells(), patchInternalField());
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix,
const Field<scalar>& weights
)
{
if (manipulatedMatrix())
{
return;
}
DynamicList<label> constraintCells(weights.size());
DynamicList<scalar> constraintomega(weights.size());
const labelUList& faceCells = patch().faceCells();
const DimensionedField<scalar, volMesh>& omega = internalField();
label nConstrainedCells = 0;
forAll(weights, facei)
{
// only set the values if the weights are > tolerance
if (weights[facei] > tolerance_)
{
++nConstrainedCells;
label celli = faceCells[facei];
constraintCells.append(celli);
constraintomega.append(omega[celli]);
}
}
if (debug)
{
Pout<< "Patch: " << patch().name()
<< ": number of constrained cells = " << nConstrainedCells
<< " out of " << patch().size()
<< endl;
}
matrix.setValues
(
constraintCells,
scalarField(constraintomega.xfer())
);
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
{
writeLocalEntries(os);
fixedValueFvPatchField<scalar>::write(os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchScalarField,
omegaWallFunctionFvPatchScalarField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //