1)Adding alphaEqn.H and alphaEqnSubCycle.H specialized version for MPPICInterFoam

2)Adapting divU in TEqn.H for compressibleInterDyMFoam and compressibleInterFoam
3)Re-instated sixDoFRigidBodyDisplacement as patch for pointFields. It allows to use a different fvDynamincMesh type
independently of the BC's

applications/solvers/multiphase/MPPICInterFoam/Make/options

@@ -2,9 +2,9 @@ interFoamPath = $(FOAM_SOLVERS)/multiphase/interFoam
 
 EXE_INC =  \
     -I. \
+    -I../VoF \
     -I./IncompressibleTwoPhaseMixtureTurbulenceModels/lnInclude \
     -I$(interFoamPath) \
-    -I../VoF \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H

@@ -0,0 +1,258 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // Set the off-centering coefficient according to ddt scheme
+    scalar ocCoeff = 0;
+    {
+        tmp<fv::ddtScheme<scalar>> tddtAlpha
+        (
+            fv::ddtScheme<scalar>::New
+            (
+                mesh,
+                mesh.ddtScheme("ddt(alpha)")
+            )
+        );
+        const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
+
+        if
+        (
+            isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
+         || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
+        )
+        {
+            ocCoeff = 0;
+        }
+        else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
+        {
+            if (nAlphaSubCycles > 1)
+            {
+                FatalErrorInFunction
+                    << "Sub-cycling is not supported "
+                       "with the CrankNicolson ddt scheme"
+                    << exit(FatalError);
+            }
+
+            if
+            (
+                alphaRestart
+             || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
+            )
+            {
+                ocCoeff =
+                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
+                   .ocCoeff();
+            }
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Only Euler and CrankNicolson ddt schemes are supported"
+                << exit(FatalError);
+        }
+    }
+
+    // Set the time blending factor, 1 for Euler
+    scalar cnCoeff = 1.0/(1.0 + ocCoeff);
+
+    // 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));
+    }
+
+    surfaceScalarField::Boundary& phicBf =
+        phic.boundaryFieldRef();
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalarField& phicp = phicBf[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    tmp<surfaceScalarField> phiCN(alphaPhic);
+
+    // Calculate the Crank-Nicolson off-centred volumetric flux
+    if (ocCoeff > 0)
+    {
+        phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime();
+    }
+
+    if (MULESCorr)
+    {
+        #include "alphaSuSp.H"
+
+        fvScalarMatrix alpha1Eqn
+        (
+            (
+                LTS
+              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1)
+              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+            )
+          + fv::gaussConvectionScheme<scalar>
+            (
+                mesh,
+                phiCN,
+                upwind<scalar>(mesh, phiCN)
+            ).fvmDiv(phiCN, alpha1)
+         - fvm::Sp(fvc::ddt(alphac) + fvc::div(phiCN), alpha1)
+         ==
+            Su + fvm::Sp(Sp + divU, alpha1)
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+            << endl;
+
+        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        alphaPhi = talphaPhiUD();
+
+        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        {
+            Info<< "Applying the previous iteration compression flux" << endl;
+            MULES::correct
+            (
+                alphac,
+                alpha1,
+                alphaPhi,
+                talphaPhiCorr0.ref(),
+                zeroField(), zeroField(),
+                1, 0
+            );
+
+            alphaPhi += talphaPhiCorr0();
+        }
+
+        // Cache the upwind-flux
+        talphaPhiCorr0 = talphaPhiUD;
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        #include "alphaSuSp.H"
+
+        surfaceScalarField phir(phic*mixture.nHatf());
+
+        tmp<surfaceScalarField> talphaPhiUn
+        (
+            fvc::flux
+            (
+                phiCN(),
+                cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
+                alphaScheme
+            )
+          + fvc::flux
+            (
+               -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+        );
+
+        if (MULESCorr)
+        {
+            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            volScalarField alpha10("alpha10", alpha1);
+
+            MULES::correct
+            (
+                alphac,
+                alpha1,
+                talphaPhiUn(),
+                talphaPhiCorr.ref(),
+                Sp,
+                (-Sp*alpha1)(),
+                1,
+                0
+            );
+
+            // Under-relax the correction for all but the 1st corrector
+            if (aCorr == 0)
+            {
+                alphaPhi += talphaPhiCorr();
+            }
+            else
+            {
+                alpha1 = 0.5*alpha1 + 0.5*alpha10;
+                alphaPhi += 0.5*talphaPhiCorr();
+            }
+        }
+        else
+        {
+            alphaPhi = talphaPhiUn;
+
+            MULES::explicitSolve
+            (
+                alphac,
+                alpha1,
+                phiCN,
+                alphaPhi,
+                Sp,
+                (Su + divU*min(alpha1(), scalar(1)))(),
+                1,
+                0
+            );
+        }
+
+        alpha2 = 1.0 - alpha1;
+
+        mixture.correct();
+    }
+
+    if (alphaApplyPrevCorr && MULESCorr)
+    {
+        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
+        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+    }
+    else
+    {
+        talphaPhiCorr0.clear();
+    }
+
+    if
+    (
+        word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::EulerDdtScheme<vector>::typeName
+    )
+    {
+        #include "rhofs.H"
+        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
+    }
+    else
+    {
+        if (ocCoeff > 0)
+        {
+            // Calculate the end-of-time-step alpha flux
+            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+        }
+
+        // Calculate the end-of-time-step mass flux
+        #include "rhofs.H"
+        rhoPhi = alphaPhi*(rho1f - rho2f) + alphaPhic*rho2f;
+    }
+
+    Info<< "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+        << endl;
+}

applications/solvers/multiphase/MPPICInterFoam/alphaEqnSubCycle.H

@@ -0,0 +1,43 @@
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    tmp<volScalarField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
+    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/compressibleInterFoam/TEqn.H

@@ -5,7 +5,7 @@
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
-            fvc::div(fvc::absolute(phi, U), p)
+            divU*p
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(

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

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -167,6 +167,12 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = alpha1*rho1 + alpha2*rho2;
+
+        // Correct p_rgh for consistency with p and the updated densities
+        p_rgh = p - rho*gh;
+        p_rgh.correctBoundaryConditions();
+
         runTime.write();
 
         Info<< "ExecutionTime = "

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) OpenCFD Ltd. 2017
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,6 +113,7 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
+            volScalarField divU(fvc::div(fvc::absolute(phi, U)));
             #include "TEqn.H"
 
             // --- Pressure corrector loop