Program Listing for File LOTWWallModelFvPatchScalarField.C
↰ Return to documentation for file (wallModels/LOTWWallModelFvPatchScalarField.C)
/*---------------------------------------------------------------------------* \
License
This file is part of libWallModelledLES.
libWallModelledLES 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.
libWallModelledLES 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 libWallModelledLES.
If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "LOTWWallModelFvPatchScalarField.H"
#include "fvPatchFieldMapper.H"
#include "addToRunTimeSelectionTable.H"
#include "codeRules.H"
#include "scalarListIOList.H"
#include "helpers.H"
#include "LawOfTheWall.H"
#include "RootFinder.H"
#include "SingleCellSampler.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::LOTWWallModelFvPatchScalarField::writeLocalEntries(Ostream& os) const
{
wallModelFvPatchScalarField::writeLocalEntries(os);
rootFinder_->write(os);
law_->write(os);
sampler_->write(os);
}
Foam::tmp<Foam::scalarField>
Foam::LOTWWallModelFvPatchScalarField::calcNut() const
{
if (debug)
{
Info<< "Updating nut for patch " << patch().name() << nl;
}
const label patchi = patch().index();
tmp<scalarField> nuw = this->nu(patchi);
const scalarListIOList & wallGradU =
sampler_->db().lookupObject<scalarListIOList>("wallGradU");
scalarField magGradU(Helpers::mag(wallGradU));
return max
(
scalar(0),
(sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw)
);
}
Foam::tmp<Foam::scalarField>
Foam::LOTWWallModelFvPatchScalarField::
calcUTau(const scalarField & magGradU) const
{
const label patchi = patch().index();
const label patchSize = patch().size();
tmp<scalarField> tnuw = this->nu(patchi);
const scalarField& nuw = tnuw();
// Turbulent viscosity
const scalarField & nutw = *this;
// Computed uTau
tmp<scalarField> tuTau(new scalarField(patchSize, 0.0));
scalarField & uTau = tuTau.ref();
// Function to give to the root finder
std::function<scalar(scalar)> value;
std::function<scalar(scalar)> derivValue;
// Grab global uTau field
volScalarField & uTauField =
const_cast<volScalarField &>
(
db().lookupObject<volScalarField>("uTauPredicted")
);
const scalarField & uTauFieldBoundary = uTauField.boundaryFieldRef()[patchi];
// Compute uTau for each face
const scalarListIOList & sampledU =
sampler_().db().lookupObject<scalarListIOList>("U");
forAll(uTau, faceI)
{
// Starting guess using utau from the last timestep or an estimate
// using the current velocity gradient and nut value.
scalar ut = 0;
if (uTauFieldBoundary[faceI] > 0)
{
ut = uTauFieldBoundary[faceI];
}
else
{
ut = sqrt((nuw[faceI] + nutw[faceI]) * magGradU[faceI]);
}
if (ut > ROOTVSMALL)
{
scalar sampledUI = mag(vector(
sampledU[faceI][0],
sampledU[faceI][1],
sampledU[faceI][2])
);
scalar nuwI = nuw[faceI];
// Initial guess for lower bound, solution corresponding to nut = 0
// Since nut is strictly positive, we cannot predict a lower stress
scalar lowerBound = sqrt(nuw[faceI]*magGradU[faceI]);
// We consider u+ >= 0.025, which in a classical TBl corresponds to
// y+ = 0.025, so very very close to the wall.
scalar upperBound = sampledUI / 0.025;
std::function<scalar(scalar)> func =
[nuwI, faceI, this](const scalar & uTau)
{
return this->law_().value(
this->sampler_(),
faceI,
uTau,
nuwI
);
};
std::function<scalar(scalar)> deriv =
[nuwI, faceI, this](const scalar & uTau)
{
return this->law_().derivative(
this->sampler_(),
faceI,
uTau,
nuwI
);
};
// Set lower bound so that we bracket the root.
lowerBound = getLowerBound(lowerBound, upperBound, func);
// Supply the functions to the root finder
const_cast<RootFinder &>(rootFinder_()).setFunction(func);
const_cast<RootFinder &>(rootFinder_()).setDerivative(deriv);
std::pair<scalar, label> sol =
rootFinder_->root(ut, lowerBound, upperBound);
// Compute root to get uTau
uTau[faceI] = max(0.0, sol.first);
}
}
// Assign computed uTau to the boundary field of the global field
uTauField.boundaryFieldRef()[patchi] == uTau;
return tuTau;
}
Foam::scalar Foam::getLowerBound(scalar initial, scalar upperBound,
std::function<scalar(scalar)> func)
{
Foam::scalar lowerBound = initial;
for (int i=0; i<10; i++)
{
if (lowerBound >= upperBound || func(lowerBound)*func(upperBound) > 0)
{
lowerBound *= 0.1;
}
else
{
break;
}
}
return lowerBound;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::LOTWWallModelFvPatchScalarField::
LOTWWallModelFvPatchScalarField
(
const fvPatch & p,
const DimensionedField<scalar, volMesh> & iF
)
:
wallModelFvPatchScalarField(p, iF),
rootFinder_(),
law_(),
sampler_(nullptr)
{
if (debug)
{
Info<< "Constructing LOTWwallModelFvPatchScalarField (lotw1) "
<< "from fvPatch and DimensionedField for patch " << patch().name()
<< nl;
}
}
Foam::LOTWWallModelFvPatchScalarField::
LOTWWallModelFvPatchScalarField
(
const LOTWWallModelFvPatchScalarField & orig,
const fvPatch & p,
const DimensionedField<scalar, volMesh> & iF,
const fvPatchFieldMapper & mapper
)
:
wallModelFvPatchScalarField(orig, p, iF, mapper),
#ifdef FOAM_AUTOPTR_HAS_CLONE_METHOD
rootFinder_(orig.rootFinder_.clone()),
law_(orig.law_.clone()),
#else
rootFinder_(orig.rootFinder_, false),
law_(orig.law_, false),
#endif
sampler_(new SingleCellSampler(orig.sampler()))
{
if (debug)
{
Info<< "Constructing LOTWWallModelFvPatchScalarField (lotw2) "
<< "from copy, fvPatch, DimensionedField, and fvPatchFieldMapper"
<< " for patch " << patch().name() << nl;
}
law_->addFieldsToSampler(sampler());
}
Foam::LOTWWallModelFvPatchScalarField::
LOTWWallModelFvPatchScalarField
(
const fvPatch & p,
const DimensionedField<scalar, volMesh> & iF,
const dictionary & dict
)
:
wallModelFvPatchScalarField(p, iF, dict),
rootFinder_(RootFinder::New(dict.subDict("RootFinder"))),
law_(LawOfTheWall::New(dict.subDict("Law"))),
sampler_
(
new SingleCellSampler
(
p,
averagingTime(),
dict.lookupOrDefault<word>("interpolationType", "cell"),
dict.lookupOrDefault<word>("sampler", "Tree"),
dict.lookupOrDefault<word>("lengthScale", "CubeRootVol"),
dict.lookupOrDefault<bool>("hIsIndex", false)
)
)
{
if (debug)
{
Info<< "Constructing LOTWWallModelFvPatchScalarField (lotw3) "
<< "from fvPatch, DimensionedField, and dictionary for patch "
<< patch().name() << nl;
}
law_->addFieldsToSampler(sampler());
if (!db().found("solverIterations"))
{
db().store
(
new volVectorField
(
IOobject
(
"solverIterations",
db().time().timeName(),
db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
patch().boundaryMesh().mesh(),
vector(0,0,0),
dimless
)
);
}
if (!db().found("wallModellingInfo"))
{
db().store
(
new volTensorField
(
IOobject
(
"wallModellingInfo",
db().time().timeName(),
db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
patch().boundaryMesh().mesh(),
tensor(0,0,0,0,0,0,0,0,0),
dimless
)
);
}
}
#ifdef FOAM_FVPATCHFIELD_NO_COPY
#else
Foam::LOTWWallModelFvPatchScalarField::
LOTWWallModelFvPatchScalarField
(
const LOTWWallModelFvPatchScalarField & orig
)
:
wallModelFvPatchScalarField(orig),
#ifdef FOAM_AUTOPTR_HAS_CLONE_METHOD
rootFinder_(orig.rootFinder_.clone()),
law_(orig.law_.clone()),
#else
rootFinder_(orig.rootFinder_, false),
law_(orig.law_, false),
#endif
sampler_(new SingleCellSampler(orig.sampler_()))
{
if (debug)
{
Info<< "Constructing LOTWWallModelFvPatchScalarField (lotw4)"
<< "from copy for patch " << patch().name() << nl;
}
law_->addFieldsToSampler(sampler());
}
#endif
Foam::LOTWWallModelFvPatchScalarField::
LOTWWallModelFvPatchScalarField
(
const LOTWWallModelFvPatchScalarField & orig,
const DimensionedField<scalar, volMesh> & iF
)
:
wallModelFvPatchScalarField(orig, iF),
#ifdef FOAM_AUTOPTR_HAS_CLONE_METHOD
rootFinder_(orig.rootFinder_.clone()),
law_(orig.law_.clone()),
#else
rootFinder_(orig.rootFinder_, false),
law_(orig.law_, false),
#endif
sampler_(new SingleCellSampler(orig.sampler_()))
{
if (debug)
{
Info<< "Constructing LOTWModelFvPatchScalarField (lotw5) "
<< "from copy and DimensionedField for patch " << patch().name()
<< nl;
}
law_->addFieldsToSampler(sampler());
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::LOTWWallModelFvPatchScalarField::write(Ostream& os) const
{
wallModelFvPatchScalarField::write(os);
}
void Foam::LOTWWallModelFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
sampler().recomputeFields();
sampler().sample();
wallModelFvPatchScalarField::updateCoeffs();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#if !defined(DOXYGEN_SHOULD_SKIP_THIS)
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
LOTWWallModelFvPatchScalarField
);
}
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //