BUG: nutWallFunction: prevent casting errors (fixes #1457, #1966)

Previously, a nutWallFunctionFvPatchScalarField ref should be
created in epsilon, k, and omega wall functions to fetch various
common wall-function coefficients necessary to carry out and complete
local operations inside these wall functions.

However, this arrangement required the use of a nut wall function,
even when unnecessary, when any of non-nut wall functions are being used.
Therefore, some users had been redundantly restrained and
obstructed with rather obscure casting-error messages.

Also, the wall-function coefficients Cmu, kappa and E have been obtained
from the specified nutWallFunction in order to ensure that each patch
possesses the same set of values for these coefficients.

Although the motivation sounds reasonable, it has also been putting redundant
restraints on users and disregarding the specifics of each wall-function.
For example, the variation of epsilon in near-wall regions is usually very
steep and non-monotonic specific - an expert user may therefore want to use
an epsilon-specific coefficient, and this was not allowed by the previous
arrangement.

This commit introduces a new class (i.e. wallFunctionCoefficients) comprising
all common wall-function coefficients and yPlus calculations.
This commit is contained in:
Kutalmis Bercin
2022-05-03 13:57:21 +01:00
parent 486be34631
commit c7836f9d83
23 changed files with 380 additions and 158 deletions

View File

@ -32,6 +32,8 @@ derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
/* Wall function BCs */
wallFunctions = derivedFvPatchFields/wallFunctions
$(wallFunctions)/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.C
nutWallFunctions = $(wallFunctions)/nutWallFunctions
$(nutWallFunctions)/nutWallFunction/nutWallFunctionFvPatchScalarField.C

View File

@ -202,8 +202,8 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
{
const label patchi = patch.index();
const nutWallFunctionFvPatchScalarField& nutw =
nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const scalarField& y = turbModel.y()[patchi];
@ -217,8 +217,10 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(nutw.Cmu());
const scalar Cmu75 = pow(nutw.Cmu(), 0.75);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75);
const scalar kappa = wallCoeffs_.kappa();
const scalar yPlusLam = wallCoeffs_.yPlusLam();
// Set epsilon and G
forAll(nutw, facei)
@ -236,13 +238,13 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
// Contribution from the inertial sublayer
const scalar epsilonLog =
w*Cmu75*pow(k[celli], 1.5)/(nutw.kappa()*y[facei]);
w*Cmu75*pow(k[celli], 1.5)/(kappa*y[facei]);
switch (blending_)
{
case blendingType::STEPWISE:
{
if (lowReCorrection_ && yPlus < nutw.yPlusLam())
if (lowReCorrection_ && yPlus < yPlusLam)
{
epsilonBlended = epsilonVis;
}
@ -285,14 +287,14 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
epsilon0[celli] += epsilonBlended;
if (!(lowReCorrection_ && yPlus < nutw.yPlusLam()))
if (!(lowReCorrection_ && yPlus < yPlusLam))
{
G0[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(nutw.kappa()*y[facei]);
/(kappa*y[facei]);
}
}
}
@ -310,6 +312,8 @@ void Foam::epsilonWallFunctionFvPatchScalarField::writeLocalEntries
{
os.writeEntry("n", n_);
}
wallCoeffs_.writeEntries(os);
}
@ -328,6 +332,7 @@ epsilonWallFunctionFvPatchScalarField
lowReCorrection_(false),
initialised_(false),
master_(-1),
wallCoeffs_(),
G_(),
epsilon_(),
cornerWeights_()
@ -349,6 +354,7 @@ epsilonWallFunctionFvPatchScalarField
lowReCorrection_(ptf.lowReCorrection_),
initialised_(false),
master_(-1),
wallCoeffs_(ptf.wallCoeffs_),
G_(),
epsilon_(),
cornerWeights_()
@ -385,6 +391,7 @@ epsilonWallFunctionFvPatchScalarField
lowReCorrection_(dict.getOrDefault("lowReCorrection", false)),
initialised_(false),
master_(-1),
wallCoeffs_(dict),
G_(),
epsilon_(),
cornerWeights_()
@ -406,6 +413,7 @@ epsilonWallFunctionFvPatchScalarField
lowReCorrection_(ewfpsf.lowReCorrection_),
initialised_(false),
master_(-1),
wallCoeffs_(ewfpsf.wallCoeffs_),
G_(),
epsilon_(),
cornerWeights_()
@ -425,6 +433,7 @@ epsilonWallFunctionFvPatchScalarField
lowReCorrection_(ewfpsf.lowReCorrection_),
initialised_(false),
master_(-1),
wallCoeffs_(ewfpsf.wallCoeffs_),
G_(),
epsilon_(),
cornerWeights_()

View File

@ -173,6 +173,7 @@ SourceFiles
#define epsilonWallFunctionFvPatchScalarField_H
#include "fixedValueFvPatchField.H"
#include "wallFunctionCoefficients.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -230,6 +231,9 @@ protected:
//- Master patch ID
label master_;
//- Wall-function coefficients
wallFunctionCoefficients wallCoeffs_;
//- Local copy of turbulence G field
scalarField G_;

View File

@ -42,6 +42,7 @@ void Foam::kLowReWallFunctionFvPatchScalarField::writeLocalEntries
os.writeEntryIfDifferent<scalar>("Ck", -0.416, Ck_);
os.writeEntryIfDifferent<scalar>("Bk", 8.366, Bk_);
os.writeEntryIfDifferent<scalar>("C", 11.0, C_);
wallCoeffs_.writeEntries(os);
}
@ -57,7 +58,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField
Ceps2_(1.9),
Ck_(-0.416),
Bk_(8.366),
C_(11.0)
C_(11.0),
wallCoeffs_()
{}
@ -73,7 +75,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField
Ceps2_(ptf.Ceps2_),
Ck_(ptf.Ck_),
Bk_(ptf.Bk_),
C_(ptf.C_)
C_(ptf.C_),
wallCoeffs_(ptf.wallCoeffs_)
{}
@ -96,7 +99,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField
),
Ck_(dict.getOrDefault<scalar>("Ck", -0.416)),
Bk_(dict.getOrDefault<scalar>("Bk", 8.366)),
C_(dict.getOrDefault<scalar>("C", 11.0))
C_(dict.getOrDefault<scalar>("C", 11.0)),
wallCoeffs_(dict)
{}
@ -109,7 +113,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField
Ceps2_(kwfpsf.Ceps2_),
Ck_(kwfpsf.Ck_),
Bk_(kwfpsf.Bk_),
C_(kwfpsf.C_)
C_(kwfpsf.C_),
wallCoeffs_(kwfpsf.wallCoeffs_)
{}
@ -123,7 +128,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField
Ceps2_(kwfpsf.Ceps2_),
Ck_(kwfpsf.Ck_),
Bk_(kwfpsf.Bk_),
C_(kwfpsf.C_)
C_(kwfpsf.C_),
wallCoeffs_(kwfpsf.wallCoeffs_)
{}
@ -147,9 +153,6 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs()
)
);
const nutWallFunctionFvPatchScalarField& nutw =
nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
const scalarField& y = turbModel.y()[patchi];
const tmp<scalarField> tnuw = turbModel.nu(patchi);
@ -158,7 +161,9 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs()
const tmp<volScalarField> tk = turbModel.k();
const volScalarField& k = tk();
const scalar Cmu25 = pow025(nutw.Cmu());
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const scalar yPlusLam = wallCoeffs_.yPlusLam();
scalarField& kw = *this;
@ -169,9 +174,9 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs()
const scalar uTau = Cmu25*sqrt(k[celli]);
const scalar yPlus = uTau*y[facei]/nuw[facei];
if (yPlus > nutw.yPlusLam())
if (yPlus > yPlusLam)
{
kw[facei] = Ck_/nutw.kappa()*log(yPlus) + Bk_;
kw[facei] = Ck_/kappa*log(yPlus) + Bk_;
}
else
{

View File

@ -101,6 +101,7 @@ SourceFiles
#define kLowReWallFunctionFvPatchScalarField_H
#include "fixedValueFvPatchField.H"
#include "wallFunctionCoefficients.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -131,6 +132,9 @@ protected:
//- C coefficient
scalar C_;
//- Wall-function coefficients
wallFunctionCoefficients wallCoeffs_;
// Protected Member Functions

View File

@ -92,6 +92,9 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau
const scalarField& nutw = *this;
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
forAll(uTaup, facei)
{
scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]);
@ -104,7 +107,7 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau
const scalar yPlus = y[facei]*ut/nuw[facei];
const scalar uTauVis = magUp[facei]/yPlus;
const scalar uTauLog =
kappa_*magUp[facei]/log(max(E_*yPlus, 1 + 1e-4));
kappa*magUp[facei]/log(max(E*yPlus, 1 + 1e-4));
const scalar utNew =
pow(pow(uTauVis, n_) + pow(uTauLog, n_), 1.0/n_);

View File

@ -56,6 +56,8 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const
// The flow velocity at the adjacent cell centre
const scalarField magUp(mag(Uw.patchInternalField() - Uw));
const scalar yPlusLam = wallCoeffs_.yPlusLam();
tmp<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus.ref();
@ -64,7 +66,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const
forAll(yPlus, facei)
{
if (yPlusLam_ < yPlus[facei])
if (yPlusLam < yPlus[facei])
{
const scalar Re = magUp[facei]*y[facei]/nuw[facei] + ROOTVSMALL;
nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
@ -95,6 +97,10 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
const scalar yPlusLam = wallCoeffs_.yPlusLam();
tmp<scalarField> tyPlus(new scalarField(patch().size(), Zero));
scalarField& yPlus = tyPlus.ref();
@ -114,9 +120,9 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
{
const scalar magUpara = magUp[facei];
const scalar Re = magUpara*y[facei]/nuw[facei];
const scalar kappaRe = kappa_*Re;
const scalar kappaRe = kappa*Re;
scalar yp = yPlusLam_;
scalar yp = yPlusLam;
const scalar ryPlusLam = 1.0/yp;
int iter = 0;
@ -155,7 +161,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
(c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
}
scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime;
scalar denom = 1.0 + log(E*yp) - G - yPlusGPrime;
if (mag(denom) > VSMALL)
{
yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
@ -178,9 +184,9 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
{
const scalar magUpara = magUp[facei];
const scalar Re = magUpara*y[facei]/nuw[facei];
const scalar kappaRe = kappa_*Re;
const scalar kappaRe = kappa*Re;
scalar yp = yPlusLam_;
scalar yp = yPlusLam;
const scalar ryPlusLam = 1.0/yp;
int iter = 0;
@ -189,7 +195,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
do
{
yPlusLast = yp;
yp = (kappaRe + yp)/(1.0 + log(E_*yp));
yp = (kappaRe + yp)/(1.0 + log(E*yp));
}
while

View File

@ -126,6 +126,9 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
const scalarField& nutw = *this;
tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
@ -146,18 +149,18 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
do
{
scalar kUu = min(kappa_*magUp[facei]/ut, 50);
scalar kUu = min(kappa*magUp[facei]/ut, 50);
scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu);
scalar f =
- ut*y[facei]/nuw[facei]
+ magUp[facei]/ut
+ 1/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu));
+ 1/E*(fkUu - 1.0/6.0*kUu*sqr(kUu));
scalar df =
y[facei]/nuw[facei]
+ magUp[facei]/sqr(ut)
+ 1/E_*kUu*fkUu/ut;
+ 1/E*kUu*fkUu/ut;
scalar uTauNew = ut + f/df;
err[facei] = mag((ut - uTauNew)/ut);

View File

@ -53,6 +53,9 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
tmp<scalarField> tyPlus = calcYPlus(magUp);
const scalarField& yPlus = tyPlus();
@ -67,7 +70,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
// Inertial sublayer contribution
const scalar nutLog =
nuw[facei]
*(yPlus[facei]*kappa_/log(max(E_*yPlus[facei], 1 + 1e-4)) - 1.0);
*(yPlus[facei]*kappa/log(max(E*yPlus[facei], 1 + 1e-4)) - 1.0);
nutw[facei] = blend(nutVis, nutLog, yPlus[facei]);
}
@ -96,14 +99,18 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
const scalar yPlusLam = wallCoeffs_.yPlusLam();
tmp<scalarField> tyPlus(new scalarField(patch().size(), Zero));
scalarField& yPlus = tyPlus.ref();
forAll(yPlus, facei)
{
const scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei];
const scalar kappaRe = kappa*magUp[facei]*y[facei]/nuw[facei];
scalar yp = yPlusLam_;
scalar yp = yPlusLam;
const scalar ryPlusLam = 1.0/yp;
int iter = 0;
@ -112,7 +119,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus
do
{
yPlusLast = yp;
yp = (kappaRe + yp)/(1.0 + log(E_*yp));
yp = (kappaRe + yp)/(1.0 + log(E*yp));
} while (mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 );
@ -219,12 +226,14 @@ Foam::nutUWallFunctionFvPatchScalarField::yPlus() const
const scalarField magUp(mag(Uw.patchInternalField() - Uw));
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar yPlusLam = wallCoeffs_.yPlusLam();
tmp<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus.ref();
forAll(yPlus, facei)
{
if (yPlusLam_ > yPlus[facei])
if (yPlusLam > yPlus[facei])
{
// viscous sublayer
yPlus[facei] =

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016, 2019 OpenFOAM Foundation
Copyright (C) 2019-2020 OpenCFD Ltd.
Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -92,10 +92,8 @@ void Foam::nutWallFunctionFvPatchScalarField::writeLocalEntries
Ostream& os
) const
{
os.writeEntry("Cmu", Cmu_);
os.writeEntry("kappa", kappa_);
os.writeEntry("E", E_);
os.writeEntryIfDifferent<word>("U", word::null, UName_);
wallCoeffs_.writeEntries(os);
}
@ -111,10 +109,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
blending_(blendingType::STEPWISE),
n_(4.0),
UName_(word::null),
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
yPlusLam_(yPlusLam(kappa_, E_))
wallCoeffs_()
{
checkType();
}
@ -132,10 +127,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
blending_(ptf.blending_),
n_(ptf.n_),
UName_(ptf.UName_),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
yPlusLam_(ptf.yPlusLam_)
wallCoeffs_(ptf.wallCoeffs_)
{
checkType();
}
@ -168,13 +160,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
)
),
UName_(dict.getOrDefault<word>("U", word::null)),
Cmu_(dict.getOrDefault<scalar>("Cmu", 0.09)),
kappa_
(
dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL))
),
E_(dict.getCheckOrDefault<scalar>("E", 9.8, scalarMinMax::ge(SMALL))),
yPlusLam_(yPlusLam(kappa_, E_))
wallCoeffs_(dict)
{
checkType();
}
@ -189,10 +175,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
blending_(wfpsf.blending_),
n_(wfpsf.n_),
UName_(wfpsf.UName_),
Cmu_(wfpsf.Cmu_),
kappa_(wfpsf.kappa_),
E_(wfpsf.E_),
yPlusLam_(wfpsf.yPlusLam_)
wallCoeffs_(wfpsf.wallCoeffs_)
{
checkType();
}
@ -208,10 +191,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField
blending_(wfpsf.blending_),
n_(wfpsf.n_),
UName_(wfpsf.UName_),
Cmu_(wfpsf.Cmu_),
kappa_(wfpsf.kappa_),
E_(wfpsf.E_),
yPlusLam_(wfpsf.yPlusLam_)
wallCoeffs_(wfpsf.wallCoeffs_)
{
checkType();
}
@ -235,23 +215,6 @@ Foam::nutWallFunctionFvPatchScalarField::nutw
}
Foam::scalar Foam::nutWallFunctionFvPatchScalarField::yPlusLam
(
const scalar kappa,
const scalar E
)
{
scalar ypl = 11.0;
for (label i = 0; i < 10; ++i)
{
ypl = log(max(E*ypl, 1.0))/kappa;
}
return ypl;
}
Foam::scalar Foam::nutWallFunctionFvPatchScalarField::blend
(
const scalar nutVis,
@ -265,7 +228,7 @@ Foam::scalar Foam::nutWallFunctionFvPatchScalarField::blend
{
case blendingType::STEPWISE:
{
if (yPlus > yPlusLam_)
if (yPlus > wallCoeffs_.yPlusLam())
{
nutw = nutLog;
}
@ -310,12 +273,6 @@ Foam::scalar Foam::nutWallFunctionFvPatchScalarField::blend
}
Foam::scalar Foam::nutWallFunctionFvPatchScalarField::yPlusLam() const
{
return yPlusLam_;
}
void Foam::nutWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())

View File

@ -167,6 +167,7 @@ SourceFiles
#define nutWallFunctionFvPatchScalarField_H
#include "fixedValueFvPatchFields.H"
#include "wallFunctionCoefficients.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -214,18 +215,8 @@ protected:
// retrieved from the turbulence model
word UName_;
//- Empirical model coefficient
scalar Cmu_;
//- von Kármán constant
scalar kappa_;
//- Wall roughness parameter
scalar E_;
//- Estimated y+ value at the intersection
//- of the viscous and inertial sublayers
scalar yPlusLam_;
//- Wall-function coefficients
wallFunctionCoefficients wallCoeffs_;
// Protected Member Functions
@ -296,24 +287,6 @@ public:
// Member Functions
//- Return Cmu
scalar Cmu() const
{
return Cmu_;
}
//- Return kappa
scalar kappa() const
{
return kappa_;
}
//- Return E
scalar E() const
{
return E_;
}
//- Return the nut patchField for the given wall patch
static const nutWallFunctionFvPatchScalarField& nutw
(
@ -321,12 +294,6 @@ public:
const label patchi
);
//- Estimate the y+ at the intersection of the two sublayers
static scalar yPlusLam(const scalar kappa, const scalar E);
//- Return the estimated y+ at the two-sublayer intersection
scalar yPlusLam() const;
//- Return the blended nut according to the chosen blending treatment
scalar blend
(

View File

@ -76,7 +76,9 @@ calcNut() const
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
tmp<scalarField> tnutw(new scalarField(*this));
scalarField& nutw = tnutw.ref();
@ -89,7 +91,7 @@ calcNut() const
const scalar yPlus = uStar*y[facei]/nuw[facei];
const scalar KsPlus = uStar*Ks_[facei]/nuw[facei];
scalar Edash = E_;
scalar Edash = E;
if (2.25 < KsPlus)
{
Edash /= fnRough(KsPlus, Cs_[facei]);
@ -105,7 +107,7 @@ calcNut() const
min
(
nuw[facei]
*(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1),
*(yPlus*kappa/log(max(Edash*yPlus, 1+1e-4)) - 1),
2*limitingNutw
), 0.5*limitingNutw
);

View File

@ -55,7 +55,9 @@ calcNut() const
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const scalar E = wallCoeffs_.E();
tmp<scalarField> tnutw(new scalarField(patch().size(), Zero));
scalarField& nutw = tnutw.ref();
@ -71,7 +73,7 @@ calcNut() const
// Inertial sublayer contribution
const scalar nutLog =
nuw[facei]*(yPlus*kappa_/log(max(E_*yPlus, 1 + 1e-4)) - 1.0);
nuw[facei]*(yPlus*kappa/log(max(E*yPlus, 1 + 1e-4)) - 1.0);
nutw[facei] = blend(nutVis, nutLog, yPlus);
}
@ -180,7 +182,8 @@ yPlus() const
const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar yPlusLam = wallCoeffs_.yPlusLam();
auto tyPlus = tmp<scalarField>::New(patch().size(), Zero);
auto& yPlus = tyPlus.ref();
@ -190,7 +193,7 @@ yPlus() const
// inertial sublayer
yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nuw[facei];
if (yPlusLam_ > yPlus[facei])
if (yPlusLam > yPlus[facei])
{
// viscous sublayer
yPlus[facei] =

View File

@ -202,8 +202,8 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate
{
const label patchi = patch.index();
const nutWallFunctionFvPatchScalarField& nutw =
nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const scalarField& y = turbModel.y()[patchi];
@ -217,7 +217,9 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(nutw.Cmu());
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const scalar yPlusLam = wallCoeffs_.yPlusLam();
// Set omega and G
forAll(nutw, facei)
@ -230,13 +232,13 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate
const scalar omegaVis = 6.0*nuw[facei]/(beta1_*sqr(y[facei]));
// Contribution from the inertial sublayer
const scalar omegaLog = sqrt(k[celli])/(Cmu25*nutw.kappa()*y[facei]);
const scalar omegaLog = sqrt(k[celli])/(Cmu25*kappa*y[facei]);
switch (blending_)
{
case blendingType::STEPWISE:
{
if (yPlus > nutw.yPlusLam())
if (yPlus > yPlusLam)
{
omega0[celli] += w*omegaLog;
}
@ -296,14 +298,14 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate
}
}
if (!(blending_ == blendingType::STEPWISE) || yPlus > nutw.yPlusLam())
if (!(blending_ == blendingType::STEPWISE) || yPlus > yPlusLam)
{
G0[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(nutw.kappa()*y[facei]);
/(kappa*y[facei]);
}
}
}
@ -321,6 +323,8 @@ void Foam::omegaWallFunctionFvPatchScalarField::writeLocalEntries
{
os.writeEntry("n", n_);
}
wallCoeffs_.writeEntries(os);
}
@ -338,6 +342,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
initialised_(false),
master_(-1),
beta1_(0.075),
wallCoeffs_(),
G_(),
omega_(),
cornerWeights_()
@ -358,6 +363,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
initialised_(false),
master_(-1),
beta1_(ptf.beta1_),
wallCoeffs_(ptf.wallCoeffs_),
G_(),
omega_(),
cornerWeights_()
@ -393,6 +399,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
initialised_(false),
master_(-1),
beta1_(dict.getOrDefault<scalar>("beta1", 0.075)),
wallCoeffs_(dict),
G_(),
omega_(),
cornerWeights_()
@ -436,6 +443,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
initialised_(false),
master_(-1),
beta1_(owfpsf.beta1_),
wallCoeffs_(owfpsf.wallCoeffs_),
G_(),
omega_(),
cornerWeights_()
@ -454,6 +462,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
initialised_(false),
master_(-1),
beta1_(owfpsf.beta1_),
wallCoeffs_(owfpsf.wallCoeffs_),
G_(),
omega_(),
cornerWeights_()

View File

@ -201,6 +201,7 @@ SourceFiles
#define omegaWallFunctionFvPatchScalarField_H
#include "fixedValueFvPatchField.H"
#include "wallFunctionCoefficients.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -264,6 +265,9 @@ protected:
//- beta1 coefficient
scalar beta1_;
//- Wall-function coefficients
wallFunctionCoefficients wallCoeffs_;
//- Local copy of turbulence G field
scalarField G_;

View File

@ -0,0 +1,96 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 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 "wallFunctionCoefficients.H"
#include "dictionary.H"
#include "MinMax.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
//- Estimate the y+ at the intersection of the two sublayers
static scalar calcYPlusLam
(
const scalar kappa,
const scalar E
)
{
scalar ypl = 11;
for (int iter = 0; iter < 10; ++iter)
{
ypl = log(max(E*ypl, scalar(1)))/kappa;
}
return ypl;
}
} // End namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::wallFunctionCoefficients::wallFunctionCoefficients()
:
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
yPlusLam_(calcYPlusLam(kappa_, E_))
{}
Foam::wallFunctionCoefficients::wallFunctionCoefficients
(
const dictionary& dict
)
:
Cmu_(dict.getOrDefault<scalar>("Cmu", 0.09)),
kappa_
(
dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL))
),
E_(dict.getCheckOrDefault<scalar>("E", 9.8, scalarMinMax::ge(SMALL))),
yPlusLam_(calcYPlusLam(kappa_, E_))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::wallFunctionCoefficients::writeEntries
(
Ostream& os
) const
{
os.writeEntryIfDifferent<scalar>("Cmu", 0.09, Cmu_);
os.writeEntryIfDifferent<scalar>("kappa", 0.41, kappa_);
os.writeEntryIfDifferent<scalar>("E", 9.8, E_);
}
// ************************************************************************* //

View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2022 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::wallFunctionCoefficients
Description
Class to describe the wall-function coefficients being
used in the wall function boundary conditions.
SourceFiles
wallFunctionCoefficients.C
\*---------------------------------------------------------------------------*/
#ifndef Foam_wallFunctionCoefficients_H
#define Foam_wallFunctionCoefficients_H
#include "scalarFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward Declarations
class dictionary;
class Ostream;
/*---------------------------------------------------------------------------*\
Class wallFunctionCoefficients Declaration
\*---------------------------------------------------------------------------*/
class wallFunctionCoefficients
{
// Private Data
//- Empirical model coefficient
scalar Cmu_;
//- von Karman constant
scalar kappa_;
//- Wall roughness parameter
scalar E_;
//- Estimated y+ value at the intersection
//- of the viscous and inertial sublayers
scalar yPlusLam_;
public:
// Constructors
//- Construct with default coefficients
wallFunctionCoefficients();
//- Construct from dictionary
explicit wallFunctionCoefficients(const dictionary& dict);
// Member Functions
// Access
//- Return the object: Cmu
scalar Cmu() const noexcept
{
return Cmu_;
}
//- Return the object: kappa
scalar kappa() const noexcept
{
return kappa_;
}
//- Return the object: E
scalar E() const noexcept
{
return E_;
}
//- Return the object: yPlusLam
scalar yPlusLam() const noexcept
{
return yPlusLam_;
}
// I-O
//- Write wall-function coefficients as dictionary entries
void writeEntries(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -45,8 +45,8 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate
{
const label patchi = patch.index();
const nutWallFunctionFvPatchScalarField& nutw =
nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const scalarField& y = turbModel.y()[patchi];
@ -60,8 +60,10 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(nutw.Cmu());
const scalar Cmu75 = pow(nutw.Cmu(), 0.75);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75);
const scalar kappa = wallCoeffs_.kappa();
const scalar yPlusLam = wallCoeffs_.yPlusLam();
const scalar t = db().time().timeOutputValue();
const scalarField z0(z0_->value(t));
@ -92,16 +94,16 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate
// (PGVB:Eq. 7, RH:Eq. 8)
scalar epsilonc =
w*Cmu75*pow(k[celli], 1.5)/(nutw.kappa()*(y[facei] + z0[facei]));
w*Cmu75*pow(k[celli], 1.5)/(kappa*(y[facei] + z0[facei]));
scalar Gc =
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(nutw.kappa()*(y[facei] + z0[facei]));
/(kappa*(y[facei] + z0[facei]));
if (lowReCorrection_ && yPlus < nutw.yPlusLam())
if (lowReCorrection_ && yPlus < yPlusLam)
{
epsilonc = w*2.0*k[celli]*nuw[facei]/sqr(y[facei] + z0[facei]);
Gc = 0;
@ -125,6 +127,8 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::writeLocalEntries
{
z0_->writeData(os);
}
wallCoeffs_.writeEntries(os);
}

View File

@ -62,6 +62,8 @@ tmp<scalarField> atmNutUWallFunctionFvPatchScalarField::calcNut() const
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar kappa = wallCoeffs_.kappa();
const scalar t = db().time().timeOutputValue();
const scalarField z0(z0_->value(t));
@ -84,7 +86,7 @@ tmp<scalarField> atmNutUWallFunctionFvPatchScalarField::calcNut() const
forAll(nutw, facei)
{
const scalar Edash = (y[facei] + z0[facei])/(z0[facei] + 1e-4);
const scalar uStar = magUpn[facei]*kappa_/log(max(Edash, 1.0 + 1e-4));
const scalar uStar = magUpn[facei]*kappa/log(max(Edash, 1.0 + 1e-4));
nutw[facei] = sqr(uStar)/max(magUpn[facei], 1e-6)*y[facei] - nuw[facei];
}

View File

@ -64,7 +64,8 @@ tmp<scalarField> atmNutWallFunctionFvPatchScalarField::calcNut() const
auto tnutw = tmp<scalarField>::New(*this);
auto& nutw = tnutw.ref();
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
const scalarField magUp(mag(Uw.patchInternalField() - Uw));
@ -95,7 +96,7 @@ tmp<scalarField> atmNutWallFunctionFvPatchScalarField::calcNut() const
const scalar Edash = (y[facei] + z0[facei])/(z0[facei] + z0Min_);
// (RH:Eq. 6)
const scalar uStarU = magUp[facei]*kappa_/log(max(Edash, 1 + SMALL));
const scalar uStarU = magUp[facei]*kappa/log(max(Edash, 1 + SMALL));
// (RH:Eq. 7)
const scalar uStarK = Cmu25*Foam::sqrt(k[celli]);

View File

@ -65,7 +65,8 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
auto tnutw = tmp<scalarField>::New(*this);
auto& nutw = tnutw.ref();
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const scalar t = db().time().timeOutputValue();
const scalarField z0(z0_->value(t));
@ -94,7 +95,7 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
const scalar yPlus = uStar*y[facei]/nuw[facei];
const scalar Edash = (y[facei] + z0[facei])/z0[facei];
nutw[facei] = nuw[facei]*(yPlus*kappa_/log(max(Edash, 1 + 1e-4)) - 1);
nutw[facei] = nuw[facei]*(yPlus*kappa/log(max(Edash, 1 + 1e-4)) - 1);
}
if (boundNut_)

View File

@ -45,8 +45,8 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate
{
const label patchi = patch.index();
const nutWallFunctionFvPatchScalarField& nutw =
nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
const tmp<scalarField> tnutw = turbModel.nut(patchi);
const scalarField& nutw = tnutw();
const scalarField& y = turbModel.y()[patchi];
@ -60,7 +60,8 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(nutw.Cmu());
const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
const scalar kappa = wallCoeffs_.kappa();
const scalar t = db().time().timeOutputValue();
const scalarField z0(z0_->value(t));
@ -88,14 +89,14 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate
const scalar w = cornerWeights[facei];
omega0[celli] +=
w*sqrt(k[celli])/(Cmu25*nutw.kappa()*(y[facei] + z0[facei]));
w*sqrt(k[celli])/(Cmu25*kappa*(y[facei] + z0[facei]));
G0[celli] +=
w
*(nutw[facei] + nuw[facei])
*magGradUw[facei]
*Cmu25*sqrt(k[celli])
/(nutw.kappa()*(y[facei] + z0[facei]));
/(kappa*(y[facei] + z0[facei]));
}
}
@ -109,6 +110,8 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::writeLocalEntries
{
z0_->writeData(os);
}
wallCoeffs_.writeEntries(os);
}

