From a92691b0dd419f47e70318ebb4788ebf8de1ad4f Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 3 Jul 2015 15:46:52 +0100
Subject: [PATCH] reactingTwoPhaseEulerFoam: Added mass-exchange contribution
 to the cell-based momentum equations

---
 .../multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H      | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H
index 424182f859..49fa1bd63c 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/UEqns.H
@@ -11,11 +11,15 @@ volScalarField Kd(fluid.Kd());
 {
     volScalarField Vm(fluid.Vm());
 
+    const volScalarField dmdt12(fluid.dmdt()*pos(fluid.dmdt()));
+    const volScalarField dmdt21(fluid.dmdt()*neg(fluid.dmdt()));
+
     {
         U1Eqn =
         (
             fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
-          - fvm::Sp(phase1.continuityError(), U1)
+          - fvm::Sp(phase1.continuityError() - dmdt21, U1)
+          - dmdt21*U2
           + MRF.DDt(alpha1*rho1, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
@@ -40,7 +44,8 @@ volScalarField Kd(fluid.Kd());
         U2Eqn =
         (
             fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
-          - fvm::Sp(phase2.continuityError(), U2)
+          - fvm::Sp(phase2.continuityError() + dmdt12, U2)
+          + dmdt12*U1
           + MRF.DDt(alpha2*rho2, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==