Add the OpenFOAM source tree

applications/solvers/multiphase/interFoam/Allwclean

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

applications/solvers/multiphase/interFoam/Allwmake

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

applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
+    interFoam
+
+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.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "fvcSmooth.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createPrghCorrTypes.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialrDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "setrDeltaT.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+
+            #define LTSSOLVE
+            #include "alphaEqnSubCycle.H"
+            #undef LTSSOLVE
+
+            mixture.correct();
+
+            turbulence->correct();
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files

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

applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options

@@ -0,0 +1,22 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/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)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling

applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H

@@ -0,0 +1,31 @@
+scalar maxDeltaT
+(
+    pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
+);
+
+volScalarField rDeltaT
+(
+    IOobject
+    (
+        "rDeltaT",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+    zeroGradientFvPatchScalarField::typeName
+);
+
+volScalarField rSubDeltaT
+(
+    IOobject
+    (
+        "rSubDeltaT",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+);

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -0,0 +1,138 @@
+{
+    const dictionary& pimpleDict = pimple.dict();
+
+    scalar maxCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    scalar maxAlphaCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    label nAlphaSpreadIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    scalar alphaSpreadDiff
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
+    );
+
+    scalar alphaSpreadMax
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
+    );
+
+    scalar alphaSpreadMin
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
+    );
+
+    label nAlphaSweepIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(mag(rhoPhi))().dimensionedInternalField()
+       /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+
+        volScalarField alpha1Bar(fvc::average(alpha1));
+
+        rDeltaT.dimensionedInternalField() = max
+        (
+            rDeltaT.dimensionedInternalField(),
+            pos(alpha1Bar.dimensionedInternalField() - alphaSpreadMin)
+           *pos(alphaSpreadMax - alpha1Bar.dimensionedInternalField())
+           *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
+           /((2*maxAlphaCo)*mesh.V())
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread
+        (
+            rDeltaT,
+            alpha1,
+            nAlphaSpreadIter,
+            alphaSpreadDiff,
+            alphaSpreadMax,
+            alphaSpreadMin
+        );
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.internalField())
+            << ", " << gMax(1/rDeltaT.internalField()) << endl;
+    }
+
+    #include "alphaControls.H"
+
+    rSubDeltaT = rDeltaT*nAlphaSubCycles;
+}

applications/solvers/multiphase/interFoam/Make/files

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

applications/solvers/multiphase/interFoam/Make/options

@@ -0,0 +1,21 @@
+EXE_INC = \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/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)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling

applications/solvers/multiphase/interFoam/UEqn.H

@@ -0,0 +1,31 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(rhoPhi, U)
+      + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }

applications/solvers/multiphase/interFoam/alphaCourantNo.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
+    alphaCourantNo
+
+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())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().internalField()
+       *fvc::surfaceSum(mag(phi))().internalField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -0,0 +1,150 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // Standard face-flux compression coefficient
+    surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
+
+    // Add the optional isotropic compression contribution
+    if (icAlpha > 0)
+    {
+        phic *= (1.0 - icAlpha);
+        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
+    }
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalarField& phicp = phic.boundaryField()[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    tmp<surfaceScalarField> tphiAlpha;
+
+    if (MULESCorr)
+    {
+        fvScalarMatrix alpha1Eqn
+        (
+            #ifdef LTSSOLVE
+            fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
+            #else
+            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+            #endif
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phi,
+                upwind<scalar>(mesh, phi)
+            ).fvmDiv(phi, alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+
+        tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
+        tphiAlpha = tmp<surfaceScalarField>
+        (
+            new surfaceScalarField(tphiAlphaUD())
+        );
+
+        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            #ifdef LTSSOLVE
+            MULES::LTScorrect(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
+            #else
+            MULES::correct(alpha1, tphiAlpha(), tphiAlphaCorr0(), 1, 0);
+            #endif
+
+            tphiAlpha() += tphiAlphaCorr0();
+        }
+
+        // Cache the upwind-flux
+        tphiAlphaCorr0 = tphiAlphaUD;
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        surfaceScalarField phir(phic*mixture.nHatf());
+
+        tmp<surfaceScalarField> tphiAlphaUn
+        (
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
+            volScalarField alpha10(alpha1);
+
+            #ifdef LTSSOLVE
+            MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
+            #else
+            MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
+            #endif
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                tphiAlpha() += tphiAlphaCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                tphiAlpha() += 0.5*tphiAlphaCorr();
+            }
+        }
+        else
+        {
+            tphiAlpha = tphiAlphaUn;
+
+            #ifdef LTSSOLVE
+            MULES::explicitLTSSolve(alpha1, phi, tphiAlpha(), 1, 0);
+            #else
+            MULES::explicitSolve(alpha1, phi, tphiAlpha(), 1, 0);
+            #endif
+        }
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        tphiAlphaCorr0 = tphiAlpha() - tphiAlphaCorr0;
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+}

applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H

@@ -0,0 +1,33 @@
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -0,0 +1,40 @@
+{
+    #include "continuityErrs.H"
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p_rgh.dimensions(), 0.0),
+        pcorrTypes
+    );
+
+    dimensionedScalar rAUf("rAUf", dimTime/rho.dimensions(), 1.0);
+
+    adjustPhi(phi, U, pcorr);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
+        );
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+    }
+
+    #include "continuityErrs.H"
+}

