settlingFoam: Changed to solving for pmh (static pressure minus hydrostatic pressure).

While p and pmh (pd in OpenFOAM-1.5.?) are equivalent on orthogonal meshes they
are not on non-orthogonal meshes and the difference is very important for
buoyancy-dominated flows such as settling.  settlingFoam is now written in terms
of pmh (static pressure minus hydrostatic pressure) which used to be called pd
but that confused too may people.

applications/solvers/multiphase/settlingFoam/UEqn.H

@@ -22,8 +22,7 @@
           ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (- ghf*fvc::snGrad(rho) - fvc::snGrad(pmh))*mesh.magSf()
-              - fvc::snGrad(p)*mesh.magSf()
             )
         );
     }

applications/solvers/multiphase/settlingFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field pmh\n" << endl;
-    volScalarField p
+    volScalarField pmh
     (
         IOobject
         (
-            "p",
+            "pmh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -119,11 +119,6 @@
     #include "compressibleCreatePhi.H"
 
 
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
     Info<< "Calculating field mul\n" << endl;
     volScalarField mul
     (
@@ -346,3 +341,48 @@
         ),
         mut + mul
     );
+
+
+    Info<< "Calculating field (g.h)f\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf());
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pmh + rho*gh
+    );
+
+    label pmhRefCell = 0;
+    scalar pmhRefValue = 0.0;
+    setRefCell
+    (
+        pmh,
+        mesh.solutionDict().subDict("PISO"),
+        pmhRefCell,
+        pmhRefValue
+    );
+
+    scalar pRefValue = 0.0;
+
+    if (pmh.needReference())
+    {
+        pRefValue = readScalar
+        (
+            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pmhRefCell)
+        );
+    }

applications/solvers/multiphase/settlingFoam/pmhEqn.H

@@ -0,0 +1,52 @@
+volScalarField rUA = 1.0/UEqn.A();
+
+surfaceScalarField rUAf
+(
+    "(rho*(1|A(U)))",
+    fvc::interpolate(rho)*fvc::interpolate(rUA)
+);
+
+U = rUA*UEqn.H();
+phi =
+    fvc::interpolate(rho)
+   *(
+       (fvc::interpolate(U) & mesh.Sf())
+     + fvc::ddtPhiCorr(rUA, rho, U, phi)
+    );
+
+surfaceScalarField phiU("phiU", phi);
+phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf();
+
+for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+{
+    fvScalarMatrix pmhEqn
+    (
+        fvm::laplacian(rUAf, pmh) == fvc::ddt(rho) + fvc::div(phi)
+    );
+
+    pmhEqn.setReference(pmhRefCell, pmhRefValue);
+    pmhEqn.solve();
+
+    if (nonOrth == nNonOrthCorr)
+    {
+        phi -= pmhEqn.flux();
+    }
+}
+
+p == pmh + rho*gh;
+
+if (pmh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pmhRefCell)
+    );
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+U.correctBoundaryConditions();

applications/solvers/multiphase/settlingFoam/settlingFoam.C

@@ -64,22 +64,24 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "calcVdj.H"
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
-
-        #include "UEqn.H"
-
-        #include "alphaEqn.H"
-
-        #include "correctViscosity.H"
-
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
         {
-            #include "pEqn.H"
+            #include "calcVdj.H"
-        }
 
-        #include "kEpsilon.H"
+            #include "UEqn.H"
+
+            #include "alphaEqn.H"
+
+            #include "correctViscosity.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pmhEqn.H"
+            }
+
+            #include "kEpsilon.H"
+        }
 
         runTime.write();
 

tutorials/multiphase/settlingFoam/ras/dahl/0/pmh

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

tutorials/multiphase/settlingFoam/ras/dahl/constant/polyMesh/boundary

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6.x                                 |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/

tutorials/multiphase/settlingFoam/ras/dahl/system/fvSchemes

@@ -39,7 +39,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),pmh) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(mut,Alpha) Gauss linear corrected;
@@ -58,7 +58,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    pmh             ;
 }
 
 

tutorials/multiphase/settlingFoam/ras/dahl/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    pmh
     {
         solver          PCG;
         preconditioner  DIC;

tutorials/multiphase/settlingFoam/ras/tank3D/0/U

@@ -57,7 +57,8 @@ boundaryField
 
     WALL6
     {
-        type            fixedValue;
+        type            translatingWallVelocity;
+        U               (-0.003 0 0);
         value           uniform (-0.003 0 0);
     }
 

tutorials/multiphase/settlingFoam/ras/tank3D/0/pmh

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      p;
+    object      pmh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -28,27 +28,32 @@ boundaryField
 
     INLE1
     {
-        type            zeroGradient;
+        type            buoyantPressure;
+        value           uniform 0;
     }
 
     OUTL9
     {
-        type            zeroGradient;
+        type            buoyantPressure;
+        value           uniform 0;
     }
 
     OUTL10
     {
-        type            zeroGradient;
+        type            buoyantPressure;
+        value           uniform 0;
     }
 
     OUTL11
     {
-        type            zeroGradient;
+        type            buoyantPressure;
+        value           uniform 0;
     }
 
     OUTL12
     {
-        type            zeroGradient;
+        type            buoyantPressure;
+        value           uniform 0;
     }
 
     WALL6

tutorials/multiphase/settlingFoam/ras/tank3D/Allclean

@@ -3,4 +3,4 @@
 # Clean time folders only
 
 rm -rf *[1-9]*
-rm -f log.* 2>/dev/null
+rm log.*

tutorials/multiphase/settlingFoam/ras/tank3D/system/controlDict

@@ -17,7 +17,7 @@ FoamFile
 
 application     settlingFoam;
 
-startFrom       startTime;
+startFrom       latestTime;
 
 startTime       0;
 

tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes

@@ -39,7 +39,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),pmh) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(mut,Alpha) Gauss linear corrected;
@@ -58,7 +58,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    pmh             ;
 }
 
 

tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    pmh
     {
         solver          PCG;
         preconditioner  DIC;
@@ -68,9 +68,17 @@ solvers
 
 PISO
 {
-    nCorrectors     2;
+    momentumPredictor yes;
+    nOuterCorrectors  1;
+    nCorrectors       2;
     nNonOrthogonalCorrectors 0;
 }
 
+relaxationFactors
+{
+    U               1;
+    k               1;
+    epsilon         1;
+}
 
 // ************************************************************************* //