diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
index 355dae922..5736d3d6f 100644
--- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
@@ -21,12 +21,21 @@
 
         p.storePrevIter();
 
+        volScalarField divPhid
+        (
+            "divPhid",
+            fvc::div(phid)
+        );
+
         for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
         {
             fvScalarMatrix pEqn
             (
                 fvm::ddt(psis, p)
               + fvm::div(phid, p)
+                // Convective flux relaxation terms
+              + fvm::SuSp(-divPhid, p)
+              + divPhid*p
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
index b45e1573d..868f74040 100644
--- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
@@ -25,12 +25,21 @@
 
         p.storePrevIter();
 
+        volScalarField divPhid
+        (
+            "divPhid",
+            fvc::div(phid)
+        );
+
         for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
         {
             fvScalarMatrix pEqn
             (
                 fvm::ddt(psis, p)
               + fvm::div(phid, p)
+                // Convective flux relaxation terms
+              + fvm::SuSp(-divPhid, p)
+              + divPhid*p
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
index f161a3c4d..0718dd8d2 100644
--- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
@@ -20,12 +20,21 @@
 
         p.storePrevIter();
 
+        volScalarField divPhid
+        (
+            "divPhid",
+            fvc::div(phid)
+        );
+
         for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
         {
             fvScalarMatrix pEqn
             (
                 fvm::ddt(psis, p)
               + fvm::div(phid, p)
+                // Convective flux relaxation terms
+              + fvm::SuSp(-divPhid, p)
+              + divPhid*p
               + fvc::div(phid2)
               - fvm::laplacian(rho*rUrelA, p)
             );