diff --git a/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C b/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C
index 43bafa0342..049df0d273 100644
--- a/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C
+++ b/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C
@@ -241,7 +241,10 @@ bool Foam::functionObjects::phaseForces::execute()
             *forceFields_[virtualMassModel::typeName] +=
                 fluid_.lookupInterfacialModel
                 <blendedVirtualMassModel>(interface).K()
-               *(otherPhase.DUDt() - phase_.DUDt());
+               *(
+                    (otherPhase.DUDt() & otherPhase.U())
+                  - (phase_.DUDt() & phase_.U())
+                );
         }
 
         if (fluid_.foundInterfacialModel<blendedLiftModel>(interface))
diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
index ac12ba5f77..11f48b6204 100644
--- a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
+++ b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
@@ -68,6 +68,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         rAs.setSize(phases.size());
     }
 
+    PtrList<volVectorField> HVms(movingPhases.size());
     PtrList<PtrList<volScalarField>> invADs;
     PtrList<PtrList<surfaceScalarField>> invADfs;
     {
@@ -109,7 +110,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             }
         }
 
-        fluid.invADs(As, invADs, invADfs);
+        fluid.invADs(As, HVms, invADs, invADfs);
     }
 
     volScalarField rho("rho", fluid.rho());
@@ -168,7 +169,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     // --- Optional momentum predictor
     if (predictMomentum)
     {
-        InfoInFunction << endl;
         PtrList<volVectorField> HbyADs;
         {
             PtrList<volVectorField> Hs(movingPhases.size());
@@ -189,6 +189,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                     )
                    *phase.U()().oldTime()
                 );
+
+                if (HVms.set(movingPhasei))
+                {
+                    Hs[movingPhasei] += HVms[movingPhasei];
+                }
             }
 
             HbyADs = invADs & Hs;
@@ -253,6 +258,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                    *phase.U()().oldTime()
                 );
 
+                if (HVms.set(movingPhasei))
+                {
+                    Hs[movingPhasei] += HVms[movingPhasei];
+                }
+
                 phiHs.set
                 (
                     movingPhasei,
diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C
index 92960da404..05e5b07f35 100644
--- a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C
+++ b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C
@@ -68,9 +68,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
         rAs.setSize(phases.size());
     }
 
+    PtrList<surfaceScalarField> HVmfs(movingPhases.size());
     PtrList<PtrList<surfaceScalarField>> invADfs;
     {
-        PtrList<surfaceScalarField> Vmfs(fluid.Vmfs());
         PtrList<surfaceScalarField> Afs(movingPhases.size());
 
         forAll(fluid.movingPhases(), movingPhasei)
@@ -110,14 +110,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                     fvc::interpolate(A)
                 )
             );
-
-            if (Vmfs.set(phase.index()))
-            {
-                Afs[movingPhasei] += Vmfs[phase.index()];
-            }
         }
 
-        invADfs = fluid.invADfs(Afs);
+        invADfs = fluid.invADfs(Afs, HVmfs);
     }
 
     volScalarField rho("rho", fluid.rho());
@@ -205,6 +200,11 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                      + fvc::flux(UEqns[phase.index()].H())
                     )
                 );
+
+                if (HVmfs.set(movingPhasei))
+                {
+                    phiHs[movingPhasei] += HVmfs[movingPhasei];
+                }
             }
 
             phiHbyADs = invADfs & phiHs;
diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index b95b33f30e..b79499374b 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -224,87 +224,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
         }
     }
 