applications/solvers/multiphase/interFoam/createFields.H

@@ -0,0 +1,135 @@
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            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;
+    immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
+
+    volScalarField& alpha1(mixture.alpha1());
+    volScalarField& alpha2(mixture.alpha2());
+
+    const dimensionedScalar& rho1 = mixture.rho1();
+    const dimensionedScalar& rho2 = mixture.rho2();
+
+
+    // Need to store rho for ddt(rho, U)
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        alpha1*rho1 + alpha2*rho2,
+        alpha1.boundaryField().types()
+    );
+    rho.oldTime();
+
+
+    // Mass flux
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    );
+
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, mixture)
+    );
+
+    #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("PIMPLE"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+
+    fv::IOoptionList fvOptions(mesh);
+
+
+    // MULES Correction
+    tmp<surfaceScalarField> tphiAlphaCorr0;

applications/solvers/multiphase/interFoam/interDyMFoam/Make/files

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

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

@@ -0,0 +1,27 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/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)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling

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

@@ -0,0 +1,65 @@
+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];
+        }
+    }
+}
+
+{
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p_rgh.dimensions(), 0.0),
+        pcorrTypes
+    );
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    #ifndef divUCorr
+    #define divUCorr
+    #endif
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) divUCorr
+        );
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+    }
+
+    #undef divUCorr
+
+    #include "continuityErrs.H"
+}

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

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
+    interDyMFoam
+
+Description
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach,
+    with optional mesh motion and mesh topology changes including adaptive
+    re-meshing.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "CMULES.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "createPrghCorrTypes.H"
+
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
+    );
+
+    #include "correctPhi.H"
+    #include "createUf.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readControls.H"
+        #include "alphaCourantNo.H"
+        #include "CourantNo.H"
+
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                    gh = g & mesh.C();
+                    ghf = g & mesh.Cf();
+                }
+
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
+
+                    #include "correctPhi.H"
+
+                    // Make the flux relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+
+                    mixture.correct();
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                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/interDyMFoam/pEqn.H

@@ -0,0 +1,88 @@
+{
+    rAU = 1.0/UEqn.A();
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
+    );
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
+    surfaceScalarField phig
+    (
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
+    }
+
+    // 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,16 @@
+#include "readTimeControls.H"
+
+bool correctPhi
+(
+    pimple.dict().lookupOrDefault<Switch>("correctPhi", true)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
+);
+
+bool moveMeshOuterCorrectors
+(
+    pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
+);

applications/solvers/multiphase/interFoam/interFoam.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
+    interFoam
+
+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.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "createPrghCorrTypes.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                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/interMixingFoam/Make/files

@@ -0,0 +1,6 @@
+incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.C
+threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.C
+interMixingFoam.C
+
+EXE = $(FOAM_APPBIN)/interMixingFoam

applications/solvers/multiphase/interFoam/interMixingFoam/Make/options

@@ -0,0 +1,26 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -IimmiscibleIncompressibleThreePhaseMixture \
+    -IincompressibleThreePhaseMixture \
+    -IthreePhaseInterfaceProperties \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseProperties/alphaContactAngle/alphaContactAngle \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -ltwoPhaseMixture \
+    -ltwoPhaseProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqns.H

@@ -0,0 +1,175 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phir
+    (
+        IOobject
+        (
+            "phir",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mixture.cAlpha()*mag(phi/mesh.magSf())*mixture.nHatf()
+    );
+
+    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    {
+        // Create the limiter to be used for all phase-fractions
+        scalarField allLambda(mesh.nFaces(), 1.0);
+
+        // Split the limiter into a surfaceScalarField
+        slicedSurfaceScalarField lambda
+        (
+            IOobject
+            (
+                "lambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimless,
+            allLambda,
+            false   // Use slices for the couples
+        );
+
+
+        // Create the complete convection flux for alpha1
+        surfaceScalarField phiAlpha1
+        (
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha3, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+        );
+
+        // Create the bounded (upwind) flux for alpha1
+        surfaceScalarField phiAlpha1BD
+        (
+            upwind<scalar>(mesh, phi).flux(alpha1)
+        );
+
+        // Calculate the flux correction for alpha1
+        phiAlpha1 -= phiAlpha1BD;
+
+        // Calculate the limiter for alpha1
+        MULES::limiter
+        (
+            allLambda,
+            1.0/runTime.deltaT().value(),
+            geometricOneField(),
+            alpha1,
+            phiAlpha1BD,
+            phiAlpha1,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Create the complete flux for alpha2
+        surfaceScalarField phiAlpha2
+        (
+            fvc::flux
+            (
+                phi,
+                alpha2,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(phir, alpha1, alpharScheme),
+                alpha2,
+                alpharScheme
+            )
+        );
+
+        // Create the bounded (upwind) flux for alpha2
+        surfaceScalarField phiAlpha2BD
+        (
+            upwind<scalar>(mesh, phi).flux(alpha2)
+        );
+
+        // Calculate the flux correction for alpha2
+        phiAlpha2 -= phiAlpha2BD;
+
+        // Further limit the limiter for alpha2
+        MULES::limiter
+        (
+            allLambda,
+            1.0/runTime.deltaT().value(),
+            geometricOneField(),
+            alpha2,
+            phiAlpha2BD,
+            phiAlpha2,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Construct the limited fluxes
+        phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
+        phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
+
+        // Solve for alpha1
+        solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
+
+        // Create the diffusion coefficients for alpha2<->alpha3
+        volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2));
+        volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3));
+
+        // Add the diffusive flux for alpha3->alpha2
+        phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
+
+        // Solve for alpha2
+        fvScalarMatrix alpha2Eqn
+        (
+            fvm::ddt(alpha2)
+          + fvc::div(phiAlpha2)
+          - fvm::laplacian(Dc23 + Dc32, alpha2)
+        );
+        alpha2Eqn.solve();
+
+        // Construct the complete mass flux
+        rhoPhi =
+              phiAlpha1*(rho1 - rho3)
+            + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
+            + phi*rho3;
+
+        alpha3 = 1.0 - alpha1 - alpha2;
+    }
+
+    Info<< "Air phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha2.weightedAverage(mesh.V()).value()
+        << "  Min(alpha2) = " << min(alpha2).value()
+        << "  Max(alpha2) = " << max(alpha2).value()
+        << endl;
+}

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H

@@ -0,0 +1,41 @@
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqns.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqns.H"
+}
+
+{
+    volScalarField rhoNew(alpha1*rho1 + alpha2*rho2 + alpha3*rho3);
+
+    //solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+    //Info<< "density error = "
+    //    << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl;
+
+    rho == rhoNew;
+}

applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H

