MichaelisMentenReactionRate: Michaelis-Menten reaction rate for enzymatic reactions

Updated chemFoam to support liquid reactions

applications/solvers/combustion/chemFoam/chemFoam.C

@@ -41,6 +41,98 @@ Description
 #include "basicSpecieMixture.H"
 #include "cellModeller.H"
 
+
+template<class ThermoType>
+scalarList W(const rhoReactionThermo& thermo)
+{
+    const PtrList<ThermoType>& specieData =
+        dynamicCast<const reactingMixture<ThermoType>>(thermo)
+       .speciesData();
+
+    scalarList W(specieData.size());
+
+    forAll(specieData, i)
+    {
+        W[i] = specieData[i].W();
+    }
+
+    return W;
+}
+
+
+template<class ThermoType>
+scalar h0
+(
+    const rhoReactionThermo& thermo,
+    const scalarList& Y,
+    const scalar p,
+    const scalar T
+)
+{
+    const PtrList<ThermoType>& specieData =
+        dynamic_cast<const reactingMixture<ThermoType>&>(thermo)
+       .speciesData();
+
+    scalar h0 = 0;
+    forAll(Y, i)
+    {
+        h0 += Y[i]*specieData[i].Hs(p, T);
+    }
+
+    return h0;
+}
+
+
+scalarList W(const rhoReactionThermo& thermo)
+{
+    if (isA<reactingMixture<gasHThermoPhysics>>(thermo))
+    {
+        return W<gasHThermoPhysics>(thermo);
+    }
+    else if (isA<reactingMixture<constFluidHThermoPhysics>>(thermo))
+    {
+        return W<constFluidHThermoPhysics>(thermo);
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Thermodynamics type " << thermo.type()
+            << " not supported by chemFoam"
+            << exit(FatalError);
+
+        return scalarList::null();
+    }
+}
+
+
+scalar h0
+(
+    const rhoReactionThermo& thermo,
+    const scalarList& Y,
+    const scalar p,
+    const scalar T
+)
+{
+    if (isA<reactingMixture<gasHThermoPhysics>>(thermo))
+    {
+        return h0<gasHThermoPhysics>(thermo, Y, p, T);
+    }
+    else if (isA<reactingMixture<constFluidHThermoPhysics>>(thermo))
+    {
+        return h0<constFluidHThermoPhysics>(thermo, Y, p, T);
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Thermodynamics type " << thermo.type()
+            << " not supported by chemFoam"
+            << exit(FatalError);
+
+        return 0;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])

applications/solvers/combustion/chemFoam/pEqn.H

@@ -5,7 +5,7 @@
         scalar invW = 0.0;
         forAll(Y, i)
         {
-            invW += Y[i][0]/specieData[i].W();
+            invW += Y[i][0]/W[i];
         }
 
         Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW;

applications/solvers/combustion/chemFoam/readInitialConditions.H

@@ -1,4 +1,4 @@
-    word constProp(initialConditions.lookup("constantProperty"));
+    const word constProp(initialConditions.lookup("constantProperty"));
     if ((constProp != "pressure") && (constProp != "volume"))
     {
         FatalError << "in initialConditions, unknown constantProperty type "
@@ -6,7 +6,7 @@
             << abort(FatalError);
     }
 
-    word fractionBasis(initialConditions.lookup("fractionBasis"));
+    const word fractionBasis(initialConditions.lookup("fractionBasis"));
     if ((fractionBasis != "mass") && (fractionBasis != "mole"))
     {
         FatalError << "in initialConditions, unknown fractionBasis type " << nl
@@ -14,20 +14,8 @@
             << fractionBasis << abort(FatalError);
     }
 
-    label nSpecie = Y.size();
-    PtrList<gasHThermoPhysics> specieData(Y.size());
-    forAll(specieData, i)
-    {
-        specieData.set
-        (
-            i,
-            new gasHThermoPhysics
-            (
-                dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
-                    (thermo).speciesData()[i]
-            )
-        );
-    }
+    const label nSpecie = Y.size();
+    const scalarList W(::W(thermo));
 
     scalarList Y0(nSpecie, 0.0);
     scalarList X0(nSpecie, 0.0);
@@ -49,12 +37,12 @@
         forAll(Y, i)
         {
             X0[i] /= mTot;
-            mw += specieData[i].W()*X0[i];
+            mw += W[i]*X0[i];
         }
 
         forAll(Y, i)
         {
-            Y0[i] = X0[i]*specieData[i].W()/mw;
+            Y0[i] = X0[i]*W[i]/mw;
         }
     }
     else  // mass fraction
@@ -73,21 +61,21 @@
         forAll(Y, i)
         {
             Y0[i] /= mTot;
-            invW += Y0[i]/specieData[i].W();
+            invW += Y0[i]/W[i];
         }
         const scalar mw = 1.0/invW;
 
         forAll(Y, i)
         {
-            X0[i] = Y0[i]*mw/specieData[i].W();
+            X0[i] = Y0[i]*mw/W[i];
         }
     }
 
-    scalar h0 = 0.0;
+    const scalar h0 = ::h0(thermo, Y0, p[0], T0);
+
     forAll(Y, i)
     {
         Y[i] = Y0[i];
-        h0 += Y0[i]*specieData[i].Hs(p[0], T0);
     }
 
     thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);