-    // Initialise Vms table if required
-    if (!Vms_.size())
-    {
-        forAllConstIter
-        (
-            virtualMassModelTable,
-            virtualMassModels_,
-            virtualMassModelIter
-        )
-        {
-            const phaseInterface& interface =
-                virtualMassModelIter()->interface();
-
-            Vms_.insert
-            (
-                interface,
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        IOobject::groupName("Vm", interface.name()),
-                        this->mesh().time().name(),
-                        this->mesh()
-                    ),
-                    this->mesh(),
-                    dimensionedScalar(virtualMassModel::dimK, 0)
-                )
-            );
-        }
-    }
-
-    // Update the virtual mass coefficients
-    forAllConstIter
-    (
-        virtualMassModelTable,
-        virtualMassModels_,
-        virtualMassModelIter
-    )
-    {
-        *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K();
-    }
-
-    // Add the virtual mass force
-    forAllConstIter(VmTable, Vms_, VmIter)
-    {
-        const volScalarField& Vm(*VmIter());
-        const phaseInterface interface(*this, VmIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            if (!phase.stationary())
-            {
-                fvVectorMatrix& eqn = *eqns[phase.name()];
-
-                const volVectorField& U = eqn.psi();
-
-                const surfaceScalarField& phi = phase.phiRef();
-                const tmp<surfaceScalarField> taphi(fvc::absolute(phi, U));
-                const surfaceScalarField& aphi(taphi());
-
-                const volScalarField VmPhase
-                (
-                    (otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-                   *Vm
-                );
-
-                eqn -=
-                    VmPhase
-                   *(
-                        fvm::ddt(U)
-                      + fvm::div(aphi, U) - fvm::Sp(fvc::div(aphi), U)
-                      - otherPhase.DUDt()
-                    )
-                  + this->MRF_.DDt(VmPhase, U - otherPhase.U());
-            }
-        }
-    }
-
     return eqnsPtr;
 }
 
@@ -373,138 +292,10 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf()
         *Kdfs_[dragModelIter.key()] = dragModelIter()->Kf();
     }
 
-    // Initialise Vms table if required
-    if (!Vms_.size())
-    {
-        forAllConstIter
-        (
-            virtualMassModelTable,
-            virtualMassModels_,
-            virtualMassModelIter
-        )
-        {
-            const phaseInterface& interface =
-                virtualMassModelIter()->interface();
-
-            Vms_.insert
-            (
-                interface,
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        IOobject::groupName("Vm", interface.name()),
-                        this->mesh().time().name(),
-                        this->mesh()
-                    ),
-                    this->mesh(),
-                    dimensionedScalar(virtualMassModel::dimK, 0)
-                )
-            );
-        }
-    }
-
-    // Update the virtual mass coefficients
-    forAllConstIter
-    (
-        virtualMassModelTable,
-        virtualMassModels_,
-        virtualMassModelIter
-    )
-    {
-        *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K();
-    }
-
-    // Create U & grad(U) fields
-    PtrList<fvVectorMatrix> UgradUs(this->phaseModels_.size());
-    forAll(this->phaseModels_, phasei)
-    {
-        const phaseModel& phase = this->phaseModels_[phasei];
-
-        if (!phase.stationary())
-        {
-            const volVectorField& U = phase.URef();
-            const surfaceScalarField& phi = phase.phiRef();
-
-            const tmp<surfaceScalarField> taphi(fvc::absolute(phi, U));
-            const surfaceScalarField& aphi(taphi());
-
-            UgradUs.set
-            (
-                phasei,
-                new fvVectorMatrix
-                (
-                    fvm::div(aphi, U) - fvm::Sp(fvc::div(aphi), U)
-                  + this->MRF().DDt(U)
-                )
-            );
-        }
-    }
-
-    // Add the virtual mass force
-    forAllConstIter(VmTable, Vms_, VmIter)
-    {
-        const volScalarField& Vm(*VmIter());
-        const phaseInterface interface(*this, VmIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            if (!phase.stationary())
-            {
-                const volScalarField VmPhase
-                (
-                    (otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-                   *Vm
-                );
-
-                *eqns[phase.name()] -=
-                    VmPhase
-                   *(
-                        UgradUs[phase.index()]
-                      - (UgradUs[otherPhase.index()] & otherPhase.U())
-                    );
-            }
-        }
-    }
-
     return eqnsPtr;
 }
 
 
