Commit 0c8fb634 authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: enforce consistent boundness on turbulence models.

- remove epsilonSmall, omegaSmall
- k0/epsilon0/omega0 become kMin/epsilonMin/omegaMin
- add qMin/zetaMin for consistency

These files still need some attention:
    dynOneEqEddy.C
    NonlinearKEShih.C
    settlingFoam

BUG: incompressible::LESModels:dynOneEqEddy::correct()
- avoid tmp field destruction for consistency with the compressible
  version

Possible TODO:
   - set kMin to zero (instead of SMALL) and introduce kSmall
     to avoid division by zero
parent 24058702
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -121,9 +121,8 @@ PDRkEpsilon::PDRkEpsilon
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
autoCreateK("k", mesh_)
),
epsilon_
(
IOobject
......@@ -134,9 +133,8 @@ PDRkEpsilon::PDRkEpsilon
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
autoCreateEpsilon("epsilon", mesh_)
),
mut_
(
IOobject
......@@ -147,9 +145,8 @@ PDRkEpsilon::PDRkEpsilon
IOobject::NO_READ,
IOobject::NO_WRITE
),
Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
autoCreateMut("mut", mesh_)
),
alphat_
(
IOobject
......@@ -163,7 +160,9 @@ PDRkEpsilon::PDRkEpsilon
autoCreateAlphat("alphat", mesh_)
)
{
mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
bound(epsilon_, epsilonMin_);
mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
alphat_ = mut_/Prt_;
......@@ -250,7 +249,7 @@ void PDRkEpsilon::correct()
if (!turbulence_)
{
// Re-calculate viscosity
mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
// Re-calculate thermal diffusivity
......@@ -302,7 +301,7 @@ void PDRkEpsilon::correct()
epsEqn().boundaryManipulate(epsilon_.boundaryField());
solve(epsEqn);
bound(epsilon_, epsilon0_);
bound(epsilon_, epsilonMin_);
// Turbulent kinetic energy equation
......@@ -320,7 +319,7 @@ void PDRkEpsilon::correct()
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
bound(k_, kMin_);
// Re-calculate viscosity
mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
......
......@@ -5,8 +5,8 @@ if (turbulence)
y.correct();
}
dimensionedScalar k0("k0", k.dimensions(), SMALL);
dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), SMALL);
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
volScalarField divU = fvc::div(phi/fvc::interpolate(rho));
......@@ -15,7 +15,7 @@ if (turbulence)
tgradU.clear();
volScalarField Gcoef =
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilon0);
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin);
#include "wallFunctions.H"
......@@ -40,7 +40,7 @@ if (turbulence)
epsEqn.relax();
epsEqn.solve();
bound(epsilon, epsilon0);
bound(epsilon, epsilonMin);
// Turbulent kinetic energy equation
......@@ -60,11 +60,13 @@ if (turbulence)
- fvm::Sp(rho*epsilon/k, k)
);
bound(k, k0);
//FIXME: why no kEqn.relax()?
bound(k, kMin);
//- Re-calculate viscosity
mut = rho*Cmu*sqr(k)/(epsilon + epsilon0);
mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
#include "wallViscosity.H"
}
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -30,11 +30,11 @@ License
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
void Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
{
scalar minVsf = min(vsf).value();
const scalar minVsf = min(vsf).value();
if (minVsf < vsf0.value())
if (minVsf < lowerBound.value())
{
Info<< "bounding " << vsf.name()
<< ", min: " << minVsf
......@@ -47,13 +47,13 @@ void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
max
(
vsf.internalField(),
fvc::average(max(vsf, vsf0))().internalField()
*pos(-vsf.internalField())
fvc::average(max(vsf, lowerBound))().internalField()
* pos(-vsf.internalField())
),
vsf0.value()
lowerBound.value()
);
vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
vsf.boundaryField() = max(vsf.boundaryField(), lowerBound.value());
}
}
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -23,11 +23,13 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
InNamespace
Foam::bound
Foam
Description
Bound the given scalar field if it has gone unbounded. Used extensively
in RAS and LES turbulence models.
Bound the given scalar field if it has gone unbounded.
Used extensively in RAS and LES turbulence models, but also of use
within solvers.
SourceFiles
bound.C
......@@ -49,7 +51,7 @@ namespace Foam
//- Bound the given scalar field if it has gone unbounded.
// Used extensively in RAS and LES turbulence models.
void bound(volScalarField& vsf, const dimensionedScalar& vsf0);
void bound(volScalarField&, const dimensionedScalar& lowerBound);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -132,7 +132,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
}
K = 0.5*tr(B_);
bound(K, k0());
bound(K, kMin_);
updateSubGridScaleFields(K);
}
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -79,11 +79,11 @@ LESModel::LESModel
printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
coeffDict_(subOrEmptyDict(type + "Coeffs")),
k0_("k0", dimVelocity*dimVelocity, SMALL),
kMin_("kMin", dimVelocity*dimVelocity, SMALL),
delta_(LESdelta::New("delta", U.mesh(), *this))
{
readIfPresent("k0", k0_);
readIfPresent("kMin", kMin_);
// Force the construction of the mesh deltaCoeffs which may be needed
// for the construction of the derived models and BCs
......@@ -167,7 +167,7 @@ bool LESModel::read()
coeffDict_ <<= *dictPtr;
}
readIfPresent("k0", k0_);
readIfPresent("kMin", kMin_);
delta_().read(*this);
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -83,7 +83,7 @@ protected:
Switch printCoeffs_;
dictionary coeffDict_;
dimensionedScalar k0_;
dimensionedScalar kMin_;
autoPtr<LESdelta> delta_;
......@@ -170,16 +170,16 @@ public:
return coeffDict_;
}
//- Return the value of k0 which k is not allowed to be less than
const dimensionedScalar& k0() const
//- Return the lower allowable limit for k (default: SMALL)
const dimensionedScalar& kMin() const
{
return k0_;
return kMin_;
}
//- Allow k0 to be changed
dimensionedScalar& k0()
//- Allow kMin to be changed
dimensionedScalar& kMin()
{
return k0_;
return kMin_;
}
//- Access function to filter width
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -74,10 +74,10 @@ dimensionedScalar dynOneEqEddy::ce_(const volSymmTensorField& D) const
pow(KK + filter_(k_), 1.5)/(2*delta()) - filter_(pow(k_, 1.5))/delta();
volScalarField ee =
2*delta()*ck_(D)*
(
2*delta()*ck_(D)
*(
filter_(sqrt(k_)*magSqr(D))
- 2*sqrt(KK + filter_(k_))*magSqr(filter_(D))
- 2*sqrt(KK + filter_(k_))*magSqr(filter_(D))
);
return average(ee*mm)/average(mm*mm);
......@@ -129,7 +129,11 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
- fvm::Sp(ce_(D)*rho()*sqrt(k_)/delta(), k_)
);
bound(k_, dimensionedScalar("0", k_.dimensions(), 1.0e-10));
//FIXME: why not this?
// kEqn.relax();
// kEqn.solve();
bound(k_, kMin_);
updateSubGridScaleFields(D);
}
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -117,7 +117,7 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
- fvm::Sp(ce_*rho()*sqrt(k_)/delta(), k_)
);
bound(k_, k0());
bound(k_, kMin_);
updateSubGridScaleFields();
}
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -107,7 +107,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
kEqn.relax();
kEqn.solve();
bound(k_, k0());
bound(k_, kMin_);
updateSubGridScaleFields();
}
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -229,7 +229,9 @@ LRR::LRR
<< exit(FatalError);
}
mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
bound(epsilon_, epsilonMin_);
mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
alphat_ = mut_/Prt_;
......@@ -324,7 +326,7 @@ void LRR::correct()
if (!turbulence_)
{
// Re-calculate viscosity
mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
// Re-calculate thermal diffusivity
......@@ -359,7 +361,7 @@ void LRR::correct()
epsEqn().boundaryManipulate(epsilon_.boundaryField());
solve(epsEqn);
bound(epsilon_, epsilon0_);
bound(epsilon_, epsilonMin_);
// Reynolds stress equation
......@@ -406,15 +408,15 @@ void LRR::correct()
R_.dimensions(),
symmTensor
(
k0_.value(), -GREAT, -GREAT,
k0_.value(), -GREAT,
k0_.value()
kMin_.value(), -GREAT, -GREAT,
kMin_.value(), -GREAT,
kMin_.value()
)
)
);
k_ = 0.5*tr(R_);
bound(k_, k0_);
bound(k_, kMin_);
// Re-calculate viscosity
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -258,7 +258,9 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
<< exit(FatalError);
}
mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
bound(epsilon_, epsilonMin_);
mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
alphat_ = mut_/Prt_;
......@@ -357,7 +359,7 @@ void LaunderGibsonRSTM::correct()
if (!turbulence_)
{
// Re-calculate viscosity
mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
// Re-calculate thermal diffusivity
......@@ -397,7 +399,7 @@ void LaunderGibsonRSTM::correct()
epsEqn().boundaryManipulate(epsilon_.boundaryField());
solve(epsEqn);
bound(epsilon_, epsilon0_);
bound(epsilon_, epsilonMin_);
// Reynolds stress equation
......@@ -453,15 +455,15 @@ void LaunderGibsonRSTM::correct()
R_.dimensions(),
symmTensor
(
k0_.value(), -GREAT, -GREAT,
k0_.value(), -GREAT,
k0_.value()
kMin_.value(), -GREAT, -GREAT,
kMin_.value(), -GREAT,
kMin_.value()
)
)
);
k_ == 0.5*tr(R_);
bound(k_, k0_);
bound(k_, kMin_);
// Re-calculate turbulent viscosity
......
......@@ -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) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -187,7 +187,9 @@ LaunderSharmaKE::LaunderSharmaKE
autoCreateAlphat("alphat", mesh_)
)
{
mut_ = rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
bound(epsilon_, epsilonMin_);
mut_ = rho_*Cmu_*fMu()*sqr(k_)/epsilon_;
mut_.correctBoundaryConditions();
alphat_ = mut_/Prt_;
......@@ -275,7 +277,7 @@ void LaunderSharmaKE::correct()
if (!turbulence_)
{
// Re-calculate viscosity
mut_ == rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
mut_ == rho_*Cmu_*fMu()*sqr(k_)/epsilon_;
// Re-calculate thermal diffusivity
alphat_ = mut_/Prt_;
......@@ -320,7 +322,7 @@ void LaunderSharmaKE::correct()
epsEqn().relax();
solve(epsEqn);
bound(epsilon_, epsilon0_);
bound(epsilon_, epsilonMin_);
// Turbulent kinetic energy equation
......@@ -338,11 +340,11 @@ void LaunderSharmaKE::correct()
kEqn().relax();
solve(kEqn);
bound(k_, k0_);
bound(k_, kMin_);
// Re-calculate viscosity
mut_ == Cmu_*fMu()*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
mut_ == Cmu_*fMu()*rho_*sqr(k_)/epsilon_;
// Re-calculate thermal diffusivity
......
......@@ -81,11 +81,9 @@ RASModel::RASModel
printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
coeffDict_(subOrEmptyDict(type + "Coeffs")),
k0_("k0", dimVelocity*dimVelocity, SMALL),
epsilon0_("epsilon0", k0_.dimensions()/dimTime, SMALL),
epsilonSmall_("epsilonSmall", epsilon0_.dimensions(), SMALL),
omega0_("omega0", dimless/dimTime, SMALL),
omegaSmall_("omegaSmall", omega0_.dimensions(), SMALL),
kMin_("kMin", dimVelocity*dimVelocity, SMALL),
epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
omegaMin_("omegaMin", dimless/dimTime, SMALL),
y_(mesh_)
{
......@@ -219,11 +217,9 @@ bool RASModel::read()
coeffDict_ <<= *dictPtr;
}
k0_.readIfPresent(*this);
epsilon0_.readIfPresent(*this);
epsilonSmall_.readIfPresent(*this);
omega0_.readIfPresent(*this);
omegaSmall_.readIfPresent(*this);
kMin_.readIfPresent(*this);
epsilonMin_.readIfPresent(*this);
omegaMin_.readIfPresent(*this);
return true;
}
......
......@@ -93,19 +93,13 @@ protected:
scalar yPlusLam_;
//- Lower limit of k
dimensionedScalar k0_;
dimensionedScalar kMin_;
//- Lower limit of epsilon
dimensionedScalar epsilon0_;
//- Small epsilon value used to avoid divide by zero
dimensionedScalar epsilonSmall_;
dimensionedScalar epsilonMin_;
//- Lower limit for omega
dimensionedScalar omega0_;
//- Small omega value used to avoid divide by zero
dimensionedScalar omegaSmall_;
dimensionedScalar omegaMin_;
//- Near wall distance boundary field
nearWallDist y_;
......@@ -185,66 +179,40 @@ public:
// Access
//- Return the lower allowable limit for k
const dimensionedScalar& k0() const
{
return k0_;
}
//- Return the lower allowable limit for epsilon
const dimensionedScalar& epsilon0() const
{
return epsilon0_;
}
//- Return the value of epsilonSmall which is added to epsilon when
// calculating mut
const dimensionedScalar& epsilonSmall() const
{
return epsilonSmall_;
}
//- Return the lower allowable limit for omega
const dimensionedScalar& omega0() const