diff --git a/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C
index 7efc1b848b3e1531b4bd0a38fc4b292662b2df50..332c7877a822290a9c9fd58d82f95d52b1452c54 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C
@@ -251,7 +251,6 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    const volScalarField::Internal unlimitedNut(Cmu_*sqr(k_())/epsilon_());
     const volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
@@ -269,7 +268,7 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
     );
     tgradU.clear();
 
-    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
+    const volScalarField::Internal G(this->GName(), nut()*GbyNu);
 
     const volScalarField::Internal eta(sqrt(mag(GbyNu))*k_/epsilon_);
     const volScalarField::Internal eta3(eta*sqr(eta));
@@ -289,7 +288,7 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
-        (C1_ - R)*alpha()*rho()*G*epsilon_()/k_()
+        (C1_ - R)*alpha()*rho()*GbyNu*Cmu_*k_()
       - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
       - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
       + epsilonSource()
@@ -312,7 +311,7 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha()*rho()*GbyNu*nut()
+        alpha()*rho()*G
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
       - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
index 33704af01f49fe7d52efbc039b8cb6425ec040a8..c77aeeaad463d37da51a977fce229dc0c2e432b5 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C
@@ -231,7 +231,6 @@ void kEpsilon<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    const volScalarField::Internal unlimitedNut(Cmu_*sqr(k_())/epsilon_());
     const volScalarField& nut = this->nut_;
 
     fv::options& fvOptions(fv::options::New(this->mesh_));
@@ -248,7 +247,7 @@ void kEpsilon<BasicTurbulenceModel>::correct()
     (
         tgradU().v() && dev(twoSymm(tgradU().v()))
     );
-    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
+    const volScalarField::Internal G(this->GName(), nut()*GbyNu);
     tgradU.clear();
 
     // Update epsilon and G at the wall
@@ -261,7 +260,7 @@ void kEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
-        C1_*alpha()*rho()*G*epsilon_()/k_()
+        C1_*alpha()*rho()*GbyNu*Cmu_*k_()
       - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
       - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
       + epsilonSource()
@@ -282,7 +281,7 @@ void kEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha()*rho()*GbyNu*nut()
+        alpha()*rho()*G
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
       - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C
index 89ceedb5bb80a54708bb0e29baa3ccd8f0514a83..2a547f53ee7197e0935627ed9e4a3a3d5b71b767 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C
@@ -191,7 +191,6 @@ void kOmega<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    const volScalarField::Internal unlimitedNut(k_()/omega_());
     const volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
@@ -207,7 +206,7 @@ void kOmega<BasicTurbulenceModel>::correct()
     (
         tgradU().v() && dev(twoSymm(tgradU().v()))
     );
-    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
+    const volScalarField::Internal G(this->GName(), nut()*GbyNu);
     tgradU.clear();
 
     // Update omega and G at the wall
@@ -220,7 +219,7 @@ void kOmega<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, omega_)
       - fvm::laplacian(alpha*rho*DomegaEff(), omega_)
      ==
-        gamma_*alpha()*rho()*G*omega_()/k_()
+        gamma_*alpha()*rho()*GbyNu
       - fvm::SuSp(((2.0/3.0)*gamma_)*alpha()*rho()*divU, omega_)
       - fvm::Sp(beta_*alpha()*rho()*omega_(), omega_)
       + fvOptions(alpha, rho, omega_)
@@ -241,7 +240,7 @@ void kOmega<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha()*rho()*GbyNu*nut()
+        alpha()*rho()*G
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
       - fvm::Sp(Cmu_*alpha()*rho()*omega_(), k_)
       + fvOptions(alpha, rho, k_)