multiphaseEulerFoam: transform to solve for p_rgh

to avoid excessive pressure/buoyancy balance errors on non-orthogonal meshes

applications/solvers/multiphase/multiphaseEulerFoam/correctPhi.H

@@ -2,7 +2,7 @@ CorrectPhi
 (
     U,
     phi,
-    p,
+    p_rgh,
     dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
     geometricZeroField(),
     pimple

applications/solvers/multiphase/multiphaseEulerFoam/createFields.H

@@ -1,11 +1,9 @@
-    #include "readGravitationalAcceleration.H"
+    Info<< "Reading field p_rgh\n" << endl;
-
+    volScalarField p_rgh
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
     (
         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
+    );

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
-    // rho1 = rho10 + psi1*p;
+    // rho1 = rho10 + psi1*p_rgh;
-    // rho2 = rho20 + psi2*p;
+    // 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::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
-    //       + fvc::div(phid1, p)
+    //       + fvc::div(phid1, p_rgh)
-    //       - fvc::Sp(fvc::div(phid1), p);
+    //       - fvc::Sp(fvc::div(phid1), p_rgh);
 
     //     pEqnComp2 =
-    //         fvc::ddt(rho2) + psi2*correction(fvm::ddt(p))
+    //         fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
-    //       + fvc::div(phid2, p)
+    //       + fvc::div(phid2, p_rgh)
-    //       - fvc::Sp(fvc::div(phid2), p);
+    //       - 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()
+               fluid.surfaceTension(phase)*mesh.magSf()
-             + (g & mesh.Sf())
+             + (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]
-                      + rAlphaAUfs[phasei]*mSfGradp/phase.rho()
+                       *(
-                    );
+                            (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;
+    // rho1 = rho10 + psi1*p_rgh;
-    // rho2 = rho20 + psi2*p;
+    // rho2 = rho20 + psi2*p_rgh;
 
-    // Dp1Dt = fvc::DDt(phi1, p);
+    // Dp1Dt = fvc::DDt(phi1, p_rgh);
-    // Dp2Dt = fvc::DDt(phi2, p);
+    // Dp2Dt = fvc::DDt(phi2, p_rgh);
 }

tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/0/p_rgh

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSchemes

@@ -55,7 +55,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh           ;
     pcorr           ;
 }
 

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;
     }

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;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSchemes

@@ -53,7 +53,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
     pcorr;
 }
 

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;
     }

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;
+        U               U;
-        phi             phi.air;
+        phi             phi;
         rho             rho;
         psi             none;
         gamma           1;

tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/controlDict

@@ -33,7 +33,7 @@ writeInterval   0.01;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  6;
 

tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes

@@ -53,7 +53,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
     pcorr;
 }
 

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;
     }

tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/0/p_rgh

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSchemes

@@ -53,7 +53,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
     pcorr;
 }
 

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;
     }