reactingTwoPhaseEulerFoam: Corrected handling of heat-transfer caused by mass-transfer

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C

@@ -288,6 +288,38 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
         }
     }
 
+    // Source term due to mass trasfer
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const volScalarField& he1(pair.phase1().thermo().he());
+        const volScalarField& he2(pair.phase2().thermo().he());
+
+        const volScalarField& K1(pair.phase1().K());
+        const volScalarField& K2(pair.phase2().K());
+
+        const volScalarField dmdt(this->dmdt(pair));
+        const volScalarField dmdt12(dmdt*pos(dmdt));
+        const volScalarField dmdt21(dmdt*neg(dmdt));
+
+        *eqns[pair.phase1().name()] +=
+            fvm::Sp(dmdt21, he1) + dmdt21*K1 - dmdt21*(he2 + K2);
+
+        *eqns[pair.phase2().name()] +=
+            dmdt12*(he1 + K1) - fvm::Sp(dmdt12, he2) - dmdt12*K2;
+    }
+
     return eqnsPtr;
 }
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -154,4 +154,12 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::divU(const volScalarField& divU)
 }
 
 
+template<class BasePhaseModel>
+const Foam::volScalarField&
+Foam::AnisothermalPhaseModel<BasePhaseModel>::K() const
+{
+    return K_;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H

@@ -89,11 +89,14 @@ public:
             //- Return true if the phase is compressible otherwise false
             virtual bool compressible() const;
 
-            //- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual const volScalarField& divU() const;
 
-            //- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
+            //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(const volScalarField& divU);
+
+            //- Return the phase kinetic energy
+            virtual const volScalarField& K() const;
 };
 
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C

@@ -314,14 +314,6 @@ Foam::MovingPhaseModel<BasePhaseModel>::continuityError() const
 }
 
 
-template<class BasePhaseModel>
-Foam::volScalarField&
-Foam::MovingPhaseModel<BasePhaseModel>::continuityError()
-{
-    return continuityError_;
-}
-
-
 template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::MovingPhaseModel<BasePhaseModel>::phi() const

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H

@@ -152,9 +152,6 @@ public:
             //- Constant access the continuity error
             virtual tmp<volScalarField> continuityError() const;
 
-            //- Access the continuity error
-            virtual volScalarField& continuityError();
-
             //- Constant access the volumetric flux
             virtual tmp<surfaceScalarField> phi() const;
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C

@@ -155,6 +155,13 @@ void Foam::phaseModel::divU(const volScalarField& divU)
 }
 
 
+const Foam::volScalarField& Foam::phaseModel::K() const
+{
+    notImplemented("Foam::phaseModel::K()");
+    return volScalarField::null();
+}
+
+
 const Foam::surfaceScalarField& Foam::phaseModel::DbyA() const
 {
     return surfaceScalarField::null();

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H

@@ -175,6 +175,9 @@ public:
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(const volScalarField& divU);
 
+            //- Return the phase kinetic energy
+            virtual const volScalarField& K() const;
+
 
         // Implicit phase pressure and dispersion support
 
@@ -218,9 +221,6 @@ public:
             //- Constant access the continuity error
             virtual tmp<volScalarField> continuityError() const = 0;
 
-            //- Access the continuity error
-            virtual volScalarField& continuityError() = 0;
-
             //- Constant access the volumetric flux
             virtual tmp<surfaceScalarField> phi() const = 0;