diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
index b96872dded90f932f81a3bcb415bb9823f468d59..67d7c82baaeeae0daf1abaed249c59f3ba51766b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
@@ -1,3 +1,6 @@
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+
     Info<< "Creating twoPhaseSystem\n" << endl;
 
     twoPhaseSystem fluid(mesh, g);
@@ -27,6 +30,13 @@
         fluid.lookup("pMin")
     );
 
+
+    Info<< "Calculating field g.h\n" << endl;
+    dimensionedScalar ghRef(g & (cmptMag(g.value())/mag(g.value()))*hRef);
+    volScalarField gh("gh", (g & mesh.C()) - ghRef);
+    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);
+
+
     rhoThermo& thermo1 = phase1.thermo();
     rhoThermo& thermo2 = phase2.thermo();
 
@@ -38,6 +48,20 @@
     volScalarField& rho2 = thermo2.rho();
     const volScalarField& psi2 = thermo2.psi();
 
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
     volVectorField U
     (
         IOobject
@@ -99,7 +123,14 @@
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+    setRefCell
+    (
+        p,
+        p_rgh,
+        pimple.dict(),
+        pRefCell,
+        pRefValue
+    );
 
     Info<< "Creating field dpdt\n" << endl;
     volScalarField dpdt
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 1c5a9e471aa9f0830a1832605e78a541775b92a6..2cefd8cf82ca1706b99a0219ef585e29c67639b7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -28,49 +28,60 @@
     mrfZones.makeAbsolute(phi2);
 
     // Phase-1 pressure flux (e.g. due to particle-particle pressure)
-    surfaceScalarField phiP1
+    surfaceScalarField phiF1
     (
-        "phiP1",
+        "phiF1",
         fvc::interpolate(rAU1*phase1.turbulence().pPrime())
        *fvc::snGrad(alpha1)*mesh.magSf()
     );
-    phiP1.boundaryField() == 0;
+    phiF1.boundaryField() == 0;
 
     // Phase-2 pressure flux (e.g. due to particle-particle pressure)
-    surfaceScalarField phiP2
+    surfaceScalarField phiF2
     (
-        "phiP2",
+        "phiF2",
         fvc::interpolate(rAU2*phase2.turbulence().pPrime())
        *fvc::snGrad(alpha2)*mesh.magSf()
     );
-    phiP2.boundaryField() == 0;
-
+    phiF2.boundaryField() == 0;
+
+    volScalarField rho("rho", fluid.rho());
+    surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    // Add the phase-1 buoyancy force
+    phiF1 +=
+        rAlphaAU1f
+       *(
+            ghSnGradRho
+          - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
+        )/fvc::interpolate(rho1);
+
+    // Add the phase-2 buoyancy force
+    phiF2 +=
+        rAlphaAU2f
+       *(
+            ghSnGradRho
+          - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
+        )/fvc::interpolate(rho2);
+
+    // Phase-1 predicted flux
     surfaceScalarField phiHbyA1
     (
         IOobject::groupName("phiHbyA", phase1.name()),
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + rAlphaAU1f*fvc::ddtCorr(U1, phi1)
+      + fvc::interpolate(rAU1*dragCoeff)*phi2
+      - phiF1
     );
 
+    // Phase-2 predicted flux
     surfaceScalarField phiHbyA2
     (
         IOobject::groupName("phiHbyA", phase2.name()),
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + rAlphaAU2f*fvc::ddtCorr(U2, phi2)
-    );
-
-    phiHbyA1 +=
-    (
-        fvc::interpolate(rAU1*dragCoeff)*phi2
-      - phiP1
-      + rAlphaAU1f*(g & mesh.Sf())
-    );
-
-    phiHbyA2 +=
-    (
-        fvc::interpolate(rAU2*dragCoeff)*phi1
-      - phiP2
-      + rAlphaAU2f*(g & mesh.Sf())
+      + fvc::interpolate(rAU2*dragCoeff)*phi1
+      - phiF2
     );
 
     mrfZones.makeRelative(phiHbyA1);
@@ -98,7 +109,7 @@
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
     (
-        p.boundaryField(),
+        p_rgh.boundaryField(),
         (
             phiHbyA.boundaryField()
           - mrfZones.relative
@@ -134,8 +145,8 @@
             )/rho1
           + (alpha1/rho1)*correction
             (
-                psi1*fvm::ddt(p)
-              + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p)
+                psi1*fvm::ddt(p_rgh)
+              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
             );
         deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
         pEqnComp1().relax();
@@ -147,8 +158,8 @@
             )/rho2
           + (alpha2/rho2)*correction
             (
-                psi2*fvm::ddt(p)
-              + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p)
+                psi2*fvm::ddt(p_rgh)
+              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
             );
         deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
         pEqnComp2().relax();
