diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C index b985c0061b..0e663c1bbc 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C @@ -514,7 +514,7 @@ void Foam::chemistryModel::jacobian } // calculate the dcdT elements numerically - const scalar delta = 1.0e-8; + const scalar delta = 1.0e-3; const scalarField dcdT0(omega(c2, T - delta, p)); const scalarField dcdT1(omega(c2, T + delta, p)); @@ -772,9 +772,6 @@ Foam::scalar Foam::chemistryModel::solve this->thermo().rho() ); - tmp thc = this->thermo().hc(); - const scalarField& hc = thc(); - const scalarField& he = this->thermo().he(); const scalarField& T = this->thermo().T(); const scalarField& p = this->thermo().p(); @@ -784,8 +781,7 @@ Foam::scalar Foam::chemistryModel::solve forAll(rho, celli) { const scalar rhoi = rho[celli]; - const scalar hi = he[celli] + hc[celli]; - const scalar pi = p[celli]; + scalar pi = p[celli]; scalar Ti = T[celli]; for (label i=0; i::solve c0[i] = c[i]; } - // initialise timing parameters - scalar t = 0; + // Initialise time progress scalar timeLeft = deltaT[celli]; - scalar tauC = this->deltaTChem_[celli]; - scalar dt = min(timeLeft, tauC); - // calculate the chemical source terms + // Calculate the chemical source terms while (timeLeft > SMALL) { - tauC = this->solve(c, Ti, pi, t, dt); - t += dt; - - // update the temperature - const scalar cTot = sum(c); - ThermoType mixture((c[0]/cTot)*specieThermo_[0]); - for (label i=1; isolve(c, Ti, pi, dt, this->deltaTChem_[celli]); timeLeft -= dt; - this->deltaTChem_[celli] = tauC; - dt = max(SMALL, min(timeLeft, tauC)); } - deltaTMin = min(tauC, deltaTMin); + + deltaTMin = min(this->deltaTChem_[celli], deltaTMin); for (label i=0; i::solve template -Foam::scalar Foam::chemistryModel::solve +void Foam::chemistryModel::solve ( scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const { notImplemented @@ -871,14 +853,12 @@ Foam::scalar Foam::chemistryModel::solve "chemistryModel::solve" "(" "scalarField&, " - "const scalar, " - "const scalar, " - "const scalar, " - "const scalar" + "scalar&, " + "scalar&, " + "scalar&, " + "scalar&" ") const" ); - - return (0); } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H index dcb092c335..9bec27b25b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H @@ -78,6 +78,8 @@ class chemistryModel protected: + typedef ThermoType thermoType; + // Private data //- Reference to the field of specie mass fractions @@ -260,13 +262,13 @@ public: scalarSquareMatrix& dfdc ) const; - virtual scalar solve + virtual void solve ( scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const; }; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C index 349ef77f8c..76c93cb0ac 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C @@ -38,7 +38,8 @@ Foam::EulerImplicit::EulerImplicit chemistrySolver(mesh), coeffsDict_(this->subDict("EulerImplicitCoeffs")), cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))), - eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")) + eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter")), + cTp_(this->nEqns()) {} @@ -52,33 +53,47 @@ Foam::EulerImplicit::~EulerImplicit() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::scalar Foam::EulerImplicit::solve +void Foam::EulerImplicit::solve ( - scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalarField& c, + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const { - scalar pf, cf, pr, cr; - label lRef, rRef; + deltaT = min(deltaT, subDeltaT); const label nSpecie = this->nSpecie(); simpleMatrix RR(nSpecie, 0, 0); - for (label i = 0; i < nSpecie; i++) + for (label i=0; ispecieThermo_[0] + ); + for (label i=1; ispecieThermo_[i]; + } + scalar ha = mixture.Ha(p, T); + + for (label i=0; ireactions(), i) { + scalar pf, cf, pr, cr; + label lRef, rRef; + scalar omegai = this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef); scalar corr = 1.0; @@ -86,11 +101,11 @@ Foam::scalar Foam::EulerImplicit::solve { if (omegai < 0.0) { - corr = 1.0/(1.0 + pr*dt); + corr = 1.0/(1.0 + pr*deltaT); } else { - corr = 1.0/(1.0 + pf*dt); + corr = 1.0/(1.0 + pf*deltaT); } } @@ -98,37 +113,46 @@ Foam::scalar Foam::EulerImplicit::solve } - for (label i = 0; i < nSpecie; i++) + for (label i=0; ispecieThermo_[0]; + for (label i=1; ispecieThermo_[i]; + } + T = mixture.THa(ha, p, T); + + // Estimate the next time step scalar tMin = GREAT; const label nEqns = this->nEqns(); - scalarField c1(nEqns, 0.0); - for (label i = 0; i < nSpecie; i++) + for (label i=0; iderivatives(0.0, c1, dcdt); + this->derivatives(0.0, cTp_, dcdt); const scalar sumC = sum(c); - for (label i = 0; i < nSpecie; i++) + for (label i=0; i::solve } } - return cTauChem_*tMin; + subDeltaT = cTauChem_*tMin; } diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H index 42f192b38f..ffc2949ffd 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H @@ -65,6 +65,9 @@ class EulerImplicit //- Equilibrium rate limiter flag (on/off) Switch eqRateLimiter_; + // Solver data + mutable scalarField cTp_; + public: @@ -85,13 +88,13 @@ public: // Member Functions //- Update the concentrations and return the chemical time - virtual scalar solve + virtual void solve ( - scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalarField& c, + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const; }; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H index 574a808383..4c0043d717 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H @@ -69,13 +69,13 @@ public: // Member Functions //- Update the concentrations and return the chemical time - virtual scalar solve + virtual void solve ( scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const = 0; }; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H index 4660dcd17a..72ab52b37f 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H @@ -31,9 +31,9 @@ License #include "chemistryModel.H" #include "noChemistrySolver.H" +#include "sequential.H" #include "EulerImplicit.H" #include "ode.H" -#include "sequential.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C index efc8c10735..be4d02093b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C @@ -48,17 +48,15 @@ Foam::noChemistrySolver::~noChemistrySolver() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::scalar Foam::noChemistrySolver::solve +void Foam::noChemistrySolver::solve ( scalarField&, - const scalar, - const scalar, - const scalar, - const scalar + scalar&, + scalar&, + scalar&, + scalar& ) const -{ - return GREAT; -} +{} // ************************************************************************* // diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H index 2a84bac414..7073df725f 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H @@ -72,13 +72,13 @@ public: // Member Functions //- Update the concentrations and return the chemical time - virtual scalar solve + virtual void solve ( - scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalarField& c, + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const; }; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C index aff5be67ca..e885653ccd 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C @@ -38,7 +38,8 @@ Foam::ode::ode coeffsDict_(this->subDict("odeCoeffs")), solverName_(coeffsDict_.lookup("solver")), odeSolver_(ODESolver::New(solverName_, *this)), - eps_(readScalar(coeffsDict_.lookup("eps"))) + eps_(readScalar(coeffsDict_.lookup("eps"))), + cTp_(this->nEqns()) {} @@ -52,44 +53,41 @@ Foam::ode::~ode() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::scalar Foam::ode::solve +void Foam::ode::solve ( scalarField& c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const { label nSpecie = this->nSpecie(); - scalarField c1(this->nEqns(), 0.0); - // copy the concentration, T and P to the total solve-vector - for (label i = 0; i < nSpecie; i++) + // Copy the concentration, T and P to the total solve-vector + for (register int i=0; isolve ( *this, - t0, - t0 + dt, - c1, + 0, + deltaT, + cTp_, eps_, - dtEst + subDeltaT ); - forAll(c, i) + for (register int i=0; i::~sequential() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -Foam::scalar Foam::sequential::solve +void Foam::sequential::solve ( - scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalarField& c, + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const { - scalar tChemInv = SMALL; + deltaT = min(deltaT, subDeltaT); - scalar pf, cf, pb, cb; - label lRef, rRef; + 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*dt; + omega /= 1.0 + pb*deltaT; } else { - omega /= 1.0 + pf*dt; + omega /= 1.0 + pf*deltaT; } } tChemInv = max(tChemInv, mag(omega)); - this->updateConcsInReactionI(i, dt, omega, p, T, c); + this->updateConcsInReactionI(i, deltaT, omega, p, T, c); } - return cTauChem_/tChemInv; + // 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 index b13273590e..c65311dc8c 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.H @@ -28,9 +28,7 @@ Description Foam::sequential SourceFiles - sequentialI.H sequential.C - sequentialIO.C \*---------------------------------------------------------------------------*/ @@ -86,13 +84,13 @@ public: // Member Functions //- Update the concentrations and return the chemical time - virtual scalar solve + virtual void solve ( - scalarField &c, - const scalar T, - const scalar p, - const scalar t0, - const scalar dt + scalarField& c, + scalar& T, + scalar& p, + scalar& deltaT, + scalar& subDeltaT ) const; }; diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H index 6eb8e9b81b..92b7b7dbff 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H @@ -34,8 +34,8 @@ Description #include "noChemistrySolver.H" #include "EulerImplicit.H" -#include "ode.H" #include "sequential.H" +#include "ode.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties index 9c00510123..2abcf4ae70 100644 --- a/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties +++ b/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties @@ -40,7 +40,6 @@ odeCoeffs { solver KRR4; eps 0.05; - scale 1; } // ************************************************************************* // diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties index 9c00510123..a5a53ae5b7 100644 --- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties +++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties @@ -17,7 +17,7 @@ FoamFile chemistryType { - chemistrySolver ode; + chemistrySolver EulerImplicit; chemistryThermo psi; } @@ -27,12 +27,12 @@ initialChemicalTimeStep 1e-07; sequentialCoeffs { - cTauChem 0.001; + cTauChem 0.05; } EulerImplicitCoeffs { - cTauChem 0.05; + cTauChem 1; equilibriumRateLimiter off; } @@ -40,7 +40,6 @@ odeCoeffs { solver KRR4; eps 0.05; - scale 1; } // ************************************************************************* // diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties index e0706fc235..c0b596bb36 100644 --- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties +++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/chemistryProperties @@ -41,7 +41,6 @@ odeCoeffs { solver SIBS; eps 0.05; - scale 1; } diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties index 5aa51f0290..137499571c 100644 --- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties +++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/chemistryProperties @@ -41,7 +41,6 @@ odeCoeffs { solver RK; eps 0.05; - scale 1; } diff --git a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties index 5aa51f0290..137499571c 100644 --- a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties +++ b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/chemistryProperties @@ -41,7 +41,6 @@ odeCoeffs { solver RK; eps 0.05; - scale 1; } diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties index 5aa51f0290..137499571c 100644 --- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties +++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/chemistryProperties @@ -41,7 +41,6 @@ odeCoeffs { solver RK; eps 0.05; - scale 1; }