zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

chemistrySolver::EulerImplicit: Updated to use the StandardChemistryModel reaction Jacobian

Browse Source
This commit is contained in:
Henry Weller
2020-07-29 19:09:40 +01:00
parent 6637ed0e9f
commit bddd829fc2
17 changed files with 70 additions and 153 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
View File

@ -153,7 +153,7 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
template<class ReactionThermo, class ThermoType>
void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
(
const scalar time,
const scalar t,
const scalarField& c,
const label li,
scalarField& dcdt
@ -186,15 +186,15 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
}
cp /= rho;
scalar dT = 0;
scalar dTdt = 0;
for (label i = 0; i < nSpecie_; i++)
{
const scalar hi = specieThermos_[i].ha(p, T);
dT += hi*dcdt[i];
dTdt += hi*dcdt[i];
}
dT /= rho*cp;
dTdt /= rho*cp;
dcdt[nSpecie_] = -dT;
dcdt[nSpecie_] = -dTdt;
// dp/dt = ...
dcdt[nSpecie_ + 1] = 0;
@ -257,6 +257,11 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
}
dTdt /= -cpMean; // K/s
dcdt[nSpecie_] = dTdt;
// dp/dt = ...
dcdt[nSpecie_ + 1] = 0;
for (label i = 0; i < nSpecie_; i++)
{
J(nSpecie_, i) = 0;

3
src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H
View File

@ -71,9 +71,6 @@ class StandardChemistryModel
protected:
typedef ThermoType thermoType;
// Protected data
//- Reference to the field of specie mass fractions

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

@ -612,8 +612,8 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
c0[i] = c[i];
phiq[i] = this->Y()[i][celli];
}
phiq[this->nSpecie()]=Ti;
phiq[this->nSpecie() + 1]=pi;
phiq[this->nSpecie()] = Ti;
phiq[this->nSpecie() + 1] = pi;
if (tabulation_->variableTimeStep())
{
phiq[this->nSpecie() + 2] = deltaT[celli];

143
src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
View File

@ -25,8 +25,6 @@ License
#include "EulerImplicit.H"
#include "addToRunTimeSelectionTable.H"
#include "simpleMatrix.H"
#include "Reaction.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -39,8 +37,10 @@ Foam::EulerImplicit<ChemistryModel>::EulerImplicit
chemistrySolver<ChemistryModel>(thermo),
coeffsDict_(this->subDict("EulerImplicitCoeffs")),
cTauChem_(coeffsDict_.lookup<scalar>("cTauChem")),
eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")),
cTp_(this->nEqns())
cTp_(this->nEqns()),
R_(this->nEqns()),
J_(this->nEqns()),
E_(this->nEqns() - 2)
{}
@ -53,41 +53,6 @@ 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)
{
const label si = R.lhs()[s].index;
const scalar sl = R.lhs()[s].stoichCoeff;
RR[si][rRef] -= sl*pr*corr;
RR[si][lRef] += sl*pf*corr;
}
forAll(R.rhs(), s)
{
const label si = R.rhs()[s].index;
const 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
(
@ -100,101 +65,71 @@ void Foam::EulerImplicit<ChemistryModel>::solve
) const
{
const label nSpecie = this->nSpecie();
simpleMatrix<scalar> RR(nSpecie, 0, 0);
for (label i=0; i<nSpecie; i++)
// Map the composition, temperature and pressure into cTp
for (int i=0; i<nSpecie; i++)
{
c[i] = max(0, c[i]);
cTp_[i] = max(0, c[i]);
}
cTp_[nSpecie] = T;
cTp_[nSpecie + 1] = p;
// Calculate the absolute enthalpy
const scalar cTot = sum(c);
typename ChemistryModel::thermoType mixture
(
(this->specieThermos_[0].W()*c[0])*this->specieThermos_[0]
);
for (label i=1; i<nSpecie; i++)
{
mixture += (this->specieThermos_[i].W()*c[i])*this->specieThermos_[i];
}
const scalar ha = mixture.Ha(p, T);
const scalar deltaTEst = min(deltaT, subDeltaT);
forAll(this->reactions(), i)
{
scalar pf, cf, pr, cr;
label lRef, rRef;
const scalar omegai
(
this->omegaI(i, p, T, c, li, pf, cf, lRef, pr, cr, rRef)
);
scalar corr = 1;
if (eqRateLimiter_)
{
if (omegai < 0)
{
corr = 1/(1 + pr*deltaTEst);
}
else
{
corr = 1/(1 + pf*deltaTEst);
}
}
updateRRInReactionI(i, pr, pf, corr, lRef, rRef, p, T, RR);
}
// Calculate the reaction rate and Jacobian
this->jacobian(0, cTp_, li, R_, J_);
// Calculate the stable/accurate time-step
scalar tMin = great;
const scalar cTot = sum(c);
for (label i=0; i<nSpecie; i++)
{
scalar d = 0;
for (label j=0; j<nSpecie; j++)
if (R_[i] < -small)
{
d -= RR(i, j)*c[j];
}
if (d < -small)
{
tMin = min(tMin, -(c[i] + small)/d);
tMin = min(tMin, -(cTp_[i] + small)/R_[i]);
}
else
{
d = max(d, small);
const scalar cm = max(cTot - c[i], 1e-5);
tMin = min(tMin, cm/d);
tMin = min
(
tMin,
max(cTot - cTp_[i], 1e-5)/max(R_[i], small)
);
}
}
subDeltaT = cTauChem_*tMin;
deltaT = min(deltaT, subDeltaT);
// Add the diagonal and source contributions from the time-derivative
// Assemble the Euler implicit matrix for the composition
scalarField& source = E_.source();
for (label i=0; i<nSpecie; i++)
{
RR(i, i) += 1/deltaT;
RR.source()[i] = c[i]/deltaT;
E_(i, i) = 1/deltaT - J_(i, i);
source[i] = R_[i] + E_(i, i)*cTp_[i];
for (label j=0; j<nSpecie; j++)
{
if (i != j)
{
E_(i, j) = -J_(i, j);
source[i] += E_(i, j)*cTp_[j];
}
}
}
// Solve for the new composition
c = RR.LUsolve();
scalarField::subField(cTp_, nSpecie) = E_.LUsolve();
// Limit the composition
// Limit the composition and transfer back into c
for (label i=0; i<nSpecie; i++)
{
c[i] = max(0, c[i]);
c[i] = max(0, cTp_[i]);
}
// Update the temperature
mixture = (this->specieThermos_[0].W()*c[0])*this->specieThermos_[0];
for (label i=1; i<nSpecie; i++)
{
mixture += (this->specieThermos_[i].W()*c[i])*this->specieThermos_[i];
}
T = mixture.THa(ha, p, T);
// Euler explicit integrate the temperature.
// Separating the integration of temperature from composition
// is significantly more stable for exothermic systems
T += deltaT*R_[nSpecie];
}

46
src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,6 +27,11 @@ Class
Description
An Euler implicit solver for chemistry
Euler implicit integration based on the reaction rate Jacobian is used to
solver for the composition and Euler explicit integration for the
temperature. Separating the integration of temperature from composition
in this manner is significantly more stable for exothermic systems
SourceFiles
EulerImplicit.C
@ -36,16 +41,13 @@ SourceFiles
#define EulerImplicit_H
#include "chemistrySolver.H"
#include "Switch.H"
#include "simpleMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template <class Type>
class simpleMatrix;
/*---------------------------------------------------------------------------*\
Class EulerImplicit Declaration
\*---------------------------------------------------------------------------*/
@ -60,18 +62,21 @@ class EulerImplicit
//- Coefficients dictionary
dictionary coeffsDict_;
//- Chemistry timescale coefficient
scalar cTauChem_;
// Model constants
//- Chemistry timescale
scalar cTauChem_;
//- Equilibrium rate limiter flag (on/off)
Switch eqRateLimiter_;
// Solver data
//- Field encapsulating the composition, temperature and pressure
mutable scalarField cTp_;
//- Reaction rate field
mutable scalarField R_;
//- Reaction Jacobian
mutable scalarSquareMatrix J_;
//- Euler implicit integration matrix for composition
mutable simpleMatrix<scalar> E_;
public:
@ -91,19 +96,6 @@ 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
(

1
tutorials/combustion/chemFoam/gri/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/combustion/chemFoam/h2/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 10;
equilibriumRateLimiter no;
}
odeCoeffs

1
tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/gas/chemistryProperties
View File

@ -33,7 +33,6 @@ odeCoeffs
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
#include "reactions"

1
tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 0.05;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/multiphase/multiphaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/constant/chemistryProperties.vapor
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs

1
tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/constant/chemistryProperties.vapor
View File

@ -27,7 +27,6 @@ initialChemicalTimeStep 1e-07;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs