diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
index 1e6809288d..f6466e271d 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
@@ -272,6 +272,7 @@ void ShihQuadraticKE::correct()
         (nut_*twoSymm(gradU) - nonlinearStress_) && gradU
     );
 
+
     // Update epsilon and G at the wall
     epsilon_.boundaryField().updateCoeffs();
 
@@ -287,9 +288,7 @@ void ShihQuadraticKE::correct()
     );
 
     epsEqn().relax();
-
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
-
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
 
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallDissipationI.H b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallDissipationI.H
deleted file mode 100644
index a519333622..0000000000
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallDissipationI.H
+++ /dev/null
@@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    wallDissipation
-
-Description
-    Set wall dissipation in the epsilon matrix
-
-\*---------------------------------------------------------------------------*/
-
-{
-    const fvPatchList& patches = mesh_.boundary();
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& p = patches[patchi];
-
-        if (isA<wallFvPatch>(p))
-        {
-            epsEqn().setValues
-            (
-                p.faceCells(),
-                epsilon_.boundaryField()[patchi].patchInternalField()
-            );
-        }
-    }
-}
-
-// ************************************************************************* //
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/qZeta/qZeta.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/qZeta/qZeta.C
index 7d62a62596..fd3461af49 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/qZeta/qZeta.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/qZeta/qZeta.C
@@ -253,7 +253,6 @@ void qZeta::correct()
 
 
     // Zeta equation
-
     tmp<fvScalarMatrix> zetaEqn
     (
         fvm::ddt(zeta_)
@@ -271,7 +270,6 @@ void qZeta::correct()
 
 
     // q equation
-
     tmp<fvScalarMatrix> qEqn
     (
         fvm::ddt(q_)
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
index 15b538cf93..5303fa1f58 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
@@ -256,15 +256,12 @@ void kEpsilon<BasicTurbulenceModel>::correct()
     );
 
     epsEqn().relax();
-
     epsEqn().boundaryManipulate(epsilon_.boundaryField());
-
     solve(epsEqn);
     bound(epsilon_, this->epsilonMin_);
 
 
     // Turbulent kinetic energy equation
-
     tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(alpha, rho, k_)
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
index 536df15f50..64c710edc5 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
@@ -83,29 +83,30 @@ void epsilonLowReWallFunctionFvPatchScalarField::calculate
     const scalarField magGradUw(mag(Uw.snGrad()));
 
     // Set epsilon and G
-    forAll(nutw, faceI)
+    forAll(nutw, facei)
     {
-        label cellI = patch.faceCells()[faceI];
+        label celli = patch.faceCells()[facei];
 
-        scalar yPlus = Cmu25*sqrt(k[cellI])*y[faceI]/nuw[faceI];
+        scalar yPlus = Cmu25*sqrt(k[celli])*y[facei]/nuw[facei];
 
-        scalar w = cornerWeights[faceI];
+        scalar w = cornerWeights[facei];
 
         if (yPlus > yPlusLam_)
         {
-            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+            epsilon[celli] = w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
         }
         else
         {
-            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
+            epsilon[celli] = w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
         }
 
-        G[cellI] =
+        // It is not clear that G should be adjusted for low-Re BCs
+        G[celli] +=
             w
-           *(nutw[faceI] + nuw[faceI])
-           *magGradUw[faceI]
-           *Cmu25*sqrt(k[cellI])
-           /(kappa_*y[faceI]);
+           *(nutw[facei] + nuw[facei])
+           *magGradUw[facei]
+           *Cmu25*sqrt(k[celli])
+           /(kappa_*y[facei]);
     }
 }
 
diff --git a/tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/0/k b/tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/0/k
index becf57028a..751ffd13b8 100644
--- a/tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/0/k
+++ b/tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/0/k
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [ 0 2 -2 0 0 0 0 ];
 
-internalField   uniform 1;
+internalField   uniform 0.1;
 
 boundaryField
 {
diff --git a/tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/fvSchemes b/tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/fvSchemes
index 0fd54ce506..a90d7abcb8 100644
--- a/tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/fvSchemes
+++ b/tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/fvSchemes
@@ -33,6 +33,7 @@ divSchemes
     div(phi,R)      bounded Gauss linear;
     div(phi,nuTilda) bounded Gauss linear;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
+    div(nonlinearStress) Gauss linear;
 }
 
 laplacianSchemes