diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C index 0e663c1bbc..6d5ad3d562 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C @@ -154,70 +154,6 @@ Foam::scalar Foam::chemistryModel::omegaI } -template -void Foam::chemistryModel::updateConcsInReactionI -( - const label index, - const scalar dt, - const scalar omeg, - const scalar p, - const scalar T, - scalarField& c -) const -{ - // update species - const Reaction& 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 -void Foam::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& RR -) const -{ - const Reaction& 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 Foam::scalar Foam::chemistryModel::omega ( @@ -513,7 +449,7 @@ void Foam::chemistryModel::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::jacobian { dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta; } - } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H index 9bec27b25b..6da538d3d9 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H @@ -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& RR - ) const; - // Chemistry model functions (overriding abstract functions in // basicChemistryModel.H) diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C index 76c93cb0ac..9c8fff4c3c 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C @@ -52,6 +52,41 @@ Foam::EulerImplicit::~EulerImplicit() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +template +void Foam::EulerImplicit::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& RR +) const +{ + const Reaction& 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 void Foam::EulerImplicit::solve ( @@ -62,8 +97,6 @@ void Foam::EulerImplicit::solve scalar& subDeltaT ) const { - deltaT = min(deltaT, subDeltaT); - const label nSpecie = this->nSpecie(); simpleMatrix RR(nSpecie, 0, 0); @@ -84,10 +117,7 @@ void Foam::EulerImplicit::solve } scalar ha = mixture.Ha(p, T); - for (label i=0; ireactions(), i) { @@ -101,24 +131,54 @@ void Foam::EulerImplicit::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::solve } T = mixture.THa(ha, p, T); - // Estimate the next time step - scalar tMin = GREAT; - const label nEqns = this->nEqns(); - + /* for (label i=0; iderivatives(0.0, cTp_, dcdt); - - const scalar sumC = sum(c); - - for (label i=0; i& RR + ) const; + //- Update the concentrations and return the chemical time virtual void solve ( diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H index 8ad6b11f6e..272b78b4b3 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H @@ -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 \ - ); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C deleted file mode 100644 index ce6316a51e..0000000000 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C +++ /dev/null @@ -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 . - -\*---------------------------------------------------------------------------*/ - -#include "sequential.H" -#include "addToRunTimeSelectionTable.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -Foam::sequential::sequential -( - const fvMesh& mesh -) -: - chemistrySolver(mesh), - coeffsDict_(this->subDict("sequentialCoeffs")), - cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))), - eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -template -Foam::sequential::~sequential() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -void Foam::sequential::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; ispecieThermo_[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; ispecieThermo_[i]; - } - T = mixture.THa(ha, p, T); - - subDeltaT = cTauChem_/tChemInv; -} - - -// ************************************************************************* // diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H deleted file mode 100644 index 8a343b6791..0000000000 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H +++ /dev/null @@ -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 . - -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 sequential -: - public chemistrySolver -{ - // 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 - -// ************************************************************************* // diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C index f5cb38648b..6752614f72 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C @@ -512,70 +512,6 @@ calculate() } -template -void Foam::pyrolysisChemistryModel:: -updateConcsInReactionI -( - const label index, - const scalar dt, - const scalar omeg, - const scalar p, - const scalar T, - scalarField& c -) const -{ - // update species - const Reaction& 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 -void Foam::pyrolysisChemistryModel:: -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& RR -) const -{ - notImplemented - ( - "void Foam::pyrolysisChemistryModel::" - "updateRRInReactionI" - ); -} - - template Foam::scalar Foam::pyrolysisChemistryModel::solve diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H index fd82b4d09a..0b757eaaf8 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H @@ -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& RR - ) const; - // Chemistry model functions diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H index 3eaa2cdf9b..cf777f549e 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H @@ -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& RR - ) const = 0; - - // Solid Chemistry model functions //- Return const access to the chemical source terms for solids diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H index 92b7b7dbff..915f77e950 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H @@ -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 \ - ); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties index 2abcf4ae70..92c8953315 100644 --- a/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties +++ b/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties @@ -25,14 +25,9 @@ chemistry on; initialChemicalTimeStep 1e-07; -sequentialCoeffs -{ - cTauChem 0.001; -} - EulerImplicitCoeffs { - cTauChem 0.05; + cTauChem 1; equilibriumRateLimiter off; } diff --git a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties index 0bd6848bf9..f4dc73c192 100644 --- a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties @@ -23,7 +23,13 @@ chemistryType chemistry on; -initialChemicalTimeStep 1e-10; +initialChemicalTimeStep 1e-6; + +EulerImplicitCoeffs +{ + cTauChem 1; + equilibriumRateLimiter off; +} odeCoeffs { diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties index a5a53ae5b7..79a1768caf 100644 --- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties +++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties @@ -25,11 +25,6 @@ chemistry on; initialChemicalTimeStep 1e-07; -sequentialCoeffs -{ - cTauChem 0.05; -} - EulerImplicitCoeffs { cTauChem 1; diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties index c0b596bb36..21cf68db59 100644 --- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties +++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties @@ -25,12 +25,6 @@ chemistry on; initialChemicalTimeStep 1e-07; -sequentialCoeffs -{ - cTauChem 0.001; - equilibriumRateLimiter off; -} - EulerImplicitCoeffs { cTauChem 0.05; diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties index 137499571c..ff1c96ea96 100644 --- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties +++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties @@ -25,12 +25,6 @@ chemistry off; initialChemicalTimeStep 1e-07; -sequentialCoeffs -{ - cTauChem 0.001; - equilibriumRateLimiter off; -} - EulerImplicitCoeffs { cTauChem 0.05; diff --git a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties index 137499571c..ff1c96ea96 100644 --- a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties +++ b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties @@ -25,12 +25,6 @@ chemistry off; initialChemicalTimeStep 1e-07; -sequentialCoeffs -{ - cTauChem 0.001; - equilibriumRateLimiter off; -} - EulerImplicitCoeffs { cTauChem 0.05; diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties index 137499571c..ff1c96ea96 100644 --- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties +++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties @@ -25,12 +25,6 @@ chemistry off; initialChemicalTimeStep 1e-07; -sequentialCoeffs -{ - cTauChem 0.001; - equilibriumRateLimiter off; -} - EulerImplicitCoeffs { cTauChem 0.05;