diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
index 546c369e5a..8aa0b8af23 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/UEqns.H
@@ -29,6 +29,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
           + MRF.DDt(alpha1*rho1, U1)
           + phase1.turbulence().divDevRhoReff(U1)
           + Vm*(UgradU1 - (UgradU2 & U2))
+          - fvOptions(alpha1, rho1, U1)
         );
         U1Eqn.relax();
         fvOptions.constrain(U1Eqn);
@@ -44,6 +45,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
           + MRF.DDt(alpha2*rho2, U2)
           + phase2.turbulence().divDevRhoReff(U2)
           + Vm*(UgradU2 - (UgradU1 & U1))
+          - fvOptions(alpha2, rho2, U2)
         );
         U2Eqn.relax();
         fvOptions.constrain(U2Eqn);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H
index 257008f7ad..0f46ae5b85 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H
@@ -3,6 +3,7 @@ Info<< "Constructing face momentum equations" << endl;
 MRF.correctBoundaryVelocity(U1);
 MRF.correctBoundaryVelocity(U2);
 MRF.correctBoundaryVelocity(U);
+Info<< "Constructing face momentum equations" << endl;
 
 fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
@@ -29,6 +30,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
           + MRF.DDt(alpha1*rho1, U1)
           + phase1.turbulence().divDevRhoReff(U1)
           + Vm*(UgradU1 - (UgradU2 & U2))
+          - fvOptions(alpha1, rho1, U1)
         );
         U1Eqn.relax();
         fvOptions.constrain(U1Eqn);
@@ -43,6 +45,7 @@ fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
           + MRF.DDt(alpha2*rho2, U2)
           + phase2.turbulence().divDevRhoReff(U2)
           + Vm*(UgradU2 - (UgradU1 & U1))
+          - fvOptions(alpha2, rho2, U2)
         );
         U2Eqn.relax();
         fvOptions.constrain(U2Eqn);