fixedFluxPressure BC: the snGrad is now pushed into the BC from pEqn.H rather than being evaluated in the BC

2025-12-28 03:37:59 +00:00 · 2013-09-11 00:10:00 +01:00
parent cda1304911
commit 08baa6eda6
77 changed files with 531 additions and 778 deletions
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@ -12,16 +12,16 @@
    surfaceScalarField rhof("rhof", fvc::interpolate(rho));

    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();

    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + Dp*fvc::ddtCorr(U, phivAbs);
+         + rhorAUf*fvc::ddtCorr(U, phivAbs);
    fvc::makeRelative(phiv, U);

-    surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
+    surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));

    phiv -= phiGradp/rhof;

@ -35,7 +35,7 @@
          + psi*correction(fvm::ddt(p))
          + fvc::div(phiv, rho)
          + fvc::div(phiGradp)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
        );

        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@ -12,15 +12,15 @@
    surfaceScalarField rhof("rhof", fvc::interpolate(rho));

    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();

    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + Dp*fvc::ddtCorr(U, phiv);
+         + rhorAUf*fvc::ddtCorr(U, phiv);

-    surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
+    surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));

    phiv -= phiGradp/rhof;

@ -32,7 +32,7 @@
          - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
          + fvc::div(phiv, rho)
          + fvc::div(phiGradp)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
        );

        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@ -46,6 +46,7 @@ Description
 #include "twoPhaseMixtureThermo.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -61,7 +62,7 @@ int main(int argc, char *argv[])
    #include "readControls.H"
    #include "initContinuityErrs.H"
    #include "createFields.H"
-    #include "createPcorrTypes.H"
+    #include "createPrghCorrTypes.H"
    #include "CourantNo.H"
    #include "setInitialDeltaT.H"

--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@ -44,6 +44,7 @@ Description
 #include "twoPhaseMixtureThermo.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -11,18 +11,27 @@
        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
-    phi = phiHbyA;

    surfaceScalarField phig
    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    tmp<fvScalarMatrix> p_rghEqnComp1;
    tmp<fvScalarMatrix> p_rghEqnComp2;

@ -70,7 +79,7 @@
        fvScalarMatrix p_rghEqnIncomp
        (
            fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p_rgh)
+          - fvm::laplacian(rAUf, p_rgh)
        );

        solve
@ -97,7 +106,7 @@
            phi = phiHbyA + p_rghEqnIncomp.flux();

            U = HbyA
-              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/Dp);
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
            U.correctBoundaryConditions();
        }
    }
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@ -46,6 +46,7 @@ Description
 #include "fvcSmooth.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@ -45,6 +45,7 @@ Description
 #include "IOMRFZoneList.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -13,23 +13,32 @@
    );
    adjustPhi(phiHbyA, U, p_rgh);
    mrfZones.makeRelative(phiHbyA);
-    phi = phiHbyA;

    surfaceScalarField phig
    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -40,7 +49,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@ -41,6 +41,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -56,7 +57,7 @@ int main(int argc, char *argv[])
    #include "createFields.H"
    #include "createUf.H"
    #include "readTimeControls.H"
-    #include "createPcorrTypes.H"
+    #include "createPrghCorrTypes.H"
    #include "correctPhi.H"
    #include "CourantNo.H"
    #include "setInitialDeltaT.H"
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -19,23 +19,31 @@
        fvc::makeAbsolute(phiHbyA, U);
    }

-    surfaceScalarField phiAbs("phiAbs", phiHbyA);
-
    surfaceScalarField phig
    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -46,7 +54,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@ -45,6 +45,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -13,23 +13,32 @@
    );

    adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;

    surfaceScalarField phig
    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -40,7 +49,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@ -47,6 +47,7 @@ Description
 #include "IOporosityModelList.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@ -50,6 +50,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -64,9 +65,9 @@ int main(int argc, char *argv[])
    pimpleControl pimple(mesh);

    #include "createFields.H"
-    #include "../interFoam/interDyMFoam/createUf.H"
+    #include "createUf.H"
    #include "readTimeControls.H"
-    #include "../interFoam/interDyMFoam/createPcorrTypes.H"
+    #include "createPcorrTypes.H"
    #include "../interFoam/interDyMFoam/correctPhi.H"
    #include "CourantNo.H"
    #include "setInitialDeltaT.H"
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -19,18 +19,26 @@
        fvc::makeAbsolute(phiHbyA, U);
    }

-    surfaceScalarField phiAbs("phiAbs", phiHbyA);
-
    surfaceScalarField phig
    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
    const volScalarField& vDotcP = vDotP[0]();
    const volScalarField& vDotvP = vDotP[1]();
@ -39,7 +47,7 @@
    {
        fvScalarMatrix p_rghEqn
        (
-            fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
          - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
        );

@ -51,7 +59,7 @@
        {
            phi = phiHbyA + p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@ -48,6 +48,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -12,18 +12,27 @@
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;

    surfaceScalarField phig
    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
    const volScalarField& vDotcP = vDotP[0]();
    const volScalarField& vDotvP = vDotP[1]();
@ -32,7 +41,7 @@
    {
        fvScalarMatrix p_rghEqn
        (
-            fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
          - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
        );

@ -44,7 +53,7 @@
        {
            phi = phiHbyA + p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@ -35,12 +35,12 @@ Description
 #include "phaseModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
-#include "pimpleControl.H"
-
 #include "singlePhaseTransportModel.H"
 #include "LESModel.H"
-
+#include "pimpleControl.H"
 #include "IOMRFZoneList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@ -97,6 +97,8 @@
          + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
        );
        mrfZones.makeRelative(phiHbyAs[phasei]);
+        mrfZones.makeRelative(phase.phi().oldTime());
+        mrfZones.makeRelative(phase.phi());

        phiHbyAs[phasei] +=
            rAlphaAUfs[phasei]
@ -157,54 +159,58 @@
        phasei++;
    }

-    // Reset phi BCs
-    phi.boundaryField() = 0;
-
-    phasei = 0;
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
-    {
-        phaseModel& phase = iter();
-
-        phase.phi().boundaryField() ==
-            (mesh.Sf().boundaryField() & phase.U().boundaryField());
-
-        mrfZones.makeRelative(phase.phi().oldTime());
-        mrfZones.makeRelative(phase.phi());
-
-        // Update phi BCs before pEqn
-        phi.boundaryField() +=
-            alphafs[phasei].boundaryField()*phase.phi().boundaryField();
-
-        phasei++;
-    }
-
-    surfaceScalarField Dp
+    surfaceScalarField rAUf
    (
        IOobject
        (
-            "Dp",
+            "rAUf",
            runTime.timeName(),
            mesh
        ),
        mesh,
-        dimensionedScalar("Dp", dimensionSet(-1, 3, 1, 0, 0), 0)
+        dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
    );

    phasei = 0;
    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
    {
        phaseModel& phase = iter();
-        Dp += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
+        rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();

        phasei++;
    }

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    {
+        surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
+        phib = 0;
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        {
+            phaseModel& phase = iter();
+
+            phib +=
+                alphafs[phasei].boundaryField()
+               *(mesh.Sf().boundaryField() & phase.U().boundaryField());
+
+            phasei++;
+        }
+
+        setSnGrad<fixedFluxPressureFvPatchScalarField>
+        (
+            p.boundaryField(),
+            (
+                phiHbyA.boundaryField() - mrfZones.relative(phib)
+            )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+        );
+    }
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqnIncomp
        (
            fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rAUf, p)
        );

        pEqnIncomp.setReference(pRefCell, pRefValue);
@ -221,7 +227,7 @@

        if (pimple.finalNonOrthogonalIter())
        {
-            surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
+            surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);

            phasei = 0;
            phi = dimensionedScalar("phi", phi.dimensions(), 0);
@ -248,7 +254,7 @@
            // );

            p.relax();
-            mSfGradp = pEqnIncomp.flux()/Dp;
+            mSfGradp = pEqnIncomp.flux()/rAUf;

            U = dimensionedVector("U", dimVelocity, vector::zero);

--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@ -37,6 +37,7 @@ Description
 #include "turbulenceModel.H"
 #include "IOMRFZoneList.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -13,20 +13,29 @@
    );
    adjustPhi(phiHbyA, U, p_rgh);
    mrfZones.makeRelative(phiHbyA);
-    phi = phiHbyA;

    surfaceScalarField phig
    (
-        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -37,7 +46,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
        }
    }
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@ -36,6 +36,7 @@ Description
 #include "multiphaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -12,23 +12,32 @@
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;

    surfaceScalarField phig
    (
        (
            mixture.surfaceTensionForce()
          - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -39,7 +48,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
        }
    }
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@ -20,8 +20,15 @@ adjustPhi(phiHbyA, U, p_gh);

 fvOptions.makeRelative(phiHbyA);

-// Update the phi BCs from U before p BCs are updated
-phi.boundaryField() = mesh.Sf().boundaryField() & U.boundaryField();
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+    p_gh.boundaryField(),
+    (
+        phiHbyA.boundaryField()
+      - (mesh.Sf().boundaryField() & U.boundaryField())
+    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+);

 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
@ -39,6 +39,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -10,23 +10,33 @@
        "phiHbyA",
        (
           (fvc::interpolate(rho*HbyA) & mesh.Sf())
-         + Dp*fvc::ddtCorr(rho, U, phi)
+         + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );
-    phi = phiHbyA;

    surfaceScalarField phig
    (
-        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*rhorAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -37,7 +47,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
            U.correctBoundaryConditions();
        }
    }
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@ -38,6 +38,7 @@ Description
 #include "plasticViscosity.H"
 #include "yieldStress.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@ -1,6 +1,6 @@
 {
    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -12,20 +12,29 @@
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;

    surfaceScalarField phig
    (
-        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
    );

    phiHbyA += phig;

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -36,7 +45,7 @@
        {
            phi = phiHbyA - p_rghEqn.flux();

-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
        }
    }
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@ -37,6 +37,7 @@ Description
 #include "incompressibleTwoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@ -80,22 +80,12 @@

    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);

-    // Update phi BCs before pEqn
-    phi.boundaryField() =
-        mrfZones.relative
-        (
-            alpha1f.boundaryField()
-           *(mesh.Sf().boundaryField() & U1.boundaryField())
-          + alpha2f.boundaryField()
-           *(mesh.Sf().boundaryField() & U2.boundaryField())
-        );
-
    HbyA1 += (1.0/rho1)*rAU1*dragCoeff*U2;
    HbyA2 += (1.0/rho2)*rAU2*dragCoeff*U1;

-    surfaceScalarField Dp
+    surfaceScalarField rAUf
    (
-        "Dp",
+        "rAUf",
        mag
        (
            alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
@ -103,6 +93,22 @@
        )
    );

+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - mrfZones.relative
+            (
+                alpha1f.boundaryField()
+               *(mesh.Sf().boundaryField() & U1.boundaryField())
+              + alpha2f.boundaryField()
+               *(mesh.Sf().boundaryField() & U2.boundaryField())
+            )
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
    tmp<fvScalarMatrix> pEqnComp1;
    tmp<fvScalarMatrix> pEqnComp2;

@ -160,7 +166,7 @@
        fvScalarMatrix pEqnIncomp
        (
            fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rAUf, p)
        );

        solve
@ -175,7 +181,7 @@

        if (pimple.finalNonOrthogonalIter())
        {
-            surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
+            surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);

            phi1.boundaryField() ==
                mrfZones.relative
@ -200,7 +206,7 @@
            );

            p.relax();
-            mSfGradp = pEqnIncomp.flux()/Dp;
+            mSfGradp = pEqnIncomp.flux()/rAUf;

            U1 = HbyA1
              + fvc::reconstruct
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@ -40,6 +40,7 @@ Description
 #include "PhaseIncompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //