Commit 851a5f1c authored by henry's avatar henry
Browse files

Many corrections and improvements.

parent 4f5f66e7
......@@ -26,7 +26,6 @@ License
#include "SpalartAllmaras.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -54,7 +53,7 @@ tmp<volScalarField> SpalartAllmaras::fv1() const
tmp<volScalarField> SpalartAllmaras::fv2() const
{
volScalarField chi = nuTilda_/nu();
return 1.0/pow3(scalar(1) + chi/Cv2_);
return 1/pow3(scalar(1) + chi/Cv2_);
}
......@@ -71,70 +70,68 @@ tmp<volScalarField> SpalartAllmaras::fv3() const
}
tmp<volScalarField> SpalartAllmaras::calcS(const volTensorField& gradU)
tmp<volScalarField> SpalartAllmaras::S(const volTensorField& gradU) const
{
return ::sqrt(2.0)*mag(skew(gradU));
return sqrt(2.0)*mag(skew(gradU));
}
tmp<volScalarField> SpalartAllmaras::calcSTilda(const volTensorField& gradU)
tmp<volScalarField> SpalartAllmaras::STilda
(
const volScalarField& S,
const volScalarField& dTilda
) const
{
return fv3()*calcS(gradU) + fv2()*nuTilda_/sqr(kappa_*dTilda_);
return fv3()*S + fv2()*nuTilda_/sqr(kappa_*dTilda);
}
tmp<volScalarField> SpalartAllmaras::r
(
const volScalarField& visc,
const volScalarField& S
const volScalarField& S,
const volScalarField& dTilda
) const
{
tmp<volScalarField> tr
return min
(
new volScalarField
(
min
(
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)
*sqr(kappa_*dTilda_)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
),
scalar(10.0)
)
)
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)
*sqr(kappa_*dTilda)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
),
scalar(10)
);
return tr;
}
tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& S) const
tmp<volScalarField> SpalartAllmaras::fw
(
const volScalarField& S,
const volScalarField& dTilda
) const
{
volScalarField r = this->r(nuTilda_, S);
volScalarField r = this->r(nuTilda_, S, dTilda);
volScalarField g = r + Cw2_*(pow6(r) - r);
return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}
void SpalartAllmaras::dTildaUpdate(const volScalarField&)
tmp<volScalarField> SpalartAllmaras::dTilda(const volScalarField&) const
{
if (mesh_.changing())
{
dTilda_ = min(CDES_*delta(), wallDist(mesh_).y());
}
return min(CDES_*delta(), y_);
}
......@@ -242,6 +239,8 @@ SpalartAllmaras::SpalartAllmaras
)
),
y_(mesh_),
nuTilda_
(
IOobject
......@@ -255,8 +254,6 @@ SpalartAllmaras::SpalartAllmaras
mesh_
),
dTilda_(min(CDES_*delta(), wallDist(mesh_).y())),
nuSgs_
(
IOobject
......@@ -277,11 +274,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
{
LESModel::correct(gradU);
const volScalarField STilda = calcSTilda(gradU);
const volScalarField S = calcS(gradU);
if (mesh_.changing())
{
y_.correct();
y_.boundaryField() = max(y_.boundaryField(), VSMALL);
}
dTildaUpdate(S);
const volScalarField S = this->S(gradU);
const volScalarField dTilda = this->dTilda(S);
const volScalarField STilda = this->STilda(S, dTilda);
fvScalarMatrix nuTildaEqn
(
......@@ -296,7 +297,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
- alphaNut_*Cb2_*magSqr(fvc::grad(nuTilda_))
==
Cb1_*STilda*nuTilda_
- fvm::Sp(Cw1_*fw(STilda)*nuTilda_/sqr(dTilda_), nuTilda_)
- fvm::Sp(Cw1_*fw(STilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_)
);
nuTildaEqn.relax();
......@@ -310,9 +311,15 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
}
tmp<volScalarField> SpalartAllmaras::k() const
{
return sqr(nuSgs()/ck_/dTilda(S(fvc::grad(U()))));
}
tmp<volScalarField> SpalartAllmaras::epsilon() const
{
return 2.0*nuEff()*magSqr(symm(fvc::grad(U())));
return 2*nuEff()*magSqr(symm(fvc::grad(U())));
}
......
......@@ -38,6 +38,7 @@ SourceFiles
#include "LESModel.H"
#include "volFields.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -86,30 +87,39 @@ protected:
// Fields
wallDist y_;
volScalarField nuTilda_;
volScalarField dTilda_;
volScalarField nuSgs_;
// Protected member functions
// Helper functions
virtual tmp<volScalarField> fv1() const;
virtual tmp<volScalarField> fv2() const;
virtual tmp<volScalarField> fv3() const;
virtual tmp<volScalarField> S(const volTensorField& gradU) const;
virtual tmp<volScalarField> fv1() const;
virtual tmp<volScalarField> fv2() const;
virtual tmp<volScalarField> fv3() const;
//-
virtual tmp<volScalarField> calcS(const volTensorField& gradU);
virtual tmp<volScalarField> calcSTilda(const volTensorField& gradU);
virtual tmp<volScalarField> r
(
const volScalarField& visc,
const volScalarField& S
) const;
virtual tmp<volScalarField> fw(const volScalarField& S) const;
virtual tmp<volScalarField> STilda
(
const volScalarField& S,
const volScalarField& dTilda
) const;
virtual tmp<volScalarField> r
(
const volScalarField& visc,
const volScalarField& S,
const volScalarField& dTilda
) const;
//- Length scale calculation
virtual void dTildaUpdate(const volScalarField& S);
virtual tmp<volScalarField> fw
(
const volScalarField& S,
const volScalarField& dTilda
) const;
//- Length scale
virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
public:
......@@ -138,10 +148,7 @@ public:
// Member Functions
//- Return SGS kinetic energy
virtual tmp<volScalarField> k() const
{
return sqr(nuSgs()/ck_/dTilda_);
}
virtual tmp<volScalarField> k() const;
//- Return sub-grid disipation rate
virtual tmp<volScalarField> epsilon() const;
......
......@@ -26,7 +26,6 @@ License
#include "SpalartAllmarasDDES.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -50,51 +49,42 @@ tmp<volScalarField> SpalartAllmarasDDES::rd
const volScalarField& S
) const
{
volScalarField d = wallDist(mesh_).y();
tmp<volScalarField> trd
return min
(
new volScalarField
(
min
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*y_)
+ dimensionedScalar
(
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*d)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
), scalar(10.0)
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
)
),
scalar(10)
);
return trd;
}
tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S)
tmp<volScalarField> SpalartAllmarasDDES::fd(const volScalarField& S) const
{
return 1.0 - tanh(pow3(8.0*rd(nuSgs_ + transport().nu(), S)));
return 1 - tanh(pow3(8*rd(nuEff(), S)));
}
void SpalartAllmarasDDES::dTildaUpdate(const volScalarField& S)
tmp<volScalarField> SpalartAllmarasDDES::dTilda(const volScalarField& S) const
{
dTilda_ =
wallDist(mesh_).y()
- fd(S)*max
(
dimensionedScalar("zero", dimLength, 0.0),
wallDist(mesh_).y() - CDES_*delta()
);
return max
(
y_
- fd(S)
*max(y_ - CDES_*delta(), dimensionedScalar("zero", dimLength, 0)),
dimensionedScalar("small", dimLength, SMALL)
);
}
......
......@@ -62,6 +62,14 @@ class SpalartAllmarasDDES
{
// Private member functions
tmp<volScalarField> fd(const volScalarField& S) const;
tmp<volScalarField> rd
(
const volScalarField& visc,
const volScalarField& S
) const;
// Disallow default bitwise copy construct and assignment
SpalartAllmarasDDES(const SpalartAllmarasDDES&);
SpalartAllmarasDDES& operator=(const SpalartAllmarasDDES&);
......@@ -71,15 +79,8 @@ protected:
// Protected member functions
tmp<volScalarField> fd(const volScalarField& S);
tmp<volScalarField> rd
(
const volScalarField& visc,
const volScalarField& S
) const;
//- Length scale calculation
virtual void dTildaUpdate(const volScalarField& S);
//- Length scale
virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
public:
......
......@@ -26,7 +26,6 @@ License
#include "SpalartAllmarasIDDES.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -46,22 +45,29 @@ addToRunTimeSelectionTable(LESModel, SpalartAllmarasIDDES, dictionary);
tmp<volScalarField> SpalartAllmarasIDDES::alpha() const
{
return
0.25
- wallDist(mesh_).y()
/dimensionedScalar("hMax", dimLength, max(cmptMax(delta())));
return max
(
0.25 - y_/dimensionedScalar("hMax", dimLength, max(cmptMax(delta()))),
scalar(-5)
);
}
tmp<volScalarField> SpalartAllmarasIDDES::ft(const volScalarField& S) const
tmp<volScalarField> SpalartAllmarasIDDES::ft
(
const volScalarField& S
) const
{
return tanh(pow3(sqr(ct_)*r(nuSgs_, S)));
return tanh(pow3(sqr(ct_)*rd(nuSgs_, S)));
}
tmp<volScalarField> SpalartAllmarasIDDES::fl(const volScalarField& S) const
tmp<volScalarField> SpalartAllmarasIDDES::fl
(
const volScalarField& S
) const
{
return tanh(pow(sqr(cl_)*r(transport().nu(), S), 10));
return tanh(pow(sqr(cl_)*rd(nu(), S), 10));
}
......@@ -71,33 +77,24 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
const volScalarField& S
) const
{
volScalarField d = wallDist(mesh_).y();
tmp<volScalarField> trd
return min
(
new volScalarField
(
min
(
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*d)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
), scalar(10.0)
)
)
visc
/(
max
(
S,
dimensionedScalar("SMALL", S.dimensions(), SMALL)
)*sqr(kappa_*y_)
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, 2 , -1, 0, 0),
ROOTVSMALL
)
),
scalar(10)
);
return trd;
}
......@@ -105,43 +102,38 @@ tmp<volScalarField> SpalartAllmarasIDDES::rd
tmp<volScalarField> SpalartAllmarasIDDES::fd(const volScalarField& S) const
{
return 1.0 - tanh(pow3(8.0*rd(nuSgs_+transport().nu(), S)));
return 1 - tanh(pow3(8*rd(nuEff(), S)));
}
void SpalartAllmarasIDDES::dTildaUpdate(const volScalarField& S)
tmp<volScalarField> SpalartAllmarasIDDES::dTilda(const volScalarField& S) const
{
volScalarField alpha = this->alpha();
volScalarField expTerm = exp(sqr(alpha));
volScalarField fHill =
2.0*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
volScalarField fStep = min(2.0*pow(expTerm, -9.0), scalar(1));
volScalarField fHyb = max(1.0 - fd(S), fStep);
volScalarField fAmp = 1.0 - max(ft(S), fl(S));
volScalarField fRestore = max(fHill - 1.0, scalar(0))*fAmp;
2*(pos(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
// volScalarField ft2 = IGNORING ft2 terms
volScalarField fStep = min(2*pow(expTerm, -9.0), scalar(1));
volScalarField fHyb = max(1 - fd(S), fStep);
volScalarField fAmp = 1 - max(ft(S), fl(S));
volScalarField fRestore = max(fHill - 1, scalar(0))*fAmp;
// IGNORING ft2 terms
volScalarField Psi = sqrt
(
min
(
scalar(100),
(1.0 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
(1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*fv2())/max(SMALL, fv1())
)
);
dTilda_ = max
return max
(
dimensionedScalar("zero", dimLength, 0.0),
fHyb*(1.0 + fRestore*Psi)*wallDist(mesh_).y()
+ (1.0 - fHyb)*CDES_*Psi*delta()
dimensionedScalar("SMALL", dimLength, SMALL),
fHyb*(1 + fRestore*Psi)*y_
+ (1 - fHyb)*CDES_*Psi*delta()
);
}
......
......@@ -69,12 +69,16 @@ class SpalartAllmarasIDDES
tmp<volScalarField> alpha() const;
tmp<volScalarField> ft(const volScalarField& S) const;
tmp<volScalarField> fl(const volScalarField& S) const;
tmp<volScalarField> rd
(
const volScalarField& visc,
const volScalarField& S
) const;
//- Delay function
tmp<volScalarField> fd(const volScalarField& S) const;
// Disallow default bitwise copy construct and assignment
SpalartAllmarasIDDES(const SpalartAllmarasIDDES&);
SpalartAllmarasIDDES& operator=(const SpalartAllmarasIDDES&);
......@@ -84,11 +88,9 @@ protected:
// Protected member functions
//- Delay function
tmp<volScalarField> fd(const volScalarField& S) const;
//- Length scale
virtual tmp<volScalarField> dTilda(const volScalarField& S) const;
//- Length scale calculation
virtual void dTildaUpdate(const volScalarField& S);
public:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment