songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

TDACChemistryModel: simplified, rationalized and automated the handling of variableTimeStep

Browse Source
This commit is contained in:
Henry Weller 2017-01-09 21:40:39 +00:00
parent d08a5193fb
commit 47bd8e13f7
13 changed files with 77 additions and 120 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/thermophysicalModels/chemistryModel/Make/options
View File

@ -7,7 +7,8 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \
-lcompressibleTransportModels \
@ -16,4 +17,5 @@ LIB_LIBS = \
-lspecie \
-lthermophysicalFunctions \
-lODE \
-lfiniteVolume
-lfiniteVolume \
-lmeshTools

12
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
View File

@ -25,6 +25,7 @@ License
#include "TDACChemistryModel.H"
#include "UniformField.H"
#include "localEulerDdtScheme.H"
#include "clockTime.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -37,6 +38,11 @@ Foam::TDACChemistryModel<CompType, ThermoType>::TDACChemistryModel
)
:
chemistryModel<CompType, ThermoType>(mesh, phaseName),
variableTimeStep_
(
mesh.time().controlDict().lookupOrDefault("adjustTimeStep", false)
|| fv::localEulerDdt::enabled(mesh)
),
timeSteps_(0),
NsDAC_(this->nSpecie_),
completeC_(this->nSpecie_, 0),
@ -598,7 +604,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
const bool reduced = mechRed_->active();
label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1:0);
label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
basicMultiComponentMixture& composition = this->thermo().composition();
@ -661,10 +667,10 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
phiq[i] = this->Y()[i][celli];
}
phiq[this->nSpecie()]=Ti;
phiq[this->nSpecie()+1]=pi;
phiq[this->nSpecie() + 1]=pi;
if (tabulation_->variableTimeStep())
{
phiq[this->nSpecie()+2] = deltaT[celli];
phiq[this->nSpecie() + 2] = deltaT[celli];
}

10
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
View File

@ -85,6 +85,8 @@ class TDACChemistryModel
{
// Private member data
bool variableTimeStep_;
label timeSteps_;
// Mechanism reduction
@ -159,11 +161,11 @@ public:
// Member Functions
inline label timeSteps() const
{
return timeSteps_;
}
//- Return true if the time-step is variable and/or non-uniform
inline bool variableTimeStep() const;
//- Return the number of chemistry evaluations (used by ISAT)
inline label timeSteps() const;
//- Create and return a TDAC log file of the given name
inline autoPtr<OFstream> logFile(const word& name) const;

18
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,22 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CompType, class ThermoType>
inline bool
Foam::TDACChemistryModel<CompType, ThermoType>::variableTimeStep() const
{
return variableTimeStep_;
}
template<class CompType, class ThermoType>
inline Foam::label
Foam::TDACChemistryModel<CompType, ThermoType>::timeSteps() const
{
return timeSteps_;
}
template<class CompType, class ThermoType>
inline Foam::autoPtr<Foam::OFstream>
Foam::TDACChemistryModel<CompType, ThermoType>::logFile(const word& name) const

23
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
View File

@ -98,30 +98,13 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
}
}
scaleFactor_[Ysize] = readScalar(scaleDict.lookup("Temperature"));
scaleFactor_[Ysize+1] = readScalar(scaleDict.lookup("Pressure"));
scaleFactor_[Ysize + 1] = readScalar(scaleDict.lookup("Pressure"));
if (this->variableTimeStep())
{
scaleFactor_[Ysize + 2] = readScalar(scaleDict.lookup("deltaT"));
}
else
{
// When the variableTimeStep option is false, if the application
// has variable time step activated, the maximum lifetime of a
// chemPoint should be 1 time step.
bool adjustTimeStep =
runTime_.controlDict().lookupOrDefault("adjustTimeStep", false);
if (chPMaxLifeTime_ > 1 && adjustTimeStep)
{
WarningInFunction
<< " variableTimeStep is not activate for ISAT while"
<< " the time step might be adjusted by the application."
<< nl
<< " This might lead to errors in the chemistry." << nl
<< " To avoid this warning either set chPMaxLifeTime to 1"
<< " or activate variableTimeStep." << endl;
}
}
}
if (this->variableTimeStep())
{
nAdditionalEqns_ = 3;
@ -433,7 +416,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
// For temperature and pressure, only unity on the diagonal
A(speciesNumber, speciesNumber) = 1;
A(speciesNumber+1, speciesNumber+1) = 1;
A(speciesNumber + 1, speciesNumber + 1) = 1;
if (this->variableTimeStep())
{
A[speciesNumber + 2][speciesNumber + 2] = 1;

66
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.C
View File

@ -28,19 +28,15 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CompType, class ThermoType>
Foam::binaryNode<CompType, ThermoType>::binaryNode
(
)
Foam::binaryNode<CompType, ThermoType>::binaryNode()
:
leafLeft_(nullptr),
leafRight_(nullptr),
nodeLeft_(nullptr),
nodeRight_(nullptr),
parent_(nullptr),
variableTimeStep_(false),
nAdditionalEqns_(0)
{
}
{}
template<class CompType, class ThermoType>
@ -56,10 +52,9 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
nodeLeft_(nullptr),
nodeRight_(nullptr),
parent_(parent),
variableTimeStep_(elementLeft->variableTimeStep()),
v_(elementLeft->completeSpaceSize(), 0)
{
if (this->variableTimeStep_)
if (elementLeft->variableTimeStep())
{
nAdditionalEqns_ = 3;
}
@ -72,41 +67,11 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
a_ = calcA(elementLeft, elementRight);
}
template<class CompType, class ThermoType>
Foam::binaryNode<CompType, ThermoType>::binaryNode
(
binaryNode<CompType, ThermoType> *bn
)
:
leafLeft_(bn->leafLeft()),
leafRight_(bn->leafRight()),
nodeLeft_(bn->nodeLeft()),
nodeRight_(bn->nodeRight()),
parent_(bn->parent()),
variableTimeStep_
(
this->coeffsDict_.lookupOrDefault("variableTimeStep", false)
),
v_(bn->v()),
a_(bn->a())
{
if (this->variableTimeStep_)
{
nAdditionalEqns_ = 3;
}
else
{
nAdditionalEqns_ = 2;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CompType, class ThermoType>
void
Foam::binaryNode<CompType, ThermoType>::calcV
void Foam::binaryNode<CompType, ThermoType>::calcV
(
chemPointISAT<CompType, ThermoType>*& elementLeft,
chemPointISAT<CompType, ThermoType>*& elementRight,
@ -116,7 +81,8 @@ Foam::binaryNode<CompType, ThermoType>::calcV
// LT is the transpose of the L matrix
scalarSquareMatrix& LT = elementLeft->LT();
bool mechReductionActive = elementLeft->chemistry().mechRed()->active();
// difference of composition in the full species domain
// Difference of composition in the full species domain
scalarField phiDif(elementRight->phi() - elementLeft->phi());
const scalarField& scaleFactor(elementLeft->scaleFactor());
scalar epsTol = elementLeft->tolerance();
@ -131,9 +97,9 @@ Foam::binaryNode<CompType, ThermoType>::calcV
if (i<elementLeft->completeSpaceSize() - nAdditionalEqns_)
{
si = elementLeft->completeToSimplifiedIndex()[i];
outOfIndexI = (si==-1);
outOfIndexI = (si == -1);
}
else// temperature and pressure
else // temperature and pressure
{
outOfIndexI = false;
const label dif =
@ -143,7 +109,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
}
if (!mechReductionActive || (mechReductionActive && !(outOfIndexI)))
{
v[i]=0.0;
v[i] = 0;
for (label j=0; j<elementLeft->completeSpaceSize(); j++)
{
label sj = j;
@ -173,7 +139,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
||(mechReductionActive && !(outOfIndexJ))
)
{
// since L is a lower triangular matrix k=0->min(i, j)
// Since L is a lower triangular matrix k=0->min(i, j)
for (label k=0; k<=min(si, sj); k++)
{
v[i] += LT(k, si)*LT(k, sj)*phiDif[j];
@ -183,8 +149,8 @@ Foam::binaryNode<CompType, ThermoType>::calcV
}
else
{
// when it is an inactive species the diagonal element of LT is
// 1/(scaleFactor*epsTol)
// When it is an inactive species the diagonal element of LT is
// 1/(scaleFactor*epsTol)
v[i] = phiDif[i]/sqr(scaleFactor[i]*epsTol);
}
}
@ -198,13 +164,13 @@ Foam::scalar Foam::binaryNode<CompType, ThermoType>::calcA
chemPointISAT<CompType, ThermoType>* elementRight
)
{
scalar a = 0.0;
scalarField phih((elementLeft->phi()+elementRight->phi())/2);
label completeSpaceSize = elementLeft->completeSpaceSize();
for (label i=0; i<completeSpaceSize; i++)
scalarField phih((elementLeft->phi() + elementRight->phi())/2);
scalar a = 0;
forAll(phih, i)
{
a += v_[i]*phih[i];
}
return a;
}

8
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.H
View File

@ -67,9 +67,6 @@ public:
//- Parent node
binaryNode<CompType, ThermoType>* parent_;
//- Switch to allow variable time step (off by default)
bool variableTimeStep_;
//- Number of equations in addition to the species eqs.
label nAdditionalEqns_;
@ -137,11 +134,6 @@ public:
chemPointISAT<CompType, ThermoType>* elementRight,
binaryNode<CompType, ThermoType>* parent
);
//- Construct from another binary node
binaryNode
(
binaryNode<CompType, ThermoType> *bn
);
// Member functions

34
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
View File

@ -227,18 +227,14 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
printProportion_(coeffsDict.lookupOrDefault("printProportion",false)),
numRetrieve_(0),
nLifeTime_(0),
variableTimeStep_
(
coeffsDict.lookupOrDefault("variableTimeStep", false)
),
completeToSimplifiedIndex_
(
completeSpaceSize - (2 + (variableTimeStep_ == 1 ? 1 : 0))
completeSpaceSize - (2 + (variableTimeStep() == 1 ? 1 : 0))
)
{
tolerance_=tolerance;
tolerance_ = tolerance;
if (this->variableTimeStep_)
if (variableTimeStep())
{
nAdditionalEqns_ = 3;
iddeltaT_ = completeSpaceSize - 1;
@ -342,12 +338,11 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
maxNumNewDim_(p.maxNumNewDim()),
numRetrieve_(0),
nLifeTime_(0),
variableTimeStep_(p.variableTimeStep()),
completeToSimplifiedIndex_(p.completeToSimplifiedIndex())
{
tolerance_ = p.tolerance();
if (this->variableTimeStep_)
if (variableTimeStep())
{
nAdditionalEqns_ = 3;
idT_ = completeSpaceSize() - 3;
@ -407,7 +402,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
temp += LT_(si, dim)*dphi[idT_];
temp += LT_(si, dim+1)*dphi[idp_];
if (variableTimeStep_)
if (variableTimeStep())
{
temp += LT_(si, dim+2)*dphi[iddeltaT_];
}
@ -426,7 +421,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
}
// Temperature
if (variableTimeStep_)
if (variableTimeStep())
{
epsTemp +=
sqr
@ -447,7 +442,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
}
// Pressure
if (variableTimeStep_)
if (variableTimeStep())
{
epsTemp +=
sqr
@ -461,7 +456,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
epsTemp += sqr(LT_(dim+1, dim+1)*dphi[idp_]);
}
if (variableTimeStep_)
if (variableTimeStep())
{
epsTemp += sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
}
@ -477,7 +472,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
propEps[idp_] =
sqr(LT_(dim+1, dim+1)*dphi[idp_]);
if (variableTimeStep_)
if (variableTimeStep())
{
propEps[iddeltaT_] =
sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
@ -572,7 +567,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
}
dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
if (variableTimeStep_)
if (variableTimeStep())
{
dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
}
@ -719,7 +714,8 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
A_(nActiveSpecies_+1, nActiveSpecies_+1)=
Avar(initNActiveSpecies+1, initNActiveSpecies+1);
if (variableTimeStep_)
if (variableTimeStep())
{
LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
@ -755,9 +751,11 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
}
phiTilde[i] += LT_(i, j)*dphi[sj];
}
phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_];
phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_];
if (variableTimeStep_)
if (variableTimeStep())
{
phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_];
}
@ -827,7 +825,7 @@ simplifiedToCompleteIndex
{
return completeSpaceSize_-nAdditionalEqns_ + 1;
}
else if (variableTimeStep_ && (i == nActiveSpecies_ + 2))
else if (variableTimeStep() && (i == nActiveSpecies_ + 2))
{
return completeSpaceSize_-nAdditionalEqns_ + 2;
}

7
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
View File

@ -198,9 +198,6 @@ class chemPointISAT
// to still live according to the maxChPLifeTime_ parameter
label nLifeTime_;
//- Switch to allow variable time step (off by default)
Switch variableTimeStep_;
List<label> completeToSimplifiedIndex_;
//- Number of equations in addition to the species eqs.
@ -397,9 +394,9 @@ public:
return nLifeTime_;
}
inline Switch variableTimeStep()
inline bool variableTimeStep() const
{
return variableTimeStep_;
return chemistry_.variableTimeStep();
}
// ISAT functions

3
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
View File

@ -40,8 +40,7 @@ Foam::chemistryTabulationMethod<CompType, ThermoType>::chemistryTabulationMethod
active_(coeffsDict_.lookupOrDefault<Switch>("active", false)),
log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
chemistry_(chemistry),
tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4)),
variableTimeStep_(coeffsDict_.lookupOrDefault("variableTimeStep", false))
tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4))
{}

4
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
View File

@ -74,8 +74,6 @@ protected:
scalar tolerance_;
//- Switch to allow variable time step (off by default)
const Switch variableTimeStep_;
public:
@ -134,7 +132,7 @@ public:
inline bool variableTimeStep()
{
return variableTimeStep_;
return chemistry_.variableTimeStep();
}
inline scalar tolerance() const

4
src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/makeChemistryTabulationMethods.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -49,7 +49,7 @@ License
add##chemistryTabulationMethods##SS##Comp##Thermo##ConstructorToTable_;
#define makeChemistryTabulationMethods(CompChemModel, Thermo) \
#define makeChemistryTabulationMethods(CompChemModel, Thermo) \
\
typedef chemistryTabulationMethod<CompChemModel, Thermo> \
chemistryTabulationMethod##CompChemModel##Thermo; \

2
tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties
View File

@ -78,8 +78,6 @@ tabulation
printNumRetrieve off;
variableTimeStep on;
// Tolerance used for retrieve and grow
tolerance 1e-3;