diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index 819e4ed4294d97bde742bbcbc354d6342b54dfc0..e9567bf7e22df20fa444f3eb1096998eb4afd397 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@@ -22,8 +22,6 @@ volScalarField rAU2
     )
 );
 
-//surfaceScalarField alpharAUf1(fvc::interpolate(alpha1*rAU1));
-//surfaceScalarField alpharAUf2(fvc::interpolate(alpha2*rAU2));
 surfaceScalarField alpharAUf1
 (
     fvc::interpolate(max(alpha1, fluid.residualAlpha(phase1))*rAU1)
@@ -94,10 +92,26 @@ while (pimple.correct())
     #include "correctContErrs.H"
 
     // Correct flux BCs to be consistent with the velocity BCs
-    phi1.boundaryField() ==
-        mrfZones.relative(mesh.Sf().boundaryField() & U1.boundaryField());
-    phi2.boundaryField() ==
-        mrfZones.relative(mesh.Sf().boundaryField() & U2.boundaryField());
+    forAll(mesh.boundary(), patchi)
+    {
+        if (isA<fixedValueFvsPatchScalarField>(phi1.boundaryField()[patchi]))
+        {
+            phi1.boundaryField()[patchi] ==
+                mrfZones.relative
+                (
+                    mesh.Sf().boundaryField() & U1.boundaryField()
+                )()[patchi];
+        }
+
+        if (isA<fixedValueFvsPatchScalarField>(phi2.boundaryField()[patchi]))
+        {
+            phi2.boundaryField()[patchi] ==
+                mrfZones.relative
+                (
+                    mesh.Sf().boundaryField() & U2.boundaryField()
+                )()[patchi];
+        }
+    }
 
     volVectorField HbyA1
     (
@@ -227,12 +241,9 @@ while (pimple.correct())
         p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
-          - mrfZones.relative
-            (
-                alphaf1.boundaryField()
-               *(mesh.Sf().boundaryField() & U1.boundaryField())
-              + alphaf2.boundaryField()
-               *(mesh.Sf().boundaryField() & U2.boundaryField())
+          - (
+                alphaf1.boundaryField()*phi1.boundaryField()
+              + alphaf2.boundaryField()*phi2.boundaryField()
             )
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index 5d39ce0d7c668fcc83eba8d192c7bc8a69eb5fb8..5f5089131102a9c6ac928cd25b4e4eeaf8aff9cf 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@@ -98,11 +98,27 @@ while (pimple.correct())
     surfaceScalarField rhof1(fvc::interpolate(rho1));
     surfaceScalarField rhof2(fvc::interpolate(rho2));
 
-    // Correct flux BCs to be consistent with the velocity BCs
-    phi1.boundaryField() ==
-        mrfZones.relative(mesh.Sf().boundaryField() & U1.boundaryField());
-    phi2.boundaryField() ==
-        mrfZones.relative(mesh.Sf().boundaryField() & U2.boundaryField());
+    // Correct fixed-flux BCs to be consistent with the velocity BCs
+    forAll(mesh.boundary(), patchi)
+    {
+        if (isA<fixedValueFvsPatchScalarField>(phi1.boundaryField()[patchi]))
+        {
+            phi1.boundaryField()[patchi] ==
+                mrfZones.relative
+                (
+                    mesh.Sf().boundaryField() & U1.boundaryField()
+                )()[patchi];
+        }
+
+        if (isA<fixedValueFvsPatchScalarField>(phi2.boundaryField()[patchi]))
+        {
+            phi2.boundaryField()[patchi] ==
+                mrfZones.relative
+                (
+                    mesh.Sf().boundaryField() & U2.boundaryField()
+                )()[patchi];
+        }
+    }
 
     surfaceScalarField alpharAUf1
     (
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
index 2b99be00b8c7c8dd053a11bbb986f75cb3d57c99..5a861911f9740d266c63ef55fc6f14bba8619da7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
@@ -30,6 +30,7 @@ License
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
+#include "fixedValueFvsPatchFields.H"
 #include "fixedValueFvPatchFields.H"
 #include "slipFvPatchFields.H"
 #include "partialSlipFvPatchFields.H"
@@ -151,7 +152,7 @@ Foam::phaseModel::phaseModel
              || isA<partialSlipFvPatchVectorField>(U_.boundaryField()[i])
             )
             {
-                phiTypes[i] = fixedValueFvPatchScalarField::typeName;
+                phiTypes[i] = fixedValueFvsPatchScalarField::typeName;
             }
         }