View File

@ -92,7 +92,8 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcUTau
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
const scalar Cmu25 = pow(Cmu_, 0.25);
const scalar Cmu25 = pow(wallCoeffs_.Cmu(), 0.25);
const scalar kappa = wallCoeffs_.kappa();
forAll(uTau, facei)
{
@ -108,7 +109,7 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcUTau
if (yPlus > yPlusCrit_)
{
scalar expTerm = exp(min(50.0, B_*mStar));
scalar powTerm = pow(yPlus, mStar/kappa_);
scalar powTerm = pow(yPlus, mStar/kappa);
factor = mStar/(expTerm*powTerm - 1.0 + ROOTVSMALL);
}
else
@ -163,8 +164,6 @@ void nutkFilmWallFunctionFvPatchScalarField::writeLocalEntries
);
os.writeEntryIfDifferent<scalar>("B", 5.5, B_);
os.writeEntryIfDifferent<scalar>("yPlusCrit", 11.05, yPlusCrit_);
os.writeEntryIfDifferent<scalar>("Cmu", 0.09, Cmu_);
os.writeEntryIfDifferent<scalar>("kappa", 0.41, kappa_);
}
@ -266,7 +265,7 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::yPlus() const
void nutkFilmWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
nutWallFunctionFvPatchScalarField::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}