interFoam: Replaced twoPhaseChangeModel with the VoFCavitation fvModel

Replacing the specific twoPhaseChangeModel with a consistent and general fvModel interface will support not just cavitation using the new VoFCavitation fvModel but also other phase-change and interface manipulation models in the future and is easier to use for case-specific and other user customisation.
2022-09-27 19:28:11 +01:00
parent f412d39174
commit d3df91a5eb
47 changed files with 883 additions and 961 deletions
--- a/applications/solvers/modules/fluid/compressibleVoF/VoFTurbulenceDamping/VoFTurbulenceDamping.C
+++ b/applications/solvers/modules/fluid/compressibleVoF/VoFTurbulenceDamping/VoFTurbulenceDamping.C
@ -68,7 +68,7 @@ Foam::fv::VoFTurbulenceDamping::VoFTurbulenceDamping
            "phaseProperties"
        )
    ),
-    interface_(refCast<const interfaceProperties>(mixture_)),
+    interface_(mixture_),
    turbulence_
    (
        mesh.lookupObject<compressibleMomentumTransportModel>
--- a/applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C
+++ b/applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C
@ -103,7 +103,7 @@ Foam::autoPtr<Foam::relativeVelocityModel> Foam::relativeVelocityModel::New
    const uniformDimensionedVectorField& g
 )
 {
-    word modelType(dict.lookup(typeName));
+    const word modelType(dict.lookup(typeName));

    Info<< "Selecting relative velocity model " << modelType << endl;

--- a/applications/solvers/multiphase/interFoam/Allwclean
+++ b/applications/solvers/multiphase/interFoam/Allwclean
@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # Run from this directory

 wclean libso incompressibleInterPhaseTransportModel
 wclean libso VoFTurbulenceDamping
+wclean libso VoFCavitation

 wclean
 wclean interMixingFoam
--- a/applications/solvers/multiphase/interFoam/Allwmake
+++ b/applications/solvers/multiphase/interFoam/Allwmake
@ -6,6 +6,7 @@ cd ${0%/*} || exit 1    # Run from this directory

 wmake $targetType incompressibleInterPhaseTransportModel
 wmake $targetType VoFTurbulenceDamping
+wmake $targetType VoFCavitation

 wmake $targetType
 wmake $targetType interMixingFoam
--- a/applications/solvers/multiphase/interFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/Make/options
@ -5,7 +5,6 @@ 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)/MomentumTransportModels/phaseIncompressible/lnInclude \
@ -16,7 +15,6 @@ EXE_INC = \

 EXE_LIBS = \
    -limmiscibleIncompressibleTwoPhaseMixture \
-    -ltwoPhaseChangeModels \
    -lmomentumTransportModels \
    -lincompressibleMomentumTransportModels \
    -lincompressibleInterPhaseTransportModel \
--- a/applications/solvers/multiphase/interFoam/UEqn.H
+++ b/applications/solvers/multiphase/interFoam/UEqn.H
@ -4,8 +4,8 @@
      + MRF.DDt(rho, U)
      + turbulence.divDevTau(rho, U)
     ==
-       phaseChange.SU(rho, rhoPhi, U)
-     + fvModels.source(rho, U)
+        // phaseChange.SU(rho, rhoPhi, U)
+       fvModels.source(rho, U)
    );

    UEqn.relax();
--- a/applications/solvers/multiphase/interFoam/VoFCavitation/Make/files
+++ b/applications/solvers/multiphase/interFoam/VoFCavitation/Make/files
@ -0,0 +1,3 @@
+VoFCavitation.C
+
+LIB = $(FOAM_LIBBIN)/libVoFCavitation
--- a/applications/solvers/multiphase/interFoam/VoFCavitation/Make/options
+++ b/applications/solvers/multiphase/interFoam/VoFCavitation/Make/options
@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(LIB_SRC)/physicalProperties/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 \
+    -I$(LIB_SRC)/twoPhaseModels/cavitation/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lphysicalProperties \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lcavitationModels \
+    -lfiniteVolume \
+    -lmeshTools
--- a/applications/solvers/multiphase/interFoam/VoFCavitation/VoFCavitation.C
+++ b/applications/solvers/multiphase/interFoam/VoFCavitation/VoFCavitation.C
@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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 "VoFCavitation.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvmSup.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace fv
+    {
+        defineTypeNameAndDebug(VoFCavitation, 0);
+
+        addToRunTimeSelectionTable
+        (
+            fvModel,
+            VoFCavitation,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::VoFCavitation::VoFCavitation
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    fvModel(sourceName, modelType, dict, mesh),
+
+    mixture_
+    (
+        mesh.lookupObjectRef<immiscibleIncompressibleTwoPhaseMixture>
+        (
+            "phaseProperties"
+        )
+    ),
+
+    cavitation_(cavitationModel::New(dict, mixture_)),
+
+    alphaName_(mixture_.alpha1().name())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::wordList Foam::fv::VoFCavitation::addSupFields() const
+{
+    return {alphaName_, "p_rgh", "U"};
+}
+
+
+void Foam::fv::VoFCavitation::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const word& fieldName
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    if (fieldName == alphaName_)
+    {
+        volScalarField::Internal alphalCoeff
+        (
+            1.0/mixture_.rho1()
+          - mixture_.alpha1()()*(1.0/mixture_.rho1() - 1.0/mixture_.rho2())
+        );
+
+        const Pair<tmp<volScalarField::Internal>> mDotAlphal
+        (
+            cavitation_->mDotAlphal()
+        );
+
+        const volScalarField::Internal vDotcAlphal(alphalCoeff*mDotAlphal[0]());
+        const volScalarField::Internal vDotvAlphal(alphalCoeff*mDotAlphal[1]());
+
+        eqn += fvm::Sp(vDotvAlphal - vDotcAlphal, eqn.psi()) + vDotcAlphal;
+    }
+}
+
+
+void Foam::fv::VoFCavitation::addSup
+(
+    const volScalarField&,
+    fvMatrix<scalar>& eqn,
+    const word& fieldName
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    if (fieldName == "p_rgh")
+    {
+        const volScalarField::Internal& rho =
+            mesh().lookupObject<volScalarField>("rho");
+
+        const volScalarField::Internal& gh =
+            mesh().lookupObject<volScalarField>("gh");
+
+        const dimensionedScalar pCoeff
+        (
+            1.0/mixture_.rho1() - 1.0/mixture_.rho2()
+        );
+
+        const Pair<tmp<volScalarField::Internal>> mDotP
+        (
+            cavitation_->mDotP()
+        );
+
+        const volScalarField::Internal vDotcP(pCoeff*mDotP[0]);
+        const volScalarField::Internal vDotvP(pCoeff*mDotP[1]);
+
+        eqn +=
+            (vDotvP - vDotcP)*(cavitation_->pSat() - rho*gh)
+          - fvm::Sp(vDotvP - vDotcP, eqn.psi());
+    }
+}
+
+
+void Foam::fv::VoFCavitation::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const word& fieldName
+) const
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    if (fieldName == "U")
+    {
+        const surfaceScalarField& rhoPhi =
+            mesh().lookupObject<surfaceScalarField>("rhoPhi");
+
+        eqn += fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), eqn.psi());
+    }
+}
+
+
+void Foam::fv::VoFCavitation::topoChange(const polyTopoChangeMap&)
+{}
+
+
+void Foam::fv::VoFCavitation::mapMesh(const polyMeshMap& map)
+{}
+
+
+void Foam::fv::VoFCavitation::distribute(const polyDistributionMap&)
+{}
+
+
+bool Foam::fv::VoFCavitation::movePoints()
+{
+    return true;
+}
+
+
+void Foam::fv::VoFCavitation::correct()
+{
+    cavitation_->correct();
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/VoFCavitation/VoFCavitation.H
+++ b/applications/solvers/multiphase/interFoam/VoFCavitation/VoFCavitation.H
@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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::fv::VoFCavitation
+
+Description
+    Cavitation fvModel
+
+Usage
+    Example usage:
+    \verbatim
+    VoFCavitation
+    {
+        type    VoFCavitation;
+
+        libs    ("libVoFCavitation.so");
+
+        model   SchnerrSauer;
+
+        KunzCoeffs
+        {
+            pSat    2300;   // Saturation pressure
+
+            UInf    20.0;
+            tInf    0.005; // L = 0.1 m
+            Cc      1000;
+            Cv      1000;
+        }
+
+        MerkleCoeffs
+        {
+            pSat    2300;   // Saturation pressure
+
+            UInf    20.0;
+            tInf    0.005;  // L = 0.1 m
+            Cc      80;
+            Cv      1e-03;
+        }
+
+        SchnerrSauerCoeffs
+        {
+            pSat    2300;   // Saturation pressure
+
+            n       1.6e+13;
+            dNuc    2.0e-06;
+            Cc      1;
+            Cv      1;
+        }
+    }
+    \endverbatim
+
+SourceFiles
+    VoFCavitation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoFCavitation_H
+#define VoFCavitation_H
+
+#include "fvModel.H"
+#include "cavitationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class immiscibleIncompressibleTwoPhaseMixture;
+class interfaceProperties;
+class incompressibleMomentumTransportModel;
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class VoFCavitation Declaration
+\*---------------------------------------------------------------------------*/
+
+class VoFCavitation
+:
+    public fvModel
+{
+    // Private Data
+
+        //- Reference to the mixture properties
+        const immiscibleIncompressibleTwoPhaseMixture& mixture_;
+
+        autoPtr<cavitationModel> cavitation_;
+
+        //- The name of the VoF phase-fraction
+        word alphaName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("VoFCavitation");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        VoFCavitation
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- Disallow default bitwise copy construction
+        VoFCavitation
+        (
+            const VoFCavitation&
+        ) = delete;
+
+
+    // Member Functions
+
+        //- Return the list of fields for which the option adds source term
+        //  to the transport equation
+        virtual wordList addSupFields() const;
+
+        //- Add implicit/explicit contributions to VoF phase-fraction equation
+        virtual void addSup
+        (
+            fvMatrix<scalar>& eqn,
+            const word& fieldName
+        ) const;
+
+        //- Add implicit/explicit contributions to p_rgh equation
+        virtual void addSup
+        (
+            const volScalarField& psi,
+            fvMatrix<scalar>& eqn,
+            const word& fieldName
+        ) const;
+
+        //- Add contribution to U equation
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvMatrix<vector>& eqn,
+            const word& fieldName
+        ) const;
+
+
+        // Mesh changes
+
+            //- Update topology using the given map
+            virtual void topoChange(const polyTopoChangeMap&);
+
+            //- Update from another mesh using the given map
+            virtual void mapMesh(const polyMeshMap&);
+
+            //- Redistribute or update using the given distribution map
+            virtual void distribute(const polyDistributionMap&);
+
+            //- Update for mesh motion
+            virtual bool movePoints();
+
+
+        //- Correct the cavitation model
+        virtual void correct();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const VoFCavitation&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/VoFTurbulenceDamping/Make/options
+++ b/applications/solvers/multiphase/interFoam/VoFTurbulenceDamping/Make/options
@ -3,6 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
    -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
    -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
@ -10,6 +11,7 @@ EXE_INC = \

 LIB_LIBS = \
    -lphysicalProperties \
+    -limmiscibleIncompressibleTwoPhaseMixture \
    -ltwoPhaseMixture \
    -linterfaceProperties \
    -lmomentumTransportModels \
--- a/applications/solvers/multiphase/interFoam/VoFTurbulenceDamping/VoFTurbulenceDamping.C
+++ b/applications/solvers/multiphase/interFoam/VoFTurbulenceDamping/VoFTurbulenceDamping.C
@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/

 #include "VoFTurbulenceDamping.H"
-#include "incompressibleTwoPhaseMixture.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
 #include "interfaceProperties.H"
 #include "incompressibleMomentumTransportModel.H"
 #include "fvMatrix.H"
@ -63,12 +63,12 @@ Foam::fv::VoFTurbulenceDamping::VoFTurbulenceDamping
    delta_("delta", dimLength, dict),
    mixture_
    (
-        mesh.lookupObject<incompressibleTwoPhaseMixture>
+        mesh.lookupObject<immiscibleIncompressibleTwoPhaseMixture>
        (
            "phaseProperties"
        )
    ),
-    interface_(refCast<const interfaceProperties>(mixture_)),
+    interface_(mixture_),
    turbulence_
    (
        mesh.lookupObject<incompressibleMomentumTransportModel>
--- a/applications/solvers/multiphase/interFoam/VoFTurbulenceDamping/VoFTurbulenceDamping.H
+++ b/applications/solvers/multiphase/interFoam/VoFTurbulenceDamping/VoFTurbulenceDamping.H
@ -81,7 +81,7 @@ SourceFiles
 namespace Foam
 {

-class incompressibleTwoPhaseMixture;
+class immiscibleIncompressibleTwoPhaseMixture;
 class interfaceProperties;
 class incompressibleMomentumTransportModel;

@ -98,7 +98,7 @@ class VoFTurbulenceDamping
 {
    // Private Data

-        //- The name of the Lagrangian phase
+        //- The name of the VoF phase
        word phaseName_;

        //- Field name
@ -109,7 +109,7 @@ class VoFTurbulenceDamping
        dimensionedScalar delta_;

        //- Reference to the mixture properties
-        const incompressibleTwoPhaseMixture& mixture_;
+        const immiscibleIncompressibleTwoPhaseMixture& mixture_;

        //- Reference to the interface properties
        const interfaceProperties& interface_;
--- a/applications/solvers/multiphase/interFoam/alphaSuSp.H
+++ b/applications/solvers/multiphase/interFoam/alphaSuSp.H
@ -1,15 +1,15 @@
-// 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;
+tmp<volScalarField::Internal> Su;
+tmp<volScalarField::Internal> Sp;

-if (Sp.valid())
+if (fvModels.addsSupToField(alpha1.name()))
 {
+    // Phase change alpha1 source
+    const fvScalarMatrix alphaSup(fvModels.source(alpha1));
+
+    Su = alphaSup.Su();
+    Sp = alphaSup.Sp();
+
    divU =
    (
        mesh.moving()
--- a/applications/solvers/multiphase/interFoam/createFieldRefs.H
+++ b/applications/solvers/multiphase/interFoam/createFieldRefs.H
@ -3,5 +3,5 @@ volScalarField& alpha2(mixture.alpha2());
 const dimensionedScalar& rho1 = mixture.rho1();
 const dimensionedScalar& rho2 = mixture.rho2();

-twoPhaseChangeModel& phaseChange = phaseChangePtr();
+//twoPhaseChangeModel& phaseChange = phaseChangePtr();
 tmp<volScalarField> rAU;
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@ -34,11 +34,6 @@ volVectorField U
 Info<< "Reading phaseProperties\n" << endl;
 immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);

-autoPtr<twoPhaseChangeModel> phaseChangePtr
-(
-    twoPhaseChangeModel::New(mixture)
-);
-

 // Need to store rho for ddt(rho, U)
 volScalarField rho
--- a/applications/solvers/multiphase/interFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/interFoam/initCorrectPhi.H
@ -1,7 +1,7 @@
 if
 (
    !runTime.restart()
- || isType<twoPhaseChangeModels::noPhaseChange>(phaseChange)
+ || !fvModels.addsSupToField(alpha1.name())
 )
 {
    if (correctPhi)
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@ -40,7 +40,6 @@ Description
 #include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "immiscibleIncompressibleTwoPhaseMixture.H"
-#include "noPhaseChange.H"
 #include "incompressibleInterPhaseTransportModel.H"
 #include "pimpleControl.H"
 #include "pressureReference.H"
@ -99,7 +98,7 @@ int main(int argc, char *argv[])
        if
        (
            correctPhi
-         && !isType<twoPhaseChangeModels::noPhaseChange>(phaseChange)
+         && fvModels.addsSupToField(alpha1.name())
         && mesh.topoChanged()
        )
        {
@ -133,7 +132,7 @@ int main(int argc, char *argv[])
                if
                (
                    correctPhi
-                 && !isType<twoPhaseChangeModels::noPhaseChange>(phaseChange)
+                 && fvModels.addsSupToField(alpha1.name())
                 && !divU.valid()
                )
                {
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@ -39,7 +39,19 @@
    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);

    // Cache the phase change pressure source
-    fvScalarMatrix Sp_rgh(phaseChange.Sp_rgh(rho, gh, p_rgh));
+    fvScalarMatrix Sp_rgh
+    (
+        fvModels.source
+        (
+            volScalarField::New
+            (
+                "1",
+                mesh,
+                dimensionedScalar(dimless/dimPressure, 1)
+            ),
+            p_rgh
+        )
+    );

    while (pimple.correctNonOrthogonal())
    {
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C
@ -263,6 +263,40 @@ DimensionedField<Type, GeoMesh>::clone() const
 }


+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    const Mesh& mesh,
+    const dimensionSet& ds,
+    const Field<Type>& field
+)
+{
+    const bool cacheTmp = mesh.thisDb().cacheTemporaryObject(name);
+
+    return tmp<DimensionedField<Type, GeoMesh>>
+    (
+        new DimensionedField<Type, GeoMesh>
+        (
+            IOobject
+            (
+                name,
+                mesh.thisDb().time().timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                cacheTmp
+            ),
+            mesh,
+            ds,
+            field
+        ),
+        cacheTmp
+    );
+}
+
+
 template<class Type, class GeoMesh>
 Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
 DimensionedField<Type, GeoMesh>::New
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H
@ -229,6 +229,17 @@ public:
        //- Clone
        tmp<DimensionedField<Type, GeoMesh>> clone() const;

+        //- Return a temporary field constructed from name, mesh,
+        //  dimensionSet and field
+        static tmp<DimensionedField<Type, GeoMesh>> New
+        (
+            const word& name,
+            const Mesh& mesh,
+            const dimensionSet&,
+            const Field<Type>&
+
+        );
+
        //- Return a temporary field constructed from name, mesh
        //  and dimensionSet
        static tmp<DimensionedField<Type, GeoMesh>> New
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@ -812,7 +812,7 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
    (
        volScalarField::New
        (
-            "A("+psi_.name()+')',
+            "A(" + psi_.name() + ')',
            psi_.mesh(),
            dimensions_/psi_.dimensions()/dimVol,
            extrapolatedCalculatedFvPatchScalarField::typeName
@ -826,6 +826,42 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
 }


+template<class Type>
+Foam::tmp<Foam::volScalarField::Internal> Foam::fvMatrix<Type>::Su() const
+{
+    tmp<volScalarField::Internal> tSu
+    (
+        volScalarField::Internal::New
+        (
+            "Su(" +psi_.name() + ')',
+            psi_.mesh(),
+            dimensions_/dimVol,
+            -source()/psi_.mesh().V()
+        )
+    );
+
+    return tSu;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::volScalarField::Internal> Foam::fvMatrix<Type>::Sp() const
+{
+    tmp<volScalarField::Internal> tSp
+    (
+        volScalarField::Internal::New
+        (
+            "Sp(" + psi_.name() + ')',
+            psi_.mesh(),
+            dimensions_/psi_.dimensions()/dimVol,
+            diag()/psi_.mesh().V()
+        )
+    );
+
+    return tSp;
+}
+
+
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::fvMatrix<Type>::H() const
@ -834,7 +870,7 @@ Foam::fvMatrix<Type>::H() const
    (
        GeometricField<Type, fvPatchField, volMesh>::New
        (
-            "H("+psi_.name()+')',
+            "H(" + psi_.name() + ')',
            psi_.mesh(),
            dimensions_/dimVol,
            extrapolatedCalculatedFvPatchScalarField::typeName
@ -940,7 +976,7 @@ flux() const
    (
        GeometricField<Type, fvsPatchField, surfaceMesh>::New
        (
-            "flux("+psi_.name()+')',
+            "flux(" + psi_.name() + ')',
            psi_.mesh(),
            dimensions()
        )
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@ -288,6 +288,14 @@ public:

        // Access

+            GeometricField<Type, fvPatchField, volMesh>& psi()
+            {
+                return const_cast<GeometricField<Type, fvPatchField, volMesh>&>
+                (
+                    psi_
+                );
+            }
+
            const GeometricField<Type, fvPatchField, volMesh>& psi() const
            {
                return psi_;
@ -435,6 +443,12 @@ public:
            //- Return the central coefficient
            tmp<volScalarField> A() const;

+            //- Return the explicit source
+            tmp<volScalarField::Internal> Su() const;
+
+            //- Return the implicit source
+            tmp<volScalarField::Internal> Sp() const;
+
            //- Return the H operation source
            tmp<GeometricField<Type, fvPatchField, volMesh>> H() const;

--- a/src/twoPhaseModels/Allwmake
+++ b/src/twoPhaseModels/Allwmake
@ -9,6 +9,6 @@ wmake $targetType interfaceProperties
 wmake $targetType twoPhaseProperties
 wmake $targetType incompressibleTwoPhaseMixture
 wmake $targetType immiscibleIncompressibleTwoPhaseMixture
-wmake $targetType twoPhaseChange
+wmake $targetType cavitation

 #------------------------------------------------------------------------------
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -30,7 +30,7 @@ License

 namespace Foam
 {
-namespace twoPhaseChangeModels
+namespace cavitationModels
 {
    defineTypeNameAndDebug(Kunz, 0);
    addToRunTimeSelectionTable(cavitationModel, Kunz, dictionary);
@ -39,17 +39,18 @@ namespace twoPhaseChangeModels

 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

-Foam::twoPhaseChangeModels::Kunz::Kunz
+Foam::cavitationModels::Kunz::Kunz
 (
+    const dictionary& dict,
    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    cavitationModel(typeName, mixture),
+    cavitationModel(dict, mixture),

-    UInf_("UInf", dimVelocity, twoPhaseChangeModelCoeffs_),
-    tInf_("tInf", dimTime, twoPhaseChangeModelCoeffs_),
-    Cc_("Cc", dimless, twoPhaseChangeModelCoeffs_),
-    Cv_("Cv", dimless, twoPhaseChangeModelCoeffs_),
+    UInf_("UInf", dimVelocity, dict),
+    tInf_("tInf", dimTime, dict),
+    Cc_("Cc", dimless, dict),
+    Cv_("Cv", dimless, dict),

    p0_("0", pSat().dimensions(), 0.0),

@ -62,18 +63,18 @@ Foam::twoPhaseChangeModels::Kunz::Kunz

 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::twoPhaseChangeModels::Kunz::mDotAlphal() const
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::cavitationModels::Kunz::mDotAlphal() const
 {
-    const volScalarField& p =
+    const volScalarField::Internal& p =
        mixture_.U().db().lookupObject<volScalarField>("p");

-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+        min(max(mixture_.alpha1()(), scalar(0)), scalar(1))
    );

-    return Pair<tmp<volScalarField>>
+    return Pair<tmp<volScalarField::Internal>>
    (
        mcCoeff_*sqr(limitedAlpha1)
       *max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),
@ -83,18 +84,18 @@ Foam::twoPhaseChangeModels::Kunz::mDotAlphal() const
 }


-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::twoPhaseChangeModels::Kunz::mDotP() const
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::cavitationModels::Kunz::mDotP() const
 {
-    const volScalarField& p =
+    const volScalarField::Internal& p =
        mixture_.U().db().lookupObject<volScalarField>("p");

-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+        min(max(mixture_.alpha1()(), scalar(0)), scalar(1))
    );

-    return Pair<tmp<volScalarField>>
+    return Pair<tmp<volScalarField::Internal>>
    (
        mcCoeff_*sqr(limitedAlpha1)*(1.0 - limitedAlpha1)
       *pos0(p - pSat())/max(p - pSat(), 0.01*pSat()),
@ -104,22 +105,18 @@ Foam::twoPhaseChangeModels::Kunz::mDotP() const
 }


-void Foam::twoPhaseChangeModels::Kunz::correct()
-{
-    cavitationModel::correct();
-}
+void Foam::cavitationModels::Kunz::correct()
+{}


-bool Foam::twoPhaseChangeModels::Kunz::read()
+bool Foam::cavitationModels::Kunz::read(const dictionary& dict)
 {
-    if (cavitationModel::read())
+    if (cavitationModel::read(dict))
    {
-        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
-
-        twoPhaseChangeModelCoeffs_.lookup("UInf") >> UInf_;
-        twoPhaseChangeModelCoeffs_.lookup("tInf") >> tInf_;
-        twoPhaseChangeModelCoeffs_.lookup("Cc") >> Cc_;
-        twoPhaseChangeModelCoeffs_.lookup("Cv") >> Cv_;
+        dict.lookup("UInf") >> UInf_;
+        dict.lookup("tInf") >> tInf_;
+        dict.lookup("Cc") >> Cc_;
+        dict.lookup("Cv") >> Cv_;

        mcCoeff_ = Cc_*mixture_.rho2()/tInf_;
        mvCoeff_ = Cv_*mixture_.rho2()/(0.5*mixture_.rho1()*sqr(UInf_)*tInf_);
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -22,7 +22,7 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Class
-    Foam::twoPhaseChangeModels::Kunz
+    Foam::cavitationModels::Kunz

 Description
    Kunz cavitation model slightly modified so that the condensation term
@ -55,7 +55,7 @@ SourceFiles

 namespace Foam
 {
-namespace twoPhaseChangeModels
+namespace cavitationModels
 {

 /*---------------------------------------------------------------------------*\
@ -90,6 +90,7 @@ public:
        //- Construct for mixture
        Kunz
        (
+            const dictionary& dict,
            const immiscibleIncompressibleTwoPhaseMixture& mixture
        );

@ -104,23 +105,23 @@ public:
        //- 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;
+        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;

        //- Return the mass condensation and vaporisation rates as coefficients
        //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;

        //- Correct the Kunz phaseChange model
        virtual void correct();

        //- Read the phaseProperties dictionary and update
-        virtual bool read();
+        virtual bool read(const dictionary& dict);
 };


 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-} // End namespace twoPhaseChangeModels
+} // End namespace cavitationModels
 } // End namespace Foam

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--- a/src/twoPhaseModels/cavitation/Make/files
+++ b/src/twoPhaseModels/cavitation/Make/files
@ -0,0 +1,7 @@
+cavitationModel/cavitationModel.C
+cavitationModel/cavitationModelNew.C
+Kunz/Kunz.C
+Merkle/Merkle.C
+SchnerrSauer/SchnerrSauer.C
+
+LIB = $(FOAM_LIBBIN)/libcavitationModels
--- a/src/twoPhaseModels/twoPhaseChange/Make/options
+++ b/src/twoPhaseModels/twoPhaseChange/Make/options
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -30,7 +30,7 @@ License

 namespace Foam
 {
-namespace twoPhaseChangeModels
+namespace cavitationModels
 {
    defineTypeNameAndDebug(Merkle, 0);
    addToRunTimeSelectionTable(cavitationModel, Merkle, dictionary);
@ -39,17 +39,18 @@ namespace twoPhaseChangeModels

 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

-Foam::twoPhaseChangeModels::Merkle::Merkle
+Foam::cavitationModels::Merkle::Merkle
 (
+    const dictionary& dict,
    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    cavitationModel(typeName, mixture),
+    cavitationModel(dict, mixture),

-    UInf_("UInf", dimVelocity, twoPhaseChangeModelCoeffs_),
-    tInf_("tInf", dimTime, twoPhaseChangeModelCoeffs_),
-    Cc_("Cc", dimless, twoPhaseChangeModelCoeffs_),
-    Cv_("Cv", dimless, twoPhaseChangeModelCoeffs_),
+    UInf_("UInf", dimVelocity, dict),
+    tInf_("tInf", dimTime, dict),
+    Cc_("Cc", dimless, dict),
+    Cv_("Cv", dimless, dict),

    p0_("0", pSat().dimensions(), 0.0),

@ -62,13 +63,13 @@ Foam::twoPhaseChangeModels::Merkle::Merkle

 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::twoPhaseChangeModels::Merkle::mDotAlphal() const
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::cavitationModels::Merkle::mDotAlphal() const
 {
-    const volScalarField& p =
+    const volScalarField::Internal& p =
        mixture_.U().db().lookupObject<volScalarField>("p");

-    return Pair<tmp<volScalarField>>
+    return Pair<tmp<volScalarField::Internal>>
    (
        mcCoeff_*max(p - pSat(), p0_),
        mvCoeff_*min(p - pSat(), p0_)
@ -76,18 +77,18 @@ Foam::twoPhaseChangeModels::Merkle::mDotAlphal() const
 }


-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::twoPhaseChangeModels::Merkle::mDotP() const
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::cavitationModels::Merkle::mDotP() const
 {
-    const volScalarField& p =
+    const volScalarField::Internal& p =
        mixture_.U().db().lookupObject<volScalarField>("p");

-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+        min(max(mixture_.alpha1()(), scalar(0)), scalar(1))
    );

-    return Pair<tmp<volScalarField>>
+    return Pair<tmp<volScalarField::Internal>>
    (
        mcCoeff_*(1.0 - limitedAlpha1)*pos0(p - pSat()),
        (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
@ -95,22 +96,18 @@ Foam::twoPhaseChangeModels::Merkle::mDotP() const
 }


-void Foam::twoPhaseChangeModels::Merkle::correct()
-{
-    cavitationModel::correct();
-}
+void Foam::cavitationModels::Merkle::correct()
+{}


-bool Foam::twoPhaseChangeModels::Merkle::read()
+bool Foam::cavitationModels::Merkle::read(const dictionary& dict)
 {
-    if (cavitationModel::read())
+    if (cavitationModel::read(dict))
    {
-        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
-
-        twoPhaseChangeModelCoeffs_.lookup("UInf") >> UInf_;
-        twoPhaseChangeModelCoeffs_.lookup("tInf") >> tInf_;
-        twoPhaseChangeModelCoeffs_.lookup("Cc") >> Cc_;
-        twoPhaseChangeModelCoeffs_.lookup("Cv") >> Cv_;
+        dict.lookup("UInf") >> UInf_;
+        dict.lookup("tInf") >> tInf_;
+        dict.lookup("Cc") >> Cc_;
+        dict.lookup("Cv") >> Cv_;

        mcCoeff_ = Cc_/(0.5*sqr(UInf_)*tInf_);
        mvCoeff_ = Cv_*mixture_.rho1()/(0.5*sqr(UInf_)*tInf_*mixture_.rho2());
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -22,7 +22,7 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Class
-    Foam::twoPhaseChangeModels::Merkle
+    Foam::cavitationModels::Merkle

 Description
    Merkle cavitation model.
@ -49,7 +49,7 @@ SourceFiles

 namespace Foam
 {
-namespace twoPhaseChangeModels
+namespace cavitationModels
 {

 /*---------------------------------------------------------------------------*\
@ -84,6 +84,7 @@ public:
        //- Construct for mixture
        Merkle
        (
+            const dictionary& dict,
            const immiscibleIncompressibleTwoPhaseMixture& mixture
        );

@ -98,23 +99,23 @@ public:
        //- 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;
+        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;

        //- Return the mass condensation and vaporisation rates as coefficients
        //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;

        //- Correct the Merkle phaseChange model
        virtual void correct();

        //- Read the phaseProperties dictionary and update
-        virtual bool read();
+        virtual bool read(const dictionary& dict);
 };


 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-} // End namespace twoPhaseChangeModels
+} // End namespace cavitationModels
 } // End namespace Foam

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -31,7 +31,7 @@ License

 namespace Foam
 {
-namespace twoPhaseChangeModels
+namespace cavitationModels
 {
    defineTypeNameAndDebug(SchnerrSauer, 0);
    addToRunTimeSelectionTable
@ -46,17 +46,18 @@ namespace twoPhaseChangeModels

 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

-Foam::twoPhaseChangeModels::SchnerrSauer::SchnerrSauer
+Foam::cavitationModels::SchnerrSauer::SchnerrSauer
 (
+    const dictionary& dict,
    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    cavitationModel(typeName, mixture),
+    cavitationModel(dict, mixture),

-    n_("n", dimless/dimVolume, twoPhaseChangeModelCoeffs_),
-    dNuc_("dNuc", dimLength, twoPhaseChangeModelCoeffs_),
-    Cc_("Cc", dimless, twoPhaseChangeModelCoeffs_),
-    Cv_("Cv", dimless, twoPhaseChangeModelCoeffs_),
+    n_("n", dimless/dimVolume, dict),
+    dNuc_("dNuc", dimLength, dict),
+    Cc_("Cc", dimless, dict),
+    Cv_("Cv", dimless, dict),

    p0_("0", pSat().dimensions(), 0.0)
 {
@ -66,10 +67,10 @@ Foam::twoPhaseChangeModels::SchnerrSauer::SchnerrSauer

 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

-Foam::tmp<Foam::volScalarField>
-Foam::twoPhaseChangeModels::SchnerrSauer::rRb
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::cavitationModels::SchnerrSauer::rRb
 (
-    const volScalarField& limitedAlpha1
+    const volScalarField::Internal& limitedAlpha1
 ) const
 {
    return pow
@ -82,25 +83,25 @@ Foam::twoPhaseChangeModels::SchnerrSauer::rRb


 Foam::dimensionedScalar
-Foam::twoPhaseChangeModels::SchnerrSauer::alphaNuc() const
+Foam::cavitationModels::SchnerrSauer::alphaNuc() const
 {
    dimensionedScalar Vnuc = n_*constant::mathematical::pi*pow3(dNuc_)/6;
    return Vnuc/(1 + Vnuc);
 }


-Foam::tmp<Foam::volScalarField>
-Foam::twoPhaseChangeModels::SchnerrSauer::pCoeff
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::cavitationModels::SchnerrSauer::pCoeff
 (
-    const volScalarField& p
+    const volScalarField::Internal& p
 ) const
 {
-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+        min(max(mixture_.alpha1()(), scalar(0)), scalar(1))
    );

-    const volScalarField rho
+    const volScalarField::Internal rho
    (
        limitedAlpha1*mixture_.rho1()
      + (scalar(1) - limitedAlpha1)*mixture_.rho2()
@ -112,20 +113,20 @@ Foam::twoPhaseChangeModels::SchnerrSauer::pCoeff
 }


-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::twoPhaseChangeModels::SchnerrSauer::mDotAlphal() const
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::cavitationModels::SchnerrSauer::mDotAlphal() const
 {
-    const volScalarField& p =
+    const volScalarField::Internal& p =
        mixture_.U().db().lookupObject<volScalarField>("p");

-    const volScalarField pCoeff(this->pCoeff(p));
+    const volScalarField::Internal pCoeff(this->pCoeff(p));

-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+        min(max(mixture_.alpha1()(), scalar(0)), scalar(1))
    );

-    return Pair<tmp<volScalarField>>
+    return Pair<tmp<volScalarField::Internal>>
    (
        Cc_*limitedAlpha1*pCoeff*max(p - pSat(), p0_),

@ -134,22 +135,22 @@ Foam::twoPhaseChangeModels::SchnerrSauer::mDotAlphal() const
 }


-Foam::Pair<Foam::tmp<Foam::volScalarField>>
-Foam::twoPhaseChangeModels::SchnerrSauer::mDotP() const
+Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
+Foam::cavitationModels::SchnerrSauer::mDotP() const
 {
-    const volScalarField& p =
+    const volScalarField::Internal& p =
        mixture_.U().db().lookupObject<volScalarField>("p");

-    const volScalarField pCoeff(this->pCoeff(p));
+    const volScalarField::Internal pCoeff(this->pCoeff(p));

-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+        min(max(mixture_.alpha1()(), scalar(0)), scalar(1))
    );

-    const volScalarField apCoeff(limitedAlpha1*pCoeff);
+    const volScalarField::Internal apCoeff(limitedAlpha1*pCoeff);

-    return Pair<tmp<volScalarField>>
+    return Pair<tmp<volScalarField::Internal>>
    (
        Cc_*(1.0 - limitedAlpha1)*pos0(p - pSat())*apCoeff,

@ -158,22 +159,18 @@ Foam::twoPhaseChangeModels::SchnerrSauer::mDotP() const
 }


-void Foam::twoPhaseChangeModels::SchnerrSauer::correct()
-{
-    cavitationModel::correct();
-}
+void Foam::cavitationModels::SchnerrSauer::correct()
+{}


-bool Foam::twoPhaseChangeModels::SchnerrSauer::read()
+bool Foam::cavitationModels::SchnerrSauer::read(const dictionary& dict)
 {
-    if (cavitationModel::read())
+    if (cavitationModel::read(dict))
    {
-        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
-
-        twoPhaseChangeModelCoeffs_.lookup("n") >> n_;
-        twoPhaseChangeModelCoeffs_.lookup("dNuc") >> dNuc_;
-        twoPhaseChangeModelCoeffs_.lookup("Cc") >> Cc_;
-        twoPhaseChangeModelCoeffs_.lookup("Cv") >> Cv_;
+        dict.lookup("n") >> n_;
+        dict.lookup("dNuc") >> dNuc_;
+        dict.lookup("Cc") >> Cc_;
+        dict.lookup("Cv") >> Cv_;

        return true;
    }
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -22,7 +22,7 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Class
-    Foam::twoPhaseChangeModels::SchnerrSauer
+    Foam::cavitationModels::SchnerrSauer

 Description
    SchnerrSauer cavitation model.
@ -49,7 +49,7 @@ SourceFiles

 namespace Foam
 {
-namespace twoPhaseChangeModels
+namespace cavitationModels
 {

 /*---------------------------------------------------------------------------*\
@ -80,10 +80,16 @@ class SchnerrSauer
        dimensionedScalar alphaNuc() const;

        //- Reciprocal bubble radius
-        tmp<volScalarField>rRb(const volScalarField& limitedAlpha1) const;
+        tmp<volScalarField::Internal>rRb
+        (
+            const volScalarField::Internal& limitedAlpha1
+        ) const;

        //- Part of the condensation and vapourisation rates
-        tmp<volScalarField> pCoeff(const volScalarField& p) const;
+        tmp<volScalarField::Internal> pCoeff
+        (
+            const volScalarField::Internal& p
+        ) const;


 public:
@ -97,6 +103,7 @@ public:
        //- Construct for mixture
        SchnerrSauer
        (
+            const dictionary& dict,
            const immiscibleIncompressibleTwoPhaseMixture& mixture
        );

@ -111,23 +118,23 @@ public:
        //- 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;
+        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;

        //- Return the mass condensation and vaporisation rates as coefficients
        //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;

        //- Correct the SchnerrSauer phaseChange model
        virtual void correct();

        //- Read the phaseProperties dictionary and update
-        virtual bool read();
+        virtual bool read(const dictionary& dict);
 };


 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-} // End namespace twoPhaseChangeModels
+} // End namespace cavitationModels
 } // End namespace Foam

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--- a/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.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) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -23,83 +23,38 @@ License

 \*---------------------------------------------------------------------------*/

-#include "twoPhaseChangeModel.H"
+#include "cavitationModel.H"
+#include "fvmSup.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

 namespace Foam
 {
-    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
-{
-    typeIOobject<IOdictionary> io
-    (
-        phaseChangePropertiesName,
-        mixture.U().mesh().time().constant(),
-        mixture.U().mesh(),
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    );
-
-    if (io.headerOk())
-    {
-        io.readOpt() = IOobject::MUST_READ_IF_MODIFIED;
-        return io;
-    }
-    else
-    {
-        io.readOpt() = IOobject::NO_READ;
-        return io;
-    }
+    defineTypeNameAndDebug(cavitationModel, 0);
+    defineRunTimeSelectionTable(cavitationModel, dictionary);
 }


 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

-Foam::twoPhaseChangeModel::twoPhaseChangeModel
+Foam::cavitationModel::cavitationModel
 (
-    const word& type,
+    const dictionary& dict,
    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 :
-    IOdictionary(createIOobject(mixture)),
    mixture_(mixture),
-    twoPhaseChangeModelCoeffs_(optionalSubDict(type + "Coeffs"))
+    pSat_("pSat", dimPressure, dict.lookup("pSat"))
 {}


 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

-void Foam::twoPhaseChangeModel::correct()
-{}
-
-
-bool Foam::twoPhaseChangeModel::read()
+bool Foam::cavitationModel::read(const dictionary& dict)
 {
-    if (regIOobject::read())
-    {
-        twoPhaseChangeModelCoeffs_ = optionalSubDict(type() + "Coeffs");
+    dict.lookup("pSat") >> pSat_;

-        return true;
-    }
-    else
-    {
-        return false;
-    }
+    return true;
 }


--- a/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.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) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -22,7 +22,7 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Class
-    Foam::twoPhaseChangeModels::cavitationModel
+    Foam::cavitationModel

 Description
    Abstract base class for cavitation models
@ -35,24 +35,27 @@ SourceFiles
 #ifndef cavitationModel_H
 #define cavitationModel_H

-#include "twoPhaseChangeModel.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "fvMatricesFwd.H"
+#include "Pair.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 namespace Foam
 {
-namespace twoPhaseChangeModels
-{

 /*---------------------------------------------------------------------------*\
                                 Class cavitationModel
 \*---------------------------------------------------------------------------*/

 class cavitationModel
-:
-    public twoPhaseChangeModel
 {
-    // Private data
+protected:
+
+    // Protected data
+
+        //- Mixture properties
+        const immiscibleIncompressibleTwoPhaseMixture& mixture_;

        //- Saturation vapour pressure
        dimensionedScalar pSat_;
@ -61,7 +64,22 @@ class cavitationModel
 public:

    //- Runtime type information
-    TypeName("cavitation");
+    TypeName("cavitationModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            cavitationModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const immiscibleIncompressibleTwoPhaseMixture& mixture
+            ),
+            (dict, mixture)
+        );


    // Constructors
@ -69,11 +87,19 @@ public:
        //- Construct for mixture
        cavitationModel
        (
-            const word& type,
+            const dictionary& dict,
            const immiscibleIncompressibleTwoPhaseMixture& mixture
        );


+    // Selector
+    static autoPtr<cavitationModel> New
+    (
+        const dictionary& dict,
+        const immiscibleIncompressibleTwoPhaseMixture& mixture
+    );
+
+
    //- Destructor
    virtual ~cavitationModel()
    {}
@ -90,44 +116,22 @@ public:
        //- 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 = 0;
+        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const = 0;

        //- Return the mass condensation and vaporisation rates as coefficients
        //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField>> mDotP() const = 0;
+        virtual Pair<tmp<volScalarField::Internal>> mDotP() const = 0;

-        //- Return the cavitation explicit and implicit sources
-        //  for the phase-fraction equation
-        virtual Pair<tmp<volScalarField::Internal>> Salpha
-        (
-            volScalarField& alpha
-        ) 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;
-
-        //- Return the cavitation source matrix for the momentum equation
-        virtual tmp<fvVectorMatrix> SU
-        (
-            const volScalarField& rho,
-            const surfaceScalarField& rhoPhi,
-            volVectorField& U
-        ) const;
+        //- Correct the cavitation model
+        virtual void correct() = 0;

        //- Read the phaseProperties dictionary and update
-        virtual bool read();
+        virtual bool read(const dictionary& dict) = 0;
 };


 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-} // End namespace twoPhaseChangeModels
 } // End namespace Foam

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--- a/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModelNew.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) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -23,45 +23,19 @@ License

 \*---------------------------------------------------------------------------*/

-#include "noPhaseChange.H"
+#include "cavitationModel.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-Foam::autoPtr<Foam::twoPhaseChangeModel> Foam::twoPhaseChangeModel::New
+Foam::autoPtr<Foam::cavitationModel> Foam::cavitationModel::New
 (
+    const dictionary& dict,
    const immiscibleIncompressibleTwoPhaseMixture& mixture
 )
 {
-    typeIOobject<IOdictionary> twoPhaseChangeModelIO
-    (
-        IOobject
-        (
-            phaseChangePropertiesName,
-            mixture.U().time().constant(),
-            mixture.U().db(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
+    const word modelType(dict.lookup("model"));

-    word modelType(twoPhaseChangeModels::noPhaseChange::typeName);
-
-    if (twoPhaseChangeModelIO.headerOk())
-    {
-        IOdictionary(twoPhaseChangeModelIO).lookup
-        (
-            twoPhaseChangeModel::typeName
-        ) >> modelType;
-    }
-    else
-    {
-        Info<< "No phase change: "
-            << twoPhaseChangeModelIO.name()
-            << " not found" << endl;
-    }
-
-    Info<< "Selecting phaseChange model " << modelType << endl;
+    Info<< "Selecting cavitation model " << modelType << endl;

    dictionaryConstructorTable::iterator cstrIter =
        dictionaryConstructorTablePtr_->find(modelType);
@ -69,14 +43,17 @@ Foam::autoPtr<Foam::twoPhaseChangeModel> Foam::twoPhaseChangeModel::New
    if (cstrIter == dictionaryConstructorTablePtr_->end())
    {
        FatalErrorInFunction
-            << "Unknown " << twoPhaseChangeModel::typeName<< " type "
+            << "Unknown cavitation model "
            << modelType << nl << nl
-            << "Valid  twoPhaseChangeModels are : " << endl
+            << "Valid  cavitationModels are : " << endl
            << dictionaryConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

-    return autoPtr<twoPhaseChangeModel>(cstrIter()(mixture));
+    return autoPtr<cavitationModel>
+    (
+        cstrIter()(dict.optionalSubDict(modelType + "Coeffs"), mixture)
+    );
 }


--- a/src/twoPhaseModels/interfaceProperties/interfaceProperties.C
+++ b/src/twoPhaseModels/interfaceProperties/interfaceProperties.C
@ -31,13 +31,6 @@ License
 #include "fvcGrad.H"
 #include "fvcSnGrad.H"

-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(interfaceProperties, 0);
-}
-

 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

@ -191,7 +184,7 @@ Foam::interfaceProperties::interfaceProperties
    (
        IOobject
        (
-            typedName("interfaceProperties"),
+            "interfaceProperties:K",
            alpha1_.time().timeName(),
            alpha1_.mesh()
        ),
--- a/src/twoPhaseModels/interfaceProperties/interfaceProperties.H
+++ b/src/twoPhaseModels/interfaceProperties/interfaceProperties.H
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -90,9 +90,6 @@ class interfaceProperties

 public:

-    TypeName("interfaceProperties");
-
-
    // Constructors

        //- Construct from volume fraction field gamma and IOdictionary
--- a/src/twoPhaseModels/twoPhaseChange/Make/files
+++ b/src/twoPhaseModels/twoPhaseChange/Make/files
@ -1,9 +0,0 @@
-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
--- a/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.C
+++ b/src/twoPhaseModels/twoPhaseChange/cavitation/cavitationModel/cavitationModel.C
@ -1,127 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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;
-    }
-}
-
-
-// ************************************************************************* //
--- a/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.C
+++ b/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.C
@ -1,129 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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();
-}
-
-
-// ************************************************************************* //
--- a/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.H
+++ b/src/twoPhaseModels/twoPhaseChange/noPhaseChange/noPhaseChange.H
@ -1,128 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 phaseProperties dictionary and update
-        virtual bool read();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace twoPhaseChangeModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
--- a/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.H
+++ b/src/twoPhaseModels/twoPhaseChange/twoPhaseChangeModel/twoPhaseChangeModel.H
@ -1,173 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 phaseProperties dictionary and update
-        virtual bool read() = 0;
-
-
-    // Member Operators
-
-        //- Disallow default bitwise assignment
-        void operator=(const twoPhaseChangeModel&) = delete;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
--- a/tutorials/multiphase/interFoam/RAS/propeller/constant/fvModels
+++ b/tutorials/multiphase/interFoam/RAS/propeller/constant/fvModels
@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      fvModels;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+VoFCavitation
+{
+    type    VoFCavitation;
+
+    libs    ("libVoFCavitation.so");
+
+    model   SchnerrSauer;
+
+    KunzCoeffs
+    {
+        pSat    2300;   // Saturation pressure
+
+        UInf    20.0;
+        tInf    0.005; // L = 0.1 m
+        Cc      1000;
+        Cv      1000;
+    }
+
+    MerkleCoeffs
+    {
+        pSat    2300;   // Saturation pressure
+
+        UInf    20.0;
+        tInf    0.005;  // L = 0.1 m
+        Cc      80;
+        Cv      1e-03;
+    }
+
+    SchnerrSauerCoeffs
+    {
+        pSat    2300;   // Saturation pressure
+
+        n       1.6e+13;
+        dNuc    2.0e-06;
+        Cc      1;
+        Cv      1;
+    }
+}
+
+
+//************************************************************************* //
--- a/tutorials/multiphase/interFoam/RAS/propeller/constant/phaseChangeProperties
+++ b/tutorials/multiphase/interFoam/RAS/propeller/constant/phaseChangeProperties
@ -1,44 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  dev
-     \\/     M anipulation  |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    format      ascii;
-    class       dictionary;
-    object      phaseChangeProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-phaseChangeModel SchnerrSauer;
-
-pSat            2300;   // Saturation pressure
-
-KunzCoeffs
-{
-    UInf            20.0;
-    tInf            0.005; // L = 0.1 m
-    Cc              1000;
-    Cv              1000;
-}
-
-MerkleCoeffs
-{
-    UInf            20.0;
-    tInf            0.005;  // L = 0.1 m
-    Cc              80;
-    Cv              1e-03;
-}
-
-SchnerrSauerCoeffs
-{
-    n               1.6e+13;
-    dNuc            2.0e-06;
-    Cc              1;
-    Cv              1;
-}
-
-// ************************************************************************* //
--- a/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/fvModels
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/fvModels
@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      fvModels;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+VoFCavitation
+{
+    type    VoFCavitation;
+
+    libs    ("libVoFCavitation.so");
+
+    model   SchnerrSauer;
+
+    KunzCoeffs
+    {
+        pSat    2300;   // Saturation pressure
+
+        UInf    20.0;
+        tInf    0.005; // L = 0.1 m
+        Cc      1000;
+        Cv      1000;
+    }
+
+    MerkleCoeffs
+    {
+        pSat    2300;   // Saturation pressure
+
+        UInf    20.0;
+        tInf    0.005;  // L = 0.1 m
+        Cc      80;
+        Cv      1e-03;
+    }
+
+    SchnerrSauerCoeffs
+    {
+        pSat    2300;   // Saturation pressure
+
+        n       1.6e+13;
+        dNuc    2.0e-06;
+        Cc      1;
+        Cv      1;
+    }
+}
+
+
+//************************************************************************* //
--- a/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/phaseChangeProperties
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/phaseChangeProperties
@ -1,44 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  dev
-     \\/     M anipulation  |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    format      ascii;
-    class       dictionary;
-    object      phaseChangeProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-phaseChangeModel SchnerrSauer;
-
-pSat            2300;   // Saturation pressure
-
-KunzCoeffs
-{
-    UInf            20.0;
-    tInf            0.005; // L = 0.1 m
-    Cc              1000;
-    Cv              1000;
-}
-
-MerkleCoeffs
-{
-    UInf            20.0;
-    tInf            0.005;  // L = 0.1 m
-    Cc              80;
-    Cv              1e-03;
-}
-
-SchnerrSauerCoeffs
-{
-    n               1.6e+13;
-    dNuc            2.0e-06;
-    Cc              1;
-    Cv              1;
-}
-
-// ************************************************************************* //
--- a/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/controlDict
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/system/controlDict
@ -15,7 +15,7 @@ FoamFile

 application             interFoam;

-startFrom               latestTime;
+startFrom               startTime;

 startTime               0;