diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 4a817a1001..822e60d19a 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -76,6 +76,7 @@
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
           + rAlphaAUfs[phasei]*(g & mesh.Sf())
         );
+
         mrfZones.relativeFlux(phiHbyAs[phasei]);
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
@@ -89,11 +90,20 @@
             ++dmIter, ++dcIter
         )
         {
-            const phaseModel *phase2Ptr = &dmIter()->phase2();
-            if (phase2Ptr == &phase)
+            const phaseModel *phase2Ptr = NULL;
+
+            if (&phase == &dmIter()->phase1())
+            {
+                phase2Ptr = &dmIter()->phase2();
+            }
+            else if (&phase == &dmIter()->phase2())
             {
                 phase2Ptr = &dmIter()->phase1();
             }
+            else
+            {
+                continue;
+            }
 
             phiHbyAs[phasei] +=
                 fvc::interpolate