diff --git a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
index d73d4b5c1d72e1a9dc72f30c55484cd467caa6d3..453ca45dd2969b8acba88d58a531f506796a017a 100644
--- a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
+++ b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
@@ -472,8 +472,8 @@ tmp<surfaceScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixFlux
     surfaceScalarField rhogEfff(fvc::interpolate(rhogEff()));
 
     return
-       fvc::interpolate(rhom_()/(alphal*rholEff() + alphag*rhogEff()*Ct2_()))
-      *(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd);
+       (alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd)
+      /(alphalf*rholEfff + alphagf*rhogEfff*fvc::interpolate(Ct2_()));
 }
 
 
@@ -516,14 +516,14 @@ tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
 template<class BasicTurbulenceModel>
 tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::kSource() const
 {
-    return fvm::Su(bubbleG(), km_());
+    return fvm::Su(bubbleG()/rhom_(), km_());
 }
 
 
 template<class BasicTurbulenceModel>
 tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::epsilonSource() const
 {
-    return fvm::Su(C3_*epsilonm_()*bubbleG()/km_(), epsilonm_());
+    return fvm::Su(C3_*epsilonm_()*bubbleG()/(rhom_()*km_()), epsilonm_());
 }
 
 
@@ -647,14 +647,14 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correct()
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
     (
-        fvm::ddt(rhom, epsilonm)
+        fvm::ddt(epsilonm)
       + fvm::div(phim, epsilonm)
-      - fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), epsilonm)
-      - fvm::laplacian(DepsilonEff(rhom*nutm), epsilonm)
+      - fvm::Sp(fvc::div(phim), epsilonm)
+      - fvm::laplacian(DepsilonEff(nutm), epsilonm)
      ==
-        C1_*rhom*Gm*epsilonm/km
-      - fvm::SuSp(((2.0/3.0)*C1_)*rhom*divUm, epsilonm)
-      - fvm::Sp(C2_*rhom*epsilonm/km, epsilonm)
+        C1_*Gm*epsilonm/km
+      - fvm::SuSp(((2.0/3.0)*C1_)*divUm, epsilonm)
+      - fvm::Sp(C2_*epsilonm/km, epsilonm)
       + epsilonSource()
     );
 
@@ -669,14 +669,14 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correct()
     // Turbulent kinetic energy equation
     tmp<fvScalarMatrix> kmEqn
     (
-        fvm::ddt(rhom, km)
+        fvm::ddt(km)
       + fvm::div(phim, km)
-      - fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), km)
-      - fvm::laplacian(DkEff(rhom*nutm), km)
+      - fvm::Sp(fvc::div(phim), km)
+      - fvm::laplacian(DkEff(nutm), km)
      ==
-        rhom*Gm
-      - fvm::SuSp((2.0/3.0)*rhom*divUm, km)
-      - fvm::Sp(rhom*epsilonm/km, km)
+        Gm
+      - fvm::SuSp((2.0/3.0)*divUm, km)
+      - fvm::Sp(epsilonm/km, km)
       + kSource()
     );