chemistryModel: Remove support for the sequential solver and rationalise EulerImplicit

This commit is contained in:
Henry
2013-10-02 08:37:55 +01:00
parent da7206c69e
commit 9d45269abc
18 changed files with 95 additions and 533 deletions

View File

@ -154,70 +154,6 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omegaI
}
template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::updateConcsInReactionI
(
const label index,
const scalar dt,
const scalar omeg,
const scalar p,
const scalar T,
scalarField& c
) const
{
// update species
const Reaction<ThermoType>& R = reactions_[index];
forAll(R.lhs(), s)
{
label si = R.lhs()[s].index;
scalar sl = R.lhs()[s].stoichCoeff;
c[si] -= dt*sl*omeg;
c[si] = max(0.0, c[si]);
}
forAll(R.rhs(), s)
{
label si = R.rhs()[s].index;
scalar sr = R.rhs()[s].stoichCoeff;
c[si] += dt*sr*omeg;
c[si] = max(0.0, c[si]);
}
}
template<class CompType, class ThermoType>
void Foam::chemistryModel<CompType, ThermoType>::updateRRInReactionI
(
const label index,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const
{
const Reaction<ThermoType>& R = reactions_[index];
forAll(R.lhs(), s)
{
label si = R.lhs()[s].index;
scalar sl = R.lhs()[s].stoichCoeff;
RR[si][rRef] -= sl*pr*corr;
RR[si][lRef] += sl*pf*corr;
}
forAll(R.rhs(), s)
{
label si = R.rhs()[s].index;
scalar sr = R.rhs()[s].stoichCoeff;
RR[si][lRef] -= sr*pf*corr;
RR[si][rRef] += sr*pr*corr;
}
}
template<class CompType, class ThermoType>
Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omega
(
@ -513,7 +449,7 @@ void Foam::chemistryModel<CompType, ThermoType>::jacobian
}
}
// calculate the dcdT elements numerically
// Calculate the dcdT elements numerically
const scalar delta = 1.0e-3;
const scalarField dcdT0(omega(c2, T - delta, p));
const scalarField dcdT1(omega(c2, T + delta, p));
@ -522,7 +458,6 @@ void Foam::chemistryModel<CompType, ThermoType>::jacobian
{
dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
}
}

View File

@ -182,32 +182,6 @@ public:
//- Calculates the reaction rates
virtual void calculate();
//- Update concentrations in reaction i given dt and reaction rate omega
// used by sequential solver
void updateConcsInReactionI
(
const label i,
const scalar dt,
const scalar omega,
const scalar p,
const scalar T,
scalarField& c
) const;
//- Update matrix RR for reaction i. Used by EulerImplicit
void updateRRInReactionI
(
const label i,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const;
// Chemistry model functions (overriding abstract functions in
// basicChemistryModel.H)

View File

@ -52,6 +52,41 @@ Foam::EulerImplicit<ChemistryModel>::~EulerImplicit()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ChemistryModel>
void Foam::EulerImplicit<ChemistryModel>::updateRRInReactionI
(
const label index,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const
{
const Reaction<typename ChemistryModel::thermoType>& R =
this->reactions_[index];
forAll(R.lhs(), s)
{
label si = R.lhs()[s].index;
scalar sl = R.lhs()[s].stoichCoeff;
RR[si][rRef] -= sl*pr*corr;
RR[si][lRef] += sl*pf*corr;
}
forAll(R.rhs(), s)
{
label si = R.rhs()[s].index;
scalar sr = R.rhs()[s].stoichCoeff;
RR[si][lRef] -= sr*pf*corr;
RR[si][rRef] += sr*pr*corr;
}
}
template<class ChemistryModel>
void Foam::EulerImplicit<ChemistryModel>::solve
(
@ -62,8 +97,6 @@ void Foam::EulerImplicit<ChemistryModel>::solve
scalar& subDeltaT
) const
{
deltaT = min(deltaT, subDeltaT);
const label nSpecie = this->nSpecie();
simpleMatrix<scalar> RR(nSpecie, 0, 0);
@ -84,10 +117,7 @@ void Foam::EulerImplicit<ChemistryModel>::solve
}
scalar ha = mixture.Ha(p, T);
for (label i=0; i<nSpecie; i++)
{
RR.source()[i] = c[i]/deltaT;
}
scalar deltaTEst = min(deltaT, subDeltaT);
forAll(this->reactions(), i)
{
@ -101,24 +131,54 @@ void Foam::EulerImplicit<ChemistryModel>::solve
{
if (omegai < 0.0)
{
corr = 1.0/(1.0 + pr*deltaT);
corr = 1.0/(1.0 + pr*deltaTEst);
}
else
{
corr = 1.0/(1.0 + pf*deltaT);
corr = 1.0/(1.0 + pf*deltaTEst);
}
}
this->updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR);
updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR);
}
// Calculate the stable/accurate time-step
scalar tMin = GREAT;
for (label i=0; i<nSpecie; i++)
{
RR[i][i] += 1.0/deltaT;
scalar d = 0;
for (label j=0; j<nSpecie; j++)
{
d -= RR[i][j]*c[j];
}
if (d < -SMALL)
{
tMin = min(tMin, -(c[i] + SMALL)/d);
}
else
{
d = max(d, SMALL);
scalar cm = max(cTot - c[i], 1.0e-5);
tMin = min(tMin, cm/d);
}
}
subDeltaT = cTauChem_*tMin;
deltaT = min(deltaT, subDeltaT);
// Add the diagonal and source contributions from the time-derivative
for (label i=0; i<nSpecie; i++)
{
RR[i][i] += 1.0/deltaT;
RR.source()[i] = c[i]/deltaT;
}
// Solve for the new composition
c = RR.LUsolve();
// Limit the composition
for (label i=0; i<nSpecie; i++)
{
c[i] = max(0.0, c[i]);
@ -133,39 +193,14 @@ void Foam::EulerImplicit<ChemistryModel>::solve
}
T = mixture.THa(ha, p, T);
// Estimate the next time step
scalar tMin = GREAT;
const label nEqns = this->nEqns();
/*
for (label i=0; i<nSpecie; i++)
{
cTp_[i] = c[i];
}
cTp_[nSpecie] = T;
cTp_[nSpecie+1] = p;
scalarField dcdt(nEqns, 0.0);
this->derivatives(0.0, cTp_, dcdt);
const scalar sumC = sum(c);
for (label i=0; i<nSpecie; i++)
{
scalar d = dcdt[i];
if (d < -SMALL)
{
tMin = min(tMin, -(c[i] + SMALL)/d);
}
else
{
d = max(d, SMALL);
scalar cm = max(sumC - c[i], 1.0e-5);
tMin = min(tMin, cm/d);
}
}
subDeltaT = cTauChem_*tMin;
*/
}

