From d425421c0ab3adc3d91f90e1b3243251be884200 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Fri, 15 Sep 2017 11:43:15 +0100
Subject: [PATCH] ENH: reactingTwoPhaseEulerFoam: Added continuity fvOption
 term to U equations

Patch contributed by Juho Peltola, VTT
Resolves bug-report https://bugs.openfoam.org/view.php?id=2699
---
 .../reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H     | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
index 8aa0b8af23..cfaaeef81b 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
@@ -30,6 +30,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
           + phase1.turbulence().divDevRhoReff(U1)
           + Vm*(UgradU1 - (UgradU2 & U2))
           - fvOptions(alpha1, rho1, U1)
+          + fvm::SuSp(fvOptions(alpha1, rho1)&rho1,U1)
         );
         U1Eqn.relax();
         fvOptions.constrain(U1Eqn);
@@ -46,6 +47,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
           + phase2.turbulence().divDevRhoReff(U2)
           + Vm*(UgradU2 - (UgradU1 & U1))
           - fvOptions(alpha2, rho2, U2)
+          + fvm::SuSp(fvOptions(alpha2, rho2)&rho2,U2)
         );
         U2Eqn.relax();
         fvOptions.constrain(U2Eqn);