diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C
index 33bec74ecc..a1cf1b82b7 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -101,51 +101,69 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
+        if (pimple.frozenFlow())
         {
-            if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
+             fluid.solve(rAUs, rAUfs);
+             fluid.correct();
+             fluid.correctContinuityError();
+
+             #include "YEqns.H"
+             #include "EEqns.H"
+             #include "pEqnComps.H"
+
+             forAll(phases, phasei)
+             {
+                 phases[phasei].divU(-pEqnComps[phasei] & p_rgh);
+             }
+        }
+        else
+        {
+            // --- Pressure-velocity PIMPLE corrector loop
+            while (pimple.loop())
             {
-                mesh.update();
-
-                if (mesh.changing())
+                if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
                 {
-                    gh = (g & mesh.C()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
+                    mesh.update();
 
-                    fluid.meshUpdate();
-
-                    if (checkMeshCourantNo)
+                    if (mesh.changing())
                     {
-                        #include "meshCourantNo.H"
+                        gh = (g & mesh.C()) - ghRef;
+                        ghf = (g & mesh.Cf()) - ghRef;
+
+                        fluid.meshUpdate();
+
+                        if (checkMeshCourantNo)
+                        {
+                            #include "meshCourantNo.H"
+                        }
                     }
                 }
-            }
 
-            fluid.solve(rAUs, rAUfs);
-            fluid.correct();
-            fluid.correctContinuityError();
+                fluid.solve(rAUs, rAUfs);
+                fluid.correct();
+                fluid.correctContinuityError();
 
-            #include "YEqns.H"
+                #include "YEqns.H"
 
-            if (faceMomentum)
-            {
-                #include "pUf/UEqns.H"
-                #include "EEqns.H"
-                #include "pUf/pEqn.H"
-            }
-            else
-            {
-                #include "pU/UEqns.H"
-                #include "EEqns.H"
-                #include "pU/pEqn.H"
-            }
+                if (faceMomentum)
+                {
+                    #include "pUf/UEqns.H"
+                    #include "EEqns.H"
+                    #include "pUf/pEqn.H"
+                }
+                else
+                {
+                    #include "pU/UEqns.H"
+                    #include "EEqns.H"
+                    #include "pU/pEqn.H"
+                }
 
-            fluid.correctKinematics();
+                fluid.correctKinematics();
 
-            if (pimple.turbCorr())
-            {
-                fluid.correctTurbulence();
+                if (pimple.turbCorr())
+                {
+                    fluid.correctTurbulence();
+                }
             }
         }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H
new file mode 100644
index 0000000000..a3cfe3c6bc
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H
@@ -0,0 +1,75 @@
+PtrList<fvScalarMatrix> pEqnComps(phases.size());
+
+{
+    PtrList<volScalarField> dmdts(fluid.dmdts());
+
+    forAll(phases, phasei)
+    {
+        phaseModel& phase = phases[phasei];
+        const volScalarField& alpha = phase;
+        volScalarField& rho = phase.thermoRef().rho();
+
+        pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
+        fvScalarMatrix& pEqnComp = pEqnComps[phasei];
+
+        // Density variation
+        if (!phase.isochoric() || !phase.pure())
+        {
+            pEqnComp +=
+                (
+                    fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
+                  - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
+                )/rho;
+        }
+
+        // Mesh dilatation correction
+        if (mesh.moving())
+        {
+            pEqnComp += fvc::div(mesh.phi())*alpha;
+        }
+
+        // Compressibility
+        if (!phase.incompressible())
+        {
+            if (pimple.transonic())
+            {
+                const surfaceScalarField phid
+                (
+                    IOobject::groupName("phid", phase.name()),
+                    fvc::interpolate(phase.thermo().psi())*phase.phi()
+                );
+
+                pEqnComp +=
+                    correction
+                    (
+                        (alpha/rho)*
+                        (
+                            phase.thermo().psi()*fvm::ddt(p_rgh)
+                          + fvm::div(phid, p_rgh)
+                          - fvm::Sp(fvc::div(phid), p_rgh)
+                        )
+                    );
+
+                pEqnComps[phasei].relax();
+            }
+            else
+            {
+                pEqnComp +=
+                    (alpha*phase.thermo().psi()/rho)
+                   *correction(fvm::ddt(p_rgh));
+            }
+        }
+
+        // Option sources
+        if (fvOptions.appliesToField(rho.name()))
+        {
+            pEqnComp -= (fvOptions(alpha, rho) & rho)/rho;
+        }
+
+        // Mass transfer
+        if (dmdts.set(phasei))
+        {
+            pEqnComp -= dmdts[phasei]/rho;
+        }
+    }
+}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pU/pEqn.H
index a123d3cc95..1c1b02df50 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pU/pEqn.H
@@ -231,77 +231,7 @@ while (pimple.correct())
     }
 
     // Compressible pressure equations
-    PtrList<fvScalarMatrix> pEqnComps(phases.size());
-    PtrList<volScalarField> dmdts(fluid.dmdts());
-    forAll(phases, phasei)
-    {
-        phaseModel& phase = phases[phasei];
-        const volScalarField& alpha = phase;
-        volScalarField& rho = phase.thermoRef().rho();
-
-        pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
-        fvScalarMatrix& pEqnComp = pEqnComps[phasei];
-
-        // Density variation
-        if (!phase.isochoric() || !phase.pure())
-        {
-            pEqnComp +=
-                (
-                    fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
-                  - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
-                )/rho;
-        }
-
-        // Mesh dilatation correction
-        if (mesh.moving())
-        {
-            pEqnComp += fvc::div(mesh.phi())*alpha;
-        }
-
-        // Compressibility
-        if (!phase.incompressible())
-        {
-            if (pimple.transonic())
-            {
-                const surfaceScalarField phid
-                (
-                    IOobject::groupName("phid", phase.name()),
-                    fvc::interpolate(phase.thermo().psi())*phase.phi()
-                );
-
-                pEqnComp +=
-                    correction
-                    (
-                        (alpha/rho)*
-                        (
-                            phase.thermo().psi()*fvm::ddt(p_rgh)
-                          + fvm::div(phid, p_rgh)
-                          - fvm::Sp(fvc::div(phid), p_rgh)
-                        )
-                    );
-
-                pEqnComps[phasei].relax();
-            }
-            else
-            {
-                pEqnComp +=
-                    (alpha*phase.thermo().psi()/rho)
-                   *correction(fvm::ddt(p_rgh));
-            }
-        }
-
-        // Option sources
-        if (fvOptions.appliesToField(rho.name()))
-        {
-            pEqnComp -= (fvOptions(alpha, rho) & rho)/rho;
-        }
-
-        // Mass transfer
-        if (dmdts.set(phasei))
-        {
-            pEqnComp -= dmdts[phasei]/rho;
-        }
-    }
+    #include "pEqnComps.H"
 
     // Cache p prior to solve for density update
     volScalarField p_rgh_0(p_rgh);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pUf/pEqn.H
index 38a1fd1cf9..74b95dfc45 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pUf/pEqn.H
@@ -222,82 +222,7 @@ while (pimple.correct())
     }
 
     // Compressible pressure equations
-    PtrList<fvScalarMatrix> pEqnComps(phases.size());
-    PtrList<volScalarField> dmdts(fluid.dmdts());
-    forAll(phases, phasei)
-    {
-        phaseModel& phase = phases[phasei];
-        const volScalarField& alpha = phase;
-        volScalarField& rho = phase.thermoRef().rho();
-
-        pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
-        fvScalarMatrix& pEqnComp = pEqnComps[phasei];
-
-        // Density variation
-        if (!phase.isochoric() || !phase.pure())
-        {
-            pEqnComp +=
-                (
-                    fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
-                  - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
-                )/rho;
-        }
-
-        // Mesh dilatation correction
-        if (mesh.moving())
-        {
-            pEqnComp += fvc::div(mesh.phi())*alpha;
-        }
-
-        // Compressibility
-        if (!phase.incompressible())
-        {
-            if (pimple.transonic())
-            {
-                const surfaceScalarField phid
-                (
-                    IOobject::groupName("phid", phase.name()),
-                    fvc::interpolate(phase.thermo().psi())*phase.phi()
-                );
-
-                pEqnComp +=
-                    correction
-                    (
-                        (alpha/rho)*
-                        (
-                            phase.thermo().psi()*fvm::ddt(p_rgh)
-                          + fvm::div(phid, p_rgh)
-                          - fvm::Sp(fvc::div(phid), p_rgh)
-                        )
-                    );
-
-                deleteDemandDrivenData
-                (
-                    pEqnComps[phasei].faceFluxCorrectionPtr()
-                );
-
-                pEqnComps[phasei].relax();
-            }
-            else
-            {
-                pEqnComp +=
-                    (alpha*phase.thermo().psi()/rho)
-                   *correction(fvm::ddt(p_rgh));
-            }
-        }
-
-        // Option sources
-        if (fvOptions.appliesToField(rho.name()))
-        {
-            pEqnComp -= (fvOptions(alpha, rho) & rho)/rho;
-        }
-
-        // Mass transfer
-        if (dmdts.set(phasei))
-        {
-            pEqnComp -= dmdts[phasei]/rho;
-        }
-    }
+    #include "pEqnComps.H"
 
     // Cache p prior to solve for density update
     volScalarField p_rgh_0(p_rgh);
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
index 9064c12792..02d1019c4e 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
@@ -60,6 +60,8 @@ PIMPLE
 
     pRefCell 0;
     pRefValue 1e5;
+
+    frozenFlow      yes;
 }
 
 relaxationFactors
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
index 9064c12792..02d1019c4e 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
@@ -60,6 +60,8 @@ PIMPLE
 
     pRefCell 0;
     pRefValue 1e5;
+
+    frozenFlow      yes;
 }
 
 relaxationFactors