View File

@ -87,6 +87,19 @@ public:
// Member Functions
void updateRRInReactionI
(
const label index,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const;
//- Update the concentrations and return the chemical time
virtual void solve
(

View File

@ -31,7 +31,6 @@ License
#include "chemistryModel.H"
#include "noChemistrySolver.H"
#include "sequential.H"
#include "EulerImplicit.H"
#include "ode.H"
@ -79,13 +78,6 @@ License
CompChemModel, \
Thermo \
); \
\
makeChemistrySolverType \
( \
sequential, \
CompChemModel, \
Thermo \
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -1,118 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 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 "sequential.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class ChemistryModel>
Foam::sequential<ChemistryModel>::sequential
(
const fvMesh& mesh
)
:
chemistrySolver<ChemistryModel>(mesh),
coeffsDict_(this->subDict("sequentialCoeffs")),
cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class ChemistryModel>
Foam::sequential<ChemistryModel>::~sequential()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ChemistryModel>
void Foam::sequential<ChemistryModel>::solve
(
scalarField& c,
scalar& T,
scalar& p,
scalar& deltaT,
scalar& subDeltaT
) const
{
deltaT = min(deltaT, subDeltaT);
const label nSpecie = this->nSpecie();
// Calculate the absolute enthalpy
scalar cTot = sum(c);
typename ChemistryModel::thermoType mixture
(
(c[0]/cTot)*this->specieThermo_[0]
);
for (label i=1; i<nSpecie; i++)
{
mixture += (c[i]/cTot)*this->specieThermo_[i];
}
scalar ha = mixture.Ha(p, T);
scalar tChemInv = SMALL;
forAll(this->reactions(), i)
{
scalar pf, cf, pb, cb;
label lRef, rRef;
scalar omega = this->omegaI(i, c, T, p, pf, cf, lRef, pb, cb, rRef);
if (eqRateLimiter_)
{
if (omega < 0.0)
{
omega /= 1.0 + pb*deltaT;
}
else
{
omega /= 1.0 + pf*deltaT;
}
}
tChemInv = max(tChemInv, mag(omega));
this->updateConcsInReactionI(i, deltaT, omega, p, T, c);
}
// Update the temperature
cTot = sum(c);
mixture = (c[0]/cTot)*this->specieThermo_[0];
for (label i=1; i<nSpecie; i++)
{
mixture += (c[i]/cTot)*this->specieThermo_[i];
}
T = mixture.THa(ha, p, T);
subDeltaT = cTauChem_/tChemInv;
}
// ************************************************************************* //

View File

@ -1,112 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 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::sequential
Description
Foam::sequential
SourceFiles
sequential.C
\*---------------------------------------------------------------------------*/
#ifndef sequential_H
#define sequential_H
#include "chemistrySolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sequential Declaration
\*---------------------------------------------------------------------------*/
template<class ChemistryModel>
class sequential
:
public chemistrySolver<ChemistryModel>
{
// Private data
//- Coefficients dictionary
dictionary coeffsDict_;
// Model constants
//- Chemistry time scale
scalar cTauChem_;
//- Equilibrium rate limiter flag (on/off)
Switch eqRateLimiter_;
public:
//- Runtime type information
TypeName("sequential");
// Constructors
//- Construct from mesh
sequential(const fvMesh& mesh);
//- Destructor
virtual ~sequential();
// Member Functions
//- Update the concentrations and return the chemical time
virtual void solve
(
scalarField& c,
scalar& T,
scalar& p,
scalar& deltaT,
scalar& subDeltaT
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "sequential.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -512,70 +512,6 @@ calculate()
}
template<class CompType, class SolidThermo, class GasThermo>
void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
updateConcsInReactionI
(
const label index,
const scalar dt,
const scalar omeg,
const scalar p,
const scalar T,
scalarField& c
) const
{
// update species
const Reaction<SolidThermo>& R = this->reactions_[index];
scalar rhoL = 0.0;
forAll(R.lhs(), s)
{
label si = R.lhs()[s].index;
rhoL = this->solidThermo_[si].rho(p, T);
c[si] -= dt*omeg;
c[si] = max(0.0, c[si]);
}
scalar sr = 0.0;
forAll(R.rhs(), s)
{
label si = R.rhs()[s].index;
const scalar rhoR = this->solidThermo_[si].rho(p, T);
sr = rhoR/rhoL;
c[si] += dt*sr*omeg;
c[si] = max(0.0, c[si]);
}
forAll(R.grhs(), g)
{
label gi = R.grhs()[g].index;
c[gi + this->nSolids_] += (1.0 - sr)*omeg*dt;
}
}
template<class CompType, class SolidThermo, class GasThermo>
void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
updateRRInReactionI
(
const label index,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const
{
notImplemented
(
"void Foam::pyrolysisChemistryModel<CompType, SolidThermo,GasThermo>::"
"updateRRInReactionI"
);
}
template<class CompType, class SolidThermo, class GasThermo>
Foam::scalar
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve

View File

@ -176,34 +176,6 @@ public:
//- Calculates the reaction rates
virtual void calculate();
//- Update concentrations in reaction i given dt and reaction rate
// omega used by sequential solver
virtual void updateConcsInReactionI
(
const label i,
const scalar dt,
const scalar omega,
const scalar p,
const scalar T,
scalarField& c
) const;
//- Update matrix RR for reaction i. Used by EulerImplicit
// (Not implemented)
virtual void updateRRInReactionI
(
const label i,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const;
// Chemistry model functions

View File

@ -175,34 +175,6 @@ public:
virtual void calculate() = 0;
//- Update concentrations in reaction i given dt and reaction rate
// omega used by sequential solver
virtual void updateConcsInReactionI
(
const label i,
const scalar dt,
const scalar omega,
const scalar p,
const scalar T,
scalarField& c
) const = 0;
//- Update matrix RR for reaction i. Used by EulerImplicit
virtual void updateRRInReactionI
(
const label i,
const scalar pr,
const scalar pf,
const scalar corr,
const label lRef,
const label rRef,
const scalar p,
const scalar T,
simpleMatrix<scalar>& RR
) const = 0;
// Solid Chemistry model functions
//- Return const access to the chemical source terms for solids

View File

@ -34,7 +34,6 @@ Description
#include "noChemistrySolver.H"
#include "EulerImplicit.H"
#include "sequential.H"
#include "ode.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -84,15 +83,7 @@ namespace Foam
SThermo, \
GThermo \
); \
\
makeSolidChemistrySolverType \
( \
sequential, \
SolidChem, \
Comp, \
SThermo, \
GThermo \
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -25,14 +25,9 @@ chemistry on;
initialChemicalTimeStep 1e-07;
sequentialCoeffs
{
cTauChem 0.001;
}
EulerImplicitCoeffs
{
cTauChem 0.05;
cTauChem 1;
equilibriumRateLimiter off;
}

View File

@ -23,7 +23,13 @@ chemistryType
chemistry on;
initialChemicalTimeStep 1e-10;
initialChemicalTimeStep 1e-6;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs
{

View File

@ -25,11 +25,6 @@ chemistry on;
initialChemicalTimeStep 1e-07;
sequentialCoeffs
{
cTauChem 0.05;
}
EulerImplicitCoeffs
{
cTauChem 1;

View File

@ -25,12 +25,6 @@ chemistry on;
initialChemicalTimeStep 1e-07;
sequentialCoeffs
{
cTauChem 0.001;
equilibriumRateLimiter off;
}
EulerImplicitCoeffs
{
cTauChem 0.05;

View File

@ -25,12 +25,6 @@ chemistry off;
initialChemicalTimeStep 1e-07;
sequentialCoeffs
{
cTauChem 0.001;
equilibriumRateLimiter off;
}
EulerImplicitCoeffs
{
cTauChem 0.05;

View File

@ -25,12 +25,6 @@ chemistry off;
initialChemicalTimeStep 1e-07;
sequentialCoeffs
{
cTauChem 0.001;
equilibriumRateLimiter off;
}
EulerImplicitCoeffs
{
cTauChem 0.05;

View File

@ -25,12 +25,6 @@ chemistry off;
initialChemicalTimeStep 1e-07;
sequentialCoeffs
{
cTauChem 0.001;
equilibriumRateLimiter off;
}
EulerImplicitCoeffs
{
cTauChem 0.05;