@@ -0,0 +1,122 @@
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            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"
+
+    immiscibleIncompressibleThreePhaseMixture mixture(U, phi);
+
+    volScalarField& alpha1(mixture.alpha1());
+    volScalarField& alpha2(mixture.alpha2());
+    volScalarField& alpha3(mixture.alpha3());
+
+    const dimensionedScalar& rho1 = mixture.rho1();
+    const dimensionedScalar& rho2 = mixture.rho2();
+    const dimensionedScalar& rho3 = mixture.rho3();
+
+    dimensionedScalar D23(mixture.lookup("D23"));
+
+    // Need to store rho for ddt(rho, U)
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        alpha1*rho1 + alpha2*rho2 + alpha3*rho3,
+        alpha1.boundaryField().types()
+    );
+    rho.oldTime();
+
+
+    // Mass flux
+    // Initialisation does not matter because rhoPhi is reset after the
+    // alpha solution before it is used in the U equation.
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        rho1*phi
+    );
+
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, mixture)
+    );
+
+
+    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("PIMPLE"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    fv::IOoptionList fvOptions(mesh);

applications/solvers/multiphase/interFoam/interMixingFoam/immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.C

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "immiscibleIncompressibleThreePhaseMixture.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::immiscibleIncompressibleThreePhaseMixture::
+immiscibleIncompressibleThreePhaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    incompressibleThreePhaseMixture(U, phi),
+    threePhaseInterfaceProperties
+    (
+        static_cast<incompressibleThreePhaseMixture&>(*this)
+    )
+{}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.H

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    Foam::immiscibleIncompressibleThreePhaseMixture
+
+Description
+    An immiscible incompressible two-phase mixture transport model
+
+SourceFiles
+    immiscibleIncompressibleThreePhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immiscibleIncompressibleThreePhaseMixture_H
+#define immiscibleIncompressibleThreePhaseMixture_H
+
+#include "incompressibleThreePhaseMixture.H"
+#include "threePhaseInterfaceProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+          Class immiscibleIncompressibleThreePhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class immiscibleIncompressibleThreePhaseMixture
+:
+    public incompressibleThreePhaseMixture,
+    public threePhaseInterfaceProperties
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        immiscibleIncompressibleThreePhaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~immiscibleIncompressibleThreePhaseMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Correct the transport and interface properties
+        virtual void correct()
+        {
+            incompressibleThreePhaseMixture::correct();
+            threePhaseInterfaceProperties::correct();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.C

@@ -0,0 +1,253 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleThreePhaseMixture.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+//- Calculate and return the laminar viscosity
+void Foam::incompressibleThreePhaseMixture::calcNu()
+{
+    nuModel1_->correct();
+    nuModel2_->correct();
+    nuModel3_->correct();
+
+    // Average kinematic viscosity calculated from dynamic viscosity
+    nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::incompressibleThreePhaseMixture::incompressibleThreePhaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "transportProperties",
+            U.time().constant(),
+            U.db(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    phase1Name_(wordList(lookup("phases"))[0]),
+    phase2Name_(wordList(lookup("phases"))[1]),
+    phase3Name_(wordList(lookup("phases"))[2]),
+
+    alpha1_
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phase1Name_),
+            U.time().timeName(),
+            U.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh()
+    ),
+
+    alpha2_
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phase2Name_),
+            U.time().timeName(),
+            U.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh()
+    ),
+
+    alpha3_
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phase3Name_),
+            U.time().timeName(),
+            U.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh()
+    ),
+
+    U_(U),
+    phi_(phi),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            U.time().timeName(),
+            U.db()
+        ),
+        U.mesh(),
+        dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0),
+        calculatedFvPatchScalarField::typeName
+    ),
+
+    nuModel1_
+    (
+        viscosityModel::New
+        (
+            "nu1",
+            subDict(phase1Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel2_
+    (
+        viscosityModel::New
+        (
+            "nu2",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel3_
+    (
+        viscosityModel::New
+        (
+            "nu3",
+            subDict(phase3Name_),
+            U,
+            phi
+        )
+    ),
+
+    rho1_(nuModel1_->viscosityProperties().lookup("rho")),
+    rho2_(nuModel2_->viscosityProperties().lookup("rho")),
+    rho3_(nuModel3_->viscosityProperties().lookup("rho"))
+{
+    alpha3_ == 1.0 - alpha1_ - alpha2_;
+    calcNu();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::incompressibleThreePhaseMixture::mu() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            "mu",
+            alpha1_*rho1_*nuModel1_->nu()
+          + alpha2_*rho2_*nuModel2_->nu()
+          + alpha3_*rho3_*nuModel3_->nu()
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::incompressibleThreePhaseMixture::muf() const
+{
+    surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
+    surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
+    surfaceScalarField alpha3f(fvc::interpolate(alpha3_));
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "mu",
+            alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+          + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+          + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::incompressibleThreePhaseMixture::nuf() const
+{
+    surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
+    surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
+    surfaceScalarField alpha3f(fvc::interpolate(alpha3_));
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "nu",
+            (
+                alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+              + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+              + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+            )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_)
+        )
+    );
+}
+
+
+bool Foam::incompressibleThreePhaseMixture::read()
+{
+    if (transportModel::read())
+    {
+        if
+        (
+            nuModel1_().read(*this)
+         && nuModel2_().read(*this)
+         && nuModel3_().read(*this)
+        )
+        {
+            nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
+            nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
+            nuModel3_->viscosityProperties().lookup("rho") >> rho3_;
+
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.H

@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    incompressibleThreePhaseMixture
+
+Description
+
+SourceFiles
+    incompressibleThreePhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef incompressibleThreePhaseMixture_H
+#define incompressibleThreePhaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "IOdictionary.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class incompressibleThreePhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class incompressibleThreePhaseMixture
+:
+    public IOdictionary,
+    public transportModel
+{
+    // Private data
+
+        word phase1Name_;
+        word phase2Name_;
+        word phase3Name_;
+
+        volScalarField alpha1_;
+        volScalarField alpha2_;
+        volScalarField alpha3_;
+
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        volScalarField nu_;
+
+        autoPtr<viscosityModel> nuModel1_;
+        autoPtr<viscosityModel> nuModel2_;
+        autoPtr<viscosityModel> nuModel3_;
+
+        dimensionedScalar rho1_;
+        dimensionedScalar rho2_;
+        dimensionedScalar rho3_;
+
+
+    // Private Member Functions
+
+        //- Calculate and return the laminar viscosity
+        void calcNu();
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        incompressibleThreePhaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    ~incompressibleThreePhaseMixture()
+    {}
+
+
+    // Member Functions
+
+        const word phase1Name() const
+        {
+            return phase1Name_;
+        }
+
+        const word phase2Name() const
+        {
+            return phase2Name_;
+        }
+
+        const word phase3Name() const
+        {
+            return phase3Name_;
+        }
+
+        const volScalarField& alpha1() const
+        {
+            return alpha1_;
+        }
+
+        volScalarField& alpha1()
+        {
+            return alpha1_;
+        }
+
+        const volScalarField& alpha2() const
+        {
+            return alpha2_;
+        }
+
+        volScalarField& alpha2()
+        {
+            return alpha2_;
+        }
+
+        const volScalarField& alpha3() const
+        {
+            return alpha3_;
+        }
+
+        volScalarField& alpha3()
+        {
+            return alpha3_;
+        }
+
+        //- Return const-access to phase1 density
+        const dimensionedScalar& rho1() const
+        {
+            return rho1_;
+        }
+
+        //- Return const-access to phase2 density
+        const dimensionedScalar& rho2() const
+        {
+            return rho2_;
+        };
+
+        //- Return const-access to phase3 density
+        const dimensionedScalar& rho3() const
+        {
+            return rho3_;
+        };
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        //- Return const-access to phase1 viscosityModel
+        const viscosityModel& nuModel1() const
+        {
+            return nuModel1_();
+        }
+
+        //- Return const-access to phase2 viscosityModel
+        const viscosityModel& nuModel2() const
+        {
+            return nuModel2_();
+        }
+
+        //- Return const-access to phase3 viscosityModel
+        const viscosityModel& nuModel3() const
+        {
+            return nuModel3_();
+        }
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const
+        {
+            return nu_;
+        }
+
+        //- Return the laminar viscosity for patch
+        tmp<scalarField> nu(const label patchi) const
+        {
+            return nu_.boundaryField()[patchi];
+        }
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Correct the laminar viscosity
+        void correct()
+        {
+            calcNu();
+        }
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
+    interMixingFoam
+
+Description
+    Solver for 3 incompressible fluids, two of which are miscible,
+    using a VOF method to capture the interface.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleThreePhaseMixture.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createPrghCorrTypes.H"
+    #include "correctPhi.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "alphaEqnsSubCycle.H"
+
+            mixture.correct();
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        #include "continuityErrs.H"
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "\n end \n";
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C

@@ -0,0 +1,237 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseInterfaceProperties.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "mathematicalConstants.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "fvcSnGrad.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void Foam::threePhaseInterfaceProperties::correctContactAngle
+(
+    surfaceVectorField::GeometricBoundaryField& nHatb
+) const
+{
+    const volScalarField::GeometricBoundaryField& alpha1 =
+        mixture_.alpha1().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha2 =
+        mixture_.alpha2().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha3 =
+        mixture_.alpha3().boundaryField();
+    const volVectorField::GeometricBoundaryField& U =
+        mixture_.U().boundaryField();
+
+    const fvMesh& mesh = mixture_.U().mesh();
+    const fvBoundaryMesh& boundary = mesh.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& a2cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha2[patchi]);
+
+            const alphaContactAngleFvPatchScalarField& a3cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha3[patchi]);
+
+            scalarField twoPhaseAlpha2(max(a2cap, scalar(0)));
+            scalarField twoPhaseAlpha3(max(a3cap, scalar(0)));
+
+            scalarField sumTwoPhaseAlpha
+            (
+                twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL
+            );
+
+            twoPhaseAlpha2 /= sumTwoPhaseAlpha;
+            twoPhaseAlpha3 /= sumTwoPhaseAlpha;
+
+            fvsPatchVectorField& nHatp = nHatb[patchi];
+
+            scalarField theta
+            (
+                convertToRad
+              * (
+                   twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
+                 + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
+                )
+            );
+
+            vectorField nf(boundary[patchi].nf());
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12(nHatp & nf);
+
+            scalarField b1(cos(theta));
+
+            scalarField b2(nHatp.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det(1.0 - a12*a12);
+
+            scalarField a((b1 - a12*b2)/det);
+            scalarField b((b2 - a12*b1)/det);
+
+            nHatp = a*nf + b*nHatp;
+
+            nHatp /= (mag(nHatp) + deltaN_.value());
+        }
+    }
+}
+
+
+void Foam::threePhaseInterfaceProperties::calculateK()
+{
+    const volScalarField& alpha1 = mixture_.alpha1();
+
+    const fvMesh& mesh = alpha1.mesh();
+    const surfaceVectorField& Sf = mesh.Sf();
+
+    // Cell gradient of alpha
+    volVectorField gradAlpha(fvc::grad(alpha1));
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
+
+    // Face unit interface normal
+    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
+    correctContactAngle(nHatfv.boundaryField());
+
+    // Face unit interface normal flux
+    nHatf_ = nHatfv & Sf;
+
+    // Simple expression for curvature
+    K_ = -fvc::div(nHatf_);
+
+    // Complex expression for curvature.
+    // Correction is formally zero but numerically non-zero.
+    //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
+    //nHat.boundaryField() = nHatfv.boundaryField();
+    //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
+(
+    const incompressibleThreePhaseMixture& mixture
+)
+:
+    mixture_(mixture),
+    cAlpha_
+    (
+        readScalar
+        (
+            mixture.U().mesh().solverDict
+            (
+                mixture_.alpha1().name()
+            ).lookup("cAlpha")
+        )
+    ),
+    sigma12_(mixture.lookup("sigma12")),
+    sigma13_(mixture.lookup("sigma13")),
+
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
+    ),
+
+    nHatf_
+    (
+        (
+            fvc::interpolate(fvc::grad(mixture.alpha1()))
+           /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
+        ) & mixture.alpha1().mesh().Sf()
+    ),
+
+    K_
+    (
+        IOobject
+        (
+            "K",
+            mixture.alpha1().time().timeName(),
+            mixture.alpha1().mesh()
+        ),
+        -fvc::div(nHatf_)
+    )
+{
+    calculateK();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::threePhaseInterfaceProperties::surfaceTensionForce() const
+{
+    return fvc::interpolate(sigmaK())*fvc::snGrad(mixture_.alpha1());
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::threePhaseInterfaceProperties::nearInterface() const
+{
+    return max
+    (
+        pos(mixture_.alpha1() - 0.01)*pos(0.99 - mixture_.alpha1()),
+        pos(mixture_.alpha2() - 0.01)*pos(0.99 - mixture_.alpha2())
+    );
+}
+
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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/>.
+
+Class
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+SourceFiles
+    threePhaseInterfaceProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseInterfaceProperties_H
+#define threePhaseInterfaceProperties_H
+
+#include "incompressibleThreePhaseMixture.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class threePhaseInterfaceProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseInterfaceProperties
+{
+    // Private data
+
+        const incompressibleThreePhaseMixture& mixture_;
+
+        //- Compression coefficient
+        scalar cAlpha_;
+
+        //- Surface tension 1-2
+        dimensionedScalar sigma12_;
+
+        //- Surface tension 1-3
+        dimensionedScalar sigma13_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        surfaceScalarField nHatf_;
+        volScalarField K_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        threePhaseInterfaceProperties(const threePhaseInterfaceProperties&);
+        void operator=(const threePhaseInterfaceProperties&);
+
+        //- Correction for the boundary condition on the unit normal nHat on
+        //  walls to produce the correct contact dynamic angle
+        //  calculated from the component of U parallel to the wall
+        void correctContactAngle
+        (
+            surfaceVectorField::GeometricBoundaryField& nHat
+        ) const;
+
+        //- Re-calculate the interface curvature
+        void calculateK();
+
+
+public:
+
+    //- Conversion factor for degrees into radians
+    static const scalar convertToRad;
+
+
+    // Constructors
+
+        //- Construct from volume fraction field alpha and IOdictionary
+        threePhaseInterfaceProperties
+        (
+            const incompressibleThreePhaseMixture& mixture
+        );
+
+
+    // Member Functions
+
+        scalar cAlpha() const
+        {
+            return cAlpha_;
+        }
+
+        const dimensionedScalar& deltaN() const
+        {
+            return deltaN_;
+        }
+
+        const surfaceScalarField& nHatf() const
+        {
+            return nHatf_;
+        }
+
+        const volScalarField& K() const
+        {
+            return K_;
+        }
+
+        tmp<volScalarField> sigma() const
+        {
+            volScalarField limitedAlpha2(max(mixture_.alpha2(), scalar(0)));
+            volScalarField limitedAlpha3(max(mixture_.alpha3(), scalar(0)));
+
+            return
+                (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)
+               /(limitedAlpha2 + limitedAlpha3 + SMALL);
+        }
+
+        tmp<volScalarField> sigmaK() const
+        {
+            return sigma()*K_;
+        }
+
+        tmp<surfaceScalarField> surfaceTensionForce() const;
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        void correct()
+        {
+            calculateK();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interFoam/pEqn.H

@@ -0,0 +1,74 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+    adjustPhi(phiHbyA, U, p_rgh);
+    fvOptions.makeRelative(phiHbyA);
+
+    surfaceScalarField phig
+    (
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #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,23 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling

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

@@ -0,0 +1,42 @@
+    surfaceScalarField muEff
+    (
+        "muEff",
+        mixture.muf()
+      + fvc::interpolate(rho*turbulence->nut())
+    );
+
+    // Calculate and cache mu for the porous media
+    volScalarField mu(mixture.mu());
+
+    fvVectorMatrix UEqn
+    (
+        //pZones.ddt(rho, U)
+        fvm::ddt(rho, U)
+      + fvm::div(rhoPhi, U)
+      - fvm::laplacian(muEff, U)
+      - (fvc::grad(U) & fvc::grad(muEff))
+    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+     ==
+        fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+    pZones.addResistance(UEqn);
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+    }

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

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

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

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     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 "CMULES.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "IOporosityModelList.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createPorousZones.H"
+    #include "readTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createPrghCorrTypes.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                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) 2011 OpenFOAM Foundation
+     \\/     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.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //