Merged all multiphase developments in OpenFOAM-1.7.x

applications/solvers/heatTransfer/buoyantSimpleFoam_old/Make/files

@@ -0,0 +1,3 @@
+buoyantSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/buoyantSimpleFoam

applications/solvers/heatTransfer/buoyantSimpleFoam_old/Make/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
+    -lfiniteVolume

applications/solvers/heatTransfer/buoyantSimpleFoam_old/UEqn.H

@@ -0,0 +1,24 @@
+    // Solve the Momentum equation
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn().relax();
+
+    solve
+    (
+        UEqn()
+     ==
+        rho*g
+      - fvc::grad(p)
+        /*
+        fvc::reconstruct
+        (
+            fvc::interpolate(rho)*(g & mesh.Sf())
+          - fvc::snGrad(p)*mesh.magSf()
+        )
+        */
+    );

applications/solvers/heatTransfer/buoyantSimpleFoam_old/buoyantSimpleFoam.C

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    buoyantSimpleFoam
+
+Description
+    Steady-state solver for buoyant, turbulent flow of compressible fluids
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "RASModel.H"
+#include "fixedGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readSIMPLEControls.H"
+
+        p.storePrevIter();
+        rho.storePrevIter();
+
+        // Pressure-velocity SIMPLE corrector
+        {
+            #include "UEqn.H"
+            #include "hEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/buoyantSimpleFoam_old/createFields.H

@@ -0,0 +1,69 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::RASModel> turbulence
+    (
+        compressible::RASModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    thermo.correct();
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell
+    (
+        p,
+        mesh.solutionDict().subDict("SIMPLE"),
+        pRefCell,
+        pRefValue
+    );
+
+
+    dimensionedScalar initialMass = fvc::domainIntegrate(rho);

applications/solvers/heatTransfer/buoyantSimpleFoam_old/hEqn.H

@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::div(phi, h)
+      - fvm::Sp(fvc::div(phi), h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
+      - p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo.correct();
+}

applications/solvers/heatTransfer/buoyantSimpleFoam_old/pEqn.H

@@ -0,0 +1,57 @@
+{
+    rho = thermo.rho();
+
+    volScalarField rUA = 1.0/UEqn().A();
+    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+
+    U = rUA*UEqn().H();
+    UEqn.clear();
+
+    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    bool closedVolume = adjustPhi(phi, U, p);
+
+    surfaceScalarField buoyancyPhi =
+        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+    phi += buoyancyPhi;
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rhorUAf, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            // For closed-volume cases adjust the pressure and density levels
+            // to obey overall mass continuity
+            if (closedVolume)
+            {
+                p += (initialMass - fvc::domainIntegrate(psi*p))
+                    /fvc::domainIntegrate(psi);
+            }
+
+            // Calculate the conservative fluxes
+            phi -= pEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U += rUA*(rho*g - fvc::grad(p));
+            //U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    rho = thermo.rho();
+    rho.relax();
+    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
+        << endl;
+}

applications/solvers/multiphase/cavitatingFoam/Make/options

@@ -7,8 +7,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude
 
 EXE_LIBS = \
-    -lincompressibleTurbulenceModel \
     -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume \

applications/solvers/multiphase/cavitatingFoam/UEqn.H

@@ -20,3 +20,5 @@
     {
         solve(UEqn == -fvc::grad(p));
     }
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C

@@ -25,7 +25,8 @@ Application
     cavitatingFoam
 
 Description
-    Transient cavitation code based on the barotropic equation of state.
+    Transient cavitation code based on the homogeneous equilibrium model
+    from which the compressibility of the liquid/vapour "mixture" is obtained.
 
     Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -52,9 +52,32 @@
         }
     }
 
-    Info<< "max-min p: " << max(p).value()
+    Info<< "Predicted p max-min : " << max(p).value()
         << " " << min(p).value() << endl;
 
+    rho == max
+    (
+        psi*p
+      + (1.0 - gamma)*rhol0
+      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+        rhoMin
+    );
+
+    #include "gammaPsi.H"
+
+    p =
+    (
+        rho
+      - (1.0 - gamma)*rhol0
+      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+    )/psi;
+
+    p.correctBoundaryConditions();
+
+    Info<< "Phase-change corrected p max-min : " << max(p).value()
+        << " " << min(p).value() << endl;
+
+    // Correct velocity
 
     U = HbyA - rUA*fvc::grad(p);
 
@@ -70,17 +93,4 @@
     U.correctBoundaryConditions();
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-
-    rho == max
-    (
-        psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
-        rhoMin
-    );
-
-    Info<< "max-min rho: " << max(rho).value()
-        << " " << min(rho).value() << endl;
-
-    #include "gammaPsi.H"
 }

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H

@@ -1,86 +0,0 @@
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    surfaceScalarField phir = phic*interface.nHatf();
-
-    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
-    {
-        volScalarField::DimensionedInternalField Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
-        );
-
-        volScalarField::DimensionedInternalField Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            divU*min(alpha1, scalar(1))
-        );
-
-        forAll(dgdt, celli)
-        {
-            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
-            {
-                Sp[celli] -= dgdt[celli]*alpha1[celli];
-                Su[celli] += dgdt[celli]*alpha1[celli];
-            }
-            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
-            {
-                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
-            }
-        }
-
-
-        surfaceScalarField phiAlpha1 =
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                alphaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, alpha2, alpharScheme),
-                alpha1,
-                alpharScheme
-            );
-
-        MULES::explicitSolve
-        (
-            geometricOneField(),
-            alpha1,
-            phi,
-            phiAlpha1,
-            Sp,
-            Su,
-            1,
-            0
-        );
-
-        surfaceScalarField rho1f = fvc::interpolate(rho1);
-        surfaceScalarField rho2f = fvc::interpolate(rho2);
-        rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
-
-        alpha2 = scalar(1) - alpha1;
-    }
-
-    Info<< "Liquid phase volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
-        << "  Min(alpha1) = " << min(alpha1).value()
-        << "  Min(alpha2) = " << min(alpha2).value()
-        << endl;
-}

applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

@@ -1,144 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Calculating field alpha1\n" << endl;
-    volScalarField alpha2("alpha2", scalar(1) - alpha1);
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "createPhi.H"
-
-
-    Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi);
-
-    dimensionedScalar rho10
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase1Name()
-        ).lookup("rho0")
-    );
-
-    dimensionedScalar rho20
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase2Name()
-        ).lookup("rho0")
-    );
-
-    dimensionedScalar psi1
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase1Name()
-        ).lookup("psi")
-    );
-
-    dimensionedScalar psi2
-    (
-        twoPhaseProperties.subDict
-        (
-            twoPhaseProperties.phase2Name()
-        ).lookup("psi")
-    );
-
-    dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
-
-    volScalarField rho1 = rho10 + psi1*p;
-    volScalarField rho2 = rho20 + psi2*p;
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        alpha1*rho1 + alpha2*rho2
-    );
-
-
-    // Mass flux
-    // Initialisation does not matter because rhoPhi is reset after the
-    // alpha1 solution before it is used in the U equation.
-    surfaceScalarField rhoPhi
-    (
-        IOobject
-        (
-            "rho*phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        fvc::interpolate(rho)*phi
-    );
-
-    volScalarField dgdt =
-        pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
-
-    // Construct interface from alpha1 distribution
-    interfaceProperties interface(alpha1, U, twoPhaseProperties);
-
-    // Construct incompressible turbulence model
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
-    );
-
-
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll(p.boundaryField(), i)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }

applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H

@@ -1,95 +0,0 @@
-{
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
-
-    tmp<fvScalarMatrix> pEqnComp;
-
-    if (transonic)
-    {
-        pEqnComp =
-            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
-    }
-    else
-    {
-        pEqnComp =
-            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
-    }
-
-
-    U = rUA*UEqn.H();
-
-    surfaceScalarField phiU
-    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
-    );
-
-    phi = phiU +
-        (
-            fvc::interpolate(interface.sigmaK())
-           *fvc::snGrad(alpha1)*mesh.magSf()
-          + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf;
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix pEqnIncomp
-        (
-            fvc::div(phi)
-          - fvm::laplacian(rUAf, p)
-        );
-
-        if
-        (
-            oCorr == nOuterCorr-1
-            && corr == nCorr-1
-            && nonOrth == nNonOrthCorr
-        )
-        {
-            solve
-            (
-                (
-                    max(alpha1, scalar(0))*(psi1/rho1)
-                  + max(alpha2, scalar(0))*(psi2/rho2)
-                )
-               *pEqnComp()
-              + pEqnIncomp,
-                mesh.solver(p.name() + "Final")
-            );
-        }
-        else
-        {
-            solve
-            (
-                (
-                    max(alpha1, scalar(0))*(psi1/rho1)
-                  + max(alpha2, scalar(0))*(psi2/rho2)
-                )
-               *pEqnComp()
-              + pEqnIncomp
-            );
-        }
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            dgdt =
-                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pEqnComp & p);
-            phi += pEqnIncomp.flux();
-        }
-    }
-
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
-    U.correctBoundaryConditions();
-
-    p.max(pMin);
-
-    rho1 = rho10 + psi1*p;
-    rho2 = rho20 + psi2*p;
-
-    Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(p) " << min(p).value() << endl;
-
-    // Make the fluxes relative to the mesh motion
-    fvc::makeRelative(phi, U);
-}

applications/solvers/multiphase/compressibleInterFoam/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean
+wclean compressibleInterDyMFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleInterFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+wmake compressibleInterDyMFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleInterFoam/Make/options

@@ -6,7 +6,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - fvc::snGrad(p)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I../interFoam \
+    -I.. \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
@@ -10,7 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
@@ -18,5 +18,4 @@ EXE_LIBS = \
     -lfiniteVolume \
     -ldynamicMesh \
     -lmeshTools \
-    -ldynamicFvMesh \
+    -ldynamicFvMesh
-    -ltopoChangerFvMesh

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/alphaEqnsSubCycle.H

@@ -27,7 +27,7 @@
             !(++alphaSubCycle).end();
         )
         {
-#           include "alphaEqns.H"
+            #include "alphaEqns.H"
             rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
         }
 
@@ -35,7 +35,7 @@
     }
     else
     {
-#       include "alphaEqns.H"
+        #include "alphaEqns.H"
     }
 
     if (oCorr == 0)

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    compressibleLesInterFoam
+    compressibleInterDyMFoam
 
 Description
     Solver for 2 compressible, isothermal immiscible fluids using a VOF
@@ -30,9 +30,10 @@ Description
     with optional mesh motion and mesh topology changes including adaptive
     re-meshing.
 
-    The momentum and other fluid properties are of the "mixture" and a
+    The momentum and other fluid properties are of the "mixture" and a single
-    single momentum equation is solved.  Turbulence modelling is generic,
+    momentum equation is solved.
-    i.e. laminar, RAS or LES may be selected.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
@@ -55,6 +56,7 @@ int main(int argc, char *argv[])
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
+    #include "createPcorrTypes.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -89,6 +91,9 @@ int main(int argc, char *argv[])
                 Info<< "Execution time for mesh.update() = "
                     << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                     << " s" << endl;
+
+                gh = g & mesh.C();
+                ghf = g & mesh.Cf();
             }
 
             if (mesh.changing() && correctPhi)

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H

@@ -42,7 +42,7 @@
 
     adjustPhi(phi, U, pcorr);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pcorrEqn
         (

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/createPcorrTypes.H

@@ -0,0 +1,13 @@
+    wordList pcorrTypes
+    (
+        p_rgh.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    for (label i=0; i<p_rgh.boundaryField().size(); i++)
+    {
+        if (p_rgh.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -0,0 +1,96 @@
+{
+    volScalarField rUA = 1.0/UEqn.A();
+    surfaceScalarField rUAf = fvc::interpolate(rUA);
+
+    tmp<fvScalarMatrix> p_rghEqnComp;
+
+    if (transonic)
+    {
+        p_rghEqnComp =
+        (
+            fvm::ddt(p_rgh)
+          + fvm::div(phi, p_rgh)
+          - fvm::Sp(fvc::div(phi), p_rgh)
+        );
+    }
+    else
+    {
+        p_rghEqnComp =
+        (
+            fvm::ddt(p_rgh)
+          + fvc::div(phi, p_rgh)
+          - fvc::Sp(fvc::div(phi), p_rgh)
+        );
+    }
+
+
+    U = rUA*UEqn.H();
+
+    surfaceScalarField phiU
+    (
+        "phiU",
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+    );
+
+    phi = phiU +
+        (
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+          - ghf*fvc::snGrad(rho)
+        )*rUAf*mesh.magSf();
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix p_rghEqnIncomp
+        (
+            fvc::div(phi)
+          - fvm::laplacian(rUAf, p_rgh)
+        );
+
+        solve
+        (
+            (
+                max(alpha1, scalar(0))*(psi1/rho1)
+              + max(alpha2, scalar(0))*(psi2/rho2)
+            )
+           *p_rghEqnComp()
+          + p_rghEqnIncomp,
+            mesh.solver
+            (
+                p_rgh.select
+                (
+                    oCorr == nOuterCorr-1
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
+                )
+            )
+        );
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            dgdt =
+                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
+               *(p_rghEqnComp & p_rgh);
+            phi += p_rghEqnIncomp.flux();
+        }
+    }
+
+    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U.correctBoundaryConditions();
+
+    p = max
+    (
+        (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
+       /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+        pMin
+    );
+
+    rho1 = rho10 + psi1*p;
+    rho2 = rho20 + psi2*p;
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H

@@ -1,5 +1,5 @@
-    #include "readPISOControls.H"
+   #include "readPISOControls.H"
-    #include "readTimeControls.H"
+   #include "readTimeControls.H"
 
     label nAlphaCorr
     (
@@ -19,9 +19,14 @@
             << exit(FatalError);
     }
 
-    const bool correctPhi =
+    bool correctPhi = true;
-        piso.lookupOrDefault("correctPhi", true);
+    if (piso.found("correctPhi"))
-
+    {
-    const bool checkMeshCourantNo =
+        correctPhi = Switch(piso.lookup("correctPhi"));
-        piso.lookupOrDefault("checkMeshCourantNo", false);
+    }
 
+    bool checkMeshCourantNo = false;
+    if (piso.found("checkMeshCourantNo"))
+    {
+        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
+    }

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -83,6 +83,28 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        max
+        (
+            (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
+           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+            pMin
+        )
+    );
+
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -2,17 +2,25 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pEqnComp;
+    tmp<fvScalarMatrix> p_rghEqnComp;
 
     if (transonic)
     {
-        pEqnComp =
+        p_rghEqnComp =
-            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
+        (
+            fvm::ddt(p_rgh)
+          + fvm::div(phi, p_rgh)
+          - fvm::Sp(fvc::div(phi), p_rgh)
+        );
     }
     else
     {
-        pEqnComp =
+        p_rghEqnComp =
-            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
+        (
+            fvm::ddt(p_rgh)
+          + fvc::div(phi, p_rgh)
+          - fvc::Sp(fvc::div(phi), p_rgh)
+        );
     }
 
 
@@ -21,22 +29,22 @@
     surfaceScalarField phiU
     (
         "phiU",
-        (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-           *fvc::snGrad(alpha1)*mesh.magSf()
+          - ghf*fvc::snGrad(rho)
-          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf*mesh.magSf();
-        )*rUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqnIncomp
+        fvScalarMatrix p_rghEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, p)
+          - fvm::laplacian(rUAf, p_rgh)
         );
 
         solve
@@ -45,27 +53,41 @@
                 max(alpha1, scalar(0))*(psi1/rho1)
               + max(alpha2, scalar(0))*(psi2/rho2)
             )
-           *pEqnComp()
+           *p_rghEqnComp()
-          + pEqnIncomp
+          + p_rghEqnIncomp,
+            mesh.solver
+            (
+                p_rgh.select
+                (
+                    oCorr == nOuterCorr-1
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
+                )
+            )
         );
 
         if (nonOrth == nNonOrthCorr)
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pEqnComp & p);
+               *(p_rghEqnComp & p_rgh);
-            phi += pEqnIncomp.flux();
+            phi += p_rghEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p.max(pMin);
+    p = max
+    (
+        (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
+       /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+        pMin
+    );
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(p) " << min(p).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
 }

applications/solvers/multiphase/interDyMFoam/correctPhi.H

@@ -1,62 +0,0 @@
-{
-    if (mesh.changing())
-    {
-        forAll(U.boundaryField(), patchi)
-        {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].initEvaluate();
-            }
-        }
-
-        forAll(U.boundaryField(), patchi)
-        {
-            if (U.boundaryField()[patchi].fixesValue())
-            {
-                U.boundaryField()[patchi].evaluate();
-
-                phi.boundaryField()[patchi] =
-                U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
-            }
-        }
-    }
-
-    #include "continuityErrs.H"
-
-    volScalarField pcorr
-    (
-        IOobject
-        (
-            "pcorr",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("pcorr", p.dimensions(), 0.0),
-        pcorrTypes
-    );
-
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
-    adjustPhi(phi, U, pcorr);
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix pcorrEqn
-        (
-            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
-        );
-
-        pcorrEqn.setReference(pRefCell, pRefValue);
-        pcorrEqn.solve();
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi -= pcorrEqn.flux();
-        }
-    }
-
-    #include "continuityErrs.H"
-}

applications/solvers/multiphase/interDyMFoam/createFields.H

@@ -1,112 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi);
-
-    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
-    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
-
-
-    // Need to store rho for ddt(rho, U)
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT
-        ),
-        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
-        alpha1.boundaryField().types()
-    );
-    rho.oldTime();
-
-
-    // Mass flux
-    // Initialisation does not matter because rhoPhi is reset after the
-    // alpha1 solution before it is used in the U equation.
-    surfaceScalarField rhoPhi
-    (
-        IOobject
-        (
-            "rho*phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        rho1*phi
-    );
-
-
-    // Construct interface from alpha1 distribution
-    interfaceProperties interface(alpha1, U, twoPhaseProperties);
-
-    // Construct incompressible turbulence model
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
-    );
-
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll(p.boundaryField(), i)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -1,52 +0,0 @@
-{
-    volScalarField rAU = 1.0/UEqn.A();
-    surfaceScalarField rAUf = fvc::interpolate(rAU);
-
-    U = rAU*UEqn.H();
-    surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
-
-    if (p.needReference())
-    {
-        fvc::makeRelative(phiU, U);
-        adjustPhi(phiU, U, p);
-        fvc::makeAbsolute(phiU, U);
-    }
-
-    phi = phiU +
-    (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
-      + fvc::interpolate(rho)*(g & mesh.Sf())
-    )*rAUf;
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::laplacian(rAUf, p) == fvc::div(phi)
-        );
-
-        pEqn.setReference(pRefCell, pRefValue);
-
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
-        {
-            pEqn.solve(mesh.solver(p.name() + "Final"));
-        }
-        else
-        {
-            pEqn.solve(mesh.solver(p.name()));
-        }
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi -= pEqn.flux();
-        }
-    }
-
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
-    #include "continuityErrs.H"
-
-    // Make the fluxes relative to the mesh motion
-    fvc::makeRelative(phi, U);
-}

applications/solvers/multiphase/interDyMFoam/readControls.H

@@ -1,9 +0,0 @@
-#   include "readTimeControls.H"
-#   include "readPISOControls.H"
-
-    const bool correctPhi =
-        piso.lookupOrDefault("correctPhi", true);
-
-    const bool checkMeshCourantNo =
-        piso.lookupOrDefault("checkMeshCourantNo", false);
-

applications/solvers/multiphase/interFoam/Allwclean

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean
+wclean interDyMFoam
+wclean MRFInterFoam
+wclean porousInterFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/interFoam/Allwmake

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+wmake interDyMFoam
+wmake MRFInterFoam
+wmake porousInterFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C

@@ -51,7 +51,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -70,6 +69,7 @@ int main(int argc, char *argv[])
         #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
@@ -79,6 +79,7 @@ int main(int argc, char *argv[])
         twoPhaseProperties.correct();
 
         #include "alphaEqnSubCycle.H"
+        #include "zonePhaseVolumes.H"
 
         #include "UEqn.H"
 
@@ -88,8 +89,6 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        #include "continuityErrs.H"
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/interFoam/MRFInterFoam/Make/files

@@ -0,0 +1,3 @@
+MRFInterFoam.C
+
+EXE = $(FOAM_APPBIN)/MRFInterFoam

applications/solvers/multiphase/interFoam/MRFInterFoam/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(FOAM_SOLVERS)/multiphase/interFoam \
+    -I.. \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
@@ -7,8 +7,9 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume

applications/solvers/multiphase/interFoam/MRFInterFoam/UEqn.H

@@ -25,10 +25,10 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - fvc::snGrad(p)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H

@@ -0,0 +1,62 @@
+{
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
+
+    U = rAU*UEqn.H();
+    surfaceScalarField phiU
+    (
+        "phiU",
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+    );
+    mrfZones.relativeFlux(phiU);
+
+    adjustPhi(phiU, U, p_rgh);
+
+    phi = phiU +
+    (
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+      - ghf*fvc::snGrad(rho)
+    )*rAUf*mesh.magSf();
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve
+        (
+            mesh.solver
+            (
+                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
+            )
+        );
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= p_rghEqn.flux();
+        }
+    }
+
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/multiphase/interFoam/MRFInterFoam/zonePhaseVolumes.H

@@ -0,0 +1,21 @@
+{
+    const scalarField& V = mesh.V();
+
+    forAll(mesh.cellZones(), czi)
+    {
+        const labelList& cellLabels = mesh.cellZones()[czi];
+
+        scalar phaseVolume = 0;
+
+        forAll(cellLabels, cli)
+        {
+            label celli = cellLabels[cli];
+            phaseVolume += alpha1[celli]*V[celli];
+        }
+
+        reduce(phaseVolume, sumOp<scalar>());
+
+        Info<< "Phase volume in zone " << mesh.cellZones()[czi].name()
+            << " = " << phaseVolume*1e6 << " ml " << endl;
+    }
+}

applications/solvers/multiphase/interFoam/Make/options

@@ -6,7 +6,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \

applications/solvers/multiphase/interFoam/UBlendingFactor.H

@@ -1,6 +0,0 @@
-    surfaceScalarField alpha1f = fvc::interpolate(alpha1);
-    surfaceScalarField UBlendingFactor
-    (
-        "UBlendingFactor",
-        sqrt(max(min(4*alpha1f*(1.0 - alpha1f), 1.0), 0.0))
-    );

applications/solvers/multiphase/interFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - fvc::snGrad(p)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interFoam/alphaCourantNo.H

@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    surfaceScalarField alphaf = fvc::interpolate(alpha1);
+
+    surfaceScalarField SfUfbyDelta =
+        pos(alphaf - 0.01)*pos(0.99 - alphaf)
+       *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
+
+    alphaCoNum = max(SfUfbyDelta/mesh.magSf())
+        .value()*runTime.deltaT().value();
+
+    meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
+        .value()*runTime.deltaT().value();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -3,13 +3,13 @@
 
     wordList pcorrTypes
     (
-        p.boundaryField().size(),
+        p_rgh.boundaryField().size(),
         zeroGradientFvPatchScalarField::typeName
     );
 
-    forAll(p.boundaryField(), i)
+    forAll (p_rgh.boundaryField(), i)
     {
-        if (p.boundaryField()[i].fixesValue())
+        if (p_rgh.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
@@ -26,11 +26,11 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", p.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p_rgh.dimensions(), 0.0),
         pcorrTypes
     );
 
-    dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
+    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
 
     adjustPhi(phi, U, pcorr);
 
@@ -38,7 +38,7 @@
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);

applications/solvers/multiphase/interFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -40,7 +40,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
 
     Info<< "Reading transportProperties\n" << endl;
@@ -83,11 +83,6 @@
     );
 
 
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties);
 
@@ -97,3 +92,54 @@
     (
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
+
+    #include "readGravitationalAcceleration.H"
+
+    /*
+    dimensionedVector g0(g);
+
+    // Read the data file and initialise the interpolation table
+    interpolationTable<vector> timeSeriesAcceleration
+    (
+        runTime.path()/runTime.caseConstant()/"acceleration.dat"
+    );
+    */
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    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,
+        mesh.solutionDict().subDict("PISO"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }

applications/solvers/multiphase/interFoam/interDyMFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I.. \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
@@ -9,7 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
@@ -19,4 +20,3 @@ EXE_LIBS = \
     -lmeshTools \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh
-

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -47,7 +47,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -62,6 +61,7 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readControls.H"
+        #include "alphaCourantNo.H"
         #include "CourantNo.H"
 
         // Make the fluxes absolute
@@ -83,6 +83,9 @@ int main(int argc, char *argv[])
             Info<< "Execution time for mesh.update() = "
                 << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
                 << " s" << endl;
+
+            gh = g & mesh.C();
+            ghf = g & mesh.Cf();
         }
 
         if (mesh.changing() && correctPhi)

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -0,0 +1,64 @@
+{
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
+
+    U = rAU*UEqn.H();
+    surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiU, U);
+        adjustPhi(phiU, U, p_rgh);
+        fvc::makeAbsolute(phiU, U);
+    }
+
+    phi = phiU +
+    (
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+      - ghf*fvc::snGrad(rho)
+    )*rAUf*mesh.magSf();
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve
+        (
+            mesh.solver
+            (
+                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
+            )
+        );
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= p_rghEqn.flux();
+        }
+    }
+
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/multiphase/interFoam/interDyMFoam/readControls.H

@@ -0,0 +1,14 @@
+#   include "readTimeControls.H"
+#   include "readPISOControls.H"
+
+    bool correctPhi = true;
+    if (piso.found("correctPhi"))
+    {
+        correctPhi = Switch(piso.lookup("correctPhi"));
+    }
+
+    bool checkMeshCourantNo = false;
+    if (piso.found("checkMeshCourantNo"))
+    {
+        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
+    }

applications/solvers/multiphase/interFoam/interFoam.C

@@ -43,6 +43,7 @@ Description
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "interpolationTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +52,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -69,6 +69,7 @@ int main(int argc, char *argv[])
         #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
@@ -87,8 +88,6 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        #include "continuityErrs.H"
-
         turbulence->correct();
 
         runTime.write();

applications/solvers/multiphase/interFoam/pEqn.H

@@ -1,49 +1,61 @@
 {
-    volScalarField rUA = 1.0/UEqn.A();
+    volScalarField rAU = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
-
-    U = rUA*UEqn.H();
 
+    U = rAU*UEqn.H();
     surfaceScalarField phiU
     (
         "phiU",
         (fvc::interpolate(U) & mesh.Sf())
-      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p);
+    adjustPhi(phiU, U, p_rgh);
 
     phi = phiU +
-        (
+    (
-            fvc::interpolate(interface.sigmaK())
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-           *fvc::snGrad(alpha1)*mesh.magSf()
+      - ghf*fvc::snGrad(rho)
-          + fvc::interpolate(rho)*(g & mesh.Sf())
+    )*rAUf*mesh.magSf();
-        )*rUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rUAf, p) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        p_rghEqn.solve
-        {
+        (
-            pEqn.solve(mesh.solver(p.name() + "Final"));
+            mesh.solver
-        }
+            (
-        else
+                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-        {
+            )
-            pEqn.solve(mesh.solver(p.name()));
+        );
-        }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pEqn.flux();
+            phi -= p_rghEqn.flux();
         }
     }
 
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
     U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
 }

applications/solvers/multiphase/interFoam/porousInterFoam/Make/files

@@ -0,0 +1,3 @@
+porousInterFoam.C
+
+EXE = $(FOAM_APPBIN)/porousInterFoam

applications/solvers/multiphase/interFoam/porousInterFoam/Make/options

@@ -0,0 +1,17 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -ltwoPhaseInterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/multiphase/interFoam/porousInterFoam/UEqn.H

@@ -5,17 +5,22 @@
       + fvc::interpolate(rho*turbulence->nut())
     );
 
+    // Calculate and cache mu for the porous media
+    volScalarField mu(twoPhaseProperties.mu());
+
     fvVectorMatrix UEqn
     (
-        fvm::ddt(rho, U)
+        pZones.ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::laplacian(muEff, U)
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
+    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();
 
+    pZones.addResistance(UEqn);
+
     if (momentumPredictor)
     {
         solve
@@ -24,10 +29,10 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - fvc::snGrad(p)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interFoam/porousInterFoam/createPorousZones.H

@@ -0,0 +1 @@
+    porousZones pZones(mesh);

applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    porousInterFoam
+
+Description
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+    For a two-fluid approach see twoPhaseEulerFoam.
+
+    Explicit handling of porous zones is included.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "interfaceProperties.H"
+#include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "porousZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createPorousZones.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        twoPhaseProperties.correct();
+
+        #include "alphaEqnSubCycle.H"
+
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/setDeltaT.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum courant Number.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    scalar maxDeltaTFact =
+        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaT().value(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
+}
+
+// ************************************************************************* //

applications/solvers/multiphase/interMixingFoam/Make/options

@@ -5,12 +5,13 @@ EXE_INC = \
     -IincompressibleThreePhaseMixture \
     -IthreePhaseInterfaceProperties \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseInterfaceProperties/alphaContactAngle/alphaContactAngle \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \

applications/solvers/multiphase/interMixingFoam/alphaCourantNo.H

@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2);
+
+    surfaceScalarField SfUfbyDelta = max
+    (
+        pos(alpha1f - 0.01)*pos(0.99 - alpha1f),
+        pos(alpha2f - 0.01)*pos(0.99 - alpha2f)
+    )*mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
+
+    alphaCoNum = max(SfUfbyDelta/mesh.magSf())
+        .value()*runTime.deltaT().value();
+
+    meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
+        .value()*runTime.deltaT().value();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/interMixingFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -73,7 +73,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
     threePhaseMixture threePhaseProperties(U, phi);
 
@@ -116,11 +116,6 @@
     );
 
 
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
     // Construct interface from alpha distribution
     threePhaseInterfaceProperties interface(threePhaseProperties);
 
@@ -130,3 +125,43 @@
     (
         incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
     );
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    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,
+        mesh.solutionDict().subDict("PISO"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }

applications/solvers/multiphase/interMixingFoam/interMixingFoam.C

@@ -62,6 +62,7 @@ int main(int argc, char *argv[])
         #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;

applications/solvers/multiphase/interPhaseChangeFoam/Make/files

@@ -1,6 +1,6 @@
 interPhaseChangeFoam.C
 phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
-phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
+phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/newPhaseChangeTwoPhaseMixture.C
 phaseChangeTwoPhaseMixtures/Kunz/Kunz.C
 phaseChangeTwoPhaseMixtures/Merkle/Merkle.C
 phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C

applications/solvers/multiphase/interPhaseChangeFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -linterfaceProperties \
+    -ltwoPhaseInterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -25,10 +25,10 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - fvc::snGrad(p)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -50,18 +50,18 @@
           + vDotcAlphal
         );
 
-        // MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+        //MULES::explicitSolve
-        // MULES::explicitSolve
+        //(
-        // (
+        //    geometricOneField(),
-        //     geometricOneField(),
+        //    alpha1,
-        //     alpha1,
+        //    phi,
-        //     phi,
+        //    phiAlpha,
-        //     phiAlpha,
+        //    Sp,
-        //     Sp,
+        //    Su,
-        //     Su,
+        //    1,
-        //     1,
+        //    0
-        //     0
+        //);
-        // );
+
         MULES::implicitSolve
         (
             geometricOneField(),

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -36,12 +36,12 @@ surfaceScalarField rhoPhi
             !(++alphaSubCycle).end();
         )
         {
-#           include "alphaEqn.H"
+            #include "alphaEqn.H"
         }
     }
     else
     {
-#       include "alphaEqn.H"
+        #include "alphaEqn.H"
     }
 
     if (nOuterCorr == 1)

applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H

@@ -1,54 +0,0 @@
-{
-    #include "continuityErrs.H"
-
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll(p.boundaryField(), i)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
-    volScalarField pcorr
-    (
-        IOobject
-        (
-            "pcorr",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("pcorr", p.dimensions(), 0.0),
-        pcorrTypes
-    );
-
-    dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
-    adjustPhi(phi, U, pcorr);
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix pcorrEqn
-        (
-            fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
-        );
-
-        pcorrEqn.setReference(pRefCell, pRefValue);
-        pcorrEqn.solve();
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi -= pcorrEqn.flux();
-        }
-    }
-
-#   include "continuityErrs.H"
-}

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -40,7 +40,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
     Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
     autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
@@ -65,12 +65,6 @@
     );
     rho.oldTime();
 
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties());
 
@@ -79,3 +73,43 @@
     (
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
     );
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    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,
+        mesh.solutionDict().subDict("PISO"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
-    #include "correctPhi.H"
+    #include "../interFoam/correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -82,9 +82,11 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        // --- Outer-corrector loop
+        // --- Pressure-velocity PIMPLE corrector loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            bool finalIter = oCorr == nOuterCorr-1;
+
             #include "UEqn.H"
 
             // --- PISO loop
@@ -92,8 +94,6 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-
-            #include "continuityErrs.H"
         }
 
         twoPhaseProperties->correct();

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -11,14 +11,13 @@
       + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p);
+    adjustPhi(phiU, U, p_rgh);
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-           *fvc::snGrad(alpha1)*mesh.magSf()
+          - ghf*fvc::snGrad(rho)
-          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf*mesh.magSf();
-        )*rUAf;
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -26,29 +25,48 @@
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvc::div(phi) - fvm::laplacian(rUAf, p)
+            fvc::div(phi) - fvm::laplacian(rUAf, p_rgh)
-          - (vDotvP - vDotcP)*pSat + fvm::Sp(vDotvP - vDotcP, p)
+          - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        p_rghEqn.setReference(pRefCell, pRefValue);
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        p_rghEqn.solve
-        {
+        (
-            pEqn.solve(mesh.solver(p.name() + "Final"));
+            mesh.solver
-        }
+            (
-        else
+                p_rgh.select
-        {
+                (
-            pEqn.solve(mesh.solver(p.name()));
+                    finalIter
-        }
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pEqn.flux();
+            phi += p_rghEqn.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
 }

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.H

@@ -9,16 +9,16 @@ License
     This file is part of OpenFOAM.
 
     OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU Generac License as published by
+    under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either  2 of the License, or
+    the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
 
-    OpenFOAM is distributed in the ho it will be useful, but WITHOUT
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the imarranty of MERCHANTABILITY or
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.he GNU General Public License
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
 
-    You should have received a copy oNU General Public License
+    You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.H

@@ -9,16 +9,16 @@ License
     This file is part of OpenFOAM.
 
     OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU Generac License as published by
+    under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either  2 of the License, or
+    the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
 
-    OpenFOAM is distributed in the ho it will be useful, but WITHOUT
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the imarranty of MERCHANTABILITY or
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.he GNU General Public License
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
 
-    You should have received a copy oNU General Public License
+    You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.H

@@ -1,5 +1,5 @@
 /*---------------------------------------------------------------------------*\
-  ========Merkle=                 |
+  =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
@@ -9,16 +9,16 @@ License
     This file is part of OpenFOAM.
 
     OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU Generac License as published by
+    under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either  2 of the License, or
+    the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
 
-    OpenFOAM is distributed in the ho it will be useful, but WITHOUT
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the imarranty of MERCHANTABILITY or
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.he GNU General Public License
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
 
-    You should have received a copy oNU General Public License
+    You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/newPhaseChangeTwoPhaseMixture.C

@@ -36,27 +36,30 @@ Foam::phaseChangeTwoPhaseMixture::New
     const word& alpha1Name
 )
 {
-    // get model name, but do not register the dictionary
+    IOdictionary transportPropertiesDict
-    const word mixtureType
     (
-        IOdictionary
+        IOobject
         (
-            IOobject
+            "transportProperties",
-            (
+            U.time().constant(),
-                "transportProperties",
+            U.db(),
-                U.time().constant(),
+            IOobject::MUST_READ,
-                U.db(),
+            IOobject::NO_WRITE,
-                IOobject::MUST_READ_IF_MODIFIED,
+            false
-                IOobject::NO_WRITE,
+        )
-                false
-            )
-        ).lookup("phaseChangeTwoPhaseMixture")
     );
 
-    Info<< "Selecting phaseChange model " << mixtureType << endl;
+    word phaseChangeTwoPhaseMixtureTypeName
+    (
+        transportPropertiesDict.lookup("phaseChangeTwoPhaseMixture")
+    );
+
+    Info<< "Selecting phaseChange model "
+        << phaseChangeTwoPhaseMixtureTypeName << endl;
 
     componentsConstructorTable::iterator cstrIter =
-        componentsConstructorTablePtr_->find(mixtureType);
+        componentsConstructorTablePtr_
+            ->find(phaseChangeTwoPhaseMixtureTypeName);
 
     if (cstrIter == componentsConstructorTablePtr_->end())
     {
@@ -64,8 +67,8 @@ Foam::phaseChangeTwoPhaseMixture::New
         (
             "phaseChangeTwoPhaseMixture::New"
         )   << "Unknown phaseChangeTwoPhaseMixture type "
-            << mixtureType << nl << nl
+            << phaseChangeTwoPhaseMixtureTypeName << endl << endl
-            << "Valid phaseChangeTwoPhaseMixture types are : " << endl
+            << "Valid  phaseChangeTwoPhaseMixtures are : " << endl
             << componentsConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H

@@ -28,7 +28,7 @@ Description
 
 SourceFiles
     phaseChangeTwoPhaseMixture.C
-    phaseChangeModelNew.C
+    newPhaseChangeModel.C
 
 \*---------------------------------------------------------------------------*/
 

applications/solvers/multiphase/multiphaseInterFoam/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean
+wclean MRFMultiphaseInterFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/multiphaseInterFoam/Allwmake

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso multiphaseMixture
+wmake
+wmake MRFMultiphaseInterFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    MRFMultiphaseInterFoam
+
+Description
+    Solver for n incompressible fluids which captures the interfaces and
+    includes surface-tension and contact-angle effects for each phase.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixture.H"
+#include "turbulenceModel.H"
+#include "MRFZones.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createMRFZones.H"
+    #include "readTimeControls.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        mixture.solve();
+        rho = mixture.rho();
+        #include "zonePhaseVolumes.H"
+
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/Make/files

@@ -0,0 +1,3 @@
+MRFMultiphaseInterFoam.C
+
+EXE = $(FOAM_APPBIN)/MRFMultiphaseInterFoam

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/Make/options

@@ -0,0 +1,19 @@
+EXE_INC = \
+    -I.. \
+    -I../../interFoam \
+    -I../../interFoam/MRFInterFoam \
+    -I../multiphaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmultiphaseInterFoam \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/UEqn.H

@@ -0,0 +1,35 @@
+    surfaceScalarField muEff
+    (
+        "muEff",
+        mixture.muf()
+      + fvc::interpolate(rho*turbulence->nut())
+    );
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(mixture.rhoPhi(), U)
+      - fvm::laplacian(muEff, U)
+      - (fvc::grad(U) & fvc::grad(muEff))
+    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+    );
+    mrfZones.addCoriolis(rho, UEqn);
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+    }

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H

@@ -0,0 +1,63 @@
+{
+    volScalarField rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
+
+    U = rAU*UEqn.H();
+
+    surfaceScalarField phiU
+    (
+        "phiU",
+        (fvc::interpolate(U) & mesh.Sf())
+      //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
+    );
+    mrfZones.relativeFlux(phiU);
+
+    adjustPhi(phiU, U, p_rgh);
+
+    phi = phiU +
+    (
+        mixture.surfaceTensionForce()
+      - ghf*fvc::snGrad(rho)
+    )*rAUf*mesh.magSf();
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve
+        (
+            mesh.solver
+            (
+                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
+            )
+        );
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= p_rghEqn.flux();
+        }
+    }
+
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/zonePhaseVolumes.H

@@ -0,0 +1,26 @@
+{
+    const scalarField& V = mesh.V();
+
+    forAll(mesh.cellZones(), czi)
+    {
+        const labelList& cellLabels = mesh.cellZones()[czi];
+
+        forAllConstIter(PtrDictionary<phase>, mixture.phases(), iter)
+        {
+            const volScalarField& alpha = iter();
+            scalar phaseVolume = 0;
+
+            forAll(cellLabels, cli)
+            {
+                label celli = cellLabels[cli];
+                phaseVolume += alpha[celli]*V[celli];
+            }
+
+            reduce(phaseVolume, sumOp<scalar>());
+
+            Info<< alpha.name()
+                << " phase volume in zone " << mesh.cellZones()[czi].name()
+                << " = " << phaseVolume*1e6 << " ml " << endl;
+        }
+    }
+}

applications/solvers/multiphase/multiphaseInterFoam/Make/files

@@ -1,6 +1,3 @@
-multiphaseMixture/phase/phase.C
-multiphaseMixture/multiphaseAlphaContactAngle/multiphaseAlphaContactAngleFvPatchScalarField.C
-multiphaseMixture/multiphaseMixture.C
 multiphaseInterFoam.C
 
 EXE = $(FOAM_APPBIN)/multiphaseInterFoam

applications/solvers/multiphase/multiphaseInterFoam/Make/options

@@ -1,14 +1,14 @@
 EXE_INC = \
     -I../interFoam \
-    -ImultiphaseMixture \
+    -ImultiphaseMixture/lnInclude \
-    -ImultiphaseMixture/phase \
-    -ImultiphaseMixture/multiphaseAlphaContactAngle \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
+    -lmultiphaseInterFoam \
     -linterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \

applications/solvers/multiphase/multiphaseInterFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
+                (
-              + (
                     mixture.surfaceTensionForce()
-                  - fvc::snGrad(p)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/multiphaseInterFoam/alphaCourantNo.H

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    surfaceScalarField SfUfbyDelta =
+        mixture.nearInterface()
+       *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
+
+    alphaCoNum = max(SfUfbyDelta/mesh.magSf())
+        .value()*runTime.deltaT().value();
+
+    meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
+        .value()*runTime.deltaT().value();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseInterFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field p\n" << endl;
+    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -45,13 +45,48 @@
     rho.oldTime();
 
 
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
     // Construct incompressible turbulence model
     autoPtr<incompressible::turbulenceModel> turbulence
     (
         incompressible::turbulenceModel::New(U, phi, mixture)
     );
+
+    #include "readGravitationalAcceleration.H"
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    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,
+        mesh.solutionDict().subDict("PISO"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C

@@ -43,7 +43,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -61,13 +60,14 @@ int main(int argc, char *argv[])
         #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "alphaCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        mixture.correct();
+        mixture.solve();
         rho = mixture.rho();
 
         #include "UEqn.H"
@@ -78,18 +78,16 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        #include "continuityErrs.H"
+        turbulence->correct();
-
-        //turbulence->correct();
 
         runTime.write();
 
-        Info<< "ExecutionTime = "
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << runTime.elapsedCpuTime()
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << " s\n" << endl << endl;
+            << nl << endl;
     }
 
-    Info<< "\n end \n";
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/Make/files

@@ -0,0 +1,5 @@
+phase/phase.C
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseMixture.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseInterFoam

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -IalphaContactAngle \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lfiniteVolume

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "multiphaseAlphaContactAngleFvPatchScalarField.H"
+#include "alphaContactAngleFvPatchScalarField.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 
@@ -32,8 +32,7 @@ License
 namespace Foam
 {
 
-multiphaseAlphaContactAngleFvPatchScalarField::interfaceThetaProps::
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
-interfaceThetaProps
 (
     Istream& is
 )
@@ -48,10 +47,10 @@ interfaceThetaProps
 Istream& operator>>
 (
     Istream& is,
-    multiphaseAlphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
 )
 {
-    is  >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
     return is;
 }
 
@@ -59,7 +58,7 @@ Istream& operator>>
 Ostream& operator<<
 (
     Ostream& os,
-    const multiphaseAlphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
 )
 {
     os  << tp.theta0_ << token::SPACE
@@ -73,8 +72,7 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-multiphaseAlphaContactAngleFvPatchScalarField::
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
-multiphaseAlphaContactAngleFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF
@@ -84,10 +82,9 @@ multiphaseAlphaContactAngleFvPatchScalarField
 {}
 
 
-multiphaseAlphaContactAngleFvPatchScalarField::
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
-multiphaseAlphaContactAngleFvPatchScalarField
 (
-    const multiphaseAlphaContactAngleFvPatchScalarField& gcpsf,
+    const alphaContactAngleFvPatchScalarField& gcpsf,
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
     const fvPatchFieldMapper& mapper
@@ -98,8 +95,7 @@ multiphaseAlphaContactAngleFvPatchScalarField
 {}
 
 
-multiphaseAlphaContactAngleFvPatchScalarField::
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
-multiphaseAlphaContactAngleFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
@@ -113,10 +109,9 @@ multiphaseAlphaContactAngleFvPatchScalarField
 }
 
 
-multiphaseAlphaContactAngleFvPatchScalarField::
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
-multiphaseAlphaContactAngleFvPatchScalarField
 (
-    const multiphaseAlphaContactAngleFvPatchScalarField& gcpsf,
+    const alphaContactAngleFvPatchScalarField& gcpsf,
     const DimensionedField<scalar, volMesh>& iF
 )
 :
@@ -127,7 +122,7 @@ multiphaseAlphaContactAngleFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void multiphaseAlphaContactAngleFvPatchScalarField::write(Ostream& os) const
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
     os.writeKeyword("thetaProperties")
@@ -138,11 +133,7 @@ void multiphaseAlphaContactAngleFvPatchScalarField::write(Ostream& os) const
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField
+makePatchTypeField(fvPatchScalarField, alphaContactAngleFvPatchScalarField);
-(
-    fvPatchScalarField,
-    multiphaseAlphaContactAngleFvPatchScalarField
-);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -22,19 +22,19 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::multiphaseAlphaContactAngleFvPatchScalarField
+    Foam::alphaContactAngleFvPatchScalarField
 
 Description
     Contact-angle boundary condition for multi-phase interface-capturing
     simulations.  Used in conjuction with multiphaseMixture.
 
 SourceFiles
-    multiphaseAlphaContactAngleFvPatchScalarField.C
+    alphaContactAngleFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef multiphaseAlphaContactAngleFvPatchScalarField_H
+#ifndef alphaContactAngleFvPatchScalarField_H
-#define multiphaseAlphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
 
 #include "zeroGradientFvPatchFields.H"
 #include "multiphaseMixture.H"
@@ -45,10 +45,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class multiphaseAlphaContactAngleFvPatch Declaration
+                           Class alphaContactAngleFvPatch Declaration
 \*---------------------------------------------------------------------------*/
 
-class multiphaseAlphaContactAngleFvPatchScalarField
+class alphaContactAngleFvPatchScalarField
 :
     public zeroGradientFvPatchScalarField
 {
@@ -132,32 +132,31 @@ private:
 public:
 
     //- Runtime type information
-    TypeName("multiphaseAlphaContactAngle");
+    TypeName("alphaContactAngle");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        multiphaseAlphaContactAngleFvPatchScalarField
+        alphaContactAngleFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        multiphaseAlphaContactAngleFvPatchScalarField
+        alphaContactAngleFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const dictionary&
         );
 
-        //- Construct by mapping given
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
-        //  multiphaseAlphaContactAngleFvPatchScalarField onto a new
+        //  onto a new patch
-        //  patch
+        alphaContactAngleFvPatchScalarField
-        multiphaseAlphaContactAngleFvPatchScalarField
         (
-            const multiphaseAlphaContactAngleFvPatchScalarField&,
+            const alphaContactAngleFvPatchScalarField&,
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const fvPatchFieldMapper&
@@ -168,14 +167,14 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new multiphaseAlphaContactAngleFvPatchScalarField(*this)
+                new alphaContactAngleFvPatchScalarField(*this)
             );
         }
 
         //- Construct as copy setting internal field reference
-        multiphaseAlphaContactAngleFvPatchScalarField
+        alphaContactAngleFvPatchScalarField
         (
-            const multiphaseAlphaContactAngleFvPatchScalarField&,
+            const alphaContactAngleFvPatchScalarField&,
             const DimensionedField<scalar, volMesh>&
         );
 
@@ -187,7 +186,7 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new multiphaseAlphaContactAngleFvPatchScalarField(*this, iF)
+                new alphaContactAngleFvPatchScalarField(*this, iF)
             );
         }
 

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -24,11 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "multiphaseMixture.H"
-#include "multiphaseAlphaContactAngleFvPatchScalarField.H"
+#include "alphaContactAngleFvPatchScalarField.H"
 #include "Time.H"
 #include "subCycle.H"
 #include "fvCFD.H"
-#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
 
@@ -237,7 +236,7 @@ Foam::multiphaseMixture::surfaceTensionForce() const
 }
 
 
-void Foam::multiphaseMixture::correct()
+void Foam::multiphaseMixture::solve()
 {
     forAllIter(PtrDictionary<phase>, phases_, iter)
     {
@@ -296,6 +295,10 @@ void Foam::multiphaseMixture::correct()
 }
 
 
+void Foam::multiphaseMixture::correct()
+{}
+
+
 Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
 (
     const volScalarField& alpha1,
@@ -351,11 +354,10 @@ void Foam::multiphaseMixture::correctContactAngle
 
     forAll(boundary, patchi)
     {
-        if (isA<multiphaseAlphaContactAngleFvPatchScalarField>(gbf[patchi]))
+        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
         {
-            const multiphaseAlphaContactAngleFvPatchScalarField& acap =
+            const alphaContactAngleFvPatchScalarField& acap =
-                refCast<const multiphaseAlphaContactAngleFvPatchScalarField>
+                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
-                (gbf[patchi]);
 
             vectorField& nHatPatch = nHatb[patchi];
 
@@ -363,7 +365,7 @@ void Foam::multiphaseMixture::correctContactAngle
                 mesh_.Sf().boundaryField()[patchi]
                /mesh_.magSf().boundaryField()[patchi];
 
-            multiphaseAlphaContactAngleFvPatchScalarField::thetaPropsTable::
+            alphaContactAngleFvPatchScalarField::thetaPropsTable::
                 const_iterator tp =
                 acap.thetaProps().find(interfacePair(alpha1, alpha2));
 
@@ -455,6 +457,34 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 }
 
 
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::nearInterface() const
+{
+    tmp<surfaceScalarField> tnearInt
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "nearInterface",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("nearInterface", dimless, 0.0)
+        )
+    );
+
+    forAllConstIter(PtrDictionary<phase>, phases_, iter)
+    {
+        surfaceScalarField alphaf = fvc::interpolate(iter());
+        tnearInt() = max(tnearInt(), pos(alphaf - 0.01)*pos(0.99 - alphaf));
+    }
+
+    return tnearInt;
+}
+
+
 void Foam::multiphaseMixture::solveAlphas
 (
     const label nAlphaCorr,
@@ -466,7 +496,7 @@ void Foam::multiphaseMixture::solveAlphas
     nSolves++;
 
     word alphaScheme("div(phi,alpha)");
-    word alphacScheme("div(phic,alpha)");
+    word alphacScheme("div(phirb,alpha)");
 
     tmp<fv::convectionScheme<scalar> > mvConvection
     (

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H

@@ -164,7 +164,7 @@ private:
         multivariateSurfaceInterpolationScheme<scalar>::fieldTable alphaTable_;
 
 
-    // Private Member Functions
+    // Private member functions
 
         void calcAlphas();
 
@@ -256,6 +256,13 @@ public:
 
         tmp<surfaceScalarField> surfaceTensionForce() const;
 
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<surfaceScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+
         //- Correct the mixture properties
         void correct();
 

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -1,47 +1,62 @@
 {
-    volScalarField rUA = 1.0/UEqn.A();
+    volScalarField rAU = 1.0/UEqn.A();
-    surfaceScalarField rUAf = fvc::interpolate(rUA);
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
 
-    U = rUA*UEqn.H();
+    U = rAU*UEqn.H();
 
     surfaceScalarField phiU
     (
         "phiU",
-        (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p);
+    adjustPhi(phiU, U, p_rgh);
 
     phi = phiU +
-        (
+    (
-            mixture.surfaceTensionForce()*mesh.magSf()
+        mixture.surfaceTensionForce()
-          + fvc::interpolate(rho)*(g & mesh.Sf())
+      - ghf*fvc::snGrad(rho)
-        )*rUAf;
+    )*rAUf*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rUAf, p) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        if (corr == nCorr-1)
+        p_rghEqn.solve
-        {
+        (
-            pEqn.solve(mesh.solver(p.name() + "Final"));
+            mesh.solver
-        }
+            (
-        else
+                p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
-        {
+            )
-            pEqn.solve(mesh.solver(p.name()));
+        );
-        }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pEqn.flux();
+            phi -= p_rghEqn.flux();
         }
     }
 
-    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
     U.correctBoundaryConditions();
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
 }

applications/solvers/multiphase/settlingFoam/UEqn.H

@@ -22,7 +22,10 @@
           ==
             fvc::reconstruct
             (
-                (- ghf*fvc::snGrad(rho) - fvc::snGrad(pmh))*mesh.magSf()
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
             )
         );
     }