diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
index 1cee6d0481..6c765ec66b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
@@ -30,6 +30,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
         U1Eqn.relax();
         U1Eqn += fvm::Sp(dragCoeff, U1);
         fvOptions.constrain(U1Eqn);
+        U1.correctBoundaryConditions();
     }
 
     {
@@ -52,5 +53,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
         U2Eqn.relax();
         U2Eqn += fvm::Sp(dragCoeff, U2);
         fvOptions.constrain(U2Eqn);
+        U2.correctBoundaryConditions();
     }
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 0924c5d178..b0a5ae9e7c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -13,6 +13,12 @@ while (pimple.correct())
     // Update continuity errors due to temperature changes
     #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());
+
     volVectorField HbyA1
     (
         IOobject::groupName("HbyA", phase1.name()),
@@ -92,9 +98,9 @@ while (pimple.correct())
     surfaceScalarField phiCorrCoeff1(pos(alpha1fBar - 0.99));
     surfaceScalarField phiCorrCoeff2(pos(0.01 - alpha1fBar));
 
-    // Set ddtPhiCorr to 0 on non-coupled boundaries
     forAll(mesh.boundary(), patchi)
     {
+        // Set ddtPhiCorr to 0 on non-coupled boundaries
         if
         (
             !mesh.boundary()[patchi].coupled()
@@ -272,18 +278,7 @@ while (pimple.correct())
                     ((phi1s + D1f*phi2s) - (phi2s + D2f*phi1s))/(1 - D1f*D2f)
                 );
 
-                phi1.boundaryField() ==
-                    mrfZones.relative
-                    (
-                        mesh.Sf().boundaryField() & U1.boundaryField()
-                    );
                 phi1 = phi + alpha2f*phir;
-
-                phi2.boundaryField() ==
-                    mrfZones.relative
-                    (
-                        mesh.Sf().boundaryField() & U2.boundaryField()
-                    );
                 phi2 = phi - alpha1f*phir;
             }
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index f9effa437b..f518de6177 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -37,6 +37,7 @@ Description
 #include "IOMRFZoneList.H"
 #include "fvIOoptionList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
+#include "fixedValueFvsPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
index a41165d8ee..93d9bb4c78 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -60,13 +60,15 @@ Foam::BlendedInterfacialModel<modelType>::BlendedInterfacialModel
     const blendingMethod& blending,
     const phasePair& pair,
     const orderedPhasePair& pair1In2,
-    const orderedPhasePair& pair2In1
+    const orderedPhasePair& pair2In1,
+    const bool correctFixedFluxBCs
 )
 :
     pair_(pair),
     pair1In2_(pair1In2),
     pair2In1_(pair2In1),
-    blending_(blending)
+    blending_(blending),
+    correctFixedFluxBCs_(correctFixedFluxBCs)
 {
     if (modelTable.found(pair_))
     {
@@ -161,7 +163,11 @@ Foam::BlendedInterfacialModel<modelType>::K() const
         x() += model2In1_->K()*f2;
     }
 
-    if (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
     {
         correctFixedFluxBCs(x());
     }
@@ -217,7 +223,11 @@ Foam::BlendedInterfacialModel<modelType>::F() const
         x() -= model2In1_->F()*f2; // note : subtraction
     }
 
-    if (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
     {
         correctFixedFluxBCs(x());
     }
@@ -272,7 +282,11 @@ Foam::BlendedInterfacialModel<modelType>::D() const
         x() += model2In1_->D()*f2;
     }
 
-    if (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
     {
         correctFixedFluxBCs(x());
     }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H
index aa3bb3df2e..f56ab025b6 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H
@@ -75,6 +75,9 @@ class BlendedInterfacialModel
         //- Blending model
         const blendingMethod& blending_;
 
+        //- If true set coefficients and forces to 0 at fixed-flux BCs
+        bool correctFixedFluxBCs_;
+
 
     // Private Member Functions
 
@@ -103,7 +106,8 @@ public:
             const blendingMethod& blending,
             const phasePair& pair,
             const orderedPhasePair& pair1In2,
-            const orderedPhasePair& pair2In1
+            const orderedPhasePair& pair2In1,
+            const bool correctFixedFluxBCs = true
         );
 
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index d7a58a3809..32ca0ace85 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -184,7 +184,8 @@ Foam::twoPhaseSystem::twoPhaseSystem
             ),
             pair_,
             pair1In2_,
