diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index c9dd621d3e2dc711a3daa7aa85611fd90b01faec..65e6e420d1257faaccba352a14e73289b6c238ad 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -601,12 +601,13 @@ void Foam::MULES::limiter
         {
             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
             const scalarField& phiCorrfPf = phiCorrBf[patchi];
+            const fvPatchScalarField& psiPf = psiBf[patchi];
 
             if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
             {
                 lambdaPf = 0;
             }
-            else
+            else if (psiPf.coupled())
             {
                 const labelList& pFaceCells =
                     mesh.boundary()[patchi].faceCells();
@@ -627,6 +628,32 @@ void Foam::MULES::limiter
                     }
                 }
             }
+            else
+            {
+                const labelList& pFaceCells =
+                    mesh.boundary()[patchi].faceCells();
+                const scalarField& phiBDPf = phiBDBf[patchi];
+
+                forAll(lambdaPf, pFacei)
+                {
+                    // Limit outlet faces only
+                    if (phiBDPf[pFacei] > 0)
+                    {
+                        label pfCelli = pFaceCells[pFacei];
+
+                        if (phiCorrfPf[pFacei] > 0.0)
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdap[pfCelli]);
+                        }
+                        else
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdam[pfCelli]);
+                        }
+                    }
+                }
+            }
         }
 
         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());