diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index 6ca75c6ed4..4285a40c0c 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
@@ -77,11 +77,39 @@ void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
 }
 
 
+template<class BasePhaseModel>
+Foam::tmp<Foam::volScalarField>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::filterPressureWork
+(
+    const tmp<volScalarField>& pressureWork
+) const
+{
+    const volScalarField& alpha = *this;
+
+    scalar pressureWorkAlphaLimit =
+        this->thermo_->lookupOrDefault("pressureWorkAlphaLimit", 0.0);
+
+    if (pressureWorkAlphaLimit > 0)
+    {
+        return
+        (
+            max(alpha - pressureWorkAlphaLimit, scalar(0))
+           /max(alpha - pressureWorkAlphaLimit, pressureWorkAlphaLimit)
+        )*pressureWork;
+    }
+    else
+    {
+        return pressureWork;
+    }
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
+    const volVectorField& U = this->U();
     const surfaceScalarField& alphaPhi = this->alphaPhi();
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi();
 
@@ -93,7 +121,8 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 
     tmp<fvScalarMatrix> tEEqn
     (
-        fvm::ddt(alpha, this->rho(), he) + fvm::div(alphaRhoPhi, he)
+        fvm::ddt(alpha, this->rho(), he)
+      + fvm::div(alphaRhoPhi, he)
       - fvm::Sp(contErr, he)
 
       + fvc::ddt(alpha, this->rho(), K_) + fvc::div(alphaRhoPhi, K_)
@@ -112,13 +141,18 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
     // Add the appropriate pressure-work term
     if (he.name() == this->thermo_->phasePropertyName("e"))
     {
-        tEEqn.ref() +=
-            fvc::ddt(alpha)*this->thermo().p()
-          + fvc::div(alphaPhi, this->thermo().p());
+        tEEqn.ref() += filterPressureWork
+        (
+            fvc::div(fvc::absolute(alphaPhi, alpha, U), this->thermo().p())
+        );
     }
     else if (this->thermo_->dpdt())
     {
-        tEEqn.ref() -= alpha*this->fluid().dpdt();
+        tEEqn.ref() -= filterPressureWork
+        (
+           fvc::ddt(alpha, this->thermo().p())
+         + alpha*(this->fluid().dpdt() - fvc::ddt(this->thermo().p()))
+        );
     }
 
     return tEEqn;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
index 7a9cbd51bf..86edbcd1f8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
@@ -58,6 +58,15 @@ class AnisothermalPhaseModel
         volScalarField K_;
 
 
+    // Private member functions
+
+        //- Optionally filter the pressure work term as the phase-fraction -> 0
+        tmp<volScalarField> filterPressureWork
+        (
+            const tmp<volScalarField>& pressureWork
+        ) const;
+
+
 public:
 
     // Constructors