@@ -160,31 +171,31 @@
                 contErr1
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
-          + (alpha1*psi1/rho1)*correction(fvm::ddt(p));
+          + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
 
         pEqnComp2 =
             (
                 contErr2
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
             )/rho2
-          + (alpha2*psi2/rho2)*correction(fvm::ddt(p));
+          + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
     }
 
     // Cache p prior to solve for density update
-    volScalarField p_0(p);
+    volScalarField p_rgh_0(p_rgh);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqnIncomp
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(rAUf, p)
+          - fvm::laplacian(rAUf, p_rgh)
         );
 
         solve
         (
             pEqnComp1() + pEqnComp2() + pEqnIncomp,
-            mesh.solver(p.select(pimple.finalInnerIter()))
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
         if (pimple.finalNonOrthogonalIter())
@@ -209,22 +220,22 @@
 
             fluid.dgdt() =
             (
-                alpha1*(pEqnComp2 & p)
-              - alpha2*(pEqnComp1 & p)
+                alpha1*(pEqnComp2 & p_rgh)
+              - alpha2*(pEqnComp1 & p_rgh)
             );
 
-            p.relax();
+            p_rgh.relax();
+
+            p = max(p_rgh + rho*gh, pMin);
+            p_rgh = p - rho*gh;
+
             mSfGradp = pEqnIncomp.flux()/rAUf;
 
             U1 = HbyA1
               + fvc::reconstruct
                 (
-                    rAlphaAU1f
-                   *(
-                        (g & mesh.Sf())
-                      + mSfGradp/fvc::interpolate(rho1)
-                    )
-                  - phiP1
+                    rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
+                  - phiF1
                 );
             U1.correctBoundaryConditions();
             fvOptions.correct(U1);
@@ -232,12 +243,8 @@
             U2 = HbyA2
               + fvc::reconstruct
                 (
-                    rAlphaAU2f
-                   *(
-                        (g & mesh.Sf())
-                      + mSfGradp/fvc::interpolate(rho2)
-                    )
-                  - phiP2
+                    rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
+                  - phiF2
                 );
             U2.correctBoundaryConditions();
             fvOptions.correct(U2);
@@ -246,11 +253,9 @@
         }
     }
 
-    p = max(p, pMin);
-
     // Update densities from change in p
-    rho1 += psi1*(p - p_0);
-    rho2 += psi2*(p - p_0);
+    rho1 += psi1*(p_rgh - p_rgh_0);
+    rho2 += psi2*(p_rgh - p_rgh_0);
 
     K1 = 0.5*magSqr(U1);
     K2 = 0.5*magSqr(U2);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index f22298377e7d0a4c4afb6e695165f39c05538e74..85eb40546558c1544c2c8b54da05e22280fbb02c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -49,7 +49,6 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRFZones.H"
     #include "createFvOptions.H"
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p
index ae2c14b7460e91d82d3717d540b6a9d33dc087e8..961771c727970e5d5259b913d2531a44f025a70f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p
@@ -22,17 +22,17 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
     outlet
     {
-        type               fixedValue;
+        type               calculated;
         value              $internalField;
     }
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..c816028aaf34844b1f6bd2ec1e3460f9dac8ac29
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/0/p_rgh
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes
index ed49e2b09976577388eb6f9d6d45542532ac0609..6a4375df385fd48ea0c5164f15edf2a5bf303e40 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSchemes
@@ -61,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution
index 75cbf243d8720bec824265818a04bde311882439..1e0e5d0f07ba546fea0039cd64b1cb167b7e6d2a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/system/fvSolution
@@ -23,7 +23,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -38,9 +38,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p
index ae2c14b7460e91d82d3717d540b6a9d33dc087e8..961771c727970e5d5259b913d2531a44f025a70f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p
@@ -22,17 +22,17 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
     outlet
     {
-        type               fixedValue;
+        type               calculated;
         value              $internalField;
     }
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..c816028aaf34844b1f6bd2ec1e3460f9dac8ac29
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/0/p_rgh
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes
index a1f49ac1dc7cccbb9eb662cf0fca8c849268064f..49b18d3935377ddc0853664069169bc609a226e0 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSchemes
@@ -63,7 +63,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution
index 000280d807247dfc487db5c9fb9101fcb72c8f99..9868241ad9985ecad147c97118919a84536b1ced 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/system/fvSolution
@@ -23,7 +23,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -38,9 +38,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p
index 08560053d57ba942a2a9776708fdbf73c01af651..b318305599bff814b53f4675657a19c68bde0061 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions          [ 1 -1 -2 0 0 0 0 ];
+dimensions          [1 -1 -2 0 0 0 0];
 
 internalField       uniform 1e5;
 
@@ -22,19 +22,19 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 
     outlet
     {
-        type               fixedValue;
+        type               calculated;
         value              $internalField;
     }
 
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..383649e640bcc4e8c6b3b877e4415a79ef0a36d5
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/p_rgh
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes
index 11ec93ac77ada393bf32133c1a1a20e053b50d55..5c1adb6f1146a8cfbb6eb66c7dc4ad01b2b91c2a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSchemes
@@ -68,7 +68,7 @@ snGradSchemes
 fluxRequired
 {
     default     no;
-    p           ;
+    p_rgh;
     alpha.particles;
 }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution
index 5e444deafe964cf6d3c44b8f88eebcde7647e0c9..0328d987376fdbb0593d5ea4bc87737efca4fa65 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution
@@ -30,7 +30,7 @@ solvers
         minIter         1;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -45,9 +45,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p
index ae2c14b7460e91d82d3717d540b6a9d33dc087e8..961771c727970e5d5259b913d2531a44f025a70f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p
@@ -22,17 +22,17 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
     outlet
     {
-        type               fixedValue;
+        type               calculated;
         value              $internalField;
     }
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..c816028aaf34844b1f6bd2ec1e3460f9dac8ac29
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/0/p_rgh
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes
index a29f58d987710dd9c5e358848819a1c4a0b2bb22..50cb2350665d4ba57e9e8f94e4df717b595c0f6a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSchemes
@@ -60,7 +60,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
index 5987e36dfae7ce07e196d31f5d5412eb2803519f..6487d35276d6bdb3ef9b921b36962e74f482cfa9 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/system/fvSolution
@@ -23,7 +23,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -38,9 +38,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p
index ae2c14b7460e91d82d3717d540b6a9d33dc087e8..961771c727970e5d5259b913d2531a44f025a70f 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p
@@ -22,17 +22,17 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
     outlet
     {
-        type               fixedValue;
+        type               calculated;
         value              $internalField;
     }
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..c816028aaf34844b1f6bd2ec1e3460f9dac8ac29
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/0/p_rgh
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
index 6b34feeb58b7a6d4bb75165fcbe7ad26904dd5e5..e17b4cee5264b88b332fc8e27de7b2eac03f4c9a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSchemes
@@ -62,7 +62,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
index 94ed79268f06e2aceb0a9c090f97afa1162b6d04..171537ad3172756755c138b73845168d27b50805 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/system/fvSolution
@@ -23,7 +23,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -38,9 +38,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p
index 08560053d57ba942a2a9776708fdbf73c01af651..b318305599bff814b53f4675657a19c68bde0061 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions          [ 1 -1 -2 0 0 0 0 ];
+dimensions          [1 -1 -2 0 0 0 0];
 
 internalField       uniform 1e5;
 
