diff --git a/applications/solvers/incompressible/simpleFoam/UEqn.H b/applications/solvers/incompressible/simpleFoam/UEqn.H
index 62c9f5536..a006813c4 100644
--- a/applications/solvers/incompressible/simpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/UEqn.H
@@ -1,16 +1,20 @@
     // Solve the momentum equation
 
-    tmp<fvVectorMatrix> UEqn
+    tmp<fvVectorMatrix> HUEqn
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
     );
 
-    UEqn().relax();
+    // Get under-relaxation factor
+    const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
 
+    // Momentum solution
     eqnResidual = solve
     (
-        UEqn() == -fvc::grad(p)
+        relax(HUEqn(), UUrf)
+     ==
+        -fvc::grad(p)
     ).initialResidual();
 
     maxResidual = max(eqnResidual, maxResidual);
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 42ebfc225..021350f09 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -1,17 +1,30 @@
     p.boundaryField().updateCoeffs();
 
-    volScalarField AU = UEqn().A();
-    U = UEqn().H()/AU;
-    UEqn.clear();
+    // Prepare clean 1/Ap without contribution from under-relaxation
+    // HJ, 26/Oct/2015
+    volScalarField rUA
+    (
+        "(1|A(U))",
+        1/HUEqn().A()
+    );
+
+    // 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();
+
+    // Global flux balance
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(1.0/AU, p) == fvc::div(phi)
+            fvm::laplacian(rUA, p) == fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -39,6 +52,8 @@
     p.relax();
 
     // Momentum corrector
-    U -= fvc::grad(p)/AU;
+    // 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();
 
diff --git a/applications/solvers/incompressible/simpleFoam/simpleFoam.C b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
index 5e80c10f7..5a6e7ef7e 100644
--- a/applications/solvers/incompressible/simpleFoam/simpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
@@ -26,6 +26,10 @@ Application
 
 Description
     Steady-state solver for incompressible, turbulent flow
+    Consistent formulation without time-step and relaxation dependence by Jasak
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
index c9329c21e..514e3240f 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
@@ -21,10 +21,7 @@ ddtSchemes
 
 gradSchemes
 {
-    default         Gauss linear;
-    grad(p)         Gauss linear;
-    grad(U)         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
+    default         cellLimited leastSquares 1;
 }
 
 divSchemes
@@ -39,8 +36,6 @@ divSchemes
 laplacianSchemes
 {
     default         Gauss linear corrected;
-//    default         Gauss linear limited 0.5;
-//    default         Gauss linear limited 0.333;
 }
 
 interpolationSchemes
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
index 61742eba9..928d503a6 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
@@ -60,7 +60,7 @@ solvers
 
 SIMPLE
 {
-    nNonOrthogonalCorrectors 0;
+    nNonOrthogonalCorrectors 2;
 }
 
 relaxationFactors