Commit 481fd3e2 authored by Andrew Heather's avatar Andrew Heather
Browse files

caching yPlusLam locally

parent be655124
......@@ -40,7 +40,7 @@ namespace RASModels
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
scalar mutRoughWallFunctionFvPatchScalarField::fnRough
(
......@@ -65,6 +65,60 @@ scalar mutRoughWallFunctionFvPatchScalarField::fnRough
}
tmp<scalarField> mutRoughWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalarField& y = rasModel.y()[patchI];
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const tmp<volScalarField> tk = rasModel.k();
const volScalarField& k = tk();
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const scalar Cmu25 = pow(Cmu_, 0.25);
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(mutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar uStar = Cmu25*sqrt(k[faceCellI]);
scalar yPlus = uStar*y[faceI]/(muw[faceI]/rhow[faceI]);
scalar KsPlus = uStar*Ks_[faceI]/(muw[faceI]/rhow[faceI]);
scalar Edash = E_;
scalar yPlusLamNew = yPlusLam_;
if (KsPlus > 2.25)
{
Edash /= fnRough(KsPlus, Cs_[faceI]);
yPlusLamNew = rasModel.yPlusLam(kappa_, Edash);
}
if (debug)
{
Info<< "yPlus = " << yPlus
<< ", KsPlus = " << KsPlus
<< ", Edash = " << Edash
<< ", yPlusLam = " << yPlusLam_
<< endl;
}
if (yPlus > yPlusLamNew)
{
mutw[faceI] =
muw[faceI]*(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1);
}
}
return tmutw;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
......@@ -158,61 +212,6 @@ void mutRoughWallFunctionFvPatchScalarField::rmap
}
tmp<scalarField> mutRoughWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const tmp<volScalarField> tk = rasModel.k();
const volScalarField& k = tk();
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const scalar Cmu25 = pow(Cmu_, 0.25);
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(mutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar uStar = Cmu25*sqrt(k[faceCellI]);
scalar yPlus = uStar*y[faceI]/(muw[faceI]/rhow[faceI]);
scalar KsPlus = uStar*Ks_[faceI]/(muw[faceI]/rhow[faceI]);
scalar Edash = E_;
scalar yPlusLamNew = yPlusLam;
if (KsPlus > 2.25)
{
Edash /= fnRough(KsPlus, Cs_[faceI]);
yPlusLamNew = rasModel.yPlusLam(kappa_, Edash);
}
if (debug)
{
Info<< "yPlus = " << yPlus
<< ", KsPlus = " << KsPlus
<< ", Edash = " << Edash
<< ", yPlusLam = " << yPlusLam
<< endl;
}
if (yPlus > yPlusLamNew)
{
mutw[faceI] =
muw[faceI]*(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1);
}
}
return tmutw;
}
void mutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
......
......@@ -39,6 +39,8 @@ namespace compressible
namespace RASModels
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<scalarField>
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
(
......@@ -48,7 +50,6 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
......@@ -74,7 +75,7 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
const scalar kappaRe = kappa_*Re;
scalar yp = yPlusLam;
scalar yp = yPlusLam_;
const scalar ryPlusLam = 1.0/yp;
int iter = 0;
......@@ -146,7 +147,7 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
const scalar kappaRe = kappa_*Re;
scalar yp = yPlusLam;
scalar yp = yPlusLam_;
const scalar ryPlusLam = 1.0/yp;
int iter = 0;
......@@ -171,6 +172,38 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
}
tmp<scalarField>
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalarField& y = rasModel.y()[patchI];
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
scalarField magUp = mag(Uw.patchInternalField() - Uw);
tmp<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(yPlus, facei)
{
if (yPlus[facei] > yPlusLam_)
{
const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei];
mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1);
}
}
return tmutw;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
......@@ -247,39 +280,6 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<scalarField>
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
scalarField magUp = mag(Uw.patchInternalField() - Uw);
tmp<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(yPlus, facei)
{
if (yPlus[facei] > yPlusLam)
{
const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei];
mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1);
}
}
return tmutw;
}
tmp<scalarField>
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
{
......
......@@ -50,7 +50,6 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI];
const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
......@@ -62,7 +61,7 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
{
scalar kappaRe = kappa_*magUp[faceI]*y[faceI]/(muw[faceI]/rhow[faceI]);
scalar yp = yPlusLam;
scalar yp = yPlusLam_;
scalar ryPlusLam = 1.0/yp;
int iter = 0;
......@@ -82,6 +81,35 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
}
tmp<scalarField>
mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
tmp<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(yPlus, faceI)
{
if (yPlus[faceI] > yPlusLam_)
{
mutw[faceI] =
muw[faceI]*(yPlus[faceI]*kappa_/log(E_*yPlus[faceI]) - 1.0);
}
}
return tmutw;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
......@@ -143,36 +171,6 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<scalarField>
mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
tmp<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(yPlus, faceI)
{
if (yPlus[faceI] > yPlusLam)
{
mutw[faceI] =
muw[faceI]*(yPlus[faceI]*kappa_/log(E_*yPlus[faceI]) - 1.0);
}
}
return tmutw;
}
tmp<scalarField>
mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
{
......
......@@ -105,6 +105,21 @@ tmp<scalarField> mutSpalartAllmarasWallFunctionFvPatchScalarField::calcUTau
}
tmp<scalarField>
mutSpalartAllmarasWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magGradU = mag(Uw.snGrad());
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
return max(0.0, rhow*sqr(calcUTau(magGradU))/magGradU - muw);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mutSpalartAllmarasWallFunctionFvPatchScalarField::
......@@ -166,21 +181,6 @@ mutSpalartAllmarasWallFunctionFvPatchScalarField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<scalarField>
mutSpalartAllmarasWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magGradU = mag(Uw.snGrad());
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
return max(0.0, rhow*sqr(calcUTau(magGradU))/magGradU - muw);
}
tmp<scalarField>
mutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const
{
......
......@@ -40,7 +40,7 @@ namespace compressible
namespace RASModels
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void mutWallFunctionFvPatchScalarField::checkType()
{
......@@ -56,6 +56,63 @@ void mutWallFunctionFvPatchScalarField::checkType()
}
scalar mutWallFunctionFvPatchScalarField::calcYPlusLam
(
const scalar kappa,
const scalar E
) const
{
scalar ypl = 11.0;
for (int i=0; i<10; i++)
{
ypl = log(E*ypl)/kappa;
}
return ypl;
}
tmp<scalarField> mutWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalarField& y = rasModel.y()[patchI];
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const tmp<volScalarField> tk = rasModel.k();
const volScalarField& k = tk();
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const scalar Cmu25 = pow(Cmu_, 0.25);
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(mutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar yPlus =
Cmu25*y[faceI]*sqrt(k[faceCellI])/(muw[faceI]/rhow[faceI]);
if (yPlus > yPlusLam_)
{
mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1);
}
}
return tmutw;
}
void mutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
{
os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
......@@ -67,7 +124,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
fixedValueFvPatchScalarField(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8)
E_(9.8),
yPlusLam_(calcYPlusLam(kappa_, E_))
{}
......@@ -82,7 +140,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
fixedValueFvPatchScalarField(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_)
E_(ptf.E_),
yPlusLam_(ptf.yPlusLam_)
{}
......@@ -96,7 +155,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
fixedValueFvPatchScalarField(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8))
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
yPlusLam_(calcYPlusLam(kappa_, E_))
{}
......@@ -108,7 +168,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
fixedValueFvPatchScalarField(wfpsf),
Cmu_(wfpsf.Cmu_),
kappa_(wfpsf.kappa_),
E_(wfpsf.E_)
E_(wfpsf.E_),
yPlusLam_(wfpsf.yPlusLam_)
{}
......@@ -121,7 +182,8 @@ mutWallFunctionFvPatchScalarField::mutWallFunctionFvPatchScalarField
fixedValueFvPatchScalarField(wfpsf, iF),
Cmu_(wfpsf.Cmu_),
kappa_(wfpsf.kappa_),
E_(wfpsf.E_)
E_(wfpsf.E_),
yPlusLam_(wfpsf.yPlusLam_)
{}
......@@ -135,40 +197,6 @@ void mutWallFunctionFvPatchScalarField::updateCoeffs()
}
tmp<scalarField> mutWallFunctionFvPatchScalarField::calcMut() const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const tmp<volScalarField> tk = rasModel.k();
const volScalarField& k = tk();
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const scalar Cmu25 = pow(Cmu_, 0.25);
tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
scalarField& mutw = tmutw();
forAll(mutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar yPlus =