@@ -22,19 +22,19 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 
     outlet
     {
-        type               fixedValue;
+        type               calculated;
         value              $internalField;
     }
 
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..383649e640bcc4e8c6b3b877e4415a79ef0a36d5
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/p_rgh
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+
+    outlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes
index 5e0f79e3ce87c3235fb9b6ffaebd61963ea90297..62ff44d8e74f83effa4c19e2d2eedc9906c575e5 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSchemes
@@ -62,7 +62,7 @@ snGradSchemes
 fluxRequired
 {
     default     no;
-    p           ;
+    p_rgh;
     alpha.particles;
 }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution
index fe8ab60ad0209e057cee40631aa04d879e292366..5ca69b0fcd987d96dfd3ebaa7b9c8dda37952b09 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/system/fvSolution
@@ -30,7 +30,7 @@ solvers
         minIter         1;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -45,9 +45,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p
index dda7241905a30c01cc8dc47140ca5a66af19a20b..ae586321875e3dc411075af2f1dbcc32e24e10c3 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p
@@ -22,20 +22,12 @@ boundaryField
 {
     outlet
     {
-        //type               fixedValue;
-        //value              $internalField;
-        type            totalPressure;
-        p0              $internalField;
-        U               U.air;
-        phi             phi.air;
-        rho             thermo:rho.air;
-        psi             none;
-        gamma           1;
-        value           $internalField;
+        type               calculated;
+        value              $internalField;
     }
     walls
     {
-        type               fixedFluxPressure;
+        type               calculated;
         value              $internalField;
     }
 }
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..118e9457ca3fdb38d954161025b220efda69e063
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/0/p_rgh
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    outlet
+    {
+        type            totalPressure;
+        p0              $internalField;
+        U               U.air;
+        phi             phi.air;
+        rho             thermo:rho.air;
+        psi             none;
+        gamma           1;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSchemes
index a29f58d987710dd9c5e358848819a1c4a0b2bb22..50cb2350665d4ba57e9e8f94e4df717b595c0f6a 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSchemes
@@ -60,7 +60,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution
index 5987e36dfae7ce07e196d31f5d5412eb2803519f..6487d35276d6bdb3ef9b921b36962e74f482cfa9 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/injection/system/fvSolution
@@ -23,7 +23,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -38,9 +38,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p
index 72f221208d63bd2bb025e53a1aca041cb7c0c742..f8b3deba004c2c4faa3b934ca21a88e85e2dcd35 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p
@@ -22,13 +22,13 @@ boundaryField
 {
     rotor
     {
-        type            fixedFluxPressure;
+        type            calculated;
         value           $internalField;
     }
 
     stator
     {
-        type            fixedFluxPressure;
+        type            calculated;
         value           $internalField;
     }
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p_rgh b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..3267f0a6f35ce0ba62a997ceb356ed75cca29bbc
--- /dev/null
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/0/p_rgh
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    rotor
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    stator
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes
index 62a486c91ac58183212d97600c933081fdbe3182..2c7bb1ae15927b07482d26a9aaf5a08f8e5ceee6 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSchemes
@@ -60,7 +60,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution
index 2519034730f768970d7dc1c81d650cc2e5df391e..59c0cbf3ab2e0bddacfb803daf62b94c6f9629ca 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/system/fvSolution
@@ -23,7 +23,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -38,9 +38,9 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }