Merge pull request #280 from minzhang0929/master

Create access to the reaction rate of each species involved in a given reaction
deepmodeling · Jun 1, 2023 · 2861dcc · 2861dcc
2 parents ee94d5b + aa0d895
commit 2861dcc
Show file tree

Hide file tree

Showing 4 changed files with 78 additions and 2 deletions.
diff --git a/src/dfCanteraMixture/CanteraMixture.C b/src/dfCanteraMixture/CanteraMixture.C
@@ -63,6 +63,7 @@ Foam::CanteraMixture::CanteraMixture
 
     CanteraSolution_=Cantera::newSolution(CanteraMechanismFile_, "");
     CanteraGas_=CanteraSolution_->thermo();
+    CanteraKinetics_=CanteraSolution_->kinetics();
     CanteraTransport_=newTransportMgr(transportModelName_, CanteraGas_.get());
 
     Y_.resize(nSpecies());

diff --git a/src/dfCanteraMixture/CanteraMixture.H b/src/dfCanteraMixture/CanteraMixture.H
@@ -39,6 +39,7 @@ SourceFiles
 #include "cantera/zerodim.h"
 #include "cantera/transport.h"
 #include "cantera/thermo/IdealGasPhase.h"
+#include "cantera/kinetics.h"
 #include "IOdictionary.H"
 #include "fvMesh.H"
 #include "word.H"
@@ -66,6 +67,7 @@ class CanteraMixture
         const string CanteraMechanismFile_;
         std::shared_ptr<Cantera::Solution> CanteraSolution_;
         std::shared_ptr<Cantera::ThermoPhase> CanteraGas_;
+        std::shared_ptr<Cantera::Kinetics> CanteraKinetics_;
         word transportModelName_;
         Cantera::Transport* CanteraTransport_;
         hashedWordList species_;
@@ -223,10 +225,13 @@ public:
 
         const hashedWordList& species() const {return species_;}
         int nSpecies() {return int(CanteraGas_->nSpecies());}
+        int nReactions() {return int(CanteraKinetics_->nReactions());} 
 
         std::shared_ptr<Cantera::ThermoPhase> CanteraGas() {return CanteraGas_;}
 
         std::shared_ptr<Cantera::Solution> CanteraSolution() {return CanteraSolution_;}
+
+        std::shared_ptr<Cantera::Kinetics> CanteraKinetics() {return CanteraKinetics_;}
 
         Cantera::Transport* CanteraTransport() {return CanteraTransport_;}
 

diff --git a/src/dfChemistryModel/dfChemistryModel.C b/src/dfChemistryModel/dfChemistryModel.C
@@ -51,14 +51,15 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
     thermo_(thermo),
     mixture_(dynamic_cast<CanteraMixture&>(thermo)),
     CanteraGas_(mixture_.CanteraGas()),
+    CanteraKinetics_(mixture_.CanteraKinetics()),
     mesh_(thermo.p().mesh()),
     chemistry_(lookup("chemistry")),
     relTol_(this->subDict("odeCoeffs").lookupOrDefault("relTol",1e-9)),
     absTol_(this->subDict("odeCoeffs").lookupOrDefault("absTol",1e-15)),
     Y_(mixture_.Y()),
     rhoD_(mixture_.nSpecies()),
     hai_(mixture_.nSpecies()),
-    hc_(mixture_.nSpecies()),
+    hc_(mixture_.nSpecies()), 
     yTemp_(mixture_.nSpecies()),
     dTemp_(mixture_.nSpecies()),
     hrtTemp_(mixture_.nSpecies()),
@@ -719,7 +720,69 @@ Foam::dfChemistryModel<ThermoType>::updateReactionRates
     return deltaTMin;
 }
 
+template<class ThermoType>
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
+Foam::dfChemistryModel<ThermoType>::calculateRR
+(
+    const label reactionI,
+    const label speciei
+) const
+{
+    tmp<volScalarField::Internal> tRR
+    (
+        volScalarField::Internal::New
+        (
+           "RR",
+           mesh_,
+           dimensionedScalar(dimMass/dimVolume/dimTime, 0)
+        )
+    );
+
+    volScalarField::Internal& RR = tRR.ref();
+
+    doublereal netRate[mixture_.nReactions()];
+    doublereal X[mixture_.nSpecies()];
+
+    forAll(rho_, celli)
+    {
+        const scalar rhoi = rho_[celli];
+        const scalar Ti = T_[celli];
+        const scalar pi = p_[celli];
+
+        for (label i=0; i<mixture_.nSpecies(); i++)
+        {
+            const scalar Yi = Y_[i][celli];
+
+            X[i] = rhoi*Yi/CanteraGas_->molecularWeight(i);
+        }
+
+	CanteraGas_->setState_TPX(Ti, pi, X);
+
+	CanteraKinetics_->getNetRatesOfProgress(netRate);
 
+	auto R = CanteraKinetics_->reaction(reactionI);
+
+	for (const auto& sp : R->reactants)
+        {
+		if (speciei == static_cast<int>(CanteraGas_->speciesIndex(sp.first)))
+		{
+			RR[celli] -= sp.second*netRate[reactionI];
+		}
+
+	}			
+	for (const auto& sp : R->products)
+	{
+		if (speciei == static_cast<int>(CanteraGas_->speciesIndex(sp.first)))
+		{
+			RR[celli] += sp.second*netRate[reactionI];
+		}
+	}					
+
+        RR[celli] *= CanteraGas_->molecularWeight(speciei);
+    }
+
+    return tRR;
+}
 
 template <class ThermoType>
 Foam::LoadBalancer

diff --git a/src/dfChemistryModel/dfChemistryModel.H b/src/dfChemistryModel/dfChemistryModel.H
@@ -93,6 +93,7 @@ public IOdictionary
         ThermoType& thermo_;
         CanteraMixture& mixture_;
         std::shared_ptr<Cantera::ThermoPhase> CanteraGas_;
+        std::shared_ptr<Cantera::Kinetics> CanteraKinetics_;
         const fvMesh& mesh_;
         Switch chemistry_;
 
@@ -107,7 +108,7 @@ public IOdictionary
         // species absolute enthalpy, [J/kg]
         PtrList<volScalarField> hai_;
         // species chemistry enthalpy, [J/kg]
-        scalarList hc_;
+        scalarList hc_;         
         // temp mass fraction
         mutable scalarList yTemp_;
         // temp mass diffusion coefficients
@@ -284,6 +285,12 @@ public:
         //- Return access to chemical source terms [kg/m^3/s]
         volScalarField::Internal& RR(const label i) {return RR_[i];}
 
+        tmp<volScalarField::Internal> calculateRR
+        (
+            const label reactionI,
+            const label speciei
+        ) const;
+
         //- Return the heat release rate [J/m/s^3]
         const volScalarField& Qdot() const
         {