Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2011-10-28 11:46:45 +01:00
159 changed files with 2996 additions and 803 deletions

View File

@ -655,6 +655,27 @@ Foam::scalar Foam::face::sweptVol
const pointField& newPoints
) const
{
if (size() == 3)
{
return
(
triPointRef
(
oldPoints[operator[](0)],
oldPoints[operator[](1)],
oldPoints[operator[](2)]
).sweptVol
(
triPointRef
(
newPoints[operator[](0)],
newPoints[operator[](1)],
newPoints[operator[](2)]
)
)
);
}
scalar sv = 0;
// Calculate the swept volume by breaking the face into triangles and

View File

@ -280,6 +280,7 @@ $(limitedSchemes)/limitWith/limitWith.C
multivariateSchemes = $(surfaceInterpolation)/multivariateSchemes
$(multivariateSchemes)/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationSchemes.C
$(multivariateSchemes)/multivariateSelectionScheme/multivariateSelectionSchemes.C
$(multivariateSchemes)/multivariateIndependentScheme/multivariateIndependentSchemes.C
$(multivariateSchemes)/upwind/multivariateUpwind.C
$(multivariateSchemes)/Gamma/multivariateGamma.C
$(multivariateSchemes)/vanLeer/multivariateVanLeer.C

View File

@ -42,8 +42,8 @@ void Foam::pimpleControl::read()
// Read solution controls
const dictionary& pimpleDict = dict();
nOuterCorr_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
nCorr_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
turbOnFinalIterOnly_ =
pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
}
@ -51,12 +51,14 @@ void Foam::pimpleControl::read()
bool Foam::pimpleControl::criteriaSatisfied()
{
if ((corr_ == 0) || residualControl_.empty() || finalIter())
// no checks on first iteration - nothing has been calculated yet
if ((corr_ == 1) || residualControl_.empty() || finalIter())
{
return false;
}
bool firstIter = corr_ == 1;
bool storeIni = this->storeInitialResiduals();
bool achieved = true;
bool checked = false; // safety that some checks were indeed performed
@ -73,7 +75,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
checked = true;
if (firstIter)
if (storeIni)
{
residualControl_[fieldI].initialResidual =
sp.first().initialResidual();
@ -83,7 +85,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
bool relCheck = false;
scalar relative = 0.0;
if (!firstIter)
if (!storeIni)
{
const scalar iniRes =
residualControl_[fieldI].initialResidual
@ -97,9 +99,10 @@ bool Foam::pimpleControl::criteriaSatisfied()
if (debug)
{
Info<< algorithmName_ << "loop statistics:" << endl;
Info<< algorithmName_ << " loop:" << endl;
Info<< " " << variableName << " iter " << corr_
Info<< " " << variableName
<< " PIMPLE iter " << corr_
<< ": ini res = "
<< residualControl_[fieldI].initialResidual
<< ", abs tol = " << residual
@ -120,25 +123,25 @@ bool Foam::pimpleControl::criteriaSatisfied()
Foam::pimpleControl::pimpleControl(fvMesh& mesh)
:
solutionControl(mesh, "PIMPLE"),
nOuterCorr_(0),
nCorr_(0),
corr_(0),
nCorrPIMPLE_(0),
nCorrPISO_(0),
corrPISO_(0),
turbOnFinalIterOnly_(true)
{
read();
if (nOuterCorr_ > 1)
if (nCorrPIMPLE_ > 1)
{
Info<< nl;
if (residualControl_.empty())
{
Info<< algorithmName_ << ": no residual control data found. "
<< "Calculations will employ " << nOuterCorr_
<< "Calculations will employ " << nCorrPIMPLE_
<< " corrector loops" << nl << endl;
}
else
{
Info<< algorithmName_ << ": max iterations = " << nOuterCorr_
Info<< algorithmName_ << ": max iterations = " << nCorrPIMPLE_
<< endl;
forAll(residualControl_, i)
{
@ -164,4 +167,60 @@ Foam::pimpleControl::~pimpleControl()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::pimpleControl::loop()
{
read();
corr_++;
if (debug)
{
Info<< algorithmName_ << " loop: corr = " << corr_ << endl;
}
if (corr_ == nCorrPIMPLE_ + 1)
{
if ((!residualControl_.empty()) && (nCorrPIMPLE_ != 1))
{
Info<< algorithmName_ << ": not converged within "
<< nCorrPIMPLE_ << " iterations" << endl;
}
corr_ = 0;
mesh_.data::remove("finalIteration");
return false;
}
bool completed = false;
if (criteriaSatisfied())
{
Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
<< endl;
completed = true;
}
else
{
if (finalIter())
{
mesh_.data::add("finalIteration", true);
}
if (corr_ <= nCorrPIMPLE_)
{
if (nCorrPIMPLE_ != 1)
{
Info<< algorithmName_ << ": iteration " << corr_ << endl;
storePrevIterFields();
}
completed = false;
}
}
return !completed;
}
// ************************************************************************* //

View File

@ -55,13 +55,13 @@ protected:
// Solution controls
//- Maximum number of PIMPLE correctors
label nOuterCorr_;
label nCorrPIMPLE_;
//- Maximum number of PISO correctors
label nCorr_;
label nCorrPISO_;
//- Current PIMPLE corrector
label corr_;
//- Current PISO corrector
label corrPISO_;
//- Flag to indicate whether to only solve turbulence on final iter
bool turbOnFinalIterOnly_;
@ -105,41 +105,35 @@ public:
// Access
//- Current corrector index
inline label corr() const;
//- Maximum number of PIMPLE correctors
inline label nOuterCorr() const;
inline label nCorrPIMPLE() const;
//- Maximum number of PISO correctors
inline label nCorr() const;
inline label nCorrPISO() const;
//- Current PISO corrector index
inline label corrPISO() const;
// Solution control
//- Loop start
inline bool start();
//- PIMPLE loop
virtual bool loop();
//- Loop loop
inline bool loop();
//- Pressure corrector loop
inline bool correct();
//- Helper function to identify when to store the intial residuals
inline bool storeInitialResiduals() const;
//- Helper function to identify final PIMPLE (outer) iteration
inline bool finalIter() const;
//- Helper function to identify final inner iteration
inline bool finalInnerIter
(
const label corr,
const label nonOrth
) const;
inline bool finalInnerIter() const;
//- Helper function to identify whether to solve for turbulence
inline bool turbCorr() const;
// Member Operators
void operator++(int);
};

View File

@ -25,89 +25,64 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::label Foam::pimpleControl::corr() const
inline Foam::label Foam::pimpleControl::nCorrPIMPLE() const
{
return corr_;
return nCorrPIMPLE_;
}
inline Foam::label Foam::pimpleControl::nOuterCorr() const
inline Foam::label Foam::pimpleControl::nCorrPISO() const
{
return nOuterCorr_;
return nCorrPISO_;
}
inline Foam::label Foam::pimpleControl::nCorr() const
inline Foam::label Foam::pimpleControl::corrPISO() const
{
return nCorr_;
return corrPISO_;
}
inline bool Foam::pimpleControl::start()
inline bool Foam::pimpleControl::correct()
{
corr_ = 0;
corrPISO_++;
return true;
}
inline bool Foam::pimpleControl::loop()
{
read();
if (criteriaSatisfied())
if (debug)
{
Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
<< endl;
return false;
Info<< algorithmName_ << " correct: corrPISO = " << corrPISO_ << endl;
}
if (corrPISO_ <= nCorrPISO_)
{
return true;
}
else
{
if (finalIter())
{
mesh_.data::add("finalIteration", true);
}
if (corr_ < nOuterCorr_)
{
if (nOuterCorr_ != 1)
{
Info<< algorithmName_ << ": iteration " << corr_ + 1 << endl;
storePrevIterFields();
}
return true;
}
else
{
if ((!residualControl_.empty()) && (nOuterCorr_ != 1))
{
Info<< algorithmName_ << ": not converged within "
<< nOuterCorr_ << " iterations" << endl;
}
return false;
}
corrPISO_ = 0;
return false;
}
}
inline bool Foam::pimpleControl::storeInitialResiduals() const
{
// start from second PIMPLE iteration
return (corr_ == 2) && (corrPISO_ == 0) && (corrNonOrtho_ == 0);
}
inline bool Foam::pimpleControl::finalIter() const
{
return corr_ == nOuterCorr_-1;
return corr_ == nCorrPIMPLE_;
}
inline bool Foam::pimpleControl::finalInnerIter
(
const label corr,
const label nonOrth
) const
inline bool Foam::pimpleControl::finalInnerIter() const
{
return
corr_ == nOuterCorr_-1
&& corr == nCorr_-1
&& nonOrth == nNonOrthCorr_;
corr_ == nCorrPIMPLE_
&& corrPISO_ == nCorrPISO_
&& corrNonOrtho_ == nNonOrthCorr_ + 1;
}
@ -117,15 +92,4 @@ inline bool Foam::pimpleControl::turbCorr() const
}
inline void Foam::pimpleControl::operator++(int)
{
if (finalIter())
{
mesh_.data::remove("finalIteration");
}
corr_++;
}
// ************************************************************************* //

View File

@ -119,4 +119,38 @@ Foam::simpleControl::~simpleControl()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::simpleControl::loop()
{
read();
Time& time = const_cast<Time&>(mesh_.time());
if (initialised_)
{
if (criteriaSatisfied())
{
Info<< nl << algorithmName_ << " solution converged in "
<< time.timeName() << " iterations" << nl << endl;
// Set to finalise calculation
time.writeAndEnd();
}
else
{
storePrevIterFields();
}
}
else
{
initialised_ = true;
storePrevIterFields();
}
return time.loop();
}
// ************************************************************************* //

View File

@ -96,7 +96,7 @@ public:
// Solution control
//- Loop loop
inline bool loop();
virtual bool loop();
};
@ -106,10 +106,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "simpleControlI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -170,7 +170,9 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
algorithmName_(algorithmName),
nNonOrthCorr_(0),
momentumPredictor_(true),
transonic_(false)
transonic_(false),
corr_(0),
corrNonOrtho_(0)
{}

View File

@ -82,6 +82,15 @@ protected:
bool transonic_;
// Evolution
//- Current corrector loop index
label corr_;
//- Current non-orthogonal corrector loop index
label corrNonOrtho_;
// Protected Member Functions
//- Read controls from fvSolution dictionary
@ -137,17 +146,35 @@ public:
//- Return the solution dictionary
inline const dictionary& dict() const;
//- Current corrector loop index
inline label corr() const;
//- Current non-orthogonal corrector index
inline label corrNonOrtho() const;
// Solution control
//- Maximum number of non-orthogonal correctors
inline label nNonOrthCorr() const;
//- Helper function to identify final non-orthogonal iteration
inline bool finalNonOrthogonalIter() const;
//- Flag to indicate to solve for momentum
inline bool momentumPredictor() const;
//- Flag to indicate to solve using transonic algorithm
inline bool transonic() const;
// Evolution
//- Main control loop
virtual bool loop() = 0;
//- Non-orthogonal corrector loop
inline bool correctNonOrthogonal();
};

View File

@ -31,12 +31,30 @@ inline const Foam::dictionary& Foam::solutionControl::dict() const
}
inline Foam::label Foam::solutionControl::corr() const
{
return corr_;
}
inline Foam::label Foam::solutionControl::corrNonOrtho() const
{
return corrNonOrtho_;
}
inline Foam::label Foam::solutionControl::nNonOrthCorr() const
{
return nNonOrthCorr_;
}
inline bool Foam::solutionControl::finalNonOrthogonalIter() const
{
return corrNonOrtho_ == nNonOrthCorr_ + 1;
}
inline bool Foam::solutionControl::momentumPredictor() const
{
return momentumPredictor_;
@ -49,4 +67,26 @@ inline bool Foam::solutionControl::transonic() const
}
inline bool Foam::solutionControl::correctNonOrthogonal()
{
corrNonOrtho_++;
if (debug)
{
Info<< algorithmName_ << " correctNonOrthogonal: corrNonOrtho = "
<< corrNonOrtho_ << endl;
}
if (corrNonOrtho_ <= nNonOrthCorr_ + 1)
{
return true;
}
else
{
corrNonOrtho_ = 0;
return false;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,57 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "multivariateIndependentScheme.H"
#include "limitedSurfaceInterpolationScheme.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "upwind.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::multivariateIndependentScheme<Type>::multivariateIndependentScheme
(
const fvMesh& mesh,
const typename multivariateSurfaceInterpolationScheme<Type>::
fieldTable& fields,
const surfaceScalarField& faceFlux,
Istream& schemeData
)
:
multivariateSurfaceInterpolationScheme<Type>
(
mesh,
fields,
faceFlux,
schemeData
),
schemes_(schemeData),
faceFlux_(faceFlux)
{}
// ************************************************************************* //

View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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::multivariateIndependentScheme
Description
Generic multi-variate discretisation scheme class for which any of the
NVD, CNVD or NVDV schemes may be selected for each variable and applied
independently.
This is equivalent to using separate "div" terms and schemes for each
variable/equation.
SourceFiles
multivariateIndependentScheme.C
\*---------------------------------------------------------------------------*/
#ifndef multivariateIndependentScheme_H
#define multivariateIndependentScheme_H
#include "multivariateSurfaceInterpolationScheme.H"
#include "limitedSurfaceInterpolationScheme.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class multivariateIndependentScheme Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class multivariateIndependentScheme
:
public multivariateSurfaceInterpolationScheme<Type>
{
// Private data
dictionary schemes_;
const surfaceScalarField& faceFlux_;
// Private Member Functions
//- Disallow default bitwise copy construct
multivariateIndependentScheme(const multivariateIndependentScheme&);
//- Disallow default bitwise assignment
void operator=(const multivariateIndependentScheme&);
public:
//- Runtime type information
TypeName("multivariateIndependent");
// Constructors
//- Construct for field, faceFlux and Istream
multivariateIndependentScheme
(
const fvMesh& mesh,
const typename multivariateSurfaceInterpolationScheme<Type>::
fieldTable& fields,
const surfaceScalarField& faceFlux,
Istream& schemeData
);
// Member Operators
tmp<surfaceInterpolationScheme<Type> > operator()
(
const GeometricField<Type, fvPatchField, volMesh>& field
) const
{
return surfaceInterpolationScheme<Type>::New
(
faceFlux_.mesh(),
faceFlux_,
schemes_.lookup(field.name())
);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "multivariateIndependentScheme.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -23,40 +23,26 @@ License
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "multivariateIndependentScheme.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline bool Foam::simpleControl::loop()
namespace Foam
{
read();
Time& time = const_cast<Time&>(mesh_.time());
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
if (initialised_)
{
if (criteriaSatisfied())
{
Info<< nl << algorithmName_ << " solution converged in "
<< time.timeName() << " iterations" << nl << endl;
defineNamedTemplateTypeNameAndDebug(multivariateIndependentScheme<scalar>, 0);
// Set to finalise calculation
time.writeAndEnd();
}
else
{
storePrevIterFields();
}
}
else
{
initialised_ = true;
storePrevIterFields();
}
multivariateSurfaceInterpolationScheme<scalar>::addIstreamConstructorToTable
<multivariateIndependentScheme<scalar> >
addMultivariateIndependentSchemeScalarConstructorToTable_;
return time.loop();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -24,7 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "VoidFraction.H"
#
// * * * * * * * * * * * * * Protectd Member Functions * * * * * * * * * * * //
template<class CloudType>

View File

@ -31,7 +31,8 @@ template<class ThermoType>
Foam::autoPtr<Foam::chemistryReader<ThermoType> >
Foam::chemistryReader<ThermoType>::New
(
const dictionary& thermoDict
const dictionary& thermoDict,
speciesTable& species
)
{
// Let the chemistry reader type default to CHEMKIN
@ -50,7 +51,7 @@ Foam::chemistryReader<ThermoType>::New
{
FatalErrorIn
(
"chemistryReader::New(const dictionary& thermoDict)"
"chemistryReader::New(const dictionary&, speciesTable&)"
) << "Unknown chemistryReader type "
<< chemistryReaderTypeName << nl << nl
<< "Valid chemistryReader types are:" << nl
@ -58,7 +59,10 @@ Foam::chemistryReader<ThermoType>::New
<< exit(FatalError);
}
return autoPtr<chemistryReader<ThermoType> >(cstrIter()(thermoDict));
return autoPtr<chemistryReader<ThermoType> >
(
cstrIter()(thermoDict, species)
);
}

View File

@ -85,16 +85,21 @@ public:
chemistryReader,
dictionary,
(
const dictionary& thermoDict
const dictionary& thermoDict,
speciesTable& species
),
(thermoDict)
(thermoDict, species)
);
// Selectors
//- Select constructed from dictionary
static autoPtr<chemistryReader> New(const dictionary& thermoDict);
static autoPtr<chemistryReader> New
(
const dictionary& thermoDict,
speciesTable& species
);
//- Destructor

View File

@ -854,29 +854,31 @@ void Foam::chemkinReader::read
Foam::chemkinReader::chemkinReader
(
const fileName& CHEMKINFileName,
speciesTable& species,
const fileName& thermoFileName
)
:
lineNo_(1),
specieNames_(10),
speciesTable_(),
speciesTable_(species),
reactions_(speciesTable_, speciesThermo_)
{
read(CHEMKINFileName, thermoFileName);
}
Foam::chemkinReader::chemkinReader(const dictionary& thermoDict)
Foam::chemkinReader::chemkinReader
(
const dictionary& thermoDict,
speciesTable& species
)
:
lineNo_(1),
specieNames_(10),
speciesTable_(),
speciesTable_(species),
reactions_(speciesTable_, speciesThermo_)
{
fileName chemkinFile
(
fileName(thermoDict.lookup("CHEMKINFile")).expand()
);
fileName chemkinFile(fileName(thermoDict.lookup("CHEMKINFile")).expand());
fileName thermoFile = fileName::null;

View File

@ -193,7 +193,7 @@ private:
HashTable<label> specieIndices_;
//- Table of species
speciesTable speciesTable_;
speciesTable& speciesTable_;
//- Specie phase
HashTable<phase> speciePhase_;
@ -318,11 +318,12 @@ public:
chemkinReader
(
const fileName& chemkinFile,
speciesTable& species,
const fileName& thermoFileName = fileName::null
);
//- Construct by getting the CHEMKIN III file name from dictionary
chemkinReader(const dictionary& thermoDict);
chemkinReader(const dictionary& thermoDict, speciesTable& species);
//- Destructor

View File

@ -27,12 +27,28 @@ License
#include "IFstream.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class ThermoType>
Foam::speciesTable& Foam::foamChemistryReader<ThermoType>::setSpecies
(
const dictionary& dict,
speciesTable& species
)
{
wordList s(dict.lookup("species"));
species.transfer(s);
return species;
}
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
template<class ThermoType>
Foam::foamChemistryReader<ThermoType>::foamChemistryReader
(
const fileName& reactionsFileName,
speciesTable& species,
const fileName& thermoFileName
)
:
@ -51,8 +67,8 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
fileName(thermoFileName).expand()
)()
),
speciesTable_(setSpecies(chemDict_, species)),
speciesThermo_(thermoDict_),
speciesTable_(chemDict_.lookup("species")),
reactions_(speciesTable_, speciesThermo_, chemDict_)
{}
@ -60,7 +76,8 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
template<class ThermoType>
Foam::foamChemistryReader<ThermoType>::foamChemistryReader
(
const dictionary& thermoDict
const dictionary& thermoDict,
speciesTable& species
)
:
chemistryReader<ThermoType>(),
@ -79,7 +96,7 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
)()
),
speciesThermo_(thermoDict_),
speciesTable_(chemDict_.lookup("species")),
speciesTable_(setSpecies(chemDict_, species)),
reactions_(speciesTable_, speciesThermo_, chemDict_)
{}

View File

@ -67,7 +67,7 @@ class foamChemistryReader
HashPtrTable<ThermoType> speciesThermo_;
//- Table of species
speciesTable speciesTable_;
speciesTable& speciesTable_;
//- List of the reactions
ReactionList<ThermoType> reactions_;
@ -75,6 +75,9 @@ class foamChemistryReader
// Private Member Functions
//- Set the species list
speciesTable& setSpecies(const dictionary& dict, speciesTable& species);
//- Disallow default bitwise copy construct
foamChemistryReader(const foamChemistryReader&);
@ -94,12 +97,17 @@ public:
foamChemistryReader
(
const fileName& reactionsFileName,
speciesTable& species,
const fileName& thermoFileName
);
//- Construct by getting the foamChemistry and thermodynamics file names
// from dictionary
foamChemistryReader(const dictionary& thermoDict);
foamChemistryReader
(
const dictionary& thermoDict,
speciesTable& species
);
//- Destructor

View File

@ -49,7 +49,7 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData
template<class ThermoType>
void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
{
// It changes Yt patches to "calculated"
// It changes Yt patches to "calculated"
volScalarField Yt("Yt", 1.0*Y_[0]);
for (label n=1; n<Y_.size(); n++)

View File

@ -35,14 +35,15 @@ Foam::reactingMixture<ThermoType>::reactingMixture
const fvMesh& mesh
)
:
speciesTable(),
autoPtr<chemistryReader<ThermoType> >
(
chemistryReader<ThermoType>::New(thermoDict)
chemistryReader<ThermoType>::New(thermoDict, *this)
),
multiComponentMixture<ThermoType>
(
thermoDict,
autoPtr<chemistryReader<ThermoType> >::operator()().species(),
*this,
autoPtr<chemistryReader<ThermoType> >::operator()().speciesThermo(),
mesh
),

View File

@ -35,6 +35,7 @@ SourceFiles
#ifndef reactingMixture_H
#define reactingMixture_H
#include "speciesTable.H"
#include "chemistryReader.H"
#include "multiComponentMixture.H"
@ -50,6 +51,7 @@ namespace Foam
template<class ThermoType>
class reactingMixture
:
public speciesTable,
public autoPtr<chemistryReader<ThermoType> >,
public multiComponentMixture<ThermoType>,
public PtrList<Reaction<ThermoType> >
@ -84,6 +86,16 @@ public:
//- Read dictionary
void read(const dictionary&);
label size() const
{
return PtrList<Reaction<ThermoType> >::size();
}
Reaction<ThermoType>& operator [] (const label i)
{
return PtrList<Reaction<ThermoType> >::operator[](i);
}
};

View File

@ -68,8 +68,8 @@ NonEquilibriumReversibleReaction
)
:
Reaction<ReactionThermo>(species, thermoDatabase, dict),
fk_(species, dict),
rk_(species, dict)
fk_(species, dict.subDict("forward")),
rk_(species, dict.subDict("reverse"))
{}
@ -136,8 +136,20 @@ void Foam::NonEquilibriumReversibleReaction<ReactionThermo, ReactionRate>::write
) const
{
Reaction<ReactionThermo>::write(os);
os << indent << "forward" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
fk_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "reverse" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
rk_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
}

View File

@ -67,10 +67,10 @@ FallOffReactionRate
const dictionary& dict
)
:
k0_(species, dict),
kInf_(species, dict),
F_(dict),
thirdBodyEfficiencies_(species, dict)
k0_(species, dict.subDict("k0")),
kInf_(species, dict.subDict("kInf")),
F_(dict.subDict("F")),
thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
{}
@ -100,10 +100,33 @@ inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::write
Ostream& os
) const
{
os << indent << "k0" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
k0_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "kInf" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
kInf_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "F" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
F_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
os << indent << "thirdBodyEfficiencies" << nl;
os << indent << token::BEGIN_BLOCK << nl;
os << incrIndent;
thirdBodyEfficiencies_.write(os);
os << decrIndent;
os << indent << token::END_BLOCK << nl;
}

View File

@ -135,6 +135,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

View File

@ -89,6 +89,16 @@ Foam::eConstThermo<EquationOfState>::New(const dictionary& dict)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState>
inline Foam::scalar Foam::eConstThermo<EquationOfState>::limit
(
const scalar T
) const
{
return T;
}
template<class EquationOfState>
inline Foam::scalar Foam::eConstThermo<EquationOfState>::cp
(

View File

@ -133,6 +133,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

View File

@ -89,6 +89,16 @@ Foam::hConstThermo<equationOfState>::New(const dictionary& dict)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState>
inline Foam::scalar Foam::hConstThermo<EquationOfState>::limit
(
const scalar T
) const
{
return T;
}
template<class equationOfState>
inline Foam::scalar Foam::hConstThermo<equationOfState>::cp
(

View File

@ -151,6 +151,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

View File

@ -82,6 +82,16 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState, int PolySize>
inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::limit
(
const scalar T
) const
{
return T;
}
template<class EquationOfState, int PolySize>
inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp
(

View File

@ -119,9 +119,6 @@ private:
//- Check that input data is valid
void checkInputData() const;
//- Check given temperature is within the range of the fitted coeffs
inline void checkT(const scalar T) const;
//- Return the coefficients corresponding to the given temperature
inline const coeffArray& coeffs(const scalar T) const;
@ -153,6 +150,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

View File

@ -52,22 +52,6 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
}
template<class EquationOfState>
inline void Foam::janafThermo<EquationOfState>::checkT(const scalar T) const
{
if (T < Tlow_ || T > Thigh_)
{
FatalErrorIn
(
"janafThermo<EquationOfState>::checkT(const scalar T) const"
) << "attempt to use janafThermo<EquationOfState>"
" out of temperature range "
<< Tlow_ << " -> " << Thigh_ << "; T = " << T
<< abort(FatalError);
}
}
template<class EquationOfState>
inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
Foam::janafThermo<EquationOfState>::coeffs
@ -75,8 +59,6 @@ Foam::janafThermo<EquationOfState>::coeffs
const scalar T
) const
{
checkT(T);
if (T < Tcommon_)
{
return lowCpCoeffs_;
@ -112,6 +94,31 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::limit
(
const scalar T
) const
{
if (T < Tlow_ || T > Thigh_)
{
WarningIn
(
"janafThermo<EquationOfState>::limit(const scalar T) const"
) << "attempt to use janafThermo<EquationOfState>"
" out of temperature range "
<< Tlow_ << " -> " << Thigh_ << "; T = " << T
<< endl;
return min(max(T, Tlow_), Thigh_);
}
else
{
return T;
}
}
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::cp
(

View File

@ -103,14 +103,15 @@ class specieThermo
// Private Member Functions
//- return the temperature corresponding to the value of the
//- Return the temperature corresponding to the value of the
// thermodynamic property f, given the function f = F(T) and dF(T)/dT
inline scalar T
(
scalar f,
scalar T0,
scalar (specieThermo::*F)(const scalar) const,
scalar (specieThermo::*dFdT)(const scalar) const
scalar (specieThermo::*dFdT)(const scalar) const,
scalar (specieThermo::*limit)(const scalar) const
) const;

View File

@ -43,7 +43,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
scalar f,
scalar T0,
scalar (specieThermo<Thermo>::*F)(const scalar) const,
scalar (specieThermo<Thermo>::*dFdT)(const scalar) const
scalar (specieThermo<Thermo>::*dFdT)(const scalar) const,
scalar (specieThermo<Thermo>::*limit)(const scalar) const
) const
{
scalar Test = T0;
@ -54,7 +55,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
do
{
Test = Tnew;
Tnew = Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test);
Tnew =
(this->*limit)(Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test));
if (iter++ > maxIter_)
{
@ -276,7 +278,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TH
const scalar T0
) const
{
return T(h, T0, &specieThermo<Thermo>::H, &specieThermo<Thermo>::Cp);
return T
(
h,
T0,
&specieThermo<Thermo>::H,
&specieThermo<Thermo>::Cp,
&specieThermo<Thermo>::limit
);
}
@ -287,7 +296,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::THs
const scalar T0
) const
{
return T(hs, T0, &specieThermo<Thermo>::Hs, &specieThermo<Thermo>::Cp);
return T
(
hs,
T0,
&specieThermo<Thermo>::Hs,
&specieThermo<Thermo>::Cp,
&specieThermo<Thermo>::limit
);
}
@ -298,7 +314,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TE
const scalar T0
) const
{
return T(e, T0, &specieThermo<Thermo>::E, &specieThermo<Thermo>::Cv);
return T
(
e,
T0,
&specieThermo<Thermo>::E,
&specieThermo<Thermo>::Cv,
&specieThermo<Thermo>::limit
);
}

View File

@ -16,6 +16,7 @@ LienCubicKELowRe/LienCubicKELowRe.C
NonlinearKEShih/NonlinearKEShih.C
LienLeschzinerLowRe/LienLeschzinerLowRe.C
LamBremhorstKE/LamBremhorstKE.C
kkLOmega/kkLOmega.C
/* Wall functions */
wallFunctions = derivedFvPatchFields/wallFunctions

View File

@ -0,0 +1,779 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "kkLOmega.H"
#include "addToRunTimeSelectionTable.H"
#include "backwardsCompatibilityWallFunctions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(kkLOmega, 0);
addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<volScalarField> kkLOmega::fv(const volScalarField& Ret) const
{
return(1.0 - exp(-sqrt(Ret)/Av_));
}
tmp<volScalarField> kkLOmega::fINT() const
{
return
(
min
(
kl_/(Cint_*(kl_ + kt_ + kMin_)),
dimensionedScalar("1.0", dimless, 1.0)
)
);
}
tmp<volScalarField> kkLOmega::fSS(const volScalarField& omega) const
{
return(exp(-sqr(Css_*nu()*omega/(kt_ + kMin_))));
}
tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
{
return(1.0/(A0_ + As_*(S/(omega_ + omegaMin_))));
}
tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& Rew) const
{
return(scalar(1.0) - exp(-sqr(max(Rew - CtsCrit_, 0.0))/Ats_));
}
tmp<volScalarField> kkLOmega::fTaul
(
const volScalarField& lambdaEff,
const volScalarField& ktL
) const
{
return
(
scalar(1.0)
- exp
(
-CtauL_*ktL
/
(
sqr
(
lambdaEff*omega_
+ dimensionedScalar
(
"ROTVSMALL",
dimLength*inv(dimTime),
ROOTVSMALL
)
)
)
)
);
}
tmp<volScalarField> kkLOmega::alphaT
(
const volScalarField& lambdaEff,
const volScalarField& fv,
const volScalarField& ktS
) const
{
return(fv*CmuStd_*sqrt(ktS)*lambdaEff);
}
tmp<volScalarField> kkLOmega::fOmega
(
const volScalarField& lambdaEff,
const volScalarField& lambdaT
) const
{
return
(
scalar(1.0)
- exp
(
-0.41
* pow4
(
lambdaEff
/ (
lambdaT
+ dimensionedScalar
(
"ROTVSMALL",
lambdaT.dimensions(),
ROOTVSMALL
)
)
)
)
);
}
tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
{
return
(
max
(
kt_/nu()
/
(
omega
+ dimensionedScalar("ROTVSMALL", omega.dimensions(), ROOTVSMALL)
)
-
CbpCrit_
,
0.0
)
);
}
tmp<volScalarField> kkLOmega::gammaNAT
(
const volScalarField& ReOmega,
const volScalarField& fNatCrit
) const
{
return
(
max
(
ReOmega
- CnatCrit_
/ (
fNatCrit + dimensionedScalar("ROTVSMALL", dimless, ROOTVSMALL)
)
,
0.0
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
kkLOmega::kkLOmega
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName,
const word& modelName
)
:
RASModel(modelName, U, phi, transport, turbulenceModelName),
A0_
(
dimensioned<scalar>::lookupOrAddToDict
(
"A0",
coeffDict_,
4.04
)
),
As_
(
dimensioned<scalar>::lookupOrAddToDict
(
"As",
coeffDict_,
2.12
)
),
Av_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Av",
coeffDict_,
6.75
)
),
Abp_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Abp",
coeffDict_,
0.6
)
),
Anat_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Anat",
coeffDict_,
200
)
),
Ats_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ats",
coeffDict_,
200
)
),
CbpCrit_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CbpCrit",
coeffDict_,
1.2
)
),
Cnc_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cnc",
coeffDict_,
0.1
)
),
CnatCrit_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CnatCrit",
coeffDict_,
1250
)
),
Cint_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cint",
coeffDict_,
0.75
)
),
CtsCrit_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CtsCrit",
coeffDict_,
1000
)
),
CrNat_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CrNat",
coeffDict_,
0.02
)
),
C11_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C11",
coeffDict_,
3.4e-6
)
),
C12_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C12",
coeffDict_,
1.0e-10
)
),
CR_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CR",
coeffDict_,
0.12
)
),
CalphaTheta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CalphaTheta",
coeffDict_,
0.035
)
),
Css_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Css",
coeffDict_,
1.5
)
),
CtauL_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CtauL",
coeffDict_,
4360
)
),
Cw1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw1",
coeffDict_,
0.44
)
),
Cw2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw2",
coeffDict_,
0.92
)
),
Cw3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw3",
coeffDict_,
0.3
)
),
CwR_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CwR",
coeffDict_,
1.5
)
),
Clambda_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Clambda",
coeffDict_,
2.495
)
),
CmuStd_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CmuStd",
coeffDict_,
0.09
)
),
Prtheta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prtheta",
coeffDict_,
0.85
)
),
Sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Sigmak",
coeffDict_,
1
)
),
Sigmaw_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Sigmaw",
coeffDict_,
1.17
)
),
kt_
(
IOobject
(
"kt",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateK("kt", mesh_)
),
omega_
(
IOobject
(
"omega",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateOmega("omega", mesh_)
),
kl_
(
IOobject
(
"kl",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateK("kl", mesh_)
),
nut_
(
IOobject
(
"nut",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateNut("nut", mesh_)
),
y_(mesh_)
{
bound(kt_, kMin_);
bound(kl_, kMin_);
bound(omega_, omegaMin_);
nut_ = kt_/(omega_ + omegaMin_);
nut_.correctBoundaryConditions();
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volSymmTensorField> kkLOmega::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*(kt_) - nut_*twoSymm(fvc::grad(U_)),
kt_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> kkLOmega::devReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-nuEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
{
return
(
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(T(fvc::grad(U))))
);
}
bool kkLOmega::read()
{
if (RASModel::read())
{
A0_.readIfPresent(coeffDict());
As_.readIfPresent(coeffDict());
Av_.readIfPresent(coeffDict());
Abp_.readIfPresent(coeffDict());
Anat_.readIfPresent(coeffDict());
Abp_.readIfPresent(coeffDict());
Ats_.readIfPresent(coeffDict());
CbpCrit_.readIfPresent(coeffDict());
Cnc_.readIfPresent(coeffDict());
CnatCrit_.readIfPresent(coeffDict());
Cint_.readIfPresent(coeffDict());
CtsCrit_.readIfPresent(coeffDict());
CrNat_.readIfPresent(coeffDict());
C11_.readIfPresent(coeffDict());
C12_.readIfPresent(coeffDict());
CR_.readIfPresent(coeffDict());
CalphaTheta_.readIfPresent(coeffDict());
Css_.readIfPresent(coeffDict());
CtauL_.readIfPresent(coeffDict());
Cw1_.readIfPresent(coeffDict());
Cw2_.readIfPresent(coeffDict());
Cw3_.readIfPresent(coeffDict());
CwR_.readIfPresent(coeffDict());
Clambda_.readIfPresent(coeffDict());
CmuStd_.readIfPresent(coeffDict());
Prtheta_.readIfPresent(coeffDict());
Sigmak_.readIfPresent(coeffDict());
Sigmaw_.readIfPresent(coeffDict());
return true;
}
else
{
return false;
}
}
void kkLOmega::correct()
{
RASModel::correct();
if (!turbulence_)
{
return;
}
if (mesh_.changing())
{
y_.correct();
y_.boundaryField() = max(y_.boundaryField(), VSMALL);
}
const volScalarField kT(kt_ + kl_);
const volScalarField lambdaT(sqrt(kT)/(omega_ + omegaMin_));
const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
const volScalarField fw
(
lambdaEff/(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL))
);
const volTensorField gradU(fvc::grad(U_));
const volScalarField omega(sqrt(2.0)*mag(skew(gradU)));
const volScalarField S2(2.0*magSqr(symm(gradU)));
const volScalarField ktS(fSS(omega)*fw*kt_);
const volScalarField nuts
(
fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
*fINT()
*Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
);
const volScalarField Pkt(nuts*S2);
const volScalarField ktL(kt_ - ktS);
const volScalarField ReOmega(sqr(y_)*omega/nu());
const volScalarField nutl
(
min
(
C11_*fTaul(lambdaEff, ktL)*omega*sqr(lambdaEff)
* sqrt(ktL)*lambdaEff/nu()
+ C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*omega
,
0.5*(kl_ + ktL)/sqrt(S2)
)
);
const volScalarField Pkl(nutl*S2);
const volScalarField alphaTEff
(
alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
);
// By pass s0urce term divided by kl_
const dimensionedScalar fwMin("SMALL", dimless, ROOTVSMALL);
const volScalarField Rbp
(
CR_*(1.0 - exp(-gammaBP(omega)()/Abp_))*omega_
/ (fw + fwMin)
);
const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
// Natural source term divided by kl_
const volScalarField Rnat
(
CrNat_*(1.0 - exp(-gammaNAT(ReOmega, fNatCrit)/Anat_))*omega
);
const volScalarField Dt(nu()*magSqr(fvc::grad(sqrt(kt_))));
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> ktEqn
(
fvm::ddt(kt_)
+ fvm::div(phi_, kt_)
- fvm::Sp(fvc::div(phi_), kt_)
- fvm::laplacian(DkEff(alphaTEff), kt_, "laplacian(alphaTEff,kt)")
==
Pkt
+ (Rbp + Rnat)*kl_
- Dt
- fvm::Sp(omega_, kt_)
);
ktEqn().relax();
ktEqn().boundaryManipulate(kt_.boundaryField());
solve(ktEqn);
bound(kt_, kMin_);
const volScalarField Dl(nu()*magSqr(fvc::grad(sqrt(kl_))));
// Laminar kinetic energy equation
tmp<fvScalarMatrix> klEqn
(
fvm::ddt(kl_)
+ fvm::div(phi_, kl_)
- fvm::Sp(fvc::div(phi_), kl_)
- fvm::laplacian(nu(), kl_, "laplacian(nu,kl)")
==
Pkl
- fvm::Sp(Rbp, kl_)
- fvm::Sp(Rnat, kl_)
- Dl
);
klEqn().relax();
klEqn().boundaryManipulate(kl_.boundaryField());
solve(klEqn);
bound(kl_, kMin_);
omega_.boundaryField().updateCoeffs();
// Turbulence specific dissipation rate equation
tmp<fvScalarMatrix> omegaEqn
(
fvm::ddt(omega_)
+ fvm::div(phi_, omega_)
- fvm::Sp(fvc::div(phi_), omega_)
- fvm::laplacian
(
DomegaEff(alphaTEff),
omega_,
"laplacian(alphaTEff,omega)"
)
==
Cw1_*Pkt*omega_/(kt_ + kMin_)
+ fvm::SuSp
(
(CwR_/(fw + fwMin) - 1.0)*kl_*(Rbp + Rnat)/(kt_ + kMin_)
, omega_
)
- fvm::Sp(Cw2_*omega_, omega_)
+ Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)/pow3(y_)
);
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
solve(omegaEqn);
bound(omega_, omegaMin_);
// Re-calculate viscosity
nut_ = nuts + nutl;
nut_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,296 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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::incompressible::RASModels::kkLOmega
Description
Low Reynolds-number k-kl-omega turbulence model for
incompressible flows.
Turbulence model described in:
\verbatim
D. Keith Walters, Davor Cokljat
"A Three-Equation Eddy-Viscosity Model for Reynold-Averaged
Navier-Stokes Simulations of Transitional Flow"
\endverbatim
The default model coefficients correspond to the following:
\verbatim
kkLOmegaCoeffs
{
A0 4.04
As 2.12
Av 6.75
Abp 0.6
Anat 200
Ats 200
CbpCrit 1.2
Cnc 0.1
CnatCrit 1250
Cint 0.75
CtsCrit 1000
CrNat 0.02
C11 3.4e-6
C12 1.0e-10
CR 0.12
CalphaTheta 0.035
Css 1.5
CtauL 4360
Cw1 0.44
Cw2 0.92
Cw3 0.3
CwR 1.5
Clambda 2.495
CmuStd 0.09
Prtheta 0.85
Sigmak 1
Sigmaw 1.17
}
\endverbatim
SourceFiles
kkLOmega.C
\*---------------------------------------------------------------------------*/
#ifndef kkLOmega_H
#define kkLOmega_H
#include "RASModel.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class kkLOmega Declaration
\*---------------------------------------------------------------------------*/
class kkLOmega
:
public RASModel
{
// Private memmber functions
tmp<volScalarField> fv(const volScalarField& Ret) const;
tmp<volScalarField> fINT() const;
tmp<volScalarField> fSS(const volScalarField& omega) const;
tmp<volScalarField> Cmu(const volScalarField& S) const;
tmp<volScalarField> BetaTS(const volScalarField& Rew) const;
tmp<volScalarField> fTaul
(
const volScalarField& lambdaEff,
const volScalarField& ktL
) const;
tmp<volScalarField> alphaT
(
const volScalarField& lambdaEff,
const volScalarField& fv,
const volScalarField& ktS
) const;
tmp<volScalarField> fOmega
(
const volScalarField& lambdaEff,
const volScalarField& lambdaT
) const;
tmp<volScalarField> gammaBP(const volScalarField& omega) const;
tmp<volScalarField> gammaNAT
(
const volScalarField& ReOmega,
const volScalarField& fNatCrit
) const;
protected:
// Protected data
// Model coefficients
dimensionedScalar A0_;
dimensionedScalar As_;
dimensionedScalar Av_;
dimensionedScalar Abp_;
dimensionedScalar Anat_;
dimensionedScalar Ats_;
dimensionedScalar CbpCrit_;
dimensionedScalar Cnc_;
dimensionedScalar CnatCrit_;
dimensionedScalar Cint_;
dimensionedScalar CtsCrit_;
dimensionedScalar CrNat_;
dimensionedScalar C11_;
dimensionedScalar C12_;
dimensionedScalar CR_;
dimensionedScalar CalphaTheta_;
dimensionedScalar Css_;
dimensionedScalar CtauL_;
dimensionedScalar Cw1_;
dimensionedScalar Cw2_;
dimensionedScalar Cw3_;
dimensionedScalar CwR_;
dimensionedScalar Clambda_;
dimensionedScalar CmuStd_;
dimensionedScalar Prtheta_;
dimensionedScalar Sigmak_;
dimensionedScalar Sigmaw_;
// Fields
volScalarField kt_;
volScalarField omega_;
volScalarField kl_;
volScalarField nut_;
wallDist y_;
public:
//- Runtime type information
TypeName("kkLOmega");
// Constructors
//- Construct from components
kkLOmega
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName = turbulenceModel::typeName,
const word& modelName = typeName
);
//- Destructor
virtual ~kkLOmega()
{}
// Member Functions
//- Return the turbulence viscosity
virtual tmp<volScalarField> nut() const
{
return nut_;
}
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& alphaT) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaT/Sigmak_ + nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& alphaT) const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaT/Sigmaw_ + nu())
);
}
//- Return the laminar kinetic energy
virtual tmp<volScalarField> kl() const
{
return kl_;
}
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return kt_;
}
//- Return the turbulence specific dissipation rate
virtual tmp<volScalarField> omega() const
{
return omega_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
mesh_.time().timeName(),
mesh_
),
kt_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const;
//- Return the effective stress tensor including the laminar stress
virtual tmp<volSymmTensorField> devReff() const;
//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
//- Read RASProperties dictionary
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //