chemistryModel: Added new option to specify the initial ODE integration time-step

In constant/chemistryProperties in addition to the specification of the initial
ODE integration time-step used at the start of the run:

    initialChemicalTimeStep 1e-12;

this time step may now also be specified for every chemistry integration by
setting the optional entry maxChemicalTimeStep, e.g.

    maxChemicalTimeStep 1e-12;

etc/controlDict

@@ -788,6 +788,7 @@ DebugSwitches
     solidBodyMotionFunction 0;
     solidBodyMotionFvMesh 0;
     solution            0;
+    specie              0;
     spectEddyVisc       0;
     sphereToCell        0;
     spherical           0;

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -191,7 +191,8 @@ void Foam::ODESolver::solve
     FatalErrorInFunction
         << "Integration steps greater than maximum " << maxSteps_ << nl
         << "    xStart = " << xStart << ", xEnd = " << xEnd
-        << ", x = " << x << ", dxDid = " << step.dxDid
+        << ", x = " << x << ", dxDid = " << step.dxDid << nl
+        << "    y = " << y
         << exit(FatalError);
 }
 

src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C

@@ -196,17 +196,17 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[lRef])
         {
             const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(0.0, c[lRef]), exp);
+            pf *= pow(max(c[lRef], 0.0), exp);
             lRef = si;
             slRef = s;
         }
         else
         {
             const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(0.0, c[si]), exp);
+            pf *= pow(max(c[si], 0.0), exp);
         }
     }
-    cf = max(0.0, c[lRef]);
+    cf = max(c[lRef], 0.0);
 
     {
         const scalar exp = R.lhs()[slRef].exponent;
@@ -238,23 +238,23 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[rRef])
         {
             const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(0.0, c[rRef]), exp);
+            pr *= pow(max(c[rRef], 0.0), exp);
             rRef = si;
             srRef = s;
         }
         else
         {
             const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(0.0, c[si]), exp);
+            pr *= pow(max(c[si], 0.0), exp);
         }
     }
-    cr = max(0.0, c[rRef]);
+    cr = max(c[rRef], 0.0);
 
     {
         const scalar exp = R.rhs()[srRef].exponent;
         if (exp < 1.0)
         {
-            if (cr>small)
+            if (cr > small)
             {
                 pr *= pow(cr, exp - 1.0);
             }
@@ -284,9 +284,9 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
     const scalar T = c[nSpecie_];
     const scalar p = c[nSpecie_ + 1];
 
-    for (label i = 0; i < nSpecie_; i++)
+    forAll(c_, i)
     {
-        c_[i] = max(0.0, c[i]);
+        c_[i] = max(c[i], 0.0);
     }
 
     omega(c_, T, p, dcdt);
@@ -366,7 +366,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
                     {
                         if (c_[si] > small)
                         {
-                            kf *= el*pow(c_[si] + vSmall, el - 1.0);
+                            kf *= el*pow(c_[si], el - 1.0);
                         }
                         else
                         {
@@ -412,7 +412,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
                     {
                         if (c_[si] > small)
                         {
-                            kr *= er*pow(c_[si] + vSmall, er - 1.0);
+                            kr *= er*pow(c_[si], er - 1.0);
                         }
                         else
                         {
@@ -741,6 +741,9 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::solve
 
             deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
 
+            this->deltaTChem_[celli] =
+                min(this->deltaTChem_[celli], this->deltaTChemMax_);
+
             for (label i=0; i<nSpecie_; i++)
             {
                 RR_[i][celli] =

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C

@@ -232,17 +232,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[lRef])
         {
             const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(0.0, c[lRef]), exp);
+            pf *= pow(max(c[lRef], 0.0), exp);
             lRef = si;
             slRef = s;
         }
         else
         {
             const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(0.0, c[si]), exp);
+            pf *= pow(max(c[si], 0.0), exp);
         }
     }
-    cf = max(0.0, c[lRef]);
+    cf = max(c[lRef], 0.0);
 
     {
         const scalar exp = R.lhs()[slRef].exponent;
@@ -274,23 +274,23 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[rRef])
         {
             const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(0.0, c[rRef]), exp);
+            pr *= pow(max(c[rRef], 0.0), exp);
             rRef = si;
             srRef = s;
         }
         else
         {
             const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(0.0, c[si]), exp);
+            pr *= pow(max(c[si], 0.0), exp);
         }
     }
-    cr = max(0.0, c[rRef]);
+    cr = max(c[rRef], 0.0);
 
     {
         const scalar exp = R.rhs()[srRef].exponent;
         if (exp < 1)
         {
-            if (cr>small)
+            if (cr > small)
             {
                 pr *= pow(cr, exp - 1);
             }
@@ -334,14 +334,14 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::derivatives
         // efficiencies
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
         }
     }
     else
     {
         for (label i=0; i<this->nSpecie(); i++)
         {
-            this->c_[i] = max(0.0, c[i]);
+            this->c_[i] = max(c[i], 0.0);
         }
     }
 
@@ -416,7 +416,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
         this->c_ = completeC_;
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
         }
     }
     else
@@ -456,7 +456,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
                         {
                             if (this->c_[si] > small)
                             {
-                                kf *= el*pow(this->c_[si] + vSmall, el - 1);
+                                kf *= el*pow(this->c_[si], el - 1);
                             }
                             else
                             {
@@ -514,7 +514,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
                         {
                             if (this->c_[si] > small)
                             {
-                                kr *= er*pow(this->c_[si] + vSmall, er - 1);
+                                kr *= er*pow(this->c_[si], er - 1);
                             }
                             else
                             {
@@ -797,6 +797,9 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
                 this->nSpecie_ = mechRed_->nSpecie();
             }
             deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
+
+            this->deltaTChem_[celli] =
+                min(this->deltaTChem_[celli], this->deltaTChemMax_);
         }
 
         // Set the RR vector (used in the solver)

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,6 +58,7 @@ Foam::basicChemistryModel::basicChemistryModel(basicThermo& thermo)
     mesh_(thermo.p().mesh()),
     chemistry_(lookup("chemistry")),
     deltaTChemIni_(readScalar(lookup("initialChemicalTimeStep"))),
+    deltaTChemMax_(lookupOrDefault("maxChemicalTimeStep", great)),
     deltaTChem_
     (
         IOobject

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,9 @@ protected:
         //- Initial chemical time step
         const scalar deltaTChemIni_;
 
+        //- Maximum chemical time step
+        const scalar deltaTChemMax_;
+
         //- Latest estimation of integration step
         volScalarField::Internal deltaTChem_;
 

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C

@@ -387,9 +387,9 @@ jacobian
                 {
                     if (exp < 1.0)
                     {
-                        if (c2[si]>small)
+                        if (c2[si] > small)
                         {
-                            kf *= exp*pow(c2[si] + vSmall, exp - 1.0);
+                            kf *= exp*pow(c2[si], exp - 1.0);
                         }
                         else
                         {

src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C

@@ -125,16 +125,7 @@ Foam::scalar Foam::ReversibleReaction
     const scalarField& c
 ) const
 {
-    scalar Kc = this->Kc(p, T);
+    return kfwd/max(this->Kc(p, T), 1e-6);
-
-    if (mag(Kc) > vSmall)
-    {
-        return kfwd/Kc;
-    }
-    else
-    {
-        return 0;
-    }
 }