Commit 6f28a119 authored by Andrew Heather's avatar Andrew Heather
Browse files

further updates to turbulence models

parent dae40256
......@@ -49,7 +49,7 @@ void DeardorffDiffStress::updateSubGridScaleFields(const volScalarField& K)
muSgs_ = ck_*rho()*sqrt(K)*delta();
muSgs_.correctBoundaryConditions();
alphaSgs_ = muSgs_/Prt();
alphaSgs_ = muSgs_/Prt_;
alphaSgs_.correctBoundaryConditions();
}
......
......@@ -60,6 +60,16 @@ GenEddyVisc::GenEddyVisc
)
),
Prt_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prt",
coeffDict_,
1.0
)
),
k_
(
IOobject
......@@ -135,6 +145,7 @@ bool GenEddyVisc::read()
if (LESModel::read())
{
ce_.readIfPresent(coeffDict());
Prt_.readIfPresent(coeffDict());
return true;
}
......
......@@ -68,7 +68,13 @@ class GenEddyVisc
protected:
// Model coefficients
dimensionedScalar ce_;
dimensionedScalar Prt_;
// Fields
volScalarField k_;
volScalarField muSgs_;
......
......@@ -64,6 +64,16 @@ GenSGSStress::GenSGSStress
)
),
Prt_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prt",
coeffDict_,
1.0
)
),
B_
(
IOobject
......
......@@ -69,7 +69,12 @@ class GenSGSStress
protected:
// Model coefficients
dimensionedScalar ce_;
dimensionedScalar Prt_;
// Fields
volSymmTensorField B_;
volScalarField muSgs_;
......
......@@ -81,45 +81,7 @@ LESModel::LESModel
k0_("k0", dimVelocity*dimVelocity, SMALL),
delta_(LESdelta::New("delta", U.mesh(), *this)),
wallFunctionDict_(subDictPtr("wallFunctionCoeffs")),
kappa_
(
dimensioned<scalar>::lookupOrAddToDict
(
"kappa",
wallFunctionDict_,
0.4187
)
),
E_
(
dimensioned<scalar>::lookupOrAddToDict
(
"E",
wallFunctionDict_,
9.0
)
),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
wallFunctionDict_,
0.07
)
),
Prt_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prt",
wallFunctionDict_,
0.85
)
)
delta_(LESdelta::New("delta", U.mesh(), *this))
{
readIfPresent("k0", k0_);
}
......@@ -201,20 +163,10 @@ bool LESModel::read()
coeffDict_ <<= *dictPtr;
}
if (const dictionary* dictPtr = subDictPtr("wallFunctionCoeffs"))
{
wallFunctionDict_ <<= *dictPtr;
}
readIfPresent("k0", k0_);
delta_().read(*this);
kappa_.readIfPresent(wallFunctionDict_);
E_.readIfPresent(wallFunctionDict_);
Cmu_.readIfPresent(wallFunctionDict_);
Prt_.readIfPresent(wallFunctionDict_);
return true;
}
else
......
......@@ -87,19 +87,6 @@ protected:
autoPtr<LESdelta> delta_;
// Wall function properties
//- Wall function dictionary
dictionary wallFunctionDict_;
dimensionedScalar kappa_;
dimensionedScalar E_;
dimensionedScalar Cmu_;
dimensionedScalar Prt_;
// Protected Member Functions
......@@ -173,55 +160,34 @@ public:
// Member Functions
//- Const access to the coefficients dictionary,
// which provides info. about choice of models,
// and all related data (particularly model coefficients).
inline const dictionary& coeffDict() const
{
return coeffDict_;
}
//- Return the value of k0 which k is not allowed to be less than
const dimensionedScalar& k0() const
{
return k0_;
}
//- Allow k0 to be changed
dimensionedScalar& k0()
{
return k0_;
}
// Access
//- Access function to filter width
inline const volScalarField& delta() const
{
return delta_();
}
//- Const access to the coefficients dictionary,
// which provides info. about choice of models,
// and all related data (particularly model coefficients).
inline const dictionary& coeffDict() const
{
return coeffDict_;
}
//- Return kappa for use in wall-functions
dimensionedScalar kappa() const
{
return kappa_;
}
//- Return the value of k0 which k is not allowed to be less than
const dimensionedScalar& k0() const
{
return k0_;
}
//- Return E for use in wall-functions
dimensionedScalar E() const
{
return E_;
}
//- Allow k0 to be changed
dimensionedScalar& k0()
{
return k0_;
}
//- Return Cmu for use in wall-functions
dimensionedScalar Cmu() const
{
return Cmu_;
}
//- Access function to filter width
inline const volScalarField& delta() const
{
return delta_();
}
//- Return turbulent Prandtl number for use in wall-functions
dimensionedScalar Prt() const
{
return Prt_;
}
//- Return the SGS turbulent kinetic energy.
virtual tmp<volScalarField> k() const = 0;
......
......@@ -57,7 +57,7 @@ void Smagorinsky::updateSubGridScaleFields(const volTensorField& gradU)
muSgs_ = ck_*rho()*delta()*sqrt(k_);
muSgs_.correctBoundaryConditions();
alphaSgs_ = muSgs_/Prt();
alphaSgs_ = muSgs_/Prt_;
alphaSgs_.correctBoundaryConditions();
}
......
......@@ -50,7 +50,7 @@ void SpalartAllmaras::updateSubGridScaleFields()
muSgs_.internalField() = rho()*fv1()*nuTilda_.internalField();
muSgs_.correctBoundaryConditions();
alphaSgs_ = muSgs_/Prt();
alphaSgs_ = muSgs_/Prt_;
alphaSgs_.correctBoundaryConditions();
}
......@@ -112,13 +112,22 @@ SpalartAllmaras::SpalartAllmaras
:
LESModel(typeName, rho, U, phi, thermoPhysicalModel),
alphaNut_
sigmaNut_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaNut",
"sigmaNut",
coeffDict_,
1.5
0.66666
)
),
Prt_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prt",
coeffDict_,
1.0
)
),
......@@ -185,7 +194,7 @@ SpalartAllmaras::SpalartAllmaras
0.4187
)
),
Cw1_(Cb1_/sqr(kappa_) + alphaNut_*(1.0 + Cb2_)),
Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
Cw2_
(
dimensioned<scalar>::lookupOrAddToDict
......@@ -299,11 +308,11 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& tgradU)
+ fvm::div(phi(), nuTilda_)
- fvm::laplacian
(
alphaNut_*(nuTilda_*rho() + mu()),
(nuTilda_*rho() + mu())/sigmaNut_,
nuTilda_,
"laplacian(DnuTildaEff,nuTilda)"
)
- alphaNut_*rho()*Cb2_*magSqr(fvc::grad(nuTilda_))
- rho()*Cb2_/sigmaNut_*magSqr(fvc::grad(nuTilda_))
==
rho()*Cb1_*Stilda*nuTilda_
- fvm::Sp(rho()*Cw1_*fw(Stilda)*nuTilda_/sqr(dTilda_), nuTilda_)
......@@ -320,10 +329,11 @@ bool SpalartAllmaras::read()
{
if (LESModel::read())
{
alphaNut_.readIfPresent(coeffDict());
sigmaNut_.readIfPresent(coeffDict());
Prt_.readIfPresent(coeffDict());
Cb1_.readIfPresent(coeffDict());
Cb2_.readIfPresent(coeffDict());
Cw1_ = Cb1_/sqr(kappa_) + alphaNut_*(1.0 + Cb2_);
Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
Cw2_.readIfPresent(coeffDict());
Cw3_.readIfPresent(coeffDict());
Cv1_.readIfPresent(coeffDict());
......
......@@ -58,18 +58,29 @@ class SpalartAllmaras
{
// Private data
dimensionedScalar alphaNut_;
// Model coefficients
dimensionedScalar Cb1_;
dimensionedScalar Cb2_;
dimensionedScalar Cv1_;
dimensionedScalar Cv2_;
dimensionedScalar CDES_;
dimensionedScalar ck_;
dimensionedScalar kappa_;
dimensionedScalar Cw1_;
dimensionedScalar Cw2_;
dimensionedScalar Cw3_;
dimensionedScalar sigmaNut_;
dimensionedScalar Prt_;
dimensionedScalar Cb1_;
dimensionedScalar Cb2_;
dimensionedScalar Cv1_;
dimensionedScalar Cv2_;
dimensionedScalar CDES_;
dimensionedScalar ck_;
dimensionedScalar kappa_;
dimensionedScalar Cw1_;
dimensionedScalar Cw2_;
dimensionedScalar Cw3_;
// Fields
volScalarField nuTilda_;
volScalarField dTilda_;
volScalarField muSgs_;
volScalarField alphaSgs_;
// Private member functions
......@@ -86,11 +97,6 @@ class SpalartAllmaras
SpalartAllmaras(const SpalartAllmaras&);
SpalartAllmaras& operator=(const SpalartAllmaras&);
volScalarField nuTilda_;
volScalarField dTilda_;
volScalarField muSgs_;
volScalarField alphaSgs_;
public:
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -75,17 +75,15 @@ scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::Psmooth
scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
(
const scalar P,
const scalar Prat,
const scalar E,
const scalar kappa
const scalar Prat
) const
{
scalar ypt = 11.0;
for (int i=0; i<maxIters_; i++)
{
scalar f = ypt - (log(E*ypt)/kappa + P)/Prat;
scalar df = 1.0 - 1.0/(ypt*kappa*Prat);
scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
scalar yptNew = ypt - f/df;
if (yptNew < VSMALL)
......@@ -115,7 +113,10 @@ alphaSgsJayatillekeWallFunctionFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(p, iF)
fixedValueFvPatchScalarField(p, iF),
Prt_(0.85),
kappa_(0.41),
E_(9.0)
{
checkType();
}
......@@ -130,7 +131,10 @@ alphaSgsJayatillekeWallFunctionFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchScalarField(ptf, p, iF, mapper)
fixedValueFvPatchScalarField(ptf, p, iF, mapper),
Prt_(ptf.Prt_),
kappa_(ptf.kappa_),
E_(ptf.E_)
{}
......@@ -142,7 +146,10 @@ alphaSgsJayatillekeWallFunctionFvPatchScalarField
const dictionary& dict
)
:
fixedValueFvPatchScalarField(p, iF, dict)
fixedValueFvPatchScalarField(p, iF, dict),
Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.0))
{
checkType();
}
......@@ -151,10 +158,13 @@ alphaSgsJayatillekeWallFunctionFvPatchScalarField
alphaSgsJayatillekeWallFunctionFvPatchScalarField::
alphaSgsJayatillekeWallFunctionFvPatchScalarField
(
const alphaSgsJayatillekeWallFunctionFvPatchScalarField& tppsf
const alphaSgsJayatillekeWallFunctionFvPatchScalarField& awfpsf
)
:
fixedValueFvPatchScalarField(tppsf)
fixedValueFvPatchScalarField(awfpsf),
Prt_(awfpsf.Prt_),
kappa_(awfpsf.kappa_),
E_(awfpsf.E_)
{
checkType();
}
......@@ -163,11 +173,14 @@ alphaSgsJayatillekeWallFunctionFvPatchScalarField
alphaSgsJayatillekeWallFunctionFvPatchScalarField::
alphaSgsJayatillekeWallFunctionFvPatchScalarField
(
const alphaSgsJayatillekeWallFunctionFvPatchScalarField& tppsf,
const alphaSgsJayatillekeWallFunctionFvPatchScalarField& awfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(tppsf, iF)
fixedValueFvPatchScalarField(awfpsf, iF),
Prt_(awfpsf.Prt_),
kappa_(awfpsf.kappa_),
E_(awfpsf.E_)
{
checkType();
}
......@@ -182,11 +195,6 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
{
const LESModel& lesModel = db().lookupObject<LESModel>("LESProperties");
// Wall function constants
const scalar E = lesModel.E().value();
const scalar kappa = lesModel.kappa().value();
const scalar Prt = lesModel.Prt().value();
// Field data
const label patchI = patch().index();
......@@ -223,18 +231,18 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
do
{
scalar kUu = min(kappa*magUp[faceI]/uTau, maxExp_);
scalar kUu = min(kappa_*magUp[faceI]/uTau, maxExp_);
scalar fkUu = exp(kUu) - 1.0 - kUu*(1.0 + 0.5*kUu);
scalar f =
- uTau/(ry[faceI]*muw[faceI]/rhow[faceI])
+ magUp[faceI]/uTau
+ 1.0/E*(fkUu - 1.0/6.0*kUu*sqr(kUu));
+ 1.0/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu));
scalar df =
- 1.0/(ry[faceI]*muw[faceI]/rhow[faceI])
- magUp[faceI]/sqr(uTau)
- 1.0/E*kUu*fkUu/uTau;
- 1.0/E_*kUu*fkUu/uTau;
scalar uTauNew = uTau - f/df;
err = mag((uTau - uTauNew)/uTau);
......@@ -248,11 +256,11 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
scalar Pr = muw[faceI]/alphaw[faceI];
// Molecular-to-turbulenbt Prandtl number ratio
scalar Prat = Pr/Prt;
scalar Prat = Pr/Prt_;
// Thermal sublayer thickness
scalar P = Psmooth(Prat);
scalar yPlusTherm = this->yPlusTherm(P, Prat, E, kappa);
scalar yPlusTherm = this->yPlusTherm(P, Prat);
// Evaluate new effective thermal diffusivity
scalar alphaEff = 0.0;
......@@ -266,11 +274,11 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
else
{
scalar A = qDot[faceI]*rhow[faceI]*uTau/ry[faceI];
scalar B = qDot[faceI]*Prt*(1.0/kappa*log(E*yPlus) + P);
scalar magUc = uTau/kappa*log(E*yPlusTherm) - mag(Uw[faceI]);
scalar B = qDot[faceI]*Prt_*(1.0/kappa_*log(E_*yPlus) + P);
scalar magUc = uTau/kappa_*log(E_*yPlusTherm) - mag(Uw[faceI]);
scalar C =
0.5*rhow[faceI]*uTau
*(Prt*sqr(magUp[faceI]) + (Pr - Prt)*sqr(magUc));
*(Prt_*sqr(magUp[faceI]) + (Pr - Prt_)*sqr(magUc));
alphaEff = A/(B + C + VSMALL);
}
......@@ -281,7 +289,7 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
{
Info<< " uTau = " << uTau << nl
<< " Pr = " << Pr << nl
<< " Prt = " << Prt << nl
<< " Prt = " << Prt_ << nl
<< " qDot = " << qDot[faceI] << nl
<< " yPlus = " << yPlus << nl
<< " yPlusTherm = " << yPlusTherm << nl
......@@ -299,6 +307,16 @@ void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
}
void alphaSgsJayatillekeWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
......