diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H
index b4c78c0ac..1b7c2d865 100644
--- a/applications/solvers/incompressible/simpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/UEqn.H
@@ -1,18 +1,13 @@
     // Solve the momentum equation
 
-    tmp<fvVectorMatrix> HUEqn
+    tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
       + turbulence->divDevReff()
     );
 
-    // Get under-relaxation factor
-    const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(U.name());
+    // Relax the equation
+    UEqn().relax();
 
     // Momentum solution
-    solve
-    (
-        relax(HUEqn(), UUrf)
-     ==
-        -fvc::grad(p)
-    );
+    solve(UEqn() == -fvc::grad(p));
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 9a1c0a4ef..6193b94ea 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -1,20 +1,11 @@
     p.boundaryField().updateCoeffs();
 
-    // Prepare clean 1/Ap without contribution from under-relaxation
-    // HJ, 26/Oct/2015
-    volScalarField rUA
-    (
-        "(1|A(U))",
-        1/HUEqn().A()
-    );
+    volScalarField rAU = 1.0/UEqn().A();
+    U = rAU*UEqn().H();
+    UEqn.clear();
 
-    // Store velocity under-relaxation point before using U for
-    // the flux precursor
-    U.storePrevIter();
-
-    U = rUA*HUEqn().H();
-    HUEqn.clear();
-    phi = fvc::interpolate(U) & mesh.Sf();
+    // Calculate under-relaxation consistent flux
+    simple.calcSteadyConsistentFlux(phi, U);
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
@@ -22,7 +13,14 @@
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUA, p) == fvc::div(phi)
+            fvm::laplacian
+            (
+                fvc::interpolate(rAU)/simple.aCoeff(),
+                p,
+                "laplacian(rAU" + p.name() + ')'
+            )
+         ==
+            fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -40,9 +38,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
-    U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
-    U.correctBoundaryConditions();
-
+    // Reconstruct consistent velocity after pressure equation
+    simple.reconstructSteadyVelocity(U, rAU, p);
diff --git a/applications/solvers/incompressible/simpleFoam/simpleFoam.C b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
index 5a90dc7a0..4b2a148be 100644
--- a/applications/solvers/incompressible/simpleFoam/simpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
@@ -25,8 +25,9 @@ Application
     simpleFoam
 
 Description
-    Steady-state solver for incompressible, turbulent flow
+    Steady-state solver for incompressible, turbulent flow.
     Consistent formulation without time-step and relaxation dependence by Jasak
+    and Tukovic.
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved