chemistryModel: Mass conservative specie sources

Chemistry integration is assumed to occur for a fixed mass of fluid at a
constant pressure. Therefore, both the density and volume change during
integration.

When computing specie source terms, it is important that these relate
exactly to the mass of fluid considered during integration in order for
mass to be conserved. That means the change in mass fractions that
results from integration needs to be multiplied by the old-time density,
as that relates to the mass in the cell's volume.

The new time (end of integration) density could be used instead to
calculate the species source terms, but then a new volume would also
have to be calculated and accounted for such that the correct total mass
for the cell is maintained. It is equivalent just to use the old-time
density, for which no additional calculation is needed.

This change additionally means that chemistry integration over a
timestep has no dependence on new-time properties. So, there is never
any reason to integrate chemistry multiple times in a time-step, as the
result will always be the same. As such, the laminar combustion model
now no longer provides the option to outer correct when integrating the
chemistry. The option is still available when not integrating (i.e.,
when chemical source terms are calculated instantaneously).

src/combustionModels/laminar/laminar.C

@@ -57,10 +57,6 @@ Foam::combustionModels::laminar::laminar
         turb,
         combustionProperties
     ),
-    outerCorrect_
-    (
-        this->coeffs().lookupOrDefault("outerCorrect", false)
-    ),
     integrateReactionRate_
     (
         this->coeffs().lookupOrDefault("integrateReactionRate", true)
@@ -69,6 +65,10 @@ Foam::combustionModels::laminar::laminar
     (
         this->coeffs().lookupOrDefault("maxIntegrationTime", vGreat)
     ),
+    outerCorrect_
+    (
+        this->coeffs().lookupOrDefault("outerCorrect", false)
+    ),
     timeIndex_(-1),
     chemistryPtr_(basicChemistryModel::New(thermo))
 {
@@ -93,7 +93,11 @@ Foam::combustionModels::laminar::~laminar()
 
 void Foam::combustionModels::laminar::correct()
 {
-    if (!outerCorrect_ && timeIndex_ == this->mesh().time().timeIndex())
+    if
+    (
+        (integrateReactionRate_ || !outerCorrect_)
+     && timeIndex_ == this->mesh().time().timeIndex()
+    )
     {
         return;
     }
@@ -146,12 +150,12 @@ bool Foam::combustionModels::laminar::read()
 {
     if (combustionModel::read())
     {
-        outerCorrect_ =
-            this->coeffs().lookupOrDefault("outerCorrect", false);
         integrateReactionRate_ =
             this->coeffs().lookupOrDefault("integrateReactionRate", true);
         maxIntegrationTime_ =
             this->coeffs().lookupOrDefault("maxIntegrationTime", vGreat);
+        outerCorrect_ =
+            this->coeffs().lookupOrDefault("outerCorrect", false);
         return true;
     }
     else

src/combustionModels/laminar/laminar.H

@@ -56,16 +56,17 @@ class laminar
 {
     // Private Data
 
-        //- Run chemistry correction on every outer iteration. Default false.
-        bool outerCorrect_;
-
         //- Integrate reaction rate over the time-step
         //  using the selected ODE solver. Default true.
         bool integrateReactionRate_;
 
         //- Maximum integration time permitted. Default vGreat.
         scalar maxIntegrationTime_;
 
+        //- Run chemistry rate calculation on every outer iteration. Default
+        //  false. Has no effect if integrating.
+        bool outerCorrect_;
+
         //- The time index of the last correction
         label timeIndex_;
 

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C

@@ -719,9 +719,6 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
         return deltaTMin;
     }
 
-    tmp<volScalarField> trhovf(this->thermo().rho());
-    const volScalarField& rhovf = trhovf();
-
     const volScalarField& rho0vf =
         this->mesh().template lookupObject<volScalarField>
         (
@@ -743,7 +740,6 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
 
     forAll(rho0vf, celli)
     {
-        const scalar rho = rhovf[celli];
         const scalar rho0 = rho0vf[celli];
 
         scalar p = p0vf[celli];
@@ -883,7 +879,7 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
         // Set the RR vector (used in the solver)
         for (label i=0; i<nSpecie_; i++)
         {
-            RR_[i][celli] = (Y_[i]*rho - Y0[i]*rho0)/deltaT[celli];
+            RR_[i][celli] = rho0*(Y_[i] - Y0[i])/deltaT[celli];
         }
 
         if (loadBalancing_)