diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H b/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H
index 1e5fa3de4b..9fdbdedcd7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H
@@ -1,41 +1,40 @@
-volScalarField::Internal Sp
+tmp<volScalarField::Internal> Sp
 (
-    IOobject
+    volScalarField::Internal::New
     (
         "Sp",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar(dgdt.dimensions(), 0)
+        mesh,
+        dimensionedScalar(dgdt.dimensions(), 0)
+    )
 );
 
-volScalarField::Internal Su
+tmp<volScalarField::Internal> Su
 (
-    IOobject
+    volScalarField::Internal::New
     (
         "Su",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar(dgdt.dimensions(), 0)
+        mesh,
+        dimensionedScalar(dgdt.dimensions(), 0)
+    )
 );
 
+volScalarField::Internal& SpRef = Sp.ref();
+volScalarField::Internal& SuRef = Su.ref();
+
 forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        SpRef[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+        SpRef[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
     }
 }
 
-volScalarField::Internal divU
+tmp<volScalarField> divU
 (
     mesh.moving()
   ? fvc::div(phiCN() + mesh.phi())
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index 8fcd8de136..3446558d10 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,11 +77,6 @@ int main(int argc, char *argv[])
     {
         #include "readDyMControls.H"
 
-        // Store divU from the previous mesh so that it can be mapped
-        // and used in correctPhi to ensure the corrected phi has the
-        // same divergence
-        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
-
         if (LTS)
         {
             #include "setRDeltaT.H"
@@ -102,6 +97,20 @@ int main(int argc, char *argv[])
         {
             if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
             {
+                // Store divU from the previous mesh so that it can be mapped
+                // and used in correctPhi to ensure the corrected phi has the
+                // same divergence
+                tmp<volScalarField> divU;
+
+                if (correctPhi)
+                {
+                    divU = volScalarField::New
+                    (
+                        "divU0",
+                        fvc::div(fvc::absolute(phi, U))
+                    );
+                }
+
                 mesh.update();
 
                 if (mesh.changing())
@@ -130,9 +139,9 @@ int main(int argc, char *argv[])
                         #include "meshCourantNo.H"
                     }
                 }
-            }
 
-            divU = fvc::div(fvc::absolute(phi, U));
+                divU.clear();
+            }
 
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterFoam/correctPhi.H
index 5f3dac0830..97720389b7 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/correctPhi.H
@@ -3,8 +3,8 @@ CorrectPhi
     U,
     phi,
     p,
-    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
-    divU,
+    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+    divU(),
     pimple,
     true
 );
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFieldRefs.H b/applications/solvers/multiphase/compressibleInterFoam/createFieldRefs.H
index 7d6360ce52..44f126f853 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFieldRefs.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFieldRefs.H
@@ -7,3 +7,5 @@ volScalarField& p = mixture.p();
 volScalarField& T = mixture.T();
 const volScalarField& psi1 = mixture.thermo1().psi();
 const volScalarField& psi2 = mixture.thermo2().psi();
+
+tmp<volScalarField> rAU;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index a9535a3e21..b568e10816 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -1,12 +1,20 @@
 {
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    if (correctPhi)
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+    volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         fvc::flux(HbyA)
-      + MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi, Uf))
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
     );
     MRF.makeRelative(phiHbyA);
 
@@ -138,7 +146,7 @@
             phi = phiHbyA + p_rghEqnIncomp.flux();
 
             U = HbyA
-              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+              + rAU()*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/interFoam/Make/options b/applications/solvers/multiphase/interFoam/Make/options
index a82a8a5682..126aa11a65 100644
--- a/applications/solvers/multiphase/interFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/Make/options
@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/twoPhaseChange/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -13,6 +14,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -limmiscibleIncompressibleTwoPhaseMixture \
+    -ltwoPhaseChangeModels \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lfiniteVolume \
diff --git a/applications/solvers/multiphase/interFoam/UEqn.H b/applications/solvers/multiphase/interFoam/UEqn.H
index d6ce1ed4df..03559daae4 100644
--- a/applications/solvers/multiphase/interFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/UEqn.H
@@ -6,7 +6,8 @@
       + MRF.DDt(rho, U)
       + turbulence->divDevTau(rho, U)
      ==
-        fvOptions(rho, U)
+       phaseChange.SU(rho, rhoPhi, U)
+     + fvOptions(rho, U)
     );
 
     UEqn.relax();
diff --git a/applications/solvers/multiphase/interFoam/alphaSuSp.H b/applications/solvers/multiphase/interFoam/alphaSuSp.H
index 0b8e05e187..feea4e84df 100644
--- a/applications/solvers/multiphase/interFoam/alphaSuSp.H
+++ b/applications/solvers/multiphase/interFoam/alphaSuSp.H
@@ -1,3 +1,19 @@
-zeroField Su;
-zeroField Sp;
-zeroField divU;
+// Phase change alpha1 source
+Pair<tmp<volScalarField::Internal>> phaseChangeS
+(
+    phaseChange.Salpha(alpha1)
+);
+const tmp<volScalarField::Internal>& Su = phaseChangeS[0];
+const tmp<volScalarField::Internal>& Sp = phaseChangeS[1];
+
+tmp<volScalarField> divU;
+
+if (Sp.valid())
+{
+    divU =
+    (
+        mesh.moving()
+      ? fvc::div(phiCN() + mesh.phi())
+      : fvc::div(phiCN())
+    );
+}
diff --git a/applications/solvers/multiphase/interFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/correctPhi.H
index 89d8dd51f2..3c4d320762 100644
--- a/applications/solvers/multiphase/interFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/correctPhi.H
@@ -1,12 +1,28 @@
-CorrectPhi
-(
-    U,
-    phi,
-    p_rgh,
-    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-    geometricZeroField(),
-    pimple,
-    true
-);
+if (divU.valid())
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+        divU(),
+        pimple,
+        true
+    );
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+        geometricZeroField(),
+        pimple,
+        true
+    );
+}
 
 #include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/interFoam/createFieldRefs.H b/applications/solvers/multiphase/interFoam/createFieldRefs.H
new file mode 100644
index 0000000000..c6874b154c
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/createFieldRefs.H
@@ -0,0 +1,2 @@
+twoPhaseChangeModel& phaseChange = phaseChangePtr();
+tmp<volScalarField> rAU;
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 7c0ca261cf..be0b5a6ff9 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -34,6 +34,11 @@ volVectorField U
 Info<< "Reading transportProperties\n" << endl;
 immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
 
+autoPtr<twoPhaseChangeModel> phaseChangePtr
+(
+    twoPhaseChangeModel::New(mixture)
+);
+
 volScalarField& alpha1(mixture.alpha1());
 volScalarField& alpha2(mixture.alpha2());
 
@@ -56,21 +61,6 @@ volScalarField rho
 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::momentumTransportModel> turbulence
 (
diff --git a/applications/solvers/multiphase/interFoam/initCorrectPhi.H b/applications/solvers/multiphase/interFoam/initCorrectPhi.H
index e53175c43d..27a3179ad4 100644
--- a/applications/solvers/multiphase/interFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/interFoam/initCorrectPhi.H
@@ -1,46 +1,49 @@
-tmp<volScalarField> rAU;
-
-if (correctPhi)
+if
+(
+    !runTime.restart()
+ || isType<twoPhaseChangeModels::noPhaseChange>(phaseChange)
+)
 {
-    rAU = new volScalarField
-    (
-        IOobject
+    if (correctPhi)
+    {
+        rAU = new volScalarField
         (
-            "rAU",
-            runTime.timeName(),
+            IOobject
+            (
+                "rAU",
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
             mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar(dimTime/dimDensity, 1)
-    );
+            dimensionedScalar(dimTime/dimDensity, 1)
+        );
 
-    CorrectPhi
-    (
-        U,
-        phi,
-        p_rgh,
-        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-        geometricZeroField(),
-        pimple,
-        false
-    );
-
-    #include "continuityErrs.H"
+        CorrectPhi
+        (
+            U,
+            phi,
+            p_rgh,
+            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+            geometricZeroField(),
+            pimple,
+            false
+        );
+    }
+    else
+    {
+        CorrectPhi
+        (
+            U,
+            phi,
+            p_rgh,
+            dimensionedScalar(dimTime/rho.dimensions(), 1),
+            geometricZeroField(),
+            pimple,
+            false
+        );
+    }
 }
-else
-{
-    CorrectPhi
-    (
-        U,
-        phi,
-        p_rgh,
-        dimensionedScalar(dimTime/rho.dimensions(), 1),
-        geometricZeroField(),
-        pimple,
-        false
-    );
 
-    #include "continuityErrs.H"
-}
+#include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 6c88a77924..cdb4b80ae8 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,7 @@ Description
 #include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "noPhaseChange.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
@@ -58,6 +59,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "createDyMControls.H"
     #include "createFields.H"
+    #include "createFieldRefs.H"
     #include "createAlphaFluxes.H"
     #include "initCorrectPhi.H"
     #include "createUfIfPresent.H"
@@ -97,6 +99,24 @@ int main(int argc, char *argv[])
         {
             if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
             {
+                // Store divU from the previous mesh so that it can be mapped
+                // and used in correctPhi to ensure the corrected phi has the
+                // same divergence
+                tmp<volScalarField> divU;
+
+                if
+                (
+                    correctPhi
+                 && !isType<twoPhaseChangeModels::noPhaseChange>(phaseChange)
+                )
+                {
+                    divU = volScalarField::New
+                    (
+                        "divU0",
+                        fvc::div(fvc::absolute(phi, U))
+                    );
+                }
+
                 mesh.update();
 
                 if (mesh.changing())
@@ -132,8 +152,22 @@ int main(int argc, char *argv[])
                         #include "meshCourantNo.H"
                     }
                 }
+
+                divU.clear();
             }
 
+            surfaceScalarField rhoPhi
+            (
+                IOobject
+                (
+                    "rhoPhi",
+                    runTime.timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar(dimMass/dimTime, 0)
+            );
+
             #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interFoam/interMixingFoam/UEqn.H
similarity index 90%
rename from applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
rename to applications/solvers/multiphase/interFoam/interMixingFoam/UEqn.H
index ac127bc731..d6ce1ed4df 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/UEqn.H
@@ -1,7 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevTau(rho, U)
      ==
         fvOptions(rho, U)
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/interMixingFoam/correctPhi.H
similarity index 85%
rename from applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
rename to applications/solvers/multiphase/interFoam/interMixingFoam/correctPhi.H
index d5e94d0ea4..89d8dd51f2 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/correctPhi.H
@@ -4,7 +4,7 @@ CorrectPhi
     phi,
     p_rgh,
     surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-    divU,
+    geometricZeroField(),
     pimple,
     true
 );
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/initCorrectPhi.H b/applications/solvers/multiphase/interFoam/interMixingFoam/initCorrectPhi.H
similarity index 76%
rename from applications/solvers/multiphase/interPhaseChangeFoam/initCorrectPhi.H
rename to applications/solvers/multiphase/interFoam/interMixingFoam/initCorrectPhi.H
index e9c5a4d87a..e53175c43d 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/initCorrectPhi.H
@@ -13,18 +13,16 @@ if (correctPhi)
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("rAU", dimTime/dimDensity, 1)
+        dimensionedScalar(dimTime/dimDensity, 1)
     );
 
-    volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
-
     CorrectPhi
     (
         U,
         phi,
         p_rgh,
         surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-        divU,
+        geometricZeroField(),
         pimple,
         false
     );
@@ -38,7 +36,7 @@ else
         U,
         phi,
         p_rgh,
-        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+        dimensionedScalar(dimTime/rho.dimensions(), 1),
         geometricZeroField(),
         pimple,
         false
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interMixingFoam/pEqn.H
similarity index 71%
rename from applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
rename to applications/solvers/multiphase/interFoam/interMixingFoam/pEqn.H
index 00572edef7..af72769388 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/pEqn.H
@@ -14,8 +14,9 @@
     (
         "phiHbyA",
         fvc::flux(HbyA)
-      + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
     );
+    MRF.makeRelative(phiHbyA);
 
     if (p_rgh.needReference())
     {
@@ -35,35 +36,33 @@
     phiHbyA += phig;
 
     // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf);
-
-    Pair<tmp<volScalarField>> vDotP = mixture.vDotP();
-    const volScalarField& vDotcP = vDotP[0]();
-    const volScalarField& vDotvP = vDotP[1]();
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-          - (vDotvP - vDotcP)*(mixture.pSat() - rho*gh)
-          + fvm::Sp(vDotvP - vDotcP, p_rgh)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, pRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
         p_rghEqn.solve();
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA + p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU()*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            p_rgh.relax();
+
+            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
     }
 
+    #include "continuityErrs.H"
+
     // Correct Uf if the mesh is moving
     fvc::correctUf(Uf, U, phi);
 
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index af72769388..723b1f0e36 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -38,11 +38,15 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
+    // Cache the phase change pressure source
+    fvScalarMatrix Sp_rgh(phaseChange.Sp_rgh(rho, gh, p_rgh));
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+         == Sp_rgh
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -51,11 +55,11 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi = phiHbyA - p_rghEqn.flux();
+            phi = phiHbyA + p_rghEqn.flux();
 
             p_rgh.relax();
 
-            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU()*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Allwclean b/applications/solvers/multiphase/interPhaseChangeFoam/Allwclean
deleted file mode 100755
index 52d19bc530..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/Allwclean
+++ /dev/null
@@ -1,7 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-wclean libso phaseChangeTwoPhaseMixtures
-wclean
-
-#------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Allwmake b/applications/solvers/multiphase/interPhaseChangeFoam/Allwmake
deleted file mode 100755
index 9655af2bd3..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/Allwmake
+++ /dev/null
@@ -1,10 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Parse arguments for library compilation
-. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-
-wmake $targetType phaseChangeTwoPhaseMixtures
-wmake $targetType
-
-#------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Make/files b/applications/solvers/multiphase/interPhaseChangeFoam/Make/files
deleted file mode 100644
index 4358afef8b..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-interPhaseChangeFoam.C
-
-EXE = $(FOAM_APPBIN)/interPhaseChangeFoam
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options
deleted file mode 100644
index 6a07d632b1..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options
+++ /dev/null
@@ -1,27 +0,0 @@
-EXE_INC = \
-    -I. \
-    -I$(LIB_SRC)/transportModels/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
-    -IphaseChangeTwoPhaseMixtures/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -lphaseChangeTwoPhaseMixtures \
-    -ltwoPhaseMixture \
-    -linterfaceProperties \
-    -ltwoPhaseProperties \
-    -lincompressibleTransportModels \
-    -lmomentumTransportModels \
-    -lincompressibleMomentumTransportModels \
-    -lfiniteVolume \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh \
-    -lfvOptions \
-    -lmeshTools
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H
deleted file mode 100644
index 04afbdd19c..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaControls.H
+++ /dev/null
@@ -1,15 +0,0 @@
-const dictionary& alphaControls = mesh.solverDict(alpha1.name());
-
-const label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
-
-const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
-
-const bool MULESCorr(alphaControls.lookupOrDefault<Switch>("MULESCorr", false));
-
-// Apply the compression correction from the previous iteration
-// Improves efficiency for steady-simulations but can only be applied
-// once the alpha field is reasonably steady, i.e. fully developed
-// const bool alphaApplyPrevCorr
-// (
-//     alphaControls.lookupOrDefault<Switch>("alphaApplyPrevCorr", false)
-// );
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
deleted file mode 100644
index 808728a604..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ /dev/null
@@ -1,113 +0,0 @@
-{
-    #include "alphaScheme.H"
-
-    Pair<tmp<volScalarField>> vDotAlphal = mixture.vDotAlphal();
-    const volScalarField& vDotcAlphal = vDotAlphal[0]();
-    const volScalarField& vDotvAlphal = vDotAlphal[1]();
-    const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
-
-    tmp<surfaceScalarField> talphaPhi;
-
-    if (MULESCorr)
-    {
-        fvScalarMatrix alpha1Eqn
-        (
-            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
-          + fv::gaussConvectionScheme<scalar>
-            (
-                mesh,
-                phi,
-                upwind<scalar>(mesh, phi)
-            ).fvmDiv(phi, alpha1)
-          - fvm::Sp(divU, alpha1)
-         ==
-            fvm::Sp(vDotvmcAlphal, alpha1)
-          + vDotcAlphal
-        );
-
-        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;
-
-        talphaPhi = alpha1Eqn.flux();
-    }
-
-    volScalarField alpha10("alpha10", alpha1);
-
-    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
-    {
-        tmp<surfaceScalarField> talphaPhiCorr
-        (
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                compressionScheme.rewind()
-            )
-        );
-
-        if (MULESCorr)
-        {
-            talphaPhiCorr.ref() -= talphaPhi();
-
-            volScalarField alpha100("alpha100", alpha10);
-            alpha10 = alpha1;
-
-            MULES::correct
-            (
-                geometricOneField(),
-                alpha1,
-                talphaPhi(),
-                talphaPhiCorr.ref(),
-                vDotvmcAlphal,
-                (
-                    divU*(alpha10 - alpha100)
-                  - vDotvmcAlphal*alpha10
-                )(),
-                oneField(),
-                zeroField()
-            );
-
-            // Under-relax the correction for all but the 1st corrector
-            if (aCorr == 0)
-            {
-                talphaPhi.ref() += talphaPhiCorr();
-            }
-            else
-            {
-                alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                talphaPhi.ref() += 0.5*talphaPhiCorr();
-            }
-        }
-        else
-        {
-            MULES::explicitSolve
-            (
-                geometricOneField(),
-                alpha1,
-                phi,
-                talphaPhiCorr.ref(),
-                vDotvmcAlphal,
-                (divU*alpha1 + vDotcAlphal)(),
-                oneField(),
-                zeroField()
-            );
-
-            talphaPhi = talphaPhiCorr;
-        }
-
-        alpha2 = 1.0 - alpha1;
-    }
-
-    rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
-
-    Info<< "Liquid phase volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
-        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
-        << endl;
-}
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
deleted file mode 100644
index eed76709b1..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
+++ /dev/null
@@ -1,33 +0,0 @@
-if (nAlphaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-    surfaceScalarField rhoPhiSum
-    (
-        IOobject
-        (
-            "rhoPhiSum",
-            runTime.timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(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;
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
deleted file mode 100644
index 6f35b4edc3..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ /dev/null
@@ -1,112 +0,0 @@
-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<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
-autoPtr<phaseChangeTwoPhaseMixture> mixturePtr
-(
-    phaseChangeTwoPhaseMixture::New(U, phi)
-);
-
-phaseChangeTwoPhaseMixture& mixture = mixturePtr();
-
-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
-);
-rho.oldTime();
-
-
-// Construct incompressible turbulence model
-autoPtr<incompressible::momentumTransportModel> turbulence
-(
-    incompressible::momentumTransportModel::New(U, phi, mixture)
-);
-
-
-#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,
-    pimple.dict(),
-    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());
-
-#include "createFvOptions.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
deleted file mode 100644
index 493525e6ad..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ /dev/null
@@ -1,175 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 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
-    interPhaseChangeFoam
-
-Description
-    Solver for 2 incompressible, isothermal immiscible fluids with phase-change
-    (e.g. cavitation).  Uses a VOF (volume of fluid) phase-fraction based
-    interface capturing approach, with optional mesh motion and mesh topology
-    changes including adaptive re-meshing.
-
-    The momentum and other fluid properties are of the "mixture" and a
-    single momentum equation is solved.
-
-    The set of phase-change models provided are designed to simulate cavitation
-    but other mechanisms of phase-change are supported within this solver
-    framework.
-
-    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "dynamicFvMesh.H"
-#include "CMULES.H"
-#include "subCycle.H"
-#include "interfaceProperties.H"
-#include "phaseChangeTwoPhaseMixture.H"
-#include "kinematicMomentumTransportModel.H"
-#include "pimpleControl.H"
-#include "fvOptions.H"
-#include "CorrectPhi.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
-    #include "initContinuityErrs.H"
-    #include "createFields.H"
-    #include "initCorrectPhi.H"
-    #include "createUfIfPresent.H"
-
-    turbulence->validate();
-
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readDyMControls.H"
-
-        // Store divU from the previous mesh so that it can be mapped
-        // and used in correctPhi to ensure the corrected phi has the
-        // same divergence
-        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
-
-        #include "CourantNo.H"
-        #include "alphaCourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
-            {
-                mesh.update();
-
-                if (mesh.changing())
-                {
-                    gh = (g & mesh.C()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
-
-                    if (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 (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-            }
-
-            divU = fvc::div(fvc::absolute(phi, U));
-
-            surfaceScalarField rhoPhi
-            (
-                IOobject
-                (
-                    "rhoPhi",
-                    runTime.timeName(),
-                    mesh
-                ),
-                mesh,
-                dimensionedScalar(dimMass/dimTime, 0)
-            );
-
-            #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;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/files b/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/files
deleted file mode 100644
index 08f86c582f..0000000000
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/files
+++ /dev/null
@@ -1,7 +0,0 @@
-phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
-phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
-Kunz/Kunz.C
-Merkle/Merkle.C
-SchnerrSauer/SchnerrSauer.C
-
-LIB = $(FOAM_LIBBIN)/libphaseChangeTwoPhaseMixtures
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/Make/options b/applications/solvers/multiphase/multiphaseInterFoam/Make/options
index 8ddd48a73c..21d4fd9a75 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/Make/options
+++ b/applications/solvers/multiphase/multiphaseInterFoam/Make/options
@@ -1,6 +1,4 @@
 EXE_INC = \
-    -I. \
-    -I../interFoam \
     -ImultiphaseMixture/lnInclude \
     -I$(LIB_SRC)/transportModels/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
new file mode 100644
index 0000000000..d6ce1ed4df
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseInterFoam/UEqn.H
@@ -0,0 +1,33 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevTau(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);
+    }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/correctPhi.H b/applications/solvers/multiphase/multiphaseInterFoam/correctPhi.H
new file mode 100644
index 0000000000..89d8dd51f2
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseInterFoam/correctPhi.H
@@ -0,0 +1,12 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+    geometricZeroField(),
+    pimple,
+    true
+);
+
+#include "continuityErrs.H"
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/initCorrectPhi.H b/applications/solvers/multiphase/multiphaseInterFoam/initCorrectPhi.H
new file mode 100644
index 0000000000..e53175c43d
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseInterFoam/initCorrectPhi.H
@@ -0,0 +1,46 @@
+tmp<volScalarField> rAU;
+
+if (correctPhi)
+{
+    rAU = new volScalarField
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(dimTime/dimDensity, 1)
+    );
+
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+        geometricZeroField(),
+        pimple,
+        false
+    );
+
+    #include "continuityErrs.H"
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        dimensionedScalar(dimTime/rho.dimensions(), 1),
+        geometricZeroField(),
+        pimple,
+        false
+    );
+
+    #include "continuityErrs.H"
+}
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
new file mode 100644
index 0000000000..af72769388
--- /dev/null
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -0,0 +1,89 @@
+{
+    if (correctPhi)
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+    volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
+    );
+    MRF.makeRelative(phiHbyA);
+
+    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 pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+    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();
+
+        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"
+
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi);
+
+    // 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;
+    }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
+}
diff --git a/bin/interPhaseChangeFoam b/bin/interPhaseChangeFoam
new file mode 100755
index 0000000000..bd748206f5
--- /dev/null
+++ b/bin/interPhaseChangeFoam
@@ -0,0 +1,54 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     | Website:  https://openfoam.org
+#   \\  /    A nd           | Copyright (C) 2021 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/>.
+#
+# Script
+#     interPhaseChangeFoam
+#
+# Description
+#     Script to inform the user that interPhaseChangeFoam has been replaced
+#     by the more general interFoam solver.
+#
+#------------------------------------------------------------------------------
+
+cat << EOF
+
+The interPhaseChangeFoam solver has solver has been replaced by the more general
+interFoam solver, which now supports phase-change using the new twoPhaseChange
+models library.
+
+To run with with phase-change create a constant/phaseChangeProperties dictionary
+containing the phase-change model specification, e.g.
+
+    phaseChangeModel SchnerrSauer;
+
+    pSat            2300;   // Saturation pressure
+
+See the following cases for an example converted from interPhaseChangeFoam:
+
+    \$FOAM_TUTORIALS/multiphase/interFoam/laminar/cavitatingBullet
+    \$FOAM_TUTORIALS/multiphase/interFoam/RAS/propeller
+
+EOF
+
+#------------------------------------------------------------------------------
diff --git a/etc/config.sh/bash_completion b/etc/config.sh/bash_completion
index d6151a0ef6..e60442c523 100644
--- a/etc/config.sh/bash_completion
+++ b/etc/config.sh/bash_completion
@@ -1986,31 +1986,6 @@ _interMixingFoam_ ()
 }
 complete -o filenames -o nospace -F _interMixingFoam_ interMixingFoam
 
-_interPhaseChangeFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvOptions -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _interPhaseChangeFoam_ interPhaseChangeFoam
-
 _kivaToFoam_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
diff --git a/src/OpenFOAM/memory/tmp/tmpI.H b/src/OpenFOAM/memory/tmp/tmpI.H
index 13de8a23c1..2a2cd32597 100644
--- a/src/OpenFOAM/memory/tmp/tmpI.H
+++ b/src/OpenFOAM/memory/tmp/tmpI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -371,12 +371,12 @@ inline void Foam::tmp<T>::operator=(const tmp<T>& t)
     {
         type_ = t.type_;
 
-        if (!t.ptr_)
-        {
-            FatalErrorInFunction
-                << "Attempted assignment to a deallocated " << typeName()
-                << abort(FatalError);
-        }
+        // if (!t.ptr_)
+        // {
+        //     FatalErrorInFunction
+        //         << "Attempted assignment to a deallocated " << typeName()
+        //         << abort(FatalError);
+        // }
 
         ptr_ = t.ptr_;
         t.ptr_ = 0;
diff --git a/src/OpenFOAM/memory/tmpNrc/tmpNrcI.H b/src/OpenFOAM/memory/tmpNrc/tmpNrcI.H
index f1d2f2ad80..cff4ccab62 100644
--- a/src/OpenFOAM/memory/tmpNrc/tmpNrcI.H
+++ b/src/OpenFOAM/memory/tmpNrc/tmpNrcI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -281,12 +281,12 @@ inline void Foam::tmpNrc<T>::operator=(const tmpNrc<T>& t)
     {
         type_ = TMP;
 
-        if (!t.ptr_)
-        {
-            FatalErrorInFunction
-                << "Attempted assignment to a deallocated " << typeName()
-                << abort(FatalError);
-        }
+        // if (!t.ptr_)
+        // {
+        //     FatalErrorInFunction
+        //         << "Attempted assignment to a deallocated " << typeName()
+        //         << abort(FatalError);
+        // }
 
         ptr_ = t.ptr_;
         t.ptr_ = 0;
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.C b/src/finiteVolume/finiteVolume/fvm/fvmSup.C
index 1fb653df8d..9ff2a5ca3b 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmSup.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -256,4 +256,17 @@ Foam::fvm::SuSp
 }
 
 
+template<class Type>
+Foam::tmp<Foam::fvMatrix<Type>>
+Foam::fvm::S
+(
+    const Pair<tmp<volScalarField::Internal>>& s,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tfvm = fvm::Sp(s[1], vf) + s[0];
+    return tfvm;
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.H b/src/finiteVolume/finiteVolume/fvm/fvmSup.H
index 879146d02b..4086143bf4 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmSup.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -150,6 +150,16 @@ namespace fvm
             const zero&,
             const GeometricField<Type, fvPatchField, volMesh>&
         );
+
+
+    // Explict and implicit sources provided as a pair
+
+        template<class Type>
+        tmp<fvMatrix<Type>> S
+        (
+            const Pair<tmp<volScalarField::Internal>>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
 }
 
 
diff --git a/src/twoPhaseModels/Allwmake b/src/twoPhaseModels/Allwmake
index e997639f22..ea5680cde4 100755
--- a/src/twoPhaseModels/Allwmake
+++ b/src/twoPhaseModels/Allwmake
@@ -9,5 +9,6 @@ wmake $targetType interfaceProperties
 wmake $targetType twoPhaseProperties
 wmake $targetType incompressibleTwoPhaseMixture
 wmake $targetType immiscibleIncompressibleTwoPhaseMixture
+wmake $targetType twoPhaseChange
 
 #------------------------------------------------------------------------------
diff --git a/src/twoPhaseModels/twoPhaseChange/Make/files b/src/twoPhaseModels/twoPhaseChange/Make/files
new file mode 100644
index 0000000000..a98f303a2b
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseChange/Make/files
@@ -0,0 +1,9 @@
+twoPhaseChangeModel/twoPhaseChangeModel.C
+twoPhaseChangeModel/twoPhaseChangeModelNew.C
+noPhaseChange/noPhaseChange.C
+cavitation/cavitationModel/cavitationModel.C
+cavitation/Kunz/Kunz.C
+cavitation/Merkle/Merkle.C
+cavitation/SchnerrSauer/SchnerrSauer.C
+
+LIB = $(FOAM_LIBBIN)/libtwoPhaseChangeModels
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options b/src/twoPhaseModels/twoPhaseChange/Make/options
similarity index 88%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options
rename to src/twoPhaseModels/twoPhaseChange/Make/options
index 90ef926582..d2a9bd3038 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Make/options
+++ b/src/twoPhaseModels/twoPhaseChange/Make/options
@@ -4,8 +4,10 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
     -limmiscibleIncompressibleTwoPhaseMixture \
+    -lmeshTools \
     -lfiniteVolume
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C b/src/twoPhaseModels/twoPhaseChange/cavitation/Kunz/Kunz.C
similarity index 56%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C
rename to src/twoPhaseModels/twoPhaseChange/cavitation/Kunz/Kunz.C
index 1f6353a0d2..55a6a14354 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/Kunz/Kunz.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,32 +30,31 @@ License
 
 namespace Foam
 {
-namespace phaseChangeTwoPhaseMixtures
+namespace twoPhaseChangeModels
 {
     defineTypeNameAndDebug(Kunz, 0);
-    addToRunTimeSelectionTable(phaseChangeTwoPhaseMixture, Kunz, components);
+    addToRunTimeSelectionTable(cavitationModel, Kunz, dictionary);
 }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::phaseChangeTwoPhaseMixtures::Kunz::Kunz
+Foam::twoPhaseChangeModels::Kunz::Kunz
 (
-    const volVectorField& U,
-    const surfaceScalarField& phi
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    phaseChangeTwoPhaseMixture(typeName, U, phi),
+    cavitationModel(typeName, mixture),
 
-    UInf_("UInf", dimVelocity, phaseChangeTwoPhaseMixtureCoeffs_),
-    tInf_("tInf", dimTime, phaseChangeTwoPhaseMixtureCoeffs_),
-    Cc_("Cc", dimless, phaseChangeTwoPhaseMixtureCoeffs_),
-    Cv_("Cv", dimless, phaseChangeTwoPhaseMixtureCoeffs_),
+    UInf_("UInf", dimVelocity, twoPhaseChangeModelCoeffs_),
+    tInf_("tInf", dimTime, twoPhaseChangeModelCoeffs_),
+    Cc_("Cc", dimless, twoPhaseChangeModelCoeffs_),
+    Cv_("Cv", dimless, twoPhaseChangeModelCoeffs_),
 
     p0_("0", pSat().dimensions(), 0.0),
 
-    mcCoeff_(Cc_*rho2()/tInf_),
-    mvCoeff_(Cv_*rho2()/(0.5*rho1()*sqr(UInf_)*tInf_))
+    mcCoeff_(Cc_*mixture_.rho2()/tInf_),
+    mvCoeff_(Cv_*mixture_.rho2()/(0.5*mixture_.rho1()*sqr(UInf_)*tInf_))
 {
     correct();
 }
@@ -64,10 +63,15 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::Kunz
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
+Foam::twoPhaseChangeModels::Kunz::mDotAlphal() const
 {
-    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
+    const volScalarField& p =
+        mixture_.U().db().lookupObject<volScalarField>("p");
+
+    const volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
 
     return Pair<tmp<volScalarField>>
     (
@@ -80,10 +84,15 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
+Foam::twoPhaseChangeModels::Kunz::mDotP() const
 {
-    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
+    const volScalarField& p =
+        mixture_.U().db().lookupObject<volScalarField>("p");
+
+    const volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
 
     return Pair<tmp<volScalarField>>
     (
@@ -95,25 +104,25 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
 }
 
 
-void Foam::phaseChangeTwoPhaseMixtures::Kunz::correct()
+void Foam::twoPhaseChangeModels::Kunz::correct()
 {
-    phaseChangeTwoPhaseMixture::correct();
+    cavitationModel::correct();
 }
 
 
-bool Foam::phaseChangeTwoPhaseMixtures::Kunz::read()
+bool Foam::twoPhaseChangeModels::Kunz::read()
 {
-    if (phaseChangeTwoPhaseMixture::read())
+    if (cavitationModel::read())
     {
-        phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
+        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
 
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("UInf") >> UInf_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("tInf") >> tInf_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cc") >> Cc_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cv") >> Cv_;
+        twoPhaseChangeModelCoeffs_.lookup("UInf") >> UInf_;
+        twoPhaseChangeModelCoeffs_.lookup("tInf") >> tInf_;
+        twoPhaseChangeModelCoeffs_.lookup("Cc") >> Cc_;
+        twoPhaseChangeModelCoeffs_.lookup("Cv") >> Cv_;
 
-        mcCoeff_ = Cc_*rho2()/tInf_;
-        mvCoeff_ = Cv_*rho2()/(0.5*rho1()*sqr(UInf_)*tInf_);
+        mcCoeff_ = Cc_*mixture_.rho2()/tInf_;
+        mvCoeff_ = Cv_*mixture_.rho2()/(0.5*mixture_.rho1()*sqr(UInf_)*tInf_);
 
         return true;
     }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.H b/src/twoPhaseModels/twoPhaseChange/cavitation/Kunz/Kunz.H
similarity index 90%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.H
rename to src/twoPhaseModels/twoPhaseChange/cavitation/Kunz/Kunz.H
index 8c0b4be206..aa2cc4f433 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.H
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/Kunz/Kunz.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::phaseChangeTwoPhaseMixtures::Kunz
+    Foam::twoPhaseChangeModels::Kunz
 
 Description
     Kunz cavitation model slightly modified so that the condensation term
@@ -49,13 +49,13 @@ SourceFiles
 #ifndef Kunz_H
 #define Kunz_H
 
-#include "phaseChangeTwoPhaseMixture.H"
+#include "cavitationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace phaseChangeTwoPhaseMixtures
+namespace twoPhaseChangeModels
 {
 
 /*---------------------------------------------------------------------------*\
@@ -64,7 +64,7 @@ namespace phaseChangeTwoPhaseMixtures
 
 class Kunz
 :
-    public phaseChangeTwoPhaseMixture
+    public cavitationModel
 {
     // Private Data
 
@@ -87,11 +87,10 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct for mixture
         Kunz
         (
-            const volVectorField& U,
-            const surfaceScalarField& phi
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
         );
 
 
@@ -121,7 +120,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace phaseChangeTwoPhaseMixtures
+} // End namespace twoPhaseChangeModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C b/src/twoPhaseModels/twoPhaseChange/cavitation/Merkle/Merkle.C
similarity index 57%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C
rename to src/twoPhaseModels/twoPhaseChange/cavitation/Merkle/Merkle.C
index 2a1b29cfd4..dc346bef4c 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/Merkle/Merkle.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,32 +30,31 @@ License
 
 namespace Foam
 {
-namespace phaseChangeTwoPhaseMixtures
+namespace twoPhaseChangeModels
 {
     defineTypeNameAndDebug(Merkle, 0);
-    addToRunTimeSelectionTable(phaseChangeTwoPhaseMixture, Merkle, components);
+    addToRunTimeSelectionTable(cavitationModel, Merkle, dictionary);
 }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::phaseChangeTwoPhaseMixtures::Merkle::Merkle
+Foam::twoPhaseChangeModels::Merkle::Merkle
 (
-    const volVectorField& U,
-    const surfaceScalarField& phi
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    phaseChangeTwoPhaseMixture(typeName, U, phi),
+    cavitationModel(typeName, mixture),
 
-    UInf_("UInf", dimVelocity, phaseChangeTwoPhaseMixtureCoeffs_),
-    tInf_("tInf", dimTime, phaseChangeTwoPhaseMixtureCoeffs_),
-    Cc_("Cc", dimless, phaseChangeTwoPhaseMixtureCoeffs_),
-    Cv_("Cv", dimless, phaseChangeTwoPhaseMixtureCoeffs_),
+    UInf_("UInf", dimVelocity, twoPhaseChangeModelCoeffs_),
+    tInf_("tInf", dimTime, twoPhaseChangeModelCoeffs_),
+    Cc_("Cc", dimless, twoPhaseChangeModelCoeffs_),
+    Cv_("Cv", dimless, twoPhaseChangeModelCoeffs_),
 
     p0_("0", pSat().dimensions(), 0.0),
 
     mcCoeff_(Cc_/(0.5*sqr(UInf_)*tInf_)),
-    mvCoeff_(Cv_*rho1()/(0.5*sqr(UInf_)*tInf_*rho2()))
+    mvCoeff_(Cv_*mixture_.rho1()/(0.5*sqr(UInf_)*tInf_*mixture_.rho2()))
 {
     correct();
 }
@@ -64,9 +63,10 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::Merkle
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotAlphal() const
+Foam::twoPhaseChangeModels::Merkle::mDotAlphal() const
 {
-    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
+    const volScalarField& p =
+        mixture_.U().db().lookupObject<volScalarField>("p");
 
     return Pair<tmp<volScalarField>>
     (
@@ -77,10 +77,15 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotAlphal() const
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const
+Foam::twoPhaseChangeModels::Merkle::mDotP() const
 {
-    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
+    const volScalarField& p =
+        mixture_.U().db().lookupObject<volScalarField>("p");
+
+    const volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
 
     return Pair<tmp<volScalarField>>
     (
@@ -90,25 +95,25 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const
 }
 
 
-void Foam::phaseChangeTwoPhaseMixtures::Merkle::correct()
+void Foam::twoPhaseChangeModels::Merkle::correct()
 {
-    phaseChangeTwoPhaseMixture::correct();
+    cavitationModel::correct();
 }
 
 
-bool Foam::phaseChangeTwoPhaseMixtures::Merkle::read()
+bool Foam::twoPhaseChangeModels::Merkle::read()
 {
-    if (phaseChangeTwoPhaseMixture::read())
+    if (cavitationModel::read())
     {
-        phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
+        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
 
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("UInf") >> UInf_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("tInf") >> tInf_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cc") >> Cc_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cv") >> Cv_;
+        twoPhaseChangeModelCoeffs_.lookup("UInf") >> UInf_;
+        twoPhaseChangeModelCoeffs_.lookup("tInf") >> tInf_;
+        twoPhaseChangeModelCoeffs_.lookup("Cc") >> Cc_;
+        twoPhaseChangeModelCoeffs_.lookup("Cv") >> Cv_;
 
         mcCoeff_ = Cc_/(0.5*sqr(UInf_)*tInf_);
-        mvCoeff_ = Cv_*rho1()/(0.5*sqr(UInf_)*tInf_*rho2());
+        mvCoeff_ = Cv_*mixture_.rho1()/(0.5*sqr(UInf_)*tInf_*mixture_.rho2());
 
         return true;
     }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.H b/src/twoPhaseModels/twoPhaseChange/cavitation/Merkle/Merkle.H
similarity index 89%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.H
rename to src/twoPhaseModels/twoPhaseChange/cavitation/Merkle/Merkle.H
index 8283763888..9d37437f63 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.H
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/Merkle/Merkle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::phaseChangeTwoPhaseMixtures::Merkle
+    Foam::twoPhaseChangeModels::Merkle
 
 Description
     Merkle cavitation model.
@@ -43,13 +43,13 @@ SourceFiles
 #ifndef Merkle_H
 #define Merkle_H
 
-#include "phaseChangeTwoPhaseMixture.H"
+#include "cavitationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace phaseChangeTwoPhaseMixtures
+namespace twoPhaseChangeModels
 {
 
 /*---------------------------------------------------------------------------*\
@@ -58,7 +58,7 @@ namespace phaseChangeTwoPhaseMixtures
 
 class Merkle
 :
-    public phaseChangeTwoPhaseMixture
+    public cavitationModel
 {
     // Private Data
 
@@ -81,11 +81,10 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct for mixture
         Merkle
         (
-            const volVectorField& U,
-            const surfaceScalarField& phi
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
         );
 
 
@@ -115,7 +114,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace phaseChangeTwoPhaseMixtures
+} // End namespace twoPhaseChangeModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C b/src/twoPhaseModels/twoPhaseChange/cavitation/SchnerrSauer/SchnerrSauer.C
similarity index 57%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C
rename to src/twoPhaseModels/twoPhaseChange/cavitation/SchnerrSauer/SchnerrSauer.C
index cfd24fae3d..ed9e77910f 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/SchnerrSauer/SchnerrSauer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,14 +31,14 @@ License
 
 namespace Foam
 {
-namespace phaseChangeTwoPhaseMixtures
+namespace twoPhaseChangeModels
 {
     defineTypeNameAndDebug(SchnerrSauer, 0);
     addToRunTimeSelectionTable
     (
-        phaseChangeTwoPhaseMixture,
+        cavitationModel,
         SchnerrSauer,
-        components
+        dictionary
     );
 }
 }
@@ -46,18 +46,17 @@ namespace phaseChangeTwoPhaseMixtures
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::SchnerrSauer
+Foam::twoPhaseChangeModels::SchnerrSauer::SchnerrSauer
 (
-    const volVectorField& U,
-    const surfaceScalarField& phi
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    phaseChangeTwoPhaseMixture(typeName, U, phi),
+    cavitationModel(typeName, mixture),
 
-    n_("n", dimless/dimVolume, phaseChangeTwoPhaseMixtureCoeffs_),
-    dNuc_("dNuc", dimLength, phaseChangeTwoPhaseMixtureCoeffs_),
-    Cc_("Cc", dimless, phaseChangeTwoPhaseMixtureCoeffs_),
-    Cv_("Cv", dimless, phaseChangeTwoPhaseMixtureCoeffs_),
+    n_("n", dimless/dimVolume, twoPhaseChangeModelCoeffs_),
+    dNuc_("dNuc", dimLength, twoPhaseChangeModelCoeffs_),
+    Cc_("Cc", dimless, twoPhaseChangeModelCoeffs_),
+    Cv_("Cv", dimless, twoPhaseChangeModelCoeffs_),
 
     p0_("0", pSat().dimensions(), 0.0)
 {
@@ -68,7 +67,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::SchnerrSauer
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::rRb
+Foam::twoPhaseChangeModels::SchnerrSauer::rRb
 (
     const volScalarField& limitedAlpha1
 ) const
@@ -83,7 +82,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::rRb
 
 
 Foam::dimensionedScalar
-Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::alphaNuc() const
+Foam::twoPhaseChangeModels::SchnerrSauer::alphaNuc() const
 {
     dimensionedScalar Vnuc = n_*constant::mathematical::pi*pow3(dNuc_)/6;
     return Vnuc/(1 + Vnuc);
@@ -91,30 +90,40 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::alphaNuc() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::pCoeff
+Foam::twoPhaseChangeModels::SchnerrSauer::pCoeff
 (
     const volScalarField& p
 ) const
 {
-    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
-    volScalarField rho
+    const volScalarField limitedAlpha1
     (
-        limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    const volScalarField rho
+    (
+        limitedAlpha1*mixture_.rho1()
+      + (scalar(1) - limitedAlpha1)*mixture_.rho2()
     );
 
     return
-        (3*rho1()*rho2())*sqrt(2/(3*rho1()))
+        (3*mixture_.rho1()*mixture_.rho2())*sqrt(2/(3*mixture_.rho1()))
        *rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const
+Foam::twoPhaseChangeModels::SchnerrSauer::mDotAlphal() const
 {
-    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField pCoeff(this->pCoeff(p));
+    const volScalarField& p =
+        mixture_.U().db().lookupObject<volScalarField>("p");
 
-    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
+    const volScalarField pCoeff(this->pCoeff(p));
+
+    const volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
 
     return Pair<tmp<volScalarField>>
     (
@@ -126,13 +135,19 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const
+Foam::twoPhaseChangeModels::SchnerrSauer::mDotP() const
 {
-    const volScalarField& p = alpha1().db().lookupObject<volScalarField>("p");
-    volScalarField pCoeff(this->pCoeff(p));
+    const volScalarField& p =
+        mixture_.U().db().lookupObject<volScalarField>("p");
 
-    volScalarField limitedAlpha1(min(max(alpha1(), scalar(0)), scalar(1)));
-    volScalarField apCoeff(limitedAlpha1*pCoeff);
+    const volScalarField pCoeff(this->pCoeff(p));
+
+    const volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    const volScalarField apCoeff(limitedAlpha1*pCoeff);
 
     return Pair<tmp<volScalarField>>
     (
@@ -143,22 +158,22 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const
 }
 
 
-void Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::correct()
+void Foam::twoPhaseChangeModels::SchnerrSauer::correct()
 {
-    phaseChangeTwoPhaseMixture::correct();
+    cavitationModel::correct();
 }
 
 
-bool Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::read()
+bool Foam::twoPhaseChangeModels::SchnerrSauer::read()
 {
-    if (phaseChangeTwoPhaseMixture::read())
+    if (cavitationModel::read())
     {
-        phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
+        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
 
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("n") >> n_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("dNuc") >> dNuc_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cc") >> Cc_;
-        phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cv") >> Cv_;
+        twoPhaseChangeModelCoeffs_.lookup("n") >> n_;
+        twoPhaseChangeModelCoeffs_.lookup("dNuc") >> dNuc_;
+        twoPhaseChangeModelCoeffs_.lookup("Cc") >> Cc_;
+        twoPhaseChangeModelCoeffs_.lookup("Cv") >> Cv_;
 
         return true;
     }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.H b/src/twoPhaseModels/twoPhaseChange/cavitation/SchnerrSauer/SchnerrSauer.H
similarity index 90%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.H
rename to src/twoPhaseModels/twoPhaseChange/cavitation/SchnerrSauer/SchnerrSauer.H
index c741048645..5d0309eccd 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.H
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/SchnerrSauer/SchnerrSauer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer
+    Foam::twoPhaseChangeModels::SchnerrSauer
 
 Description
     SchnerrSauer cavitation model.
@@ -43,13 +43,13 @@ SourceFiles
 #ifndef SchnerrSauer_H
 #define SchnerrSauer_H
 
-#include "phaseChangeTwoPhaseMixture.H"
+#include "cavitationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace phaseChangeTwoPhaseMixtures
+namespace twoPhaseChangeModels
 {
 
 /*---------------------------------------------------------------------------*\
@@ -58,7 +58,7 @@ namespace phaseChangeTwoPhaseMixtures
 
 class SchnerrSauer
 :
-    public phaseChangeTwoPhaseMixture
+    public cavitationModel
 {
     // Private Data
 
@@ -94,11 +94,10 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct for mixture
         SchnerrSauer
         (
-            const volVectorField& U,
-            const surfaceScalarField& phi
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
         );
 
 
@@ -128,7 +127,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace phaseChangeTwoPhaseMixtures
+} // End namespace twoPhaseChangeModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.C b/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.C
new file mode 100644
index 0000000000..a64bb4e948
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.C
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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 "cavitationModel.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace twoPhaseChangeModels
+{
+    defineTypeNameAndDebug(cavitationModel, 0);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::twoPhaseChangeModels::cavitationModel::cavitationModel
+(
+    const word& type,
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
+)
+:
+    twoPhaseChangeModel(type, mixture),
+    pSat_("pSat", dimPressure, lookup("pSat"))
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::twoPhaseChangeModels::cavitationModel::Salpha
+(
+    volScalarField& alpha
+) const
+{
+    volScalarField::Internal alphalCoeff
+    (
+        1.0/mixture_.rho1()
+      - mixture_.alpha1()()*(1.0/mixture_.rho1() - 1.0/mixture_.rho2())
+    );
+    Pair<tmp<volScalarField>> mDotAlphal = this->mDotAlphal();
+
+    const volScalarField::Internal vDotcAlphal(alphalCoeff*mDotAlphal[0]());
+    const volScalarField::Internal vDotvAlphal(alphalCoeff*mDotAlphal[1]());
+
+    return Pair<tmp<volScalarField::Internal>>
+    (
+        1.0*vDotcAlphal,
+        vDotvAlphal - vDotcAlphal
+    );
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::twoPhaseChangeModels::cavitationModel::Sp_rgh
+(
+    const volScalarField& rho,
+    const volScalarField& gh,
+    volScalarField& p_rgh
+) const
+{
+    dimensionedScalar pCoeff(1.0/mixture_.rho1() - 1.0/mixture_.rho2());
+    Pair<tmp<volScalarField>> mDotP = this->mDotP();
+
+    const volScalarField vDotcP = pCoeff*mDotP[0];
+    const volScalarField vDotvP = pCoeff*mDotP[1];
+
+    return
+        (vDotvP - vDotcP)*(pSat() - rho*gh)
+      - fvm::Sp(vDotvP - vDotcP, p_rgh);
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::twoPhaseChangeModels::cavitationModel::SU
+(
+    const volScalarField& rho,
+    const surfaceScalarField& rhoPhi,
+    volVectorField& U
+) const
+{
+    return fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U);
+}
+
+
+bool Foam::twoPhaseChangeModels::cavitationModel::read()
+{
+    if (twoPhaseChangeModel::read())
+    {
+        lookup("pSat") >> pSat_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H b/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.H
similarity index 53%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
rename to src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.H
index b68e41274b..90aceef663 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,45 +22,37 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::phaseChangeTwoPhaseMixture
+    Foam::twoPhaseChangeModels::cavitationModel
 
 Description
+    Abstract base class for cavitation models
 
 SourceFiles
-    phaseChangeTwoPhaseMixture.C
-    phaseChangeModelNew.C
+    cavitationModel.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef phaseChangeTwoPhaseMixture_H
-#define phaseChangeTwoPhaseMixture_H
+#ifndef cavitationModel_H
+#define cavitationModel_H
 
-#include "immiscibleIncompressibleTwoPhaseMixture.H"
-#include "typeInfo.H"
-#include "runTimeSelectionTables.H"
-#include "volFields.H"
-#include "dimensionedScalar.H"
-#include "autoPtr.H"
-#include "Pair.H"
+#include "twoPhaseChangeModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace twoPhaseChangeModels
+{
 
 /*---------------------------------------------------------------------------*\
-                 Class phaseChangeTwoPhaseMixture Declaration
+                                 Class cavitationModel
 \*---------------------------------------------------------------------------*/
 
-class phaseChangeTwoPhaseMixture
+class cavitationModel
 :
-    public immiscibleIncompressibleTwoPhaseMixture
+    public twoPhaseChangeModel
 {
-protected:
-
-    // Protected data
-
-        dictionary phaseChangeTwoPhaseMixtureCoeffs_;
+    // Private data
 
         //- Saturation vapour pressure
         dimensionedScalar pSat_;
@@ -69,50 +61,21 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("phaseChangeTwoPhaseMixture");
-
-
-    // Declare run-time constructor selection table
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            phaseChangeTwoPhaseMixture,
-            components,
-            (
-                const volVectorField& U,
-                const surfaceScalarField& phi
-            ),
-            (U, phi)
-        );
+    TypeName("cavitation");
 
 
     // Constructors
 
-        //- Construct from components
-        phaseChangeTwoPhaseMixture
+        //- Construct for mixture
+        cavitationModel
         (
             const word& type,
-            const volVectorField& U,
-            const surfaceScalarField& phi
-        );
-
-        //- Disallow default bitwise copy construction
-        phaseChangeTwoPhaseMixture(const phaseChangeTwoPhaseMixture&);
-
-
-    // Selectors
-
-        //- Return a reference to the selected phaseChange model
-        static autoPtr<phaseChangeTwoPhaseMixture> New
-        (
-            const volVectorField& U,
-            const surfaceScalarField& phi
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
         );
 
 
     //- Destructor
-    virtual ~phaseChangeTwoPhaseMixture()
+    virtual ~cavitationModel()
     {}
 
 
@@ -133,31 +96,38 @@ public:
         //  to multiply (p - pSat)
         virtual Pair<tmp<volScalarField>> mDotP() const = 0;
 
-        //- Return the volumetric condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        Pair<tmp<volScalarField>> vDotAlphal() const;
+        //- Return the cavitation explicit and implicit sources
+        //  for the phase-fraction equation
+        virtual Pair<tmp<volScalarField::Internal>> Salpha
+        (
+            volScalarField& alpha
+        ) const;
 
-        //- Return the volumetric condensation and vaporisation rates as
-        //  coefficients to multiply (p - pSat)
-        Pair<tmp<volScalarField>> vDotP() const;
+        //- Return the cavitation source matrix
+        //  for the p_rgh pressure equation
+        virtual tmp<fvScalarMatrix> Sp_rgh
+        (
+            const volScalarField& rho,
+            const volScalarField& gh,
+            volScalarField& p_rgh
+        ) const;
 
-        //- Correct the phaseChange model
-        virtual void correct() = 0;
+        //- Return the cavitation source matrix for the momentum equation
+        virtual tmp<fvVectorMatrix> SU
+        (
+            const volScalarField& rho,
+            const surfaceScalarField& rhoPhi,
+            volVectorField& U
+        ) const;
 
         //- Read the transportProperties dictionary and update
-        virtual bool read() = 0;
-
-
-    // Member Operators
-
-        //- Disallow default bitwise assignment
-        void operator=(const phaseChangeTwoPhaseMixture&) = delete;
+        virtual bool read();
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace twoPhaseChangeModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.C b/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.C
new file mode 100644
index 0000000000..8a3a9d1c8c
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.C
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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 "noPhaseChange.H"
+#include "fvScalarMatrix.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace twoPhaseChangeModels
+{
+    defineTypeNameAndDebug(noPhaseChange, 0);
+    addToRunTimeSelectionTable(twoPhaseChangeModel, noPhaseChange, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::twoPhaseChangeModels::noPhaseChange::noPhaseChange
+(
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
+)
+:
+    twoPhaseChangeModel(typeName, mixture)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+Foam::twoPhaseChangeModels::noPhaseChange::mDotAlphal() const
+{
+    return Pair<tmp<volScalarField>>
+    (
+        volScalarField::null(),
+        volScalarField::null()
+    );
+}
+
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+Foam::twoPhaseChangeModels::noPhaseChange::mDotP() const
+{
+    return Pair<tmp<volScalarField>>
+    (
+        volScalarField::null(),
+        volScalarField::null()
+    );
+}
+
+
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::twoPhaseChangeModels::noPhaseChange::Salpha
+(
+    volScalarField& alpha
+) const
+{
+    return Pair<tmp<volScalarField::Internal>>
+    (
+        tmp<volScalarField::Internal>(nullptr),
+        tmp<volScalarField::Internal>(nullptr)
+    );
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::twoPhaseChangeModels::noPhaseChange::Sp_rgh
+(
+    const volScalarField& rho,
+    const volScalarField& gh,
+    volScalarField& p_rgh
+) const
+{
+    return tmp<fvScalarMatrix>(new fvScalarMatrix(p_rgh, dimVolume/dimTime));
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::twoPhaseChangeModels::noPhaseChange::SU
+(
+    const volScalarField& rho,
+    const surfaceScalarField& rhoPhi,
+    volVectorField& U
+) const
+{
+    return tmp<fvVectorMatrix>
+    (
+        new fvVectorMatrix(U, dimMass*dimVelocity/dimTime)
+    );
+}
+
+
+void Foam::twoPhaseChangeModels::noPhaseChange::correct()
+{
+    twoPhaseChangeModel::correct();
+}
+
+
+bool Foam::twoPhaseChangeModels::noPhaseChange::read()
+{
+    return twoPhaseChangeModel::read();
+}
+
+
+// ************************************************************************* //
diff --git a/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.H b/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.H
new file mode 100644
index 0000000000..9d8bc3861f
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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::twoPhaseChangeModels::noPhaseChange
+
+Description
+    Dummy phase change model for 'no phase change'
+
+SourceFiles
+    noPhaseChange.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noPhaseChange_H
+#define noPhaseChange_H
+
+#include "twoPhaseChangeModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace twoPhaseChangeModels
+{
+
+/*---------------------------------------------------------------------------*\
+                                 Class noPhaseChange
+\*---------------------------------------------------------------------------*/
+
+class noPhaseChange
+:
+    public twoPhaseChangeModel
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct for mixture
+        noPhaseChange
+        (
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
+        );
+
+
+    //- Destructor
+    virtual ~noPhaseChange()
+    {}
+
+
+    // Member Functions
+
+        //- Return the mass condensation and vaporisation rates as a
+        //  coefficient to multiply (1 - alphal) for the condensation rate
+        //  and a coefficient to multiply  alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField>> mDotAlphal() const;
+
+        //- Return the mass condensation and vaporisation rates as coefficients
+        //  to multiply (p - pSat)
+        virtual Pair<tmp<volScalarField>> mDotP() const;
+
+        //- Return unallocated explicit and implicit sources
+        //  for the phase-fraction equation
+        virtual Pair<tmp<volScalarField::Internal>> Salpha
+        (
+            volScalarField&
+        ) const;
+
+        //- Return the zero source matrix
+        //  for the p_rgh pressure equation
+        virtual tmp<fvScalarMatrix> Sp_rgh
+        (
+            const volScalarField& rho,
+            const volScalarField& gh,
+            volScalarField& p_rgh
+        ) const;
+
+        //- Return the zero source matrix for the momentum equation
+        virtual tmp<fvVectorMatrix> SU
+        (
+            const volScalarField& rho,
+            const surfaceScalarField& rhoPhi,
+            volVectorField& U
+        ) const;
+
+        //- Correct the noPhaseChange phaseChange model
+        virtual void correct();
+
+        //- Read the transportProperties dictionary and update
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace twoPhaseChangeModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.C
similarity index 51%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
rename to src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.C
index f4c7b22fa8..b5a5df8cec 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
+++ b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,69 +23,76 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "phaseChangeTwoPhaseMixture.H"
+#include "twoPhaseChangeModel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(phaseChangeTwoPhaseMixture, 0);
-    defineRunTimeSelectionTable(phaseChangeTwoPhaseMixture, components);
+    defineTypeNameAndDebug(twoPhaseChangeModel, 0);
+    defineRunTimeSelectionTable(twoPhaseChangeModel, dictionary);
 }
 
+const Foam::word Foam::twoPhaseChangeModel::phaseChangePropertiesName
+(
+    "phaseChangeProperties"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::IOobject Foam::twoPhaseChangeModel::createIOobject
+(
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
+) const
+{
+    IOobject io
+    (
+        phaseChangePropertiesName,
+        mixture.U().mesh().time().constant(),
+        mixture.U().mesh(),
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (io.typeHeaderOk<IOdictionary>(true))
+    {
+        io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
+        return io;
+    }
+    else
+    {
+        io.readOpt() = IOobject::NO_READ;
+        return io;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::phaseChangeTwoPhaseMixture::phaseChangeTwoPhaseMixture
+Foam::twoPhaseChangeModel::twoPhaseChangeModel
 (
     const word& type,
-    const volVectorField& U,
-    const surfaceScalarField& phi
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    immiscibleIncompressibleTwoPhaseMixture(U, phi),
-    phaseChangeTwoPhaseMixtureCoeffs_(optionalSubDict(type + "Coeffs")),
-    pSat_("pSat", dimPressure, lookup("pSat"))
+    IOdictionary(createIOobject(mixture)),
+    mixture_(mixture),
+    twoPhaseChangeModelCoeffs_(optionalSubDict(type + "Coeffs"))
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixture::vDotAlphal() const
+void Foam::twoPhaseChangeModel::correct()
+{}
+
+
+bool Foam::twoPhaseChangeModel::read()
 {
-    volScalarField alphalCoeff(1.0/rho1() - alpha1()*(1.0/rho1() - 1.0/rho2()));
-    Pair<tmp<volScalarField>> mDotAlphal = this->mDotAlphal();
-
-    return Pair<tmp<volScalarField>>
-    (
-        alphalCoeff*mDotAlphal[0],
-        alphalCoeff*mDotAlphal[1]
-    );
-}
-
-
-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::phaseChangeTwoPhaseMixture::vDotP() const
-{
-    dimensionedScalar pCoeff(1.0/rho1() - 1.0/rho2());
-    Pair<tmp<volScalarField>> mDotP = this->mDotP();
-
-    return Pair<tmp<volScalarField>>(pCoeff*mDotP[0], pCoeff*mDotP[1]);
-}
-
-
-void Foam::phaseChangeTwoPhaseMixture::correct()
-{
-    immiscibleIncompressibleTwoPhaseMixture::correct();
-}
-
-
-bool Foam::phaseChangeTwoPhaseMixture::read()
-{
-    if (immiscibleIncompressibleTwoPhaseMixture::read())
+    if (regIOobject::read())
     {
-        phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
-        lookup("pSat") >> pSat_;
+        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
 
         return true;
     }
diff --git a/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.H b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.H
new file mode 100644
index 0000000000..bfe464d33e
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.H
@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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::twoPhaseChangeModel
+
+Description
+
+SourceFiles
+    twoPhaseChangeModel.C
+    phaseChangeModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef twoPhaseChangeModel_H
+#define twoPhaseChangeModel_H
+
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "fvMatricesFwd.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class twoPhaseChangeModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class twoPhaseChangeModel
+:
+    public IOdictionary
+{
+    // Private Member Functions
+
+        //- Construct the base IO object
+        IOobject createIOobject
+        (
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
+        ) const;
+
+
+protected:
+
+    // Protected data
+
+        //- Reference to the two-phase mixture
+        const immiscibleIncompressibleTwoPhaseMixture& mixture_;
+
+        //- Model coefficient dictionary
+        dictionary twoPhaseChangeModelCoeffs_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("phaseChangeModel");
+
+    //- Default phaseChangeProperties dictionary name
+    static const word phaseChangePropertiesName;
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            twoPhaseChangeModel,
+            dictionary,
+            (
+                const immiscibleIncompressibleTwoPhaseMixture& mixture
+            ),
+            (mixture)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        twoPhaseChangeModel
+        (
+            const word& type,
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
+        );
+
+        //- Disallow default bitwise copy construction
+        twoPhaseChangeModel(const twoPhaseChangeModel&);
+
+
+    // Selectors
+
+        //- Return a reference to the selected phaseChange model
+        static autoPtr<twoPhaseChangeModel> New
+        (
+            const immiscibleIncompressibleTwoPhaseMixture& mixture
+        );
+
+
+    //- Destructor
+    virtual ~twoPhaseChangeModel()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phase-change explicit and implicit sources
+        //  for the phase-fraction equation
+        virtual Pair<tmp<volScalarField::Internal>> Salpha
+        (
+            volScalarField& alpha
+        ) const = 0;
+
+        //- Return the phase-change source matrix
+        //  for the p_rgh pressure equation
+        virtual tmp<fvScalarMatrix> Sp_rgh
+        (
+            const volScalarField& rho,
+            const volScalarField& gh,
+            volScalarField& p_rgh
+        ) const = 0;
+
+        //- Return the phase-change source matrix for the momentum equation
+        virtual tmp<fvVectorMatrix> SU
+        (
+            const volScalarField& rho,
+            const surfaceScalarField& rhoPhi,
+            volVectorField& U
+        ) const = 0;
+
+        //- Correct the phaseChange model
+        virtual void correct() = 0;
+
+        //- Read the transportProperties dictionary and update
+        virtual bool read() = 0;
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const twoPhaseChangeModel&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModelNew.C
similarity index 53%
rename from applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
rename to src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModelNew.C
index 013b6cf25c..56c27d36f8 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
+++ b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModelNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,54 +23,60 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "phaseChangeTwoPhaseMixture.H"
-#include "incompressibleTwoPhaseMixture.H"
+#include "noPhaseChange.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::autoPtr<Foam::phaseChangeTwoPhaseMixture>
-Foam::phaseChangeTwoPhaseMixture::New
+Foam::autoPtr<Foam::twoPhaseChangeModel> Foam::twoPhaseChangeModel::New
 (
-    const volVectorField& U,
-    const surfaceScalarField& phi
+    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 {
-    IOdictionary transportPropertiesDict
+    IOobject twoPhaseChangeModelIO
     (
         IOobject
         (
-            "transportProperties",
-            U.time().constant(),
-            U.db(),
+            phaseChangePropertiesName,
+            mixture.U().time().constant(),
+            mixture.U().db(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE,
             false
         )
     );
 
-    word phaseChangeTwoPhaseMixtureTypeName
-    (
-        transportPropertiesDict.lookup("phaseChangeTwoPhaseMixture")
-    );
+    word modelType(twoPhaseChangeModels::noPhaseChange::typeName);
 
-    Info<< "Selecting phaseChange model "
-        << phaseChangeTwoPhaseMixtureTypeName << endl;
+    if (twoPhaseChangeModelIO.typeHeaderOk<IOdictionary>(false))
+    {
+        IOdictionary(twoPhaseChangeModelIO).lookup
+        (
+            twoPhaseChangeModel::typeName
+        ) >> modelType;
+    }
+    else
+    {
+        Info<< "No phase change: "
+            << twoPhaseChangeModelIO.name()
+            << " not found" << endl;
+    }
 
-    componentsConstructorTable::iterator cstrIter =
-        componentsConstructorTablePtr_
-            ->find(phaseChangeTwoPhaseMixtureTypeName);
+    Info<< "Selecting phaseChange model " << modelType << endl;
 
-    if (cstrIter == componentsConstructorTablePtr_->end())
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorInFunction
-            << "Unknown phaseChangeTwoPhaseMixture type "
-            << phaseChangeTwoPhaseMixtureTypeName << endl << endl
-            << "Valid  phaseChangeTwoPhaseMixtures are : " << endl
-            << componentsConstructorTablePtr_->sortedToc()
+            << "Unknown " << twoPhaseChangeModel::typeName<< " type "
+            << modelType << nl << nl
+            << "Valid  twoPhaseChangeModels are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
 
-    return autoPtr<phaseChangeTwoPhaseMixture>(cstrIter()(U, phi));
+    return autoPtr<twoPhaseChangeModel>(cstrIter()(mixture));
 }
 
 
diff --git a/src/twoPhaseModels/twoPhaseMixture/Make/files b/src/twoPhaseModels/twoPhaseMixture/Make/files
index c537f522c0..6cb3ae7197 100644
--- a/src/twoPhaseModels/twoPhaseMixture/Make/files
+++ b/src/twoPhaseModels/twoPhaseMixture/Make/files
@@ -1,6 +1,7 @@
 twoPhaseMixture/twoPhaseMixture.C
 
 interfaceCompression/interfaceCompression.C
+noInterfaceCompression/noInterfaceCompression.C
 
 MPLIC/MPLICface.C
 MPLIC/MPLICcellStorage.C
diff --git a/src/twoPhaseModels/twoPhaseMixture/VoF/alphaEqn.H b/src/twoPhaseModels/twoPhaseMixture/VoF/alphaEqn.H
index 42b7aa5f73..b368fc87b7 100644
--- a/src/twoPhaseModels/twoPhaseMixture/VoF/alphaEqn.H
+++ b/src/twoPhaseModels/twoPhaseMixture/VoF/alphaEqn.H
@@ -66,10 +66,10 @@
         );
     }
 
+    #include "alphaSuSp.H"
+
     if (MULESCorr)
     {
-        #include "alphaSuSp.H"
-
         fvScalarMatrix alpha1Eqn
         (
             (
@@ -83,12 +83,13 @@
                 phiCN,
                 upwind<scalar>(mesh, phiCN)
             ).fvmDiv(phiCN, alpha1)
-       // - fvm::Sp(fvc::ddt(dimensionedScalar(dimless, 1), mesh)
-       //           + fvc::div(phiCN), alpha1)
-         ==
-            Su + fvm::Sp(Sp + divU, alpha1)
         );
 
+        if (divU.valid())
+        {
+            alpha1Eqn -= Su() + fvm::Sp(Sp() + divU(), alpha1);
+        }
+
         alpha1Eqn.solve();
 
         Info<< "Phase-1 volume fraction = "
@@ -126,8 +127,6 @@
 
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
     {
-        #include "alphaSuSp.H"
-
         // Split operator
         tmp<surfaceScalarField> talphaPhi1Un
         (
@@ -144,17 +143,32 @@
             tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
             volScalarField alpha10("alpha10", alpha1);
 
-            MULES::correct
-            (
-                geometricOneField(),
-                alpha1,
-                talphaPhi1Un(),
-                talphaPhi1Corr.ref(),
-                Sp,
-                (-Sp*alpha1)(),
-                oneField(),
-                zeroField()
-            );
+            if (divU.valid())
+            {
+                MULES::correct
+                (
+                    geometricOneField(),
+                    alpha1,
+                    talphaPhi1Un(),
+                    talphaPhi1Corr.ref(),
+                    Sp(),
+                    (-Sp()*alpha1)(),
+                    oneField(),
+                    zeroField()
+                );
+            }
+            else
+            {
+                MULES::correct
+                (
+                    geometricOneField(),
+                    alpha1,
+                    talphaPhi1Un(),
+                    talphaPhi1Corr.ref(),
+                    oneField(),
+                    zeroField()
+                );
+            }
 
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
@@ -171,17 +185,32 @@
         {
             alphaPhi10 = talphaPhi1Un;
 
-            MULES::explicitSolve
-            (
-                geometricOneField(),
-                alpha1,
-                phiCN,
-                alphaPhi10,
-                Sp,
-                (Su + divU*min(alpha1(), scalar(1)))(),
-                oneField(),
-                zeroField()
-            );
+            if (divU.valid())
+            {
+                MULES::explicitSolve
+                (
+                    geometricOneField(),
+                    alpha1,
+                    phiCN,
+                    alphaPhi10,
+                    Sp(),
+                    (Su() + divU()*min(alpha1(), scalar(1)))(),
+                    oneField(),
+                    zeroField()
+                );
+            }
+            else
+            {
+                MULES::explicitSolve
+                (
+                    geometricOneField(),
+                    alpha1,
+                    phiCN,
+                    alphaPhi10,
+                    oneField(),
+                    zeroField()
+                );
+            }
         }
 
         alpha2 = 1.0 - alpha1;
diff --git a/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H b/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H
index c7f8837709..76b2b37def 100644
--- a/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H
+++ b/src/twoPhaseModels/twoPhaseMixture/VoF/alphaScheme.H
@@ -1,6 +1,7 @@
 static const wordHashSet compressionSchemes
 {
     "interfaceCompression",
+    "noInterfaceCompression",
     "PLIC",
     "PLICU",
     "MPLIC",
diff --git a/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H b/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H
index db549c9738..5748069779 100644
--- a/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H
+++ b/src/twoPhaseModels/twoPhaseMixture/VoF/setRDeltaT.H
@@ -59,8 +59,8 @@
     rDeltaT.ref() = max
     (
         1/dimensionedScalar(dimTime, maxDeltaT),
-        fvc::surfaceSum(mag(rhoPhi))()()
-       /((2*maxCo)*mesh.V()*rho())
+        fvc::surfaceSum(mag(phi))()()
+       /((2*maxCo)*mesh.V())
     );
 
     if (maxAlphaCo < maxCo)
diff --git a/src/twoPhaseModels/twoPhaseMixture/noInterfaceCompression/noInterfaceCompression.C b/src/twoPhaseModels/twoPhaseMixture/noInterfaceCompression/noInterfaceCompression.C
new file mode 100644
index 0000000000..61b1277112
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseMixture/noInterfaceCompression/noInterfaceCompression.C
@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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 "noInterfaceCompression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(noInterfaceCompressionNew, 0);
+
+    surfaceInterpolationScheme<scalar>::
+        addMeshFluxConstructorToTable<noInterfaceCompressionNew>
+        addnoInterfaceCompressionScalarMeshFluxConstructorToTable_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/twoPhaseModels/twoPhaseMixture/noInterfaceCompression/noInterfaceCompression.H b/src/twoPhaseModels/twoPhaseMixture/noInterfaceCompression/noInterfaceCompression.H
new file mode 100644
index 0000000000..d84c490ee5
--- /dev/null
+++ b/src/twoPhaseModels/twoPhaseMixture/noInterfaceCompression/noInterfaceCompression.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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::noInterfaceCompression
+
+Description
+    Wrapper scheme to allow VoF solvers to run efficiently without interface
+    compression, e.g. for cavitation simulations.
+
+    Example:
+    \verbatim
+    divSchemes
+    {
+        .
+        .
+        div(phi,alpha)     Gauss noInterfaceCompression vanLeer;
+        .
+        .
+    }
+    \endverbatim
+
+See also
+    Foam::interfaceCompression
+    Foam::vanLeer
+    Foam::vanAlbada
+
+SourceFiles
+    noInterfaceCompression.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noInterfaceCompression_H
+#define noInterfaceCompression_H
+
+#include "surfaceInterpolationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class noInterfaceCompression Declaration
+\*---------------------------------------------------------------------------*/
+
+class noInterfaceCompressionNew
+:
+    public surfaceInterpolationScheme<scalar>
+{
+    // Private member data
+
+        const surfaceScalarField& phi_;
+
+        //- Base scheme to which the compression is applied
+        tmp<surfaceInterpolationScheme<scalar>> tScheme_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("noInterfaceCompression");
+
+
+    // Constructors
+
+        //- Construct from faceFlux and Istream
+        noInterfaceCompressionNew
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<scalar>(mesh),
+            phi_(faceFlux),
+            tScheme_
+            (
+                surfaceInterpolationScheme<scalar>::New(mesh, faceFlux, is)
+            )
+        {}
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors
+        virtual tmp<surfaceScalarField> weights
+        (
+            const volScalarField&
+        ) const
+        {
+            NotImplemented;
+
+            return tmp<surfaceScalarField>(nullptr);
+        }
+
+        //- Return the face-interpolate of the given cell field
+        virtual tmp<surfaceScalarField> interpolate
+        (
+            const volScalarField& vf
+        ) const
+        {
+            return tScheme_().interpolate(vf);
+        }
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const noInterfaceCompressionNew&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleFoam/RAS/propeller/Allrun b/tutorials/incompressible/pimpleFoam/RAS/propeller/Allrun
index 6b31605e57..667cbaf269 100755
--- a/tutorials/incompressible/pimpleFoam/RAS/propeller/Allrun
+++ b/tutorials/incompressible/pimpleFoam/RAS/propeller/Allrun
@@ -5,7 +5,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
 # Mesh if necessary
-if ! cloneMesh ../../../../multiphase/interPhaseChangeFoam/propeller .
+if ! cloneMesh ../../../../multiphase/interFoam/RAS/propeller .
 then
     ./Allmesh
 fi
diff --git a/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/alpha.water.orig b/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/alpha.water.orig
index 8629aca4dd..cbb792affb 100644
--- a/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/alpha.water.orig
+++ b/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/alpha.water.orig
@@ -25,28 +25,34 @@ boundaryField
     {
         type            zeroGradient;
     }
+
     rightWall
     {
         type            zeroGradient;
     }
+
     lowerWall
     {
         type            zeroGradient;
     }
+
     atmosphere
     {
         type            inletOutlet;
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     defaultFaces
     {
         type            empty;
     }
+
     porous_half0
     {
         type            cyclic;
     }
+
     porous_half1
     {
         type            cyclic;
diff --git a/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/p_rgh b/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/p_rgh
index 0f7fb6d620..ae30622c15 100644
--- a/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/p_rgh
+++ b/tutorials/multiphase/interFoam/RAS/damBreakPorousBaffle/0/p_rgh
@@ -64,13 +64,7 @@ boundaryField
     }
     porous_half1
     {
-        type            porousBafflePressure;
-        patchType       cyclic;
-        D               10;
-        I               5;
-        length          0.15;
-        jump            uniform 0;
-        value           uniform 0;
+        type            cyclic;
     }
 }
 
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/0/U b/tutorials/multiphase/interFoam/RAS/propeller/0/U
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/0/U
rename to tutorials/multiphase/interFoam/RAS/propeller/0/U
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/0/alpha.water b/tutorials/multiphase/interFoam/RAS/propeller/0/alpha.water
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/0/alpha.water
rename to tutorials/multiphase/interFoam/RAS/propeller/0/alpha.water
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/0/epsilon b/tutorials/multiphase/interFoam/RAS/propeller/0/epsilon
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/0/epsilon
rename to tutorials/multiphase/interFoam/RAS/propeller/0/epsilon
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/0/k b/tutorials/multiphase/interFoam/RAS/propeller/0/k
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/0/k
rename to tutorials/multiphase/interFoam/RAS/propeller/0/k
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/0/nut b/tutorials/multiphase/interFoam/RAS/propeller/0/nut
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/0/nut
rename to tutorials/multiphase/interFoam/RAS/propeller/0/nut
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/0/p_rgh b/tutorials/multiphase/interFoam/RAS/propeller/0/p_rgh
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/0/p_rgh
rename to tutorials/multiphase/interFoam/RAS/propeller/0/p_rgh
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/Allclean b/tutorials/multiphase/interFoam/RAS/propeller/Allclean
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/Allclean
rename to tutorials/multiphase/interFoam/RAS/propeller/Allclean
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/Allmesh b/tutorials/multiphase/interFoam/RAS/propeller/Allmesh
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/Allmesh
rename to tutorials/multiphase/interFoam/RAS/propeller/Allmesh
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/Allrun b/tutorials/multiphase/interFoam/RAS/propeller/Allrun
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/Allrun
rename to tutorials/multiphase/interFoam/RAS/propeller/Allrun
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/constant/dynamicMeshDict b/tutorials/multiphase/interFoam/RAS/propeller/constant/dynamicMeshDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/constant/dynamicMeshDict
rename to tutorials/multiphase/interFoam/RAS/propeller/constant/dynamicMeshDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/constant/g b/tutorials/multiphase/interFoam/RAS/propeller/constant/g
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/constant/g
rename to tutorials/multiphase/interFoam/RAS/propeller/constant/g
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/constant/momentumTransport b/tutorials/multiphase/interFoam/RAS/propeller/constant/momentumTransport
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/constant/momentumTransport
rename to tutorials/multiphase/interFoam/RAS/propeller/constant/momentumTransport
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/constant/transportProperties b/tutorials/multiphase/interFoam/RAS/propeller/constant/phaseChangeProperties
similarity index 75%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/constant/transportProperties
rename to tutorials/multiphase/interFoam/RAS/propeller/constant/phaseChangeProperties
index 838850c107..50012e31aa 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/propeller/constant/transportProperties
+++ b/tutorials/multiphase/interFoam/RAS/propeller/constant/phaseChangeProperties
@@ -10,33 +10,14 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
-    object      transportProperties;
+    object      phaseChangeProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phases (water vapour);
-
-phaseChangeTwoPhaseMixture SchnerrSauer;
+phaseChangeModel SchnerrSauer;
 
 pSat            2300;   // Saturation pressure
 
-sigma           0.07;
-
-water
-{
-    transportModel  Newtonian;
-    nu              9e-07;
-    rho             1000;
-}
-
-vapour
-{
-    transportModel Newtonian;
-    nu              4.273e-04;
-    rho             0.02308;
-}
-
 KunzCoeffs
 {
     UInf            20.0;
diff --git a/tutorials/multiphase/interFoam/RAS/propeller/constant/transportProperties b/tutorials/multiphase/interFoam/RAS/propeller/constant/transportProperties
new file mode 100644
index 0000000000..9f9784f2d1
--- /dev/null
+++ b/tutorials/multiphase/interFoam/RAS/propeller/constant/transportProperties
@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water vapour);
+
+pSat            2300;   // Saturation pressure
+
+sigma           0.07;
+
+water
+{
+    transportModel  Newtonian;
+    nu              9e-07;
+    rho             1000;
+}
+
+vapour
+{
+    transportModel Newtonian;
+    nu              4.273e-04;
+    rho             0.02308;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/constant/triSurface/README b/tutorials/multiphase/interFoam/RAS/propeller/constant/triSurface/README
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/constant/triSurface/README
rename to tutorials/multiphase/interFoam/RAS/propeller/constant/triSurface/README
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/blockMeshDict b/tutorials/multiphase/interFoam/RAS/propeller/system/blockMeshDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/blockMeshDict
rename to tutorials/multiphase/interFoam/RAS/propeller/system/blockMeshDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/controlDict b/tutorials/multiphase/interFoam/RAS/propeller/system/controlDict
similarity index 96%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/controlDict
rename to tutorials/multiphase/interFoam/RAS/propeller/system/controlDict
index da214a14b2..51b7e5a4ed 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/controlDict
+++ b/tutorials/multiphase/interFoam/RAS/propeller/system/controlDict
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     interPhaseChangeFoam;
+application     interFoam;
 
 startFrom       startTime;
 
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/createBafflesDict b/tutorials/multiphase/interFoam/RAS/propeller/system/createBafflesDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/createBafflesDict
rename to tutorials/multiphase/interFoam/RAS/propeller/system/createBafflesDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/decomposeParDict b/tutorials/multiphase/interFoam/RAS/propeller/system/decomposeParDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/decomposeParDict
rename to tutorials/multiphase/interFoam/RAS/propeller/system/decomposeParDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/fvSchemes b/tutorials/multiphase/interFoam/RAS/propeller/system/fvSchemes
similarity index 95%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/fvSchemes
rename to tutorials/multiphase/interFoam/RAS/propeller/system/fvSchemes
index 94807b787d..db93e285e9 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/RAS/propeller/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
 {
     default             none;
 
-    div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
+    div(phi,alpha)      Gauss noInterfaceCompression vanLeer;
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(phi,k)          Gauss upwind;
     div(phi,epsilon)    Gauss upwind;
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/fvSolution b/tutorials/multiphase/interFoam/RAS/propeller/system/fvSolution
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/fvSolution
rename to tutorials/multiphase/interFoam/RAS/propeller/system/fvSolution
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/snappyHexMeshDict b/tutorials/multiphase/interFoam/RAS/propeller/system/snappyHexMeshDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/snappyHexMeshDict
rename to tutorials/multiphase/interFoam/RAS/propeller/system/snappyHexMeshDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/propeller/system/surfaceFeaturesDict b/tutorials/multiphase/interFoam/RAS/propeller/system/surfaceFeaturesDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/propeller/system/surfaceFeaturesDict
rename to tutorials/multiphase/interFoam/RAS/propeller/system/surfaceFeaturesDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/Phi b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/Phi
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/Phi
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/Phi
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/U.orig b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/U.orig
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/U.orig
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/U.orig
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/alpha.water b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/alpha.water
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/alpha.water
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/alpha.water
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/p_rgh.orig b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/p_rgh.orig
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/0/p_rgh.orig
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/0/p_rgh.orig
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allclean b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/Allclean
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allclean
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/Allclean
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allrun b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/Allrun
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/Allrun
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/Allrun
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/g b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/g
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/g
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/g
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/momentumTransport b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/momentumTransport
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/momentumTransport
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/momentumTransport
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/transportProperties b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/phaseChangeProperties
similarity index 68%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/transportProperties
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/phaseChangeProperties
index dcc15846a1..50012e31aa 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/transportProperties
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/phaseChangeProperties
@@ -10,38 +10,20 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      transportProperties;
+    object      phaseChangeProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phases (water vapour);
-
-phaseChangeTwoPhaseMixture SchnerrSauer;
+phaseChangeModel SchnerrSauer;
 
 pSat            2300;   // Saturation pressure
 
-sigma           0.07;
-
-water
-{
-    transportModel  Newtonian;
-    nu              9e-07;
-    rho             1000;
-}
-
-vapour
-{
-    transportModel  Newtonian;
-    nu              4.273e-04;
-    rho             0.02308;
-}
-
 KunzCoeffs
 {
-    UInf            U20.0;
-    tInf            t0.005; // L = 0.1 m
-    Cc              C1000;
-    Cv              C1000;
+    UInf            20.0;
+    tInf            0.005; // L = 0.1 m
+    Cc              1000;
+    Cv              1000;
 }
 
 MerkleCoeffs
@@ -60,5 +42,4 @@ SchnerrSauerCoeffs
     Cv              1;
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/transportProperties b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/transportProperties
new file mode 100644
index 0000000000..f5a9155e97
--- /dev/null
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/transportProperties
@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water vapour);
+
+sigma           0.07;
+
+water
+{
+    transportModel  Newtonian;
+    nu              9e-07;
+    rho             1000;
+}
+
+vapour
+{
+    transportModel  Newtonian;
+    nu              4.273e-04;
+    rho             0.02308;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/triSurface/bullet.stl.gz b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/triSurface/bullet.stl.gz
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/constant/triSurface/bullet.stl.gz
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/triSurface/bullet.stl.gz
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/blockMeshDict b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/blockMeshDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/blockMeshDict
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/blockMeshDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/controlDict
similarity index 96%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/controlDict
index 77de46b40b..90e5ffcdc4 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/controlDict
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application             interPhaseChangeFoam;
+application             interFoam;
 
 startFrom               latestTime;
 
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/decomposeParDict b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/decomposeParDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/decomposeParDict
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/decomposeParDict
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/fvSchemes
similarity index 95%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/fvSchemes
index adb8bf8afe..f275dec749 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/fvSchemes
@@ -30,7 +30,7 @@ divSchemes
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(phi,omega)      Gauss linearUpwind grad(omega);
     div(phi,k)          Gauss linearUpwind grad(k);
-    div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
+    div(phi,alpha)      Gauss interfaceCompression vanLeer 1;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
 
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/fvSolution
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/fvSolution
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/snappyHexMeshDict b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/snappyHexMeshDict
similarity index 100%
rename from tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/snappyHexMeshDict
rename to tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/snappyHexMeshDict