diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
index 3911fa339d..b64277749e 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H
@@ -21,6 +21,10 @@
     (
         fvm::ddt(alpha1, he1)
       + fvm::div(alphaPhi1, he1)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), he1)
+
       - fvm::laplacian(k1, he1)
      ==
         heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
@@ -33,6 +37,10 @@
     (
         fvm::ddt(alpha2, he2)
       + fvm::div(alphaPhi2, he2)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), he2)
+
       - fvm::laplacian(k2, he2)
      ==
         heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
index a243db6300..ce5c622f9a 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H
@@ -31,6 +31,9 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
             fvm::ddt(alpha1, U1)
           + fvm::div(alphaPhi1, U1)
 
+            // Compressibity correction
+          - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
+
           + Cvm*rho2*alpha1*alpha2/rho1*
             (
                 fvm::ddt(U1)
@@ -61,6 +64,9 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
             fvm::ddt(alpha2, U2)
           + fvm::div(alphaPhi2, U2)
 
+            // Compressibity correction
+          - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
+
           + Cvm*rho2*alpha1*alpha2/rho2*
             (
                 fvm::ddt(U2)
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
index bc9a07b0a8..b74e1b18d7 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H
@@ -16,6 +16,10 @@ if (turbulence)
     (
         fvm::ddt(alpha2, epsilon)
       + fvm::div(alphaPhi2, epsilon)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), epsilon)
+
       - fvm::laplacian
         (
             alpha1Eps*nuEff2, epsilon,
@@ -41,6 +45,10 @@ if (turbulence)
     (
         fvm::ddt(alpha2, k)
       + fvm::div(alphaPhi2, k)
+
+        // Compressibity correction
+      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), k)
+
       - fvm::laplacian
         (
             alpha1k*nuEff2, k,
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
index 37a7524ee6..1851801037 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/bubbleColumn/system/fvSchemes
@@ -18,8 +18,6 @@ FoamFile
 ddtSchemes
 {
     default         Euler;
-
-    "ddt\(alpha.*,.*\)"  bounded Euler;
 }
 
 gradSchemes
@@ -34,15 +32,15 @@ divSchemes
     div(phi,alphaair)       Gauss vanLeer;
     div(phir,alphaair)      Gauss vanLeer;
 
-    "div\(alphaPhi.*,U.*\)" bounded Gauss limitedLinearV 1;
+    "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
     "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
     "div\(\(alpha.*Rc\)\)"  Gauss linear;
     "div\(phid.*,p\)"       Gauss upwind;
 
-    "div\(alphaPhi.*,h.*\)" bounded Gauss limitedLinear 1;
+    "div\(alphaPhi.*,h.*\)" Gauss limitedLinear 1;
     "div\(phi.*,K.*\)"      Gauss limitedLinear 1;
 
-    "div\(alphaPhi.*,(k|epsilon)\)"  bounded Gauss limitedLinear 1;
+    "div\(alphaPhi.*,(k|epsilon)\)"  Gauss limitedLinear 1;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
index 72cf3cb54b..8799476cc7 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/fluidisedBed/system/fvSchemes
@@ -18,8 +18,6 @@ FoamFile
 ddtSchemes
 {
     default             Euler;
-
-    "ddt\(alpha.,.*\)"  bounded Euler;
 }
 
 gradSchemes
@@ -34,16 +32,16 @@ divSchemes
     div(phi,alpha1)         Gauss vanLeer;
     div(phir,alpha1)        Gauss vanLeer;
 
-    "div\(alphaPhi.,U.\)"   bounded Gauss limitedLinearV 1;
+    "div\(alphaPhi.,U.\)"   Gauss limitedLinearV 1;
     "div\(phi.,U.\)"        Gauss limitedLinearV 1;
     "div\(\(alpha.*Rc\)\)"  Gauss linear;
     "div\(phid.,p\)"        Gauss upwind;
 
-    "div\(alphaPhi.,h.\)"   bounded Gauss limitedLinear 1;
+    "div\(alphaPhi.,h.\)"   Gauss limitedLinear 1;
     "div\(phi.,K.\)"        Gauss limitedLinear 1;
 
-    div(alphaPhi2,k)        bounded Gauss limitedLinear 1;
-    div(alphaPhi2,epsilon)  bounded Gauss limitedLinear 1;
+    div(alphaPhi2,k)        Gauss limitedLinear 1;
+    div(alphaPhi2,epsilon)  Gauss limitedLinear 1;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
index e1b193c952..c88d130244 100644
--- a/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/compressibleTwoPhaseEulerFoam/mixerVessel2D/system/fvSchemes
@@ -18,8 +18,6 @@ FoamFile
 ddtSchemes
 {
     default         Euler;
-
-    "ddt\(alpha.*,.*\)"  bounded Euler;
 }
 
 gradSchemes
@@ -34,15 +32,15 @@ divSchemes
     div(phi,alpha1)         Gauss vanLeer;
     div(phir,alpha1)        Gauss vanLeer;
 
-    "div\(alphaPhi.,U.\)"   bounded Gauss limitedLinearV 1;
+    "div\(alphaPhi.,U.\)"   Gauss limitedLinearV 1;
     "div\(phi.,U.\)"        Gauss limitedLinearV 1;
     "div\(\(alpha.*Rc\)\)"  Gauss linear;
     "div\(phid.,p\)"        Gauss linear;
 
-    "div\(alphaPhi.,h.\)"   bounded Gauss limitedLinear 1;
+    "div\(alphaPhi.,h.\)"   Gauss limitedLinear 1;
     "div\(phi.,K.\)"        Gauss linear;
 
-    "div\(alphaPhi.,(k|epsilon)\)"  bounded Gauss limitedLinear 1;
+    "div\(alphaPhi.,(k|epsilon)\)"  Gauss limitedLinear 1;
 }
 
 laplacianSchemes
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
index 56e0a545c1..bf47f69643 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
+++ b/tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/polyMesh/boundary
@@ -38,6 +38,7 @@ FoamFile
     defaultFaces
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3750;
         startFace       3850;
     }