Merge pull request #304 from pkuLmq/master

add calculation of net production rates in dfChemistryModel

Author: OpenFOAMFans

src/dfChemistryModel/dfChemistryModel.C

@@ -65,6 +65,7 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
     hrtTemp_(mixture_.nSpecies()),
     cTemp_(mixture_.nSpecies()),
     RR_(mixture_.nSpecies()),
+    wrate_(mixture_.nSpecies()),
     alpha_(const_cast<volScalarField&>(thermo.alpha())),
     T_(thermo.T()),
     p_(thermo.p()),
@@ -231,7 +232,26 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
             )
         );
     }
-
+    forAll(wrate_, fieldi)
+    {
+        wrate_.set
+        (
+            fieldi,
+            new volScalarField::Internal
+            (
+                IOobject
+                (
+                    "wrate." + Y_[fieldi].name(),
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar(dimMass/dimVolume/dimTime, 0)
+            )
+        );
+    }
     forAll(rhoD_, i)
     {
         rhoD_.set
@@ -727,6 +747,34 @@ Foam::dfChemistryModel<ThermoType>::updateReactionRates
     return deltaTMin;
 }
 
+template<class ThermoType>
+void Foam::dfChemistryModel<ThermoType>::calculateW()
+{
+    doublereal Yi[mixture_.nSpecies()];
+    doublereal twrate[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++)
+        {
+            Yi[i] = Y_[i][celli];
+        }
+
+        CanteraGas_->setState_TPY(Ti, pi, Yi);
+
+        CanteraKinetics_->getNetProductionRates(twrate);
+
+        for (label i=0; i<mixture_.nSpecies(); i++)
+        {
+            wrate_[i][celli] = twrate[i]*CanteraGas_->molecularWeight(i);
+        }
+    }
+}
+
 template<class ThermoType>
 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
 Foam::dfChemistryModel<ThermoType>::calculateRR

src/dfChemistryModel/dfChemistryModel.H

@@ -119,6 +119,8 @@ public IOdictionary
         mutable scalarList cTemp_;
         // mass change rate, [kg/m^3/s]
         PtrList<volScalarField::Internal> RR_;
+        // net production rates, [kg/m^3/s]
+        PtrList<volScalarField::Internal> wrate_;
         hashedWordList species_;
         volScalarField& alpha_;
         volScalarField& T_;
@@ -285,6 +287,10 @@ public:
         //- Return access to chemical source terms [kg/m^3/s]
         volScalarField::Internal& RR(const label i) {return RR_[i];}
 
+        void calculateW();
+
+        volScalarField::Internal& wrate(const label i) {return wrate_[i];}
+
         tmp<volScalarField::Internal> calculateRR
         (
             const label reactionI,