Merge branch 'feature-MPPICInterFoam' into 'develop'

Feature mppic inter foam

New MPPICInterFoam solver. Add MPPIC cloud to a VOF approach. Particles volume are considered into transport Eq fluxes.
    Solves 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.

Solver:
/applications/solvers/multiphase/MPPICInterFoam
Tutorial:
/tutorials/multiphase/MPPICInterFoam/twoPhasePachuka

See merge request !41

applications/solvers/multiphase/MPPICInterFoam/Allwclean

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

applications/solvers/multiphase/MPPICInterFoam/Allwmake

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

applications/solvers/multiphase/MPPICInterFoam/CompressibleTwoPhaseMixtureTurbulenceModels/CompressibleTwoPhaseMixtureTurbulenceModels.C

@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PhaseCompressibleTurbulenceModel.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
+
+#include "laminar.H"
+#include "turbulentTransportModel.H"
+#include "LESModel.H"
+
+makeBaseTurbulenceModel
+(
+    volScalarField,
+    volScalarField,
+    compressibleTurbulenceModel,
+    PhaseCompressibleTurbulenceModel,
+    immiscibleIncompressibleTwoPhaseMixture
+);
+
+#define makeRASModel(Type)                                                     \
+    makeTemplatedTurbulenceModel                                               \
+    (immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel, RAS, Type)
+
+#define makeLESModel(Type)                                                     \
+    makeTemplatedTurbulenceModel                                               \
+    (immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel, LES, Type)
+
+#include "kEpsilon.H"
+makeRASModel(kEpsilon);
+
+#include "Smagorinsky.H"
+makeLESModel(Smagorinsky);
+
+#include "kEqn.H"
+makeLESModel(kEqn);
+
+#include "kOmega.H"
+makeRASModel(kOmega);
+
+// ************************************************************************* //

applications/solvers/multiphase/MPPICInterFoam/CompressibleTwoPhaseMixtureTurbulenceModels/Make/files

@@ -0,0 +1,3 @@
+CompressibleTwoPhaseMixtureTurbulenceModels.C
+
+LIB = $(FOAM_LIBBIN)/libCompressibleTwoPhaseMixtureTurbulenceModels

applications/solvers/multiphase/MPPICInterFoam/CompressibleTwoPhaseMixtureTurbulenceModels/Make/options

@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude

applications/solvers/multiphase/MPPICInterFoam/MPPICInterFoam.C

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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
+    MPPICInterFoam
+
+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.
+
+    It includes MRF and MPPIC clouds
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "subCycle.H"
+
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "PhaseCompressibleTurbulenceModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+#include "basicKinematicMPPICCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createMRF.H"
+    #include "createFvOptions.H"
+    #include "correctPhi.H"
+
+    turbulence->validate();
+
+
+    #include "readTimeControls.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;
+
+        Info<< "Evolving " << kinematicCloud.name() << endl;
+
+        kinematicCloud.evolve();
+
+        // Update continuous phase volume fraction field
+        alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+        alphac.correctBoundaryConditions();
+
+         Info<< "Continous phase-1 volume fraction = "
+            << alphac.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alphac) = " << min(alphac).value()
+            << "  Max(alphac) = " << max(alphac).value()
+            << endl;
+
+        alphacf = fvc::interpolate(alphac);
+        alphaRhoPhic = alphacf*rhoPhi;
+        alphaPhic = alphacf*phi;
+        alphacRho = alphac*rho;
+
+        fvVectorMatrix cloudSU(kinematicCloud.SU(U));
+        volVectorField cloudVolSUSu
+        (
+            IOobject
+            (
+                "cloudVolSUSu",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedVector
+            (
+                "0",
+                cloudSU.dimensions()/dimVolume,
+                vector::zero
+            ),
+            zeroGradientFvPatchVectorField::typeName
+        );
+
+        cloudVolSUSu.internalField() = -cloudSU.source()/mesh.V();
+        cloudVolSUSu.correctBoundaryConditions();
+
+        cloudSU.source() = vector::zero;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            #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/MPPICInterFoam/Make/files

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

applications/solvers/multiphase/MPPICInterFoam/Make/options

@@ -0,0 +1,42 @@
+interFoamPath = $(FOAM_SOLVERS)/multiphase/interFoam
+
+EXE_INC =  \
+    -I./IncompressibleTwoPhaseMixtureTurbulenceModels/lnInclude \
+    -I$(interFoamPath) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lthermophysicalFunctions \
+    -lspecie \
+    -lincompressibleTransportModels \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lturbulenceModels \
+    -lsampling \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lfvOptions \
+    -lCompressibleTwoPhaseMixtureTurbulenceModels

applications/solvers/multiphase/MPPICInterFoam/UEqn.H

@@ -0,0 +1,50 @@
+MRF.correctBoundaryVelocity(U);
+
+fvVectorMatrix UEqn
+(
+    fvm::ddt(alphacRho, U)
+  + MRF.DDt(alphacRho, U)
+  - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
+  + fvm::div(rhoPhi, U)
+  + turbulence->divDevRhoReff(U)
+  ==
+    fvOptions(rho, U)
+  + cloudSU
+);
+
+UEqn.relax();
+fvOptions.constrain(UEqn);
+
+volScalarField rAUc(1.0/UEqn.A());
+surfaceScalarField rAUcf(fvc::interpolate(rAUc));
+
+
+surfaceScalarField phicForces
+(
+    (fvc::interpolate(rAUc*cloudVolSUSu) & mesh.Sf())
+);
+
+
+if (pimple.momentumPredictor())
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            phicForces/rAUcf
+          +
+            (
+                fvc::interpolate
+                (
+                    mixture.sigmaK()
+                )*fvc::snGrad(alpha1)
+                - ghf*fvc::snGrad(rho)
+                - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        )
+    );
+
+    fvOptions.correct(U);
+}

applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H

@@ -0,0 +1,163 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // Standard face-flux compression coefficient
+    surfaceScalarField phic
+    (
+        mixture.cAlpha()*mag(alphaPhic/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
+        (
+            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1)
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                alphaPhic,
+                upwind<scalar>(mesh, alphaPhic)
+            ).fvmDiv(alphaPhic, alpha1)
+          - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), 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());
+        alphaPhi = tphiAlphaUD();
+
+        if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+
+            MULES::correct
+            (
+                alphac,
+                alpha1,
+                alphaPhi,
+                tphiAlphaCorr0.ref(),
+                zeroField(), zeroField(),
+                1, 0
+            );
+
+            alphaPhi += 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
+            (
+                alphaPhic,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - alphaPhi);
+            volScalarField alpha10("alpha10", alpha1);
+
+            //MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
+
+            MULES::correct
+            (
+                alphac,
+                alpha1,
+                tphiAlphaUn(),
+                tphiAlphaCorr.ref(),
+                zeroField(), zeroField(),
+                1, 0
+            );
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                alphaPhi += tphiAlphaCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                alphaPhi += 0.5*tphiAlphaCorr();
+            }
+        }
+        else
+        {
+            alphaPhi = tphiAlphaUn;
+
+            MULES::explicitSolve
+            (
+                alphac,
+                alpha1,
+                alphaPhic,
+                alphaPhi,
+                zeroField(), zeroField(),
+                1, 0
+            );
+
+        }
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    rhoPhi = alphaPhi*(rho1 - rho2) + alphaPhic*rho2;
+
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        tphiAlphaCorr0 = alphaPhi - 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/MPPICInterFoam/alphaEqnSubCycle.H

@@ -0,0 +1,34 @@
+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;
+mu = mixture.mu();

applications/solvers/multiphase/MPPICInterFoam/continuityErrs.H

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 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
+    continuityErrs
+
+Description
+    Calculates and prints the continuity errors.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    volScalarField contErr(fvc::ddt(alphac) + fvc::div(alphacf*phi));
+
+    scalar sumLocalContErr = runTime.deltaTValue()*
+        mag(contErr)().weightedAverage(mesh.V()).value();
+
+    scalar globalContErr = runTime.deltaTValue()*
+        contErr.weightedAverage(mesh.V()).value();
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
+
+// ************************************************************************* //

applications/solvers/multiphase/MPPICInterFoam/correctPhi.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/multiphase/MPPICInterFoam/createFields.H

@@ -0,0 +1,219 @@
+    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,
+            IOobject::AUTO_WRITE
+        ),
+        alpha1*rho1 + alpha2*rho2
+    );
+    rho.oldTime();
+
+    // Need to store mu as incompressibleTwoPhaseMixture does not store it
+    volScalarField mu
+    (
+        IOobject
+        (
+            "mu",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        mixture.mu(),
+        calculatedFvPatchScalarField::typeName
+    );
+
+
+    // Mass flux
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    );
+
+
+    #include "readGravitationalAcceleration.H"
+
+    #include "readhRef.H"
+    #include "gh.H"
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rho*gh
+    );
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell
+    (
+        p,
+        p_rgh,
+        mesh.solutionDict().subDict("PIMPLE"),
+        pRefCell,
+        pRefValue
+    );
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    mesh.setFluxRequired(p_rgh.name());
+    mesh.setFluxRequired(alpha1.name());
+
+    // MULES flux from previous time-step
+    surfaceScalarField alphaPhi
+    (
+        IOobject
+        (
+            "alphaPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        phi*fvc::interpolate(alpha1)
+    );
+
+
+    tmp<surfaceScalarField> tphiAlphaCorr0;
+
+    // alphac must be constructed before the cloud
+    // so that the drag-models can find it
+    volScalarField alphac
+    (
+        IOobject
+        (
+            "alphac",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimless, 0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+    alphac.oldTime();
+
+    volScalarField alphacRho(alphac*rho);
+    alphacRho.oldTime();
+
+    Info<< "Constructing kinematicCloud " << endl;
+    basicKinematicMPPICCloud kinematicCloud
+    (
+        "kinematicCloud",
+        rho,
+        U,
+        mu,
+        g
+    );
+
+    // Particle fraction upper limit
+    scalar alphacMin
+    (
+        1.0
+      - readScalar
+        (
+            kinematicCloud.particleProperties().subDict("constantProperties")
+           .lookup("alphaMax")
+        )
+    );
+
+    // Update alphac from the particle locations
+    alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
+    alphac.correctBoundaryConditions();
+
+    surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
+
+    // Phase mass flux
+    surfaceScalarField alphaRhoPhic("alphaRhoPhic", alphacf*rhoPhi);
+
+    // Volumetric phase flux
+    surfaceScalarField alphaPhic("alphaPhic", alphacf*phi);
+
+    autoPtr
+    <
+        PhaseCompressibleTurbulenceModel
+        <
+            immiscibleIncompressibleTwoPhaseMixture
+        >
+    >turbulence
+    (
+        PhaseCompressibleTurbulenceModel
+        <
+            immiscibleIncompressibleTwoPhaseMixture
+        >::New
+        (
+            alphac,
+            rho,
+            U,
+            alphaRhoPhic,
+            rhoPhi,
+            mixture
+        )
+    );

applications/solvers/multiphase/MPPICInterFoam/pEqn.H

@@ -0,0 +1,73 @@
+{
+    volVectorField HbyA(constrainHbyA(rAUc*UEqn.H(), U, p_rgh));
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + alphacf*fvc::interpolate(rho*rAUc)*fvc::ddtCorr(U, phi)
+    );
+
+    MRF.makeRelative(phiHbyA);
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    surfaceScalarField phig
+    (
+        phicForces +
+        (
+            fvc::interpolate(mixture.sigmaK())*fvc::snGrad(alpha1)
+          - ghf*fvc::snGrad(rho)
+        )*rAUcf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, U, phiHbyA, rAUcf, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        surfaceScalarField Dp("Dp", alphacf*rAUcf);
+
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(Dp, p_rgh)
+         ==
+            fvc::ddt(alphac) + fvc::div(alphacf*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()/alphacf;
+
+            p_rgh.relax();
+
+            U =
+                HbyA
+              + rAUc*fvc::reconstruct((phig - p_rghEqn.flux()/alphacf)/rAUcf);
+
+            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;
+    }
+}