diff --git a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
index da03517ee04da2cbf86745e076160da77ff93a59..8ab48d7b58b80f92f9b90920f1067dfe3c5a6cb6 100644
--- a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
+++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
@@ -87,7 +87,6 @@ int main(int argc, char *argv[])
             #include "alphaEqnSubCycle.H"
 
             twoPhaseProperties.correct();
-            Info<< average(twoPhaseProperties.mu()) << endl;
 
             #include "UEqn.H"
 
diff --git a/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
index 66216b87ce31a92283211024cc7dd3aad15c376d..ddcc8cf21067c9fd26262392e06f645d734bca4c 100644
--- a/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
+++ b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
@@ -10,8 +10,6 @@ if (turbulence)
     dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
     dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
 
-    volScalarField divU(fvc::div(phi));
-
     tmp<volTensorField> tgradU = fvc::grad(U);
     volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
@@ -21,7 +19,7 @@ if (turbulence)
         Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
     );
 
-    volScalarField muc(twoPhaseProperties.nucModel().nu()*rho2);
+    volScalarField mul(twoPhaseProperties.mu());
 
     #include "wallFunctions.H"
 
@@ -32,12 +30,12 @@ if (turbulence)
       + fvm::div(rhoPhi, epsilon)
       - fvm::laplacian
         (
-            mut/sigmaEps + muc, epsilon,
+            mut/sigmaEps + mul, epsilon,
             "laplacian(DepsilonEff,epsilon)"
         )
      ==
         C1*G*epsilon/(k + kMin)
-      - fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
+      - fvm::SuSp(C1*(1.0 - C3)*Gcoef, epsilon)
       - fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
     );
 
@@ -56,12 +54,12 @@ if (turbulence)
       + fvm::div(rhoPhi, k)
       - fvm::laplacian
         (
-            mut/sigmak + muc, k,
+            mut/sigmak + mul, k,
             "laplacian(DkEff,k)"
         )
      ==
         G
-      - fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
+      - fvm::SuSp(Gcoef, k)
       - fvm::Sp(rho*epsilon/(k + kMin), k)
     );
 
@@ -75,6 +73,10 @@ if (turbulence)
     mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
 
     #include "wallViscosity.H"
-}
 
-muEff = mut + twoPhaseProperties.mu();
+    muEff = mut + mul;
+}
+else
+{
+    muEff = mut + twoPhaseProperties.mu();
+}
diff --git a/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
index b9ff84817a36907a0a274d6fe99d81b6a53186bb..2a064b47f1b104f21f6fb30dcf87c309cdd967eb 100644
--- a/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
+++ b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
@@ -33,7 +33,7 @@
         if (isA<wallFvPatch>(curPatch))
         {
             const scalarField& mutw = mut.boundaryField()[patchi];
-            const scalarField& mucw = muc.boundaryField()[patchi];
+            const scalarField& mulw = mul.boundaryField()[patchi];
 
             scalarField magFaceGradU
             (
@@ -55,7 +55,7 @@
                    /(kappa_*y[patchi][facei]);
 
                 G[faceCelli] +=
-                    (mutw[facei] + mucw[facei])
+                    (mutw[facei] + mulw[facei])
                    *magFaceGradU[facei]
                    *Cmu25*::sqrt(k[faceCelli])
                    /(kappa_*y[patchi][facei]);
diff --git a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
index 2b54f6c2a85fcdceb09514441532437694dee30f..d73e0a5d1d7d6f4fc29ba807563d95fe112ee36c 100644
--- a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
+++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
@@ -12,7 +12,7 @@
         if (isA<wallFvPatch>(curPatch))
         {
             scalarField& mutw = mut.boundaryField()[patchi];
-            const scalarField& mucw = muc.boundaryField()[patchi];
+            const scalarField& mulw = mul.boundaryField()[patchi];
 
             forAll(curPatch, facei)
             {
@@ -20,12 +20,12 @@
 
                 scalar yPlus =
                     Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
-                   /(mucw[facei]/rho2.value());
+                   /(mulw[facei]/rho2.value());
 
                 if (yPlus > 11.6)
                 {
                     mutw[facei] =
-                        mucw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                        mulw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1);
                 }
                 else
                 {