diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
index f9dba26fd0..c72c82c678 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
@@ -47,10 +47,10 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     PtrList<surfaceScalarField> alphafs(phases.size());
     forAll(phases, phasei)
     {
-        phaseModel& phase = phases[phasei];
+        const phaseModel& phase = phases[phasei];
         const volScalarField& alpha = phase;
 
-        alphafs.set(phasei, fvc::interpolate(alpha).ptr());
+        alphafs.set(phasei, fvc::interpolate(max(alpha, scalar(0))).ptr());
         alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
     }
 
@@ -65,7 +65,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 
         forAll(fluid.movingPhases(), movingPhasei)
         {
-            phaseModel& phase = fluid.movingPhases()[movingPhasei];
+            const phaseModel& phase = fluid.movingPhases()[movingPhasei];
             const volScalarField& alpha = phase;
 
             volScalarField AU
@@ -103,7 +103,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     PtrList<surfaceScalarField> alpharAUfs(phases.size());
     forAll(phases, phasei)
     {
-        phaseModel& phase = phases[phasei];
+        const phaseModel& phase = phases[phasei];
         const volScalarField& alpha = phase;
 
         alpharAUfs.set
@@ -145,7 +145,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 
             forAll(phases, phasei)
             {
-                phaseModel& phase = phases[phasei];
+                const phaseModel& phase = phases[phasei];
 
                 phigFs.set
                 (
@@ -171,7 +171,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 
             forAll(fluid.movingPhases(), movingPhasei)
             {
-                phaseModel& phase = fluid.movingPhases()[movingPhasei];
+                const phaseModel& phase = fluid.movingPhases()[movingPhasei];
                 const volScalarField& alpha = phase;
                 const label phasei = phase.index();
 
@@ -260,17 +260,14 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                 mesh
             ),
             mesh,
-            dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), 0)
+            dimensionedScalar(dimTime/dimDensity, 0)
         );
 
         forAll(phases, phasei)
         {
             rAUf += alphafs[phasei]*alpharAUfs[phasei];
-            // rAUf += alphafs[phasei]*alphafs[phasei]*rAUfs[phasei];
         }
 
-        rAUf = mag(rAUf);
-
         // Update the fixedFluxPressure BCs to ensure flux consistency
         {
             surfaceScalarField::Boundary phib
@@ -282,7 +279,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 
             forAll(phases, phasei)
             {
-                phaseModel& phase = phases[phasei];
+                const phaseModel& phase = phases[phasei];
                 phib +=
                     alphafs[phasei].boundaryField()
                    *phase.phi()().boundaryField();
diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 2a2f797881..f8e71e6ad3 100644
--- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -195,6 +195,31 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
     )
     {
         *Kds_[dragModelIter.key()] = dragModelIter()->K();
+
+        // Zero-gradient the drag coefficient to boundaries with fixed velocity
+        const phaseInterface& interface = dragModelIter()->interface();
+        volScalarField& K = *Kds_[dragModelIter.key()];
+
+        forAll(K.boundaryField(), patchi)
+        {
+            if
+            (
+                (
+                    !interface.phase1().stationary()
+                 && interface.phase1().U()()
+                   .boundaryField()[patchi].fixesValue()
+                )
+             && (
+                    !interface.phase2().stationary()
+                 && interface.phase2().U()()
+                   .boundaryField()[patchi].fixesValue()
+                )
+            )
+            {
+                K.boundaryFieldRef()[patchi] =
+                    K.boundaryField()[patchi].patchInternalField();
+            }
+        }
     }
 
     // Initialise Vms table if required