diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
index 32cedb0fe0..cf70a99e18 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
@@ -31,5 +31,6 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
 
         UEqns[phasei].relax();
         fvOptions.constrain(UEqns[phasei]);
+        fvOptions.correct(U);
     }
 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H
index 9252630429..99bc95cfd6 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/UEqns.H
@@ -3,61 +3,36 @@ Info<< "Constructing momentum equations" << endl;
 fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
 
-volScalarField Kd(fluid.Kd());
-
 {
-    volScalarField Vm(fluid.Vm());
+    autoPtr<phaseSystem::momentumTransferTable>
+        momentumTransferPtr(fluid.momentumTransfer());
 
-    const volScalarField dmdt12(posPart(fluid.dmdt()));
-    const volScalarField dmdt21(negPart(fluid.dmdt()));
+    phaseSystem::momentumTransferTable&
+        momentumTransfer(momentumTransferPtr());
 
     {
         U1Eqn =
         (
-            fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
-          - fvm::Sp(phase1.continuityError() + dmdt21, U1) + dmdt21*U2
-          + MRF.DDt(alpha1*rho1, U1)
-          + phase1.turbulence().divDevRhoReff(U1)
+            phase1.UEqn()
          ==
-          - Vm
-           *(
-                fvm::ddt(U1)
-              + fvm::div(phi1, U1)
-              - fvm::Sp(fvc::div(phi1), U1)
-              - phase2.DUDt()
-            )
-          - MRF.DDt(Vm, U1 - U2)
+           *momentumTransfer[phase1.name()]
           + fvOptions(alpha1, rho1, U1)
         );
         U1Eqn.relax();
-        U1Eqn += fvm::Sp(Kd, U1);
         fvOptions.constrain(U1Eqn);
-        U1.correctBoundaryConditions();
         fvOptions.correct(U1);
     }
 
     {
         U2Eqn =
         (
-            fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
-          - fvm::Sp(phase2.continuityError() - dmdt12, U2) - dmdt12*U1
-          + MRF.DDt(alpha2*rho2, U2)
-          + phase2.turbulence().divDevRhoReff(U2)
+            phase2.UEqn()
          ==
-          - Vm
-           *(
-                fvm::ddt(U2)
-              + fvm::div(phi2, U2)
-              - fvm::Sp(fvc::div(phi2), U2)
-              - phase1.DUDt()
-            )
-          - MRF.DDt(Vm, U2 - U1)
+           *momentumTransfer[phase2.name()]
           + fvOptions(alpha2, rho2, U2)
         );
         U2Eqn.relax();
-        U2Eqn += fvm::Sp(Kd, U2);
         fvOptions.constrain(U2Eqn);
-        U2.correctBoundaryConditions();
         fvOptions.correct(U2);
     }
 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 1af167cc40..f269a79760 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -31,6 +31,8 @@ surfaceScalarField alpharAUf2
     fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2)
 );
 
+volScalarField Kd(fluid.Kd());
+
 // Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
 tmp<surfaceScalarField> phiF1;
 tmp<surfaceScalarField> phiF2;