reactingEulerFoam: Removed continuity error from species equations

This fixes a consistency issue in the interface-composition method, and
also seems to improve stability/convergence of the pimple iteration in
the presence of significant mass transfer.

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -167,6 +167,10 @@ addMassTransferMomentumTransfer(phaseSystem::momentumTransferTable& eqns) const
             continue;
         }
 
+        // Note that the phase UEqn contains a continuity error term, which
+        // implicitly adds a mass transfer term of fvm::Sp(dmdt, U). These
+        // additions do not include this term.
+
         const volScalarField dmdt(this->dmdt(pair));
 
         if (!pair.phase1().stationary())

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C

@@ -129,6 +129,10 @@ heatTransfer() const
         const volScalarField K1(phase1.K());
         const volScalarField K2(phase2.K());
 
+        // Note that the phase heEqn contains a continuity error term, which
+        // implicitly adds a mass transfer term of fvm::Sp(dmdt, he). These
+        // additions do not include this term.
+
         const volScalarField dmdt(this->dmdt(pair));
         const volScalarField dmdt21(posPart(dmdt));
         const volScalarField dmdt12(negPart(dmdt));

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C

@@ -161,6 +161,9 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::massTransfer() const
         const phaseModel& phase = pair.phase1();
         const phaseModel& otherPhase = pair.phase2();
 
+        // Note that the phase YiEqn does not contain a continuity error term,
+        // so these additions represent the entire mass transfer
+
         const volScalarField dmdt(this->rDmdt(pair));
         const volScalarField dmdt12(negPart(dmdt));
         const volScalarField dmdt21(posPart(dmdt));
@@ -181,11 +184,11 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::massTransfer() const
 
             *eqns[name] +=
                 dmdt21*eqns[otherName]->psi()
-              - fvm::Sp(dmdt21, eqns[name]->psi());
+              + fvm::Sp(dmdt12, eqns[name]->psi());
 
             *eqns[otherName] -=
                 dmdt12*eqns[name]->psi()
-              - fvm::Sp(dmdt12, eqns[otherName]->psi());
+              + fvm::Sp(dmdt21, eqns[otherName]->psi());
         }
 
     }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C

@@ -192,6 +192,9 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::massTransfer() const
         const phaseModel& phase = pair.phase1();
         const phaseModel& otherPhase = pair.phase2();
 
+        // Note that the phase YiEqn does not contain a continuity error term,
+        // so these additions represent the entire mass transfer
+
         const volScalarField dmdt(this->pDmdt(pair));
         const volScalarField dmdt12(negPart(dmdt));
         const volScalarField dmdt21(posPart(dmdt));
@@ -212,11 +215,11 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::massTransfer() const
 
             *eqns[name] +=
                 dmdt21*eqns[otherName]->psi()
-              - fvm::Sp(dmdt21, eqns[name]->psi());
+              + fvm::Sp(dmdt12, eqns[name]->psi());
 
             *eqns[otherName] -=
                 dmdt12*eqns[name]->psi()
-              - fvm::Sp(dmdt12, eqns[otherName]->psi());
+              + fvm::Sp(dmdt21, eqns[otherName]->psi());
         }
     }
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -283,6 +283,11 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
 
         forAll(Yi, i)
         {
+            if (Yi[i].member() != volatile_)
+            {
+                continue;
+            }
+
             const word name
             (
                 IOobject::groupName(volatile_, phase.name())
@@ -293,18 +298,15 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
                 IOobject::groupName(volatile_, otherPhase.name())
             );
 
+            // Note that the phase YiEqn does not contain a continuity error
+            // term, so these additions represent the entire mass transfer
+
             const volScalarField dmdt(this->iDmdt(pair) + this->wDmdt(pair));
 
-            *eqns[name] -= fvm::Sp(dmdt, eqns[name]->psi());
-            *eqns[otherName] += fvm::Sp(dmdt, eqns[otherName]->psi());
-
-            if (Yi[i].member() == volatile_)
-            {
             *eqns[name] += dmdt;
             *eqns[otherName] -= dmdt;
         }
     }
-    }
 
     return eqnsPtr;
 }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C

@@ -155,7 +155,6 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn(volScalarField& Yi)
     (
         fvm::ddt(alpha, rho, Yi)
       + fvm::div(alphaRhoPhi, Yi, "div(" + alphaRhoPhi.name() + ",Yi)")
-      - fvm::Sp(this->continuityError(), Yi)
 
       - fvm::laplacian
         (