diff --git a/applications/solvers/incompressible/simpleSRFFoam/UEqn.H b/applications/solvers/incompressible/simpleSRFFoam/UEqn.H
index a10560319..10fcbf895 100644
--- a/applications/solvers/incompressible/simpleSRFFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleSRFFoam/UEqn.H
@@ -1,15 +1,11 @@
     // Solve the momentum equation
 
-    tmp<fvVectorMatrix> HUrelEqn
+    tmp<fvVectorMatrix> UrelEqn
     (
         fvm::div(phi, Urel)
       + turbulence->divDevReff()
       + SRF->Su()
     );
 
-    // Get under-relaxation factor
-    const scalar UUrf =
-        mesh.solutionDict().equationRelaxationFactor(Urel.name());
-
     // Momentum solution
-    solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p));
+    solve(relax(UrelEqn()) == -fvc::grad(p));
diff --git a/applications/solvers/incompressible/simpleSRFFoam/pEqn.H b/applications/solvers/incompressible/simpleSRFFoam/pEqn.H
index 04e6813d2..f506264ff 100644
--- a/applications/solvers/incompressible/simpleSRFFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleSRFFoam/pEqn.H
@@ -1,36 +1,32 @@
     p.boundaryField().updateCoeffs();
 
-    // Prepare clean 1/Ap without contribution from under-relaxation
-    // HJ, 26/Oct/2015
-    volScalarField rUA
-    (
-        "(1|A(Urel))",
-        1/HUrelEqn().A()
-    );
+    volScalarField rAUrel = 1.0/UrelEqn().A();
+    Urel = rAUrel*UrelEqn().H();
+    UrelEqn.clear();
 
-    // Store velocity under-relaxation point before using U for the flux
-    // precursor
-    Urel.storePrevIter();
-
-    Urel = rUA*HUrelEqn().H();
-    HUrelEqn.clear();
-    phi = fvc::interpolate(Urel) & mesh.Sf();
-
-    // Global flux balance
+    // Calculate under-relaxation consistent flux
+    simple.calcSteadyConsistentFlux(phi, Urel);
     adjustPhi(phi, Urel, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUA, p) == fvc::div(phi)
+            fvm::laplacian
+            (
+                fvc::interpolate(rAUrel)/simple.aCoeff(),
+                p,
+                "laplacian(rAU," + p.name() + ')'
+            )
+          ==
+            fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }
@@ -41,8 +37,5 @@
     // Explicitly relax pressure for momentum corrector
     p.relax();
 
-    // Momentum corrector
-    // Note: since under-relaxation does not change aU, H/a in U can be
-    // re-used.  HJ, 22/Jan/2016
-    Urel = UUrf*(Urel - rUA*fvc::grad(p)) + (1 - UUrf)*Urel.prevIter();
-    Urel.correctBoundaryConditions();
+    // Reconstruct consistent velocity after pressure equation
+    simple.reconstructSteadyVelocity(Urel, rAUrel, p);
diff --git a/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C b/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C
index c1106ad72..ec77f71f1 100644
--- a/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C
+++ b/applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C
@@ -65,50 +65,8 @@ int main(int argc, char *argv[])
 
         // Pressure-velocity SIMPLE corrector
         {
-            // Momentum predictor
-            tmp<fvVectorMatrix> UrelEqn
-            (
-                fvm::div(phi, Urel)
-              + turbulence->divDevReff()
-              + SRF->Su()
-            );
-
-            UrelEqn().relax();
-
-            solve(UrelEqn() == -fvc::grad(p));
-
-            p.boundaryField().updateCoeffs();
-            volScalarField AUrel = UrelEqn().A();
-            Urel = UrelEqn().H()/AUrel;
-            UrelEqn.clear();
-            phi = fvc::interpolate(Urel) & mesh.Sf();
-            adjustPhi(phi, Urel, p);
-
-            // Non-orthogonal pressure corrector loop
-            while (simple.correctNonOrthogonal())
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::laplacian(1.0/AUrel, p) == fvc::div(phi)
-                );
-
-                pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
-
-                if (simple.finalNonOrthogonalIter())
-                {
-                    phi -= pEqn.flux();
-                }
-            }
-
-#           include "continuityErrs.H"
-
-            // Explicitly relax pressure for momentum corrector
-            p.relax();
-
-            // Momentum corrector
-            Urel -= fvc::grad(p)/AUrel;
-            Urel.correctBoundaryConditions();
+#           include "UEqn.H"
+#           include "pEqn.H"
         }
 
         turbulence->correct();