From a7f0de9aa99242f541358b18905649f188e7a5f0 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 19 Mar 2015 21:40:41 +0000
Subject: [PATCH] multiphaseEulerFoam: transform to solve for p_rgh to avoid
 excessive pressure/buoyancy balance errors on non-orthogonal meshes

---
 .../multiphaseEulerFoam/correctPhi.H          |  2 +-
 .../multiphaseEulerFoam/createFields.H        | 47 +++++++++++-----
 .../multiphase/multiphaseEulerFoam/pEqn.H     | 56 +++++++++++--------
 .../bubbleColumn/0/{p => p_rgh}               |  2 +-
 .../bubbleColumn/system/fvSchemes             |  2 +-
 .../bubbleColumn/system/fvSolution            |  8 +--
 .../damBreak4phase/0.org/{p => p_rgh}         |  2 +-
 .../damBreak4phase/system/fvSchemes           |  2 +-
 .../damBreak4phase/system/fvSolution          |  6 +-
 .../damBreak4phaseFine/0.org/{p => p_rgh}     |  6 +-
 .../damBreak4phaseFine/system/controlDict     |  2 +-
 .../damBreak4phaseFine/system/fvSchemes       |  2 +-
 .../damBreak4phaseFine/system/fvSolution      |  6 +-
 .../mixerVessel2D/0/{p => p_rgh}              |  2 +-
 .../mixerVessel2D/system/fvSchemes            |  2 +-
 .../mixerVessel2D/system/fvSolution           |  6 +-
 16 files changed, 92 insertions(+), 61 deletions(-)
 rename tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/{p => p_rgh} (98%)
 rename tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/{p => p_rgh} (98%)
 rename tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/{p => p_rgh} (94%)
 rename tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/{p => p_rgh} (98%)

diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
index 9767adda1dd..9afcd58a66e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H
@@ -2,7 +2,7 @@ CorrectPhi
 (
     U,
     phi,
-    p,
+    p_rgh,
     dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
     geometricZeroField(),
     pimple
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
index 5087d815952..249673a9b2d 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
@@ -1,11 +1,9 @@
-    #include "readGravitationalAcceleration.H"
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -65,12 +63,6 @@
         fluid.lookupOrDefault<scalar>("maxSlamVelocity", GREAT)
     );
 
-    // dimensionedScalar pMin
-    // (
-    //     "pMin",
-    //     dimPressure,
-    //     fluid.lookup("pMin")
-    // );
 
     volScalarField rho
     (
@@ -85,12 +77,39 @@
         fluid.rho()
     );
 
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
 
     // Construct incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
     (
         incompressible::turbulenceModel::New(U, phi, fluid)
     );
+
+
+    #include "readGravitationalAcceleration.H"
+    #include "readhRef.H"
+    #include "gh.H"
+
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rho*gh
+    );
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell
+    (
+        p,
+        p_rgh,
+        pimple.dict(),
+        pRefCell,
+        pRefValue
+    );
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index e8112b300aa..5305e7387e6 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    // rho1 = rho10 + psi1*p;
-    // rho2 = rho20 + psi2*p;
+    // rho1 = rho10 + psi1*p_rgh;
+    // rho2 = rho20 + psi2*p_rgh;
 
     // tmp<fvScalarMatrix> pEqnComp1;
     // tmp<fvScalarMatrix> pEqnComp2;
@@ -14,14 +14,14 @@
     //     surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
 
     //     pEqnComp1 =
-    //         fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
-    //       + fvc::div(phid1, p)
-    //       - fvc::Sp(fvc::div(phid1), p);
+    //         fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+    //       + fvc::div(phid1, p_rgh)
+    //       - fvc::Sp(fvc::div(phid1), p_rgh);
 
     //     pEqnComp2 =
-    //         fvc::ddt(rho2) + psi2*correction(fvm::ddt(p))
-    //       + fvc::div(phid2, p)
-    //       - fvc::Sp(fvc::div(phid2), p);
+    //         fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+    //       + fvc::div(phid2, p_rgh)
+    //       - fvc::Sp(fvc::div(phid2), p_rgh);
     // }
 
     PtrList<surfaceScalarField> alphafs(fluid.phases().size());
@@ -58,6 +58,9 @@
         dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
     );
 
+    volScalarField rho("rho", fluid.rho());
+    surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -104,9 +107,10 @@
         phiHbyAs[phasei] +=
             rAlphaAUfs[phasei]
            *(
-               fluid.surfaceTension(phase)*mesh.magSf()/phase.rho()
-             + (g & mesh.Sf())
-            );
+               fluid.surfaceTension(phase)*mesh.magSf()
+             + (phase.rho() - fvc::interpolate(rho))*(g & mesh.Sf())
+             - ghSnGradRho
+            )/phase.rho();
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
             fluid.dragModels().begin();
@@ -199,7 +203,7 @@
 
         setSnGrad<fixedFluxPressureFvPatchScalarField>
         (
-            p.boundaryField(),
+            p_rgh.boundaryField(),
             (
                 phiHbyA.boundaryField() - mrfZones.relative(phib)
             )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
@@ -211,7 +215,7 @@
         fvScalarMatrix pEqnIncomp
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(rAUf, p)
+          - fvm::laplacian(rAUf, p_rgh)
         );
 
         pEqnIncomp.setReference(pRefCell, pRefValue);
@@ -223,7 +227,7 @@
             //   + (alpha2/rho2)*pEqnComp2()
             // ) +
             pEqnIncomp,
-            mesh.solver(p.select(pimple.finalInnerIter()))
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
         if (pimple.finalNonOrthogonalIter())
@@ -254,7 +258,10 @@
             //   - pos(alpha1)*(pEqnComp1 & p)/rho1
             // );
 
-            p.relax();
+            p_rgh.relax();
+
+            p = p_rgh + rho*gh;
+
             mSfGradp = pEqnIncomp.flux()/rAUf;
 
             U = dimensionedVector("U", dimVelocity, vector::zero);
@@ -269,9 +276,14 @@
                     HbyAs[phasei]
                   + fvc::reconstruct
                     (
-                        rAlphaAUfs[phasei]*(g & mesh.Sf())
-                      + rAlphaAUfs[phasei]*mSfGradp/phase.rho()
-                    );
+                        rAlphaAUfs[phasei]
+                       *(
+                            (phase.rho() - fvc::interpolate(rho))
+                           *(g & mesh.Sf())
+                          - ghSnGradRho
+                          + mSfGradp
+                        )
+                    )/phase.rho();
 
                 //phase.U() = fvc::reconstruct(phase.phi());
                 phase.U().correctBoundaryConditions();
@@ -287,9 +299,9 @@
 
     #include "continuityErrs.H"
 
-    // rho1 = rho10 + psi1*p;
-    // rho2 = rho20 + psi2*p;
+    // rho1 = rho10 + psi1*p_rgh;
+    // rho2 = rho20 + psi2*p_rgh;
 
-    // Dp1Dt = fvc::DDt(phi1, p);
-    // Dp2Dt = fvc::DDt(phi2, p);
+    // Dp1Dt = fvc::DDt(phi1, p_rgh);
+    // Dp2Dt = fvc::DDt(phi2, p_rgh);
 }
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p_rgh
similarity index 98%
rename from tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p
rename to tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p_rgh
index 278db348b8d..73510043119 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p_rgh
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSchemes b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSchemes
index feb67860ae3..a9532d13760 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSchemes
@@ -55,7 +55,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh           ;
     pcorr           ;
 }
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSolution
index 74edd57bbf5..6ed886cf3fb 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSolution
@@ -22,7 +22,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         smoother        DIC;
@@ -37,16 +37,16 @@ solvers
         relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         tolerance       1e-9;
         relTol          0;
     }
 
     pcorr
     {
-        $p;
+        $p_rgh;
         tolerance       1e-5;
         relTol          0;
     }
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/p b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/p_rgh
similarity index 98%
rename from tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/p
rename to tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/p_rgh
index a4f1d992e52..22671e295d2 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/p
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/0.org/p_rgh
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSchemes b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSchemes
index 5a4bb79d962..95a8eb87f13 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSchemes
@@ -53,7 +53,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
     pcorr;
 }
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSolution
index 94423be84d3..827873c725c 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSolution
@@ -22,7 +22,7 @@ solvers
         nAlphaSubCycles 3;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         tolerance       1e-7;
@@ -37,7 +37,7 @@ solvers
         mergeLevels     1;
     }
 
-    pFinal
+    p_rghFinal
     {
         solver          PCG;
         preconditioner
@@ -62,7 +62,7 @@ solvers
 
     pcorr
     {
-        $pFinal;
+        $p_rghFinal;
         tolerance       1e-5;
         relTol          0;
     }
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p_rgh
similarity index 94%
rename from tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p
rename to tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p_rgh
index a4f1d992e52..cfabea4db1c 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/0.org/p_rgh
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,8 +42,8 @@ boundaryField
     {
         type            totalPressure;
         p0              uniform 0;
-        U               U.air;
-        phi             phi.air;
+        U               U;
+        phi             phi;
         rho             rho;
         psi             none;
         gamma           1;
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict
index 90224412879..9020d64f8cd 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict
@@ -33,7 +33,7 @@ writeInterval   0.01;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  6;
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes
index 5a4bb79d962..95a8eb87f13 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes
@@ -53,7 +53,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
     pcorr;
 }
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution
index 94423be84d3..827873c725c 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSolution
@@ -22,7 +22,7 @@ solvers
         nAlphaSubCycles 3;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         tolerance       1e-7;
@@ -37,7 +37,7 @@ solvers
         mergeLevels     1;
     }
 
-    pFinal
+    p_rghFinal
     {
         solver          PCG;
         preconditioner
@@ -62,7 +62,7 @@ solvers
 
     pcorr
     {
-        $pFinal;
+        $p_rghFinal;
         tolerance       1e-5;
         relTol          0;
     }
diff --git a/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p b/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p_rgh
similarity index 98%
rename from tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p
rename to tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p_rgh
index 32d92ebc363..1e324882729 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p
+++ b/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p_rgh
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSchemes b/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSchemes
index 5a4bb79d962..95a8eb87f13 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSchemes
@@ -53,7 +53,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
     pcorr;
 }
 
diff --git a/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSolution
index 4c352dce5f8..de5c1c53ada 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSolution
@@ -22,7 +22,7 @@ solvers
         nAlphaSubCycles 2;
     }
 
-    p
+    p_rgh
     {
         solver          GAMG;
         tolerance       1e-7;
@@ -37,7 +37,7 @@ solvers
         mergeLevels     1;
     }
 
-    pFinal
+    p_rghFinal
     {
         solver          PCG;
         preconditioner
@@ -62,7 +62,7 @@ solvers
 
     pcorr
     {
-        $pFinal;
+        $p_rghFinal;
         tolerance       1e-5;
         relTol          0;
     }
-- 
GitLab