-template<class BasePhaseSystem>
-Foam::PtrList<Foam::surfaceScalarField>
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vmfs() const
-{
-    PtrList<surfaceScalarField> Vmfs(this->phaseModels_.size());
-
-    // Add the implicit part of the virtual mass force
-    forAllConstIter(VmTable, Vms_, VmIter)
-    {
-        const volScalarField& Vm(*VmIter());
-        const phaseInterface interface(*this, VmIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            const volScalarField VmPhase
-            (
-                (otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-               *Vm
-            );
-
-            addField(phase, "Vmf", byDt(fvc::interpolate(VmPhase)), Vmfs);
-        }
-    }
-
-    return Vmfs;
-}
-
-
 template<class BasePhaseSystem>
 Foam::PtrList<Foam::surfaceScalarField>
 Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
@@ -633,44 +424,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Ffs() const
 {
     PtrList<surfaceScalarField> Ffs(this->phaseModels_.size());
 
-    // Add the explicit part of the virtual mass force
-    forAllConstIter(VmTable, Vms_, VmIter)
-    {
-        const volScalarField& Vm(*VmIter());
-        const phaseInterface interface(*this, VmIter.key());
-
-        forAllConstIter(phaseInterface, interface, iter)
-        {
-            const phaseModel& phase = iter();
-            const phaseModel& otherPhase = iter.otherPhase();
-
-            const volScalarField VmPhase
-            (
-                (otherPhase/max(otherPhase, otherPhase.residualAlpha()))
-               *Vm
-            );
-
-            addField
-            (
-                phase,
-                "Ff",
-               -fvc::interpolate(VmPhase)
-               *(
-                   byDt
-                   (
-                       fvc::absolute
-                       (
-                           this->MRF().absolute(iter().phi()().oldTime()),
-                           iter().U()
-                       )
-                   )
-                 + otherPhase.DUDtf()
-                ),
-                Ffs
-            );
-        }
-    }
-
     // Add the lift force
     forAllConstIter
     (
@@ -871,6 +624,7 @@ template<class BasePhaseSystem>
 void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADs
 (
     const PtrList<volScalarField>& As,
+    PtrList<volVectorField>& HVms,
     PtrList<PtrList<volScalarField>>& invADs,
     PtrList<PtrList<surfaceScalarField>>& invADfs
 ) const
@@ -962,6 +716,94 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADs
         }
     }
 
+    // Cache the phase acceleration As and Hs
+    PtrList<volScalarField> ADUDts(movingPhases.size());
+    PtrList<volVectorField> HDUDts(movingPhases.size());
+
+    forAllConstIter
+    (
+        virtualMassModelTable,
+        virtualMassModels_,
+        VmIter
+    )
+    {
+        const phaseInterface& interface = VmIter()->interface();
+
+        forAllConstIter(phaseInterface, interface, iter)
+        {
+            const phaseModel& phase = iter();
+            const label i = movingPhases[phase.index()];
+
+            if (i != -1 && !ADUDts.set(i))
+            {
+                const fvVectorMatrix DUDt(phase.DUDt());
+                ADUDts.set(i, DUDt.A());
+                HDUDts.set(i, DUDt.H());
+            }
+        }
+    }
+
+    // Add the virtual mass contributions
+    forAllConstIter
+    (
+        virtualMassModelTable,
+        virtualMassModels_,
+        VmIter
+    )
+    {
+        const phaseInterface& interface = VmIter()->interface();
+        const volScalarField Vm(VmIter()->K());
+
+        forAllConstIter(phaseInterface, interface, iter)
+        {
+            const phaseModel& phase = iter();
+            const phaseModel& otherPhase = iter.otherPhase();
+
+            const label i = movingPhases[phase.index()];
+
+            if (i != -1)
+            {
+                const volScalarField VmPhase
+                (
+                    (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*Vm
+                );
+
+                {
+                    const volScalarField AVm(VmPhase*ADUDts[i]);
+
+                    invADs[i][i] += AVm;
+                    invADfs[i][i] += fvc::interpolate(AVm);
+
+                    addField
+                    (
+                        i,
+                        "HVm",
+                        VmPhase*HDUDts[i],
+                        HVms
+                    );
+                }
+
+                const label j = movingPhases[otherPhase.index()];
+
+                if (j != -1)
+                {
+                    const volScalarField AVm(VmPhase*ADUDts[j]);
+
+                    invADs[i][j] -= AVm;
+                    invADfs[i][j] -= fvc::interpolate(AVm);
+
+                    addField
+                    (
+                        i,
+                        "HVm",
+                        -VmPhase*HDUDts[j],
+                        HVms
+                    );
+                }
+            }
+        }
+    }
+
     MomentumTransferPhaseSystem<BasePhaseSystem>::invADs(invADs);
     MomentumTransferPhaseSystem<BasePhaseSystem>::invADs(invADfs);
 }
@@ -971,17 +813,18 @@ template<class BasePhaseSystem>
 Foam::PtrList<Foam::PtrList<Foam::surfaceScalarField>>
 Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADfs
 (
-    const PtrList<surfaceScalarField>& As
+    const PtrList<surfaceScalarField>& Afs,
+    PtrList<surfaceScalarField>& HVmfs
 ) const
 {
-    const label n = As.size();
+    const label n = Afs.size();
 
     PtrList<PtrList<surfaceScalarField>> invADfs(n);
 
     forAll(invADfs, i)
     {
         invADfs.set(i, new PtrList<surfaceScalarField>(n));
-        invADfs[i].set(i, As[i].clone());
+        invADfs[i].set(i, Afs[i].clone());
     }
 
     labelList movingPhases(this->phases().size(), -1);
@@ -1040,13 +883,119 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::invADfs
                     (
                         "0",
                         this->mesh(),
-                        dimensionedScalar(As[0].dimensions(), 0)
+                        dimensionedScalar(Afs[0].dimensions(), 0)
                     )
                 );
             }
         }
     }
 
+    // Cache the phase acceleration Afs and Hs
+    PtrList<volScalarField> AUgradUs(movingPhases.size());
+    PtrList<volVectorField> HUgradUs(movingPhases.size());
+
+    forAllConstIter
+    (
+        virtualMassModelTable,
+        virtualMassModels_,
+        VmIter
+    )
+    {
+        const phaseInterface& interface = VmIter()->interface();
+
+        forAllConstIter(phaseInterface, interface, iter)
+        {
+            const phaseModel& phase = iter();
+            const label i = movingPhases[phase.index()];
+
+            if (i != -1 && !AUgradUs.set(i))
+            {
+                const fvVectorMatrix UgradU(phase.UgradU());
+                AUgradUs.set(i, UgradU.A());
+                HUgradUs.set(i, UgradU.H());
+            }
+        }
+    }
+
+    // Add the virtual mass contributions
+    forAllConstIter
+    (
+        virtualMassModelTable,
+        virtualMassModels_,
+        VmIter
+    )
+    {
+        const phaseInterface& interface = VmIter()->interface();
+        const volScalarField Vm(VmIter()->K());
+
+        forAllConstIter(phaseInterface, interface, iter)
+        {
+            const phaseModel& phase = iter();
+            const phaseModel& otherPhase = iter.otherPhase();
+
+            const label i = movingPhases[phase.index()];
+
+            if (i != -1)
+            {
+                const volScalarField VmPhase
+                (
+                    (otherPhase/max(otherPhase, otherPhase.residualAlpha()))
+                   *Vm
+                );
+
+                const surfaceScalarField VmPhasef(fvc::interpolate(VmPhase));
+
+                invADfs[i][i] +=
+                    byDt(VmPhasef) + fvc::interpolate(VmPhase*AUgradUs[i]);
+
+                addField
+                (
+                    i,
+                    "HVmf",
+                    VmPhasef
+                   *byDt
+                    (
+                        fvc::absolute
+                        (
+                            this->MRF().absolute(phase.phi()().oldTime()),
+                            phase.U()
+                        )
+                    )
+                  + fvc::flux(VmPhase*HUgradUs[i]),
+                    HVmfs
+                );
+
+                const label j = movingPhases[otherPhase.index()];
+
+                if (j != -1)
+                {
+                    invADfs[i][j] -=
+                        byDt(VmPhasef) + fvc::interpolate(VmPhase*AUgradUs[j]);
+
+                    addField
+                    (
+                        i,
+                        "HVmf",
+                       -VmPhasef
+                       *byDt
+                        (
+                            fvc::absolute
+                            (
+                                this->MRF().absolute
+                                (
+                                    otherPhase.phi()().oldTime()
+                                ),
+                                otherPhase.U()
+                            )
+                        )
+                      - fvc::flux(VmPhase*HUgradUs[j]),
+                        HVmfs
+                    );
+                }
+            }
+        }
+    }
+
     invADs(invADfs);
 
     return invADfs;
@@ -1194,12 +1143,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::ddtCorrs() const
             const phaseModel& phase = this->movingPhases()[movingPhasei];
             const label phasei = phase.index();
 
-            VmDdtCoeffs.set
-            (
-                phasei,
-                fvm::ddt(phase.U()())().A()
-            );
-
             VmDdtCorrs.set
             (
                 phasei,
@@ -1212,10 +1155,15 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::ddtCorrs() const
             );
         }
 
-        forAllConstIter(VmTable, Vms_, VmIter)
+        forAllConstIter
+        (
+            virtualMassModelTable,
+            virtualMassModels_,
+            VmIter
+        )
         {
-            const volScalarField& Vm(*VmIter());
-            const phaseInterface interface(*this, VmIter.key());
+            const phaseInterface& interface = VmIter()->interface();
+            const volScalarField Vm(VmIter()->K());
 
             forAllConstIter(phaseInterface, interface, iter)
             {
@@ -1235,19 +1183,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::ddtCorrs() const
                     phase,
                     "ddtCorr",
                     fvc::interpolate(VmPhase)
-                   *(
-                        VmDdtCorrs[phasei]
-                      + (
-                            fvc::interpolate(VmDdtCoeffs[otherPhasei])
-                           *(
-                               otherPhase.Uf().valid()
-                             ? (this->mesh_.Sf() & otherPhase.Uf()())()
-                             : otherPhase.phi()()
-                            )
-                          - fvc::flux(VmDdtCoeffs[otherPhasei]*otherPhase.U())
-                        )
-                      - VmDdtCorrs[otherPhase.index()]
-                    ),
+                   *(VmDdtCorrs[phasei] - VmDdtCorrs[otherPhasei]),
                     ddtCorrs
                 );
             }
diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index 8c35852810..9333d011f3 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
@@ -79,13 +79,6 @@ class MomentumTransferPhaseSystem
             phaseInterfaceKey::hash
         > KdfTable;
 
-        typedef HashPtrTable
-        <
-            volScalarField,
-            phaseInterfaceKey,
-            phaseInterfaceKey::hash
-        > VmTable;
-
         typedef HashTable
         <
             autoPtr<blendedDragModel>,
@@ -130,9 +123,6 @@ class MomentumTransferPhaseSystem
         //- Face drag coefficients
         KdfTable Kdfs_;
 
-        //- Virtual mass coefficients
-        VmTable Vms_;
-
 
         // Sub Models
 
@@ -202,10 +192,6 @@ public:
         //- As momentumTransfer, but for the face-based algorithm
         virtual autoPtr<phaseSystem::momentumTransferTable> momentumTransferf();
 
-        //- Return implicit force coefficients on the faces, for the face-based
-        //  algorithm.
-        virtual PtrList<surfaceScalarField> Vmfs() const;
-
         //- Return the explicit force fluxes for the cell-based algorithm, that
         //  do not depend on phase mass/volume fluxes, and can therefore be
         //  evaluated outside the corrector loop. This includes things like
@@ -219,6 +205,7 @@ public:
         virtual void invADs
         (
             const PtrList<volScalarField>& As,
+            PtrList<volVectorField>& HVms,
             PtrList<PtrList<volScalarField>>& invADs,
             PtrList<PtrList<surfaceScalarField>>& invADfs
         ) const;
@@ -226,7 +213,8 @@ public:
         //- Return the inverse of (Afs + Dfs)
         virtual PtrList<PtrList<surfaceScalarField>> invADfs
         (
-            const PtrList<surfaceScalarField>& Afs
+            const PtrList<surfaceScalarField>& Afs,
+            PtrList<surfaceScalarField>& HVmfs
         ) const;
 
         //- Returns true if the phase pressure is treated implicitly
diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
index a75563de45..049c339224 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
+++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
@@ -162,8 +162,6 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
         dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), 0)
     ),
     Uf_(nullptr),
-    DUDt_(nullptr),
-    DUDtf_(nullptr),
     divU_(nullptr),
     momentumTransport_
     (
@@ -252,18 +250,6 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctKinematics()
 {
     BasePhaseModel::correctKinematics();
 
-    if (DUDt_.valid())
-    {
-        DUDt_.clear();
-        DUDt();
-    }
-
-    if (DUDtf_.valid())
-    {
-        DUDtf_.clear();
-        DUDtf();
-    }
-
     if (K_.valid())
     {
         K_.ref() = 0.5*magSqr(this->U());
@@ -509,40 +495,23 @@ Foam::MovingPhaseModel<BasePhaseModel>::alphaRhoPhiRef() const
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::volVectorField>
-Foam::MovingPhaseModel<BasePhaseModel>::DUDt() const
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::MovingPhaseModel<BasePhaseModel>::UgradU() const
 {
-    if (!DUDt_.valid())
-    {
-        const tmp<surfaceScalarField> taphi(fvc::absolute(phi_, U_));
-        const surfaceScalarField& aphi(taphi());
-        DUDt_ =
-            new volVectorField
-            (
-                IOobject::groupName("DUDt", this->name()),
-                fvc::ddt(U_) + fvc::div(aphi, U_) - fvc::div(aphi)*U_
-            );
-    }
+    const tmp<surfaceScalarField> taphi(fvc::absolute(phi_, U_));
+    const surfaceScalarField& aphi(taphi());
 
-    return tmp<volVectorField>(DUDt_());
+    return
+        fvm::div(aphi, U_) - fvm::Sp(fvc::div(aphi), U_)
+      + this->fluid().MRF().DDt(U_);
 }
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::surfaceScalarField>
-Foam::MovingPhaseModel<BasePhaseModel>::DUDtf() const
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::MovingPhaseModel<BasePhaseModel>::DUDt() const
 {
-    if (!DUDtf_.valid())
-    {
-        DUDtf_ =
-            new surfaceScalarField
-            (
-                IOobject::groupName("DUDtf", this->name()),
-                byDt(phi_ - phi_.oldTime())
-            );
-    }
-
-    return tmp<surfaceScalarField>(DUDtf_());
+    return fvm::ddt(U_) + UgradU();
 }
 
 
diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
index 29431aceb6..b1e2463087 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
+++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
@@ -112,12 +112,6 @@ protected:
         //- Face velocity field
         autoPtr<surfaceVectorField> Uf_;
 
-        //- Lagrangian acceleration field (needed for virtual-mass)
-        mutable tmp<volVectorField> DUDt_;
-
-        //- Lagrangian acceleration field on the faces (needed for virtual-mass)
-        mutable tmp<surfaceScalarField> DUDtf_;
-
         //- Dilatation rate
         autoPtr<volScalarField> divU_;
 
@@ -253,11 +247,11 @@ public:
             //- Access the mass flux of the phase
             virtual const surfaceScalarField& alphaRhoPhiRef() const;
 
-            //- Return the substantive acceleration
-            virtual tmp<volVectorField> DUDt() const;
+            //- Return the velocity transport matrix
+            virtual tmp<fvVectorMatrix> UgradU() const;
 
-            //- Return the substantive acceleration on the faces
-            virtual tmp<surfaceScalarField> DUDtf() const;
+            //- Return the substantive acceleration matrix
+            virtual tmp<fvVectorMatrix> DUDt() const;
 
             //- Return the continuity error
             virtual tmp<volScalarField> continuityError() const;
diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
index f5f21a7759..fddd2054aa 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
+++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
@@ -267,28 +267,26 @@ Foam::StationaryPhaseModel<BasePhaseModel>::alphaRhoPhiRef() const
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::volVectorField>
-Foam::StationaryPhaseModel<BasePhaseModel>::DUDt() const
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::StationaryPhaseModel<BasePhaseModel>::UgradU() const
 {
-    return volVectorField::New
-    (
-        IOobject::groupName("DUDt", this->name()),
-        this->mesh(),
-        dimensionedVector(dimVelocity/dimTime, Zero)
-    );
+    FatalErrorInFunction
+        << "Cannot calculate DUDt of a stationary phase"
+        << exit(FatalError);
+
+    return tmp<fvVectorMatrix>(nullptr);
 }
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::surfaceScalarField>
-Foam::StationaryPhaseModel<BasePhaseModel>::DUDtf() const
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::StationaryPhaseModel<BasePhaseModel>::DUDt() const
 {
-    return surfaceScalarField::New
-    (
-        IOobject::groupName("DUDtf", this->name()),
-        this->mesh(),
-        dimensionedScalar(dimVolume/sqr(dimTime), 0)
-    );
+    FatalErrorInFunction
+        << "Cannot calculate DUDt of a stationary phase"
+        << exit(FatalError);
+
+    return tmp<fvVectorMatrix>(nullptr);
 }
 
 
diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
index b23aa2ce00..3d4b7ded34 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
+++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H
@@ -136,11 +136,11 @@ public:
             //- Access the mass flux of the phase
             virtual const surfaceScalarField& alphaRhoPhiRef() const;
 
-            //- Return the substantive acceleration
-            virtual tmp<volVectorField> DUDt() const;
+            //- Return the velocity transport matrix
+            virtual tmp<fvVectorMatrix> UgradU() const;
 
-            //- Return the substantive acceleration on the faces
-            virtual tmp<surfaceScalarField> DUDtf() const;
+            //- Return the substantive acceleration matrix
+            virtual tmp<fvVectorMatrix> DUDt() const;
 
             //- Return the continuity error
             virtual tmp<volScalarField> continuityError() const;
diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H
index f42991ab34..4a355e9f95 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -359,11 +359,11 @@ public:
             //- Access the mass flux of the phase
             virtual const surfaceScalarField& alphaRhoPhiRef() const = 0;
 
-            //- Return the substantive acceleration
-            virtual tmp<volVectorField> DUDt() const = 0;
+            //- Return the velocity transport matrix
+            virtual tmp<fvVectorMatrix> UgradU() const = 0;
 
-            //- Return the substantive acceleration on the faces
-            virtual tmp<surfaceScalarField> DUDtf() const = 0;
+            //- Return the substantive acceleration matrix
+            virtual tmp<fvVectorMatrix> DUDt() const = 0;
 
             //- Return the continuity error
             virtual tmp<volScalarField> continuityError() const = 0;
diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H
index 4a87aaa0a4..c26213db0b 100644
--- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H
@@ -521,10 +521,6 @@ public:
             //  algorithm
             virtual autoPtr<momentumTransferTable> momentumTransferf() = 0;
 
-            //- Return the implicit force coefficients for the face-based
-            //  algorithm
-            virtual PtrList<surfaceScalarField> Vmfs() const = 0;
-
             //- Return the force fluxes for the cell-based algorithm
             virtual PtrList<surfaceScalarField> Fs() const = 0;
 
@@ -535,6 +531,7 @@ public:
             virtual void invADs
             (
                 const PtrList<volScalarField>& As,
+                PtrList<volVectorField>& HVms,
                 PtrList<PtrList<volScalarField>>& invADs,
                 PtrList<PtrList<surfaceScalarField>>& invADfs
             ) const = 0;
@@ -542,7 +539,8 @@ public:
             //- Return the inverse of (Afs + Dfs)
             virtual PtrList<PtrList<surfaceScalarField>> invADfs
             (
-                const PtrList<surfaceScalarField>& Afs
+                const PtrList<surfaceScalarField>& Afs,
+                PtrList<surfaceScalarField>& HVmfs
             ) const = 0;
 
             //- Returns true if the phase pressure is treated implicitly