-            pair2In1_
+            pair2In1_,
+            false // Do not zero drag coefficent at fixed-flux BCs
         )
     );
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/T.water b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/T.water
index d19f29f191..3c62df75d0 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/T.water
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/T.water
@@ -28,7 +28,7 @@ boundaryField
     {
         type               inletOutlet;
         phi                phi.water;
-        inletValue         uniform 300;
+        inletValue         $internalField;
         value              $internalField;
     }
     inlet
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution
index 1e0e5d0f07..de1aecacf9 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution
@@ -35,7 +35,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/T.water b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/T.water
index d19f29f191..3c62df75d0 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/T.water
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/T.water
@@ -28,7 +28,7 @@ boundaryField
     {
         type               inletOutlet;
         phi                phi.water;
-        inletValue         uniform 300;
+        inletValue         $internalField;
         value              $internalField;
     }
     inlet
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution
index 9868241ad9..f6f72076f0 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution
@@ -35,7 +35,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution
index 0328d98737..49d9bfea71 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution
@@ -42,7 +42,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/T.water b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/T.water
index 1cfd38f926..e592c501d5 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/T.water
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/T.water
@@ -28,7 +28,7 @@ boundaryField
     {
         type               inletOutlet;
         phi                phi.water;
-        inletValue         uniform 300;
+        inletValue         $internalField;
         value              $internalField;
     }
     inlet
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
index 6487d35276..350b0c0c1f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
@@ -35,7 +35,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water
index d19f29f191..3c62df75d0 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/T.water
@@ -28,7 +28,7 @@ boundaryField
     {
         type               inletOutlet;
         phi                phi.water;
-        inletValue         uniform 300;
+        inletValue         $internalField;
         value              $internalField;
     }
     inlet
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
index 171537ad31..28572b1e0a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
@@ -35,7 +35,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution
index 5ca69b0fcd..d5d97565bd 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution
@@ -42,7 +42,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/constant/fvOptions b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/constant/fvOptions
index 5bcb5c0dbb..f231c3d605 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/constant/fvOptions
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/constant/fvOptions
@@ -23,7 +23,7 @@ injector1
     selectionMode   points;
     points
     (
-        (0.075 0.1 0.05)
+        (0.075 0.2 0.05)
     );
 }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution
index 6487d35276..350b0c0c1f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution
@@ -35,7 +35,7 @@ solvers
         agglomerator    faceAreaPair;
         mergeLevels     1;
         tolerance       1e-8;
-        relTol          0.01;
+        relTol          0;
     }
 
     p_rghFinal
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution
index 59c0cbf3ab..a95a952401 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution
@@ -25,17 +25,21 @@ solvers
 
     p_rgh
     {
-        solver          GAMG;
-        smoother        DIC;
-        nPreSweeps      0;
-        nPostSweeps     2;
-        nFinestSweeps   2;
-        cacheAgglomeration true;
-        nCellsInCoarsestLevel 10;
-        agglomerator    faceAreaPair;
-        mergeLevels     1;
-        tolerance       1e-7;
-        relTol          0.01;
+        solver          PCG;
+        preconditioner
+        {
+            preconditioner  GAMG;
+            smoother        GaussSeidel;
+            nPreSweeps      0;
+            nPostSweeps     2;
+            nFinestSweeps   2;
+            cacheAgglomeration true;
+            nCellsInCoarsestLevel 10;
+            agglomerator    faceAreaPair;
+            mergeLevels     1;
+        }
+        tolerance       1e-6;
+        relTol          0;
     }
 
     p_rghFinal