diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
index bac588c73895f5eac7fff922f101cd7edb72b631..3f659d44ec4cfbe0b2a65cd35545f137870c7f66 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
@@ -76,34 +76,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField bp(pow(beta, -2.65));
     volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
 
-    volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     // Wen and Yu (1966)
-    tmp<volScalarField> tKWenYu(0.75*Cds*phase2_.rho()*Ur*bp/d);
-    volScalarField& KWenYu = tKWenYu();
-
-    // Ergun
-    forAll (beta, cellj)
-    {
-        if (beta[cellj] <= 0.8)
-        {
-            KWenYu[cellj] =
-                150.0*alpha_[cellj]*phase2_.nu().value()*phase2_.rho().value()
-               /sqr(beta[cellj]*d[cellj])
-              + 1.75*phase2_.rho().value()*Ur[cellj]
-               /(beta[cellj]*d[cellj]);
-        }
-    }
-
-    return tKWenYu;
+    return
+    (
+        pos(beta - 0.8)
+       *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+      + neg(beta - 0.8)
+       *(
+           150.0*alpha_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+         + 1.75*phase2_.rho()*Ur/(beta*d)
+        )
+    );
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
index 99b5cacdf2fdc150df0039ce77247d99c628abb5..0a624c6b0a2a8a7bb8e8f4ba552b83eb13b3b967 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
@@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
index 79259cbe05858bdb9138023972809325a2df5d01..f68e5379d2403e64ead7aadf0ac369d5bfa5ed8d 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
@@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
index 4f971b4bb7639a29991a074ff7fc48d75f6a5ecb..541868cf870c581b0a50e8bf44e29e7593e11bc8 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
@@ -73,15 +73,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
 {
     volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
     volScalarField A(pow(beta, 4.14));
-    volScalarField B(0.8*pow(beta, 1.28));
-
-    forAll (beta, celli)
-    {
-        if (beta[celli] > 0.85)
-        {
-            B[celli] = pow(beta[celli], 2.65);
-        }
-    }
+    volScalarField B
+    (
+        neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+      + pos(beta - 0.85)*(pow(beta, 2.65))
+    );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
index f1126125bb1be131e9f60f34bf83158510683a92..bf279a3990a0d065d2a680896bb78e2212ef9f0b 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
@@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 4320bbad354f0329c8b9cbc253c92c0e681e85c8..d458f07b754218cb44de74b2cbe483b68626f188 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
+    muff.correctBoundaryConditions();
+
     return tmuf;
 }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
index e34f219ebc0cbefa3972d31399b26ff852b52fbf..b52b766de9238653fdc699a75c6baa14e914cbc6 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
@@ -75,34 +75,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField bp(pow(beta, -2.65));
     volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
 
-    volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     // Wen and Yu (1966)
-    tmp<volScalarField> tKWenYu = 0.75*Cds*phase2_.rho()*Ur*bp/d;
-    volScalarField& KWenYu = tKWenYu();
-
-    // Ergun
-    forAll (beta, cellj)
-    {
-        if (beta[cellj] <= 0.8)
-        {
-            KWenYu[cellj] =
-                150.0*phase1_[cellj]*phase2_.nu().value()*phase2_.rho().value()
-               /sqr(beta[cellj]*d[cellj])
-              + 1.75*phase2_.rho().value()*Ur[cellj]
-               /(beta[cellj]*d[cellj]);
-        }
-    }
-
-    return tKWenYu;
+    return
+    (
+        pos(beta - 0.8)
+       *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+      + neg(beta - 0.8)
+       *(
+           150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+         + 1.75*phase2_.rho()*Ur/(beta*d)
+        )
+    );
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
index ecd25eb386fe6523e6266993e4df67072d1b34f3..25de6b1dc7001cd66e7094996655a3913dca7bed 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
@@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
index 3259b5d0e08b708a80e1f4b66c42322248abefa5..d4d77a2bb35df65363b95393b57ffd8d52d9c19e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
@@ -71,15 +71,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
index 736be885e0fdaa596b12b592989100c2870d5305..07f85b38a1f8e04809689ea6e1d94e13255552a8 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
@@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
 {
     volScalarField beta(max(phase2_, scalar(1.0e-6)));
     volScalarField A(pow(beta, 4.14));
-    volScalarField B(0.8*pow(beta, 1.28));
-
-    forAll (beta, celli)
-    {
-        if (beta[celli] > 0.85)
-        {
-            B[celli] = pow(beta[celli], 2.65);
-        }
-    }
+    volScalarField B
+    (
+        neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+      + pos(beta - 0.85)*(pow(beta, 2.65))
+    );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
index 6cfc119c1c4ed8e2e79ae7e7332df7a17eba06dc..6319276408b935a09c5102ed726af4fb54004c9c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
@@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 4320bbad354f0329c8b9cbc253c92c0e681e85c8..d458f07b754218cb44de74b2cbe483b68626f188 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
+    muff.correctBoundaryConditions();
+
     return tmuf;
 }