Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

applications/legacy/combustion/PDRFoam/Make/files

@@ -27,6 +27,9 @@ PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
 
 laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
 
+derivedFvPatchFields/activeBaffleVelocity/activeBaffleVelocityFvPatchVectorField.C
+derivedFvPatchFields/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
+
 PDRFoam.C
 
 EXE = $(FOAM_APPBIN)/PDRFoam

applications/modules/multiphaseEuler/phaseSystem/Make/options

@@ -28,6 +28,7 @@ LIB_LIBS = \
     -lcompressibleTwoPhases \
     -lmultiphaseProperties \
     -lphaseCompressibleMomentumTransportModels \
+    -lphaseFluidThermophysicalTransportModels \
     -lphaseFluidMulticomponentThermophysicalTransportModels \
     -lphaseSolidThermophysicalTransportModels \
     -lcombustionModels

applications/modules/multiphaseEuler/populationBalance/Make/options

@@ -20,4 +20,5 @@ EXE_INC = \
     -I$(LIB_SRC)/combustionModels/lnInclude
 
 LIB_LIBS = \
-    -lphaseSystem
+    -lphaseSystem \
+    -lmultiphaseEulerInterfacialModels

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/LehrMilliesMewes/LehrMilliesMewes.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LehrMilliesMewes.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/Liao/Liao.C

@@ -24,8 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Liao.H"
-#include "addToRunTimeSelectionTable.H"
 #include "fvcGrad.H"
+#include "phaseCompressibleMomentumTransportModel.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/Liao/LiaoBase.C

@@ -25,6 +25,7 @@ License
 
 #include "LiaoBase.H"
 #include "fvcGrad.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C

@@ -24,8 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LuoSvendsen.H"
-#include "addToRunTimeSelectionTable.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "linearInterpolationWeights.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

applications/modules/multiphaseEuler/populationBalance/breakupModels/Kusters/Kusters.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Kusters.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/breakupModels/Laakkonen/Laakkonen.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Laakkonen.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/AdachiStuartFokkink/AdachiStuartFokkink.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "AdachiStuartFokkink.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/CoulaloglouTavlarides/CoulaloglouTavlarides.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "CoulaloglouTavlarides.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/LehrMilliesMewesCoalescence/LehrMilliesMewesCoalescence.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LehrMilliesMewesCoalescence.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/LiaoCoalescence/LiaoCoalescence.C

@@ -24,10 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LiaoCoalescence.H"
-#include "addToRunTimeSelectionTable.H"
+#include "phaseSystem.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "fvcGrad.H"
 #include "dragModel.H"
-#include "phaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/Luo/Luo.C

@@ -24,9 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Luo.H"
-#include "addToRunTimeSelectionTable.H"
-#include "dispersedVirtualMassModel.H"
 #include "phaseSystem.H"
+#include "phaseCompressibleMomentumTransportModel.H"
+#include "dispersedVirtualMassModel.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/PrinceBlanch/PrinceBlanch.C

@@ -24,8 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PrinceBlanch.H"
-#include "addToRunTimeSelectionTable.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "fvcGrad.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/turbulentShear/turbulentShear.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "turbulentShear.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

applications/modules/multiphaseEuler/populationBalance/diameterModels/velocityGroup/velocityGroup.C

@@ -143,6 +143,7 @@ Foam::diameterModels::velocityGroup::velocityGroup
 :
     diameterModel(diameterProperties, phase),
     popBalName_(diameterProperties.lookup("populationBalance")),
+    popBalPtr_(nullptr),
     sizeGroups_
     (
         diameterProperties.lookup("sizeGroups"),
@@ -174,6 +175,19 @@ Foam::diameterModels::velocityGroup::~velocityGroup()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::diameterModels::populationBalanceModel&
+Foam::diameterModels::velocityGroup::popBal() const
+{
+    if (popBalPtr_ == nullptr)
+    {
+        popBalPtr_ =
+            &phase().mesh().lookupObject<populationBalanceModel>(popBalName_);
+    }
+
+    return *popBalPtr_;
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::diameterModels::velocityGroup::d() const
 {
     return d_;
@@ -207,8 +221,7 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::velocityGroup::Av() const
 
 void Foam::diameterModels::velocityGroup::correct()
 {
-    const populationBalanceModel& popBal =
-        phase().mesh().lookupObject<populationBalanceModel>(popBalName_);
+    const populationBalanceModel& popBal = this->popBal();
 
     if (!popBal.solveOnFinalIterOnly() || popBal.fluid().pimple().finalIter())
     {

applications/modules/multiphaseEuler/populationBalance/diameterModels/velocityGroup/velocityGroup.H

@@ -76,6 +76,7 @@ namespace diameterModels
 
 // Forward declaration of classes
 class sizeGroup;
+class populationBalanceModel;
 
 /*---------------------------------------------------------------------------*\
                         Class velocityGroup Declaration
@@ -90,6 +91,9 @@ class velocityGroup
         //- Name of the populationBalance this velocityGroup belongs to
         word popBalName_;
 
+        //- Pointer to the populationBalance this velocityGroup belongs to
+        mutable const populationBalanceModel* popBalPtr_;
+
         //- sizeGroups belonging to this velocityGroup
         PtrList<sizeGroup> sizeGroups_;
 
@@ -133,6 +137,9 @@ public:
         //- Return name of populationBalance this velocityGroup belongs to
         inline const word& popBalName() const;
 
+        //- Return the populationBalance this velocityGroup belongs to
+        const populationBalanceModel& popBal() const;
+
         //- Return sizeGroups belonging to this velocityGroup
         inline const PtrList<sizeGroup>& sizeGroups() const;
 

applications/modules/multiphaseEuler/populationBalance/nucleationModels/reactionDriven/reactionDriven.C

@@ -110,7 +110,12 @@ Foam::diameterModels::nucleationModels::reactionDriven::addToNucleationRate
         interface.index(velGroup_.phase()) == 0 ? +1 : -1;
 
     nucleationRate +=
-        popBal_.eta(i, pi/6*pow3(dNuc_))*dmidtfSign*dmidtf/rho/fi.x();
+        popBal_.eta
+        (
+            i,
+            populationBalanceModel::etaBoundsHandling::extrapolate,
+            pi/6*pow3(dNuc_)
+        )*dmidtfSign*dmidtf/rho/fi.x();
 }
 
 

applications/modules/multiphaseEuler/populationBalance/populationBalanceModel/populationBalanceModel.C

@@ -25,6 +25,7 @@ License
 
 #include "populationBalanceModel.H"
 #include "phaseSystem.H"
+#include "phaseCompressibleMomentumTransportModel.H"
 #include "coalescenceModel.H"
 #include "breakupModel.H"
 #include "binaryBreakupModel.H"
@@ -154,7 +155,7 @@ void Foam::diameterModels::populationBalanceModel::birthByCoalescence
 
     for (label i = j; i < sizeGroups().size(); i++)
     {
-        Eta = eta(i, v);
+        Eta = eta(i, etaBoundsHandling::extrapolate, v);
 
         if (Eta.value() == 0) continue;
 
@@ -328,7 +329,7 @@ void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
 
     for (label k = 0; k <= j; k++)
     {
-        Eta = eta(k, v);
+        Eta = eta(k, etaBoundsHandling::extrapolate, v);
 
         if (Eta.value() == 0) continue;
 
@@ -632,6 +633,53 @@ bool Foam::diameterModels::populationBalanceModel::updateSources()
 }
 
 
+template<class EtaType, class VType>
+EtaType Foam::diameterModels::populationBalanceModel::eta
+(
+    const label i,
+    const etaBoundsHandling ebh,
+    const VType& v
+) const
+{
+    const label n = sizeGroups().size();
+
+    static const dimensionedScalar rootVSmallV(dimVolume, rootVSmall);
+    static const dimensionedScalar rootVGreatV(dimVolume, rootVGreat);
+    static const dimensionedScalar z(dimless, scalar(0));
+
+    const dimensionedScalar& x0 =
+        i > 0 ? sizeGroups()[i - 1].x() : rootVSmallV;
+
+    const dimensionedScalar& xi = sizeGroups()[i].x();
+
+    switch (ebh)
+    {
+        case etaBoundsHandling::extrapolate:
+        {
+            const dimensionedScalar& x1 =
+                i < n - 1 ? sizeGroups()[i + 1].x() : rootVSmallV;
+
+            return max(min((v - x0)/(xi - x0), (x1 - v)/(x1 - xi)), z);
+        }
+
+        case etaBoundsHandling::clamp:
+        {
+            const dimensionedScalar& x1 =
+                i < n - 1 ? sizeGroups()[i + 1].x() : rootVGreatV;
+
+            const VType vClip
+            (
+                min(max(v, sizeGroups().first().x()), sizeGroups().last().x())
+            );
+
+            return max(min((vClip - x0)/(xi - x0), (x1 - vClip)/(x1 - xi)), z);
+        }
+    }
+
+    return NaN*v/v;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::diameterModels::populationBalanceModel::populationBalanceModel
@@ -1009,43 +1057,26 @@ bool Foam::diameterModels::populationBalanceModel::writeData(Ostream& os) const
 }
 
 
-const Foam::dimensionedScalar
-Foam::diameterModels::populationBalanceModel::eta
+Foam::dimensionedScalar Foam::diameterModels::populationBalanceModel::eta
 (
     const label i,
+    const etaBoundsHandling ebh,
     const dimensionedScalar& v
 ) const
 {
-    const dimensionedScalar& x0 = sizeGroups()[0].x();
-    const dimensionedScalar& xi = sizeGroups()[i].x();
-    const dimensionedScalar& xm = sizeGroups().last().x();
-    dimensionedScalar lowerBoundary(x0);
-    dimensionedScalar upperBoundary(xm);
+    return eta<dimensionedScalar>(i, ebh, v);
+}
 
-    if (i != 0) lowerBoundary = sizeGroups()[i-1].x();
 
-    if (i != sizeGroups().size() - 1) upperBoundary = sizeGroups()[i+1].x();
-
-    if ((i == 0 && v < x0) || (i == sizeGroups().size() - 1 && v > xm))
-    {
-        return v/xi;
-    }
-    else if (v < lowerBoundary || v > upperBoundary)
-    {
-        return 0;
-    }
-    else if (v.value() == xi.value())
-    {
-        return 1;
-    }
-    else if (v > xi)
-    {
-        return (upperBoundary - v)/(upperBoundary - xi);
-    }
-    else
-    {
-        return (v - lowerBoundary)/(xi - lowerBoundary);
-    }
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::diameterModels::populationBalanceModel::eta
+(
+    const label i,
+    const etaBoundsHandling ebh,
+    const volScalarField::Internal& v
+) const
+{
+    return eta<tmp<volScalarField::Internal>>(i, ebh, v);
 }
 
 

applications/modules/multiphaseEuler/populationBalance/populationBalanceModel/populationBalanceModel.H

@@ -206,7 +206,6 @@ SourceFiles
 
 #include "sizeGroup.H"
 #include "phaseSystem.H"
-#include "phaseCompressibleMomentumTransportModel.H"
 #include "HashPtrTable.H"
 #include "Pair.H"
 
@@ -214,6 +213,9 @@ SourceFiles
 
 namespace Foam
 {
+
+class phaseCompressibleMomentumTransportModel;
+
 namespace diameterModels
 {
 
@@ -369,6 +371,25 @@ public:
                 }
         };
 
+    // Public Enumerations
+
+        //- Enumeration for specifying how to allocate out of bounds of the
+        //  population balance size-groups.
+        //
+        //  "Clamp" will create an out-of-bounds allocation coefficient of
+        //  exactly one. This partitions unity across all size-space.
+        //
+        //  "Extrapolate" will create an out-of-bounds allocation coefficient
+        //  in proportion to the ratio between the given size and the nearest
+        //  size-group size. This does not partition unity outside the range of
+        //  the size-groups.
+        //
+        enum class etaBoundsHandling
+        {
+            clamp,
+            extrapolate
+        };
+
 
 private:
 
@@ -509,6 +530,16 @@ private:
         //- Return the interval at which the sources are updated
         inline label sourceUpdateInterval() const;
 
+        //- Calculate and return the allocation coefficient
+        template<class EtaType, class VType>
+        EtaType eta
+        (
+            const label i,
+            const etaBoundsHandling ebh,
+            const VType& v
+        ) const;
+
+
 public:
 
     //- Runtime type information
@@ -602,12 +633,21 @@ public:
         inline const volVectorField& U() const;
 
         //- Return allocation coefficient
-        const dimensionedScalar eta
+        dimensionedScalar eta
         (
             const label i,
+            const etaBoundsHandling ebh,
             const dimensionedScalar& v
         ) const;
 
+        //- Return allocation coefficient
+        tmp<volScalarField::Internal> eta
+        (
+            const label i,
+            const etaBoundsHandling ebh,
+            const volScalarField::Internal& v
+        ) const;
+
         //- Return the surface tension coefficient between a given dispersed
         //  and the continuous phase
         const tmp<volScalarField> sigmaWithContinuousPhase
@@ -616,7 +656,7 @@ public:
         ) const;
 
         //- Return reference to momentumTransport model of the continuous phase
-        const phaseCompressible::momentumTransportModel&
+        const phaseCompressibleMomentumTransportModel&
             continuousTurbulence() const;
 
         //- Solve the population balance equation

applications/modules/multiphaseEuler/thermophysicalTransportModels/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C

@@ -120,8 +120,7 @@ void Foam::diameterModels::nucleationModels::wallBoiling::precompute()
 }
 
 
-void
-Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
+void Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
 (
     volScalarField& nucleationRate,
     const label i
@@ -144,35 +143,32 @@ Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
 
     forAll(alphatBf, patchi)
     {
-        if
-        (
-            isA<alphatWallBoilingWallFunction>(alphatBf[patchi])
-        )
+        if (!isA<alphatWallBoilingWallFunction>(alphatBf[patchi])) continue;
+
+        const alphatWallBoilingWallFunction& alphatw =
+            refCast<const alphatWallBoilingWallFunction>(alphatBf[patchi]);
+
+        const scalarField& dmdt = alphatw.dmdtf();
+        const scalarField& dDep = alphatw.dDeparture();
+
+        const labelList& faceCells = alphatw.patch().faceCells();
+
+        dimensionedScalar unitLength("unitLength", dimLength, 1);
+
+        forAll(alphatw, facei)
         {
-            const alphatWallBoilingWallFunction& alphatw =
-                refCast<const alphatWallBoilingWallFunction>(alphatBf[patchi]);
-
-            const scalarField& dmdt = alphatw.dmdtf();
-            const scalarField& dDep = alphatw.dDeparture();
-
-            const labelList& faceCells = alphatw.patch().faceCells();
-
-            dimensionedScalar unitLength("unitLength", dimLength, 1);
-
-            forAll(alphatw, facei)
+            if (dmdt[facei] > small)
             {
-                if (dmdt[facei] > small)
-                {
-                    const label faceCelli = faceCells[facei];
+                const label faceCelli = faceCells[facei];
 
-                    nucleationRate[faceCelli] +=
-                        popBal_.eta
-                        (
-                            i,
-                            fi.x()/pow3(fi.dSph())*pow3(dDep[facei]*unitLength)
-                        ).value()
-                       *dmdt[facei]/rho[faceCelli]/fi.x().value();
-                }
+                nucleationRate[faceCelli] +=
+                    popBal_.eta
+                    (
+                        i,
+                        populationBalanceModel::etaBoundsHandling::extrapolate,
+                        fi.x()/pow3(fi.dSph())*pow3(dDep[facei]*unitLength)
+                    ).value()
+                   *dmdt[facei]/rho[faceCelli]/fi.x().value();
             }
         }
     }

etc/caseDicts/postProcessing/graphs/graphLayerAverage

@@ -10,18 +10,22 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-patches         (<patchNames>);  // Patches from which layers extrude
+patches         (<patchNames>);       // Patches from which layers extrude
 
-zones           ();              // Zones from which layers extrude
+zones           ();                   // Zones from which layers extrude
 
-axis            distance;        // The independent variable of the graph. Can
-                                 // be "x", "y", "z", "xyz" (all coordinates
-                                 // written out), or "distance" (from the start
-                                 // point).
+axis            distance;             // The independent variable of the graph.
+                                      // Can be "x", "y", "z", "xyz" (all
+                                      // coordinates written out), or
+                                      // "distance" (from the start point).
 
-symmetric       false;           // Are the layers symmetric about the centre?
+symmetric       false;                // Are the layers symmetric about the
+                                      // centre?
 
-fields          (<fieldsNames>); // Fields to plot
+fields          (<fieldsNames>);      // Fields to average
+
+//weightField   <weightFieldName>;    // Field or fields with which to weight
+//weightFields  (<weightFieldNames>); // the average
 
 #includeEtc "caseDicts/postProcessing/graphs/graphLayerAverage.cfg"
 

src/finiteVolume/Make/files

@@ -166,8 +166,6 @@ $(constraintFvPatchFields)/nonConformalError/nonConformalErrorFvPatchFields.C
 $(constraintFvPatchFields)/nonConformalProcessorCyclic/nonConformalProcessorCyclicFvPatchFields.C
 
 derivedFvPatchFields = $(fvPatchFields)/derived
-$(derivedFvPatchFields)/activeBaffleVelocity/activeBaffleVelocityFvPatchVectorField.C
-$(derivedFvPatchFields)/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/advective/advectiveFvPatchFields.C
 $(derivedFvPatchFields)/codedFixedValue/codedFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/codedMixed/codedMixedFvPatchFields.C

src/functionObjects/field/histogram/histogram.C

@@ -162,7 +162,7 @@ bool Foam::functionObjects::histogram::write()
                 file_.baseTimeDir(),
                 typeName,
                 coordSet(true, fieldName_, xBin),
-                field.name(),
+                "v/vTotal",
                 volFrac
             );
         }

src/functionObjects/field/layerAverage/layerAverage.C

@@ -126,10 +126,10 @@ void Foam::functionObjects::layerAverage::calcLayers()
     }
 
     // Sum number of entries per layer
-    layerCount_ = sum(scalarField(mesh_.nCells(), 1));
+    layerVolume_ = sum<scalar>(mesh_.V());
 
     // Average the cell centres
-    layerCentre_ = sum(mesh_.cellCentres())/layerCount_;
+    layerCentre_ = sum<vector>(mesh_.V()*mesh_.C())/layerVolume_;
 
     // If symmetric, keep only half of the coordinates
     if (symmetric_)
@@ -139,6 +139,34 @@ void Foam::functionObjects::layerAverage::calcLayers()
 }
 
 
+Foam::tmp<Foam::VolInternalField<Foam::scalar>>
+Foam::functionObjects::layerAverage::weight() const
+{
+    if (weightFields_.empty())
+    {
+        return tmp<VolInternalField<scalar>>();
+    }
+
+    tmp<VolInternalField<scalar>> tresult =
+        VolInternalField<scalar>::New
+        (
+            "weight",
+            mesh_,
+            dimensionedScalar(dimless, scalar(1))
+        );
+
+    forAll(weightFields_, i)
+    {
+        const VolInternalField<scalar>& weightField =
+            mesh_.lookupObject<VolInternalField<scalar>>(weightFields_[i]);
+
+        tresult.ref() *= weightField;
+    }
+
+    return tresult;
+}
+
+
 template<>
 Foam::vector
 Foam::functionObjects::layerAverage::symmetricCoeff<Foam::vector>() const
@@ -245,6 +273,13 @@ bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
 
     fields_ = dict.lookup<wordList>("fields");
 
+    weightFields_ =
+        dict.found("weightFields")
+      ? dict.lookup<wordList>("weightFields")
+      : dict.found("weightField")
+      ? wordList(1, dict.lookup<word>("weightField"))
+      : wordList();
+
     formatter_ = setWriter::New(dict.lookup("setFormat"), dict);
 
     calcLayers();
@@ -255,7 +290,9 @@ bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
 
 Foam::wordList Foam::functionObjects::layerAverage::fields() const
 {
-    return fields_;
+    wordList result(fields_);
+    result.append(weightFields_);
+    return result;
 }
 
 
@@ -267,6 +304,13 @@ bool Foam::functionObjects::layerAverage::execute()
 
 bool Foam::functionObjects::layerAverage::write()
 {
+    // Get the weights
+    tmp<VolInternalField<scalar>> weight(this->weight());
+    tmp<Field<scalar>> layerWeight
+    (
+        weight.valid() ? sum<scalar>(mesh_.V()*weight) : tmp<Field<scalar>>()
+    );
+
     // Create list of available fields
     wordList fieldNames;
     forAll(fields_, fieldi)
@@ -275,7 +319,7 @@ bool Foam::functionObjects::layerAverage::write()
         (
             false
             #define FoundTypeField(Type, nullArg) \
-              || foundObject<VolField<Type>>(fields_[fieldi])
+              || foundObject<VolInternalField<Type>>(fields_[fieldi])
             FOR_ALL_FIELD_TYPES(FoundTypeField)
             #undef FoundTypeField
         )
@@ -295,16 +339,18 @@ bool Foam::functionObjects::layerAverage::write()
     #undef DeclareTypeValueSets
     forAll(fieldNames, fieldi)
     {
+        const word& fieldName = fieldNames[fieldi];
+
         #define CollapseTypeFields(Type, nullArg)                           \
-            if (mesh_.foundObject<VolField<Type>>(fieldNames[fieldi]))      \
+            if (mesh_.foundObject<VolInternalField<Type>>(fieldName))       \
             {                                                               \
-                const VolField<Type>& field =                               \
-                    mesh_.lookupObject<VolField<Type>>(fieldNames[fieldi]); \
+                const VolInternalField<Type>& field =                       \
+                    mesh_.lookupObject<VolInternalField<Type>>(fieldName);  \
                                                                             \
                 Type##ValueSets.set                                         \
                 (                                                           \
                     fieldi,                                                 \
-                    average(field.primitiveField())                         \
+                    average<Type>(weight, layerWeight, field)               \
                 );                                                          \
             }
         FOR_ALL_FIELD_TYPES(CollapseTypeFields);

src/functionObjects/field/layerAverage/layerAverage.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,9 @@ Usage
         axis          | Component of the position to plot against | yes |
         symmetric     | Is the geometry symmetric around the centre layer? \
                                                    | no         | false
-        field         | Fields to average and plot | yes        |
+        fields        | Fields to average and plot | yes        |
+        weightField   | Field with which to weight the average | no | none
+        weightFields  | Fields with which to weight the average | no | ()
     \endtable
 
 SourceFiles
@@ -74,6 +76,7 @@ SourceFiles
 #include "fvMeshFunctionObject.H"
 #include "setWriter.H"
 #include "boolList.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -113,15 +116,18 @@ class layerAverage
         //- Per cell the global layer
         labelList cellLayer_;
 
-        //- Per global layer the number of cells
-        scalarField layerCount_;
+        //- Per global layer the volume
+        scalarField layerVolume_;
 
         //- The average centre of each layer
         pointField layerCentre_;
 
-        //- Fields to sample
+        //- Fields to average
         wordList fields_;
 
+        //- Fields with which to weight the averages
+        wordList weightFields_;
+
         //- Set formatter
         autoPtr<setWriter> formatter_;
 
@@ -131,17 +137,26 @@ class layerAverage
         //- Create the layer information, the sort map, and the scalar axis
         void calcLayers();
 
+        //- Calculate and return the weight field, or a null pointer if there
+        //  are no weight fields
+        tmp<VolInternalField<scalar>> weight() const;
+
         //- Return the coefficient to multiply onto symmetric values
         template<class T>
         T symmetricCoeff() const;
 
         //- Sum field per layer
         template<class T>
-        Field<T> sum(const Field<T>& cellField) const;
+        tmp<Field<T>> sum(const VolInternalField<T>& cellField) const;
 
         //- Average a field per layer
         template<class T>
-        Field<T> average(const Field<T>& cellField) const;
+        tmp<Field<T>> average
+        (
+            const tmp<VolInternalField<scalar>>& cellWeight,
+            const tmp<Field<scalar>>& layerWeight,
+            const VolInternalField<T>& cellField
+        ) const;
 
 
 public:

src/functionObjects/field/layerAverage/layerAverageTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "layerAverage.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -35,12 +36,13 @@ T Foam::functionObjects::layerAverage::symmetricCoeff() const
 
 
 template<class T>
-Foam::Field<T> Foam::functionObjects::layerAverage::sum
+Foam::tmp<Foam::Field<T>> Foam::functionObjects::layerAverage::sum
 (
-    const Field<T>& cellField
+    const VolInternalField<T>& cellField
 ) const
 {
-    Field<T> layerField(nLayers_, Zero);
+    tmp<Field<T>> tlayerField(new Field<T>(nLayers_, Zero));
+    Field<T>& layerField = tlayerField.ref();
 
     forAll(cellLayer_, celli)
     {
@@ -53,22 +55,30 @@ Foam::Field<T> Foam::functionObjects::layerAverage::sum
     Pstream::listCombineGather(layerField, plusEqOp<T>());
     Pstream::listCombineScatter(layerField);
 
-    return layerField;
+    return tlayerField;
 }
 
 
 template<class T>
-Foam::Field<T> Foam::functionObjects::layerAverage::average
+Foam::tmp<Foam::Field<T>> Foam::functionObjects::layerAverage::average
 (
-    const Field<T>& cellField
+    const tmp<VolInternalField<scalar>>& cellWeight,
+    const tmp<Field<scalar>>& layerWeight,
+    const VolInternalField<T>& cellField
 ) const
 {
-    // Sum and average
-    Field<T> layerField(sum(cellField)/layerCount_);
+    tmp<Field<T>> tlayerField
+    (
+        cellWeight.valid()
+      ? sum<T>(mesh_.V()*cellWeight*cellField)/layerWeight
+      : sum<T>(mesh_.V()*cellField)/layerVolume_
+    );
 
     // Handle symmetry
     if (symmetric_)
     {
+        Field<T>& layerField = tlayerField.ref();
+
         const T coeff = symmetricCoeff<T>();
 
         for (label i=0; i<nLayers_/2; i++)
@@ -82,7 +92,7 @@ Foam::Field<T> Foam::functionObjects::layerAverage::average
         layerField.setSize(nLayers_/2);
     }
 
-    return layerField;
+    return tlayerField;
 }
 
 

src/fvModels/Make/files

@@ -32,6 +32,10 @@ derived/zeroDimensionalMassSource/zeroDimensionalMassSourceBase.C
 derived/zeroDimensionalMassSource/zeroDimensionalMassSource.C
 derived/massTransfer/massTransferBase.C
 derived/massTransfer/coefficientMassTransfer.C
+derived/phaseChange/phaseChangeBase.C
+derived/phaseChange/singleComponentPhaseChangeBase.C
+derived/phaseChange/multicomponentPhaseChangeBase.C
+derived/phaseChange/coefficientPhaseChange.C
 
 interRegion/interRegionModel/interRegionModel.C
 interRegion/interRegionPorosityForce/interRegionPorosityForce.C

src/fvModels/Make/options

@@ -4,8 +4,9 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \

src/fvModels/derived/massTransfer/coefficientMassTransfer.C

@@ -104,7 +104,7 @@ void Foam::fv::coefficientMassTransfer::addSup
     }
     else
     {
-        massTransferBase::addSupType(alpha, eqn);
+        massTransferBase::addSup(alpha, eqn);
     }
 }
 
@@ -123,10 +123,7 @@ void Foam::fv::coefficientMassTransfer::addSup
         const volScalarField& alpha1 =
             mesh().lookupObject<volScalarField>(alphaNames().first());
 
-        const volScalarField::Internal SByAlpha1
-        (
-            C_*mag(fvc::grad(alpha1))
-        );
+        const volScalarField::Internal SByAlpha1(C_*mag(fvc::grad(alpha1)));
 
         if (i == 0)
         {
@@ -139,7 +136,7 @@ void Foam::fv::coefficientMassTransfer::addSup
     }
     else
     {
-        massTransferBase::addSupType(alpha, rho, eqn);
+        massTransferBase::addSup(alpha, rho, eqn);
     }
 }
 

src/fvModels/derived/phaseChange/ThermoRefPair.H

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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::ThermoRefPair
+
+Description
+    Class containing a pair of thermo references. Handles down-casting to more
+    specific thermo types by constructing one pair from another (e.g.,
+    constructing a multicomponentThermo reference pair from a basicThermo
+    pair). Tracks validity of the references.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ThermoRefPair_H
+#define ThermoRefPair_H
+
+#include "physicalProperties.H"
+#include "objectRegistry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class ThermoRefPair Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ThermoType>
+class ThermoRefPair
+{
+    // Private Member Data
+
+        //- Validity flags
+        const Pair<bool> valid_;
+
+        //- The first thermo
+        const ThermoType& first_;
+
+        //- The second thermo
+        const ThermoType& second_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from a database and phase names
+        ThermoRefPair
+        (
+            const objectRegistry& db,
+            const Pair<word>& phaseNames
+        )
+        :
+            valid_(true, true),
+            first_
+            (
+                db.lookupObject<ThermoType>
+                (
+                    IOobject::groupName
+                    (
+                        physicalProperties::typeName,
+                        phaseNames.first()
+                    )
+                )
+            ),
+            second_
+            (
+                db.lookupObject<ThermoType>
+                (
+                    IOobject::groupName
+                    (
+                        physicalProperties::typeName,
+                        phaseNames.second()
+                    )
+                )
+            )
+        {}
+
+        //- Construct by casting a more primitive thermo type
+        template<class BasicThermoType>
+        ThermoRefPair(const ThermoRefPair<BasicThermoType> basicThermos)
+        :
+            valid_
+            (
+                isA<ThermoType>(basicThermos.first()),
+                isA<ThermoType>(basicThermos.second())
+            ),
+            first_
+            (
+                valid_.first()
+              ? refCast<const ThermoType>(basicThermos.first())
+              : NullObjectRef<ThermoType>()
+            ),
+            second_
+            (
+                valid_.second()
+              ? refCast<const ThermoType>(basicThermos.second())
+              : NullObjectRef<ThermoType>()
+            )
+        {}
+
+
+    // Member Functions
+
+        //- Access the validity flags
+        const Pair<bool>& valid() const
+        {
+            return valid_;
+        }
+
+        //- Access the first thermo
+        const ThermoType& first() const
+        {
+            return first_;
+        }
+
+        //- Access the second thermo
+        const ThermoType& second() const
+        {
+            return second_;
+        }
+
+
+    // Member Operators
+
+        //- Access a thermo by index
+        const ThermoType& operator[](const label i) const
+        {
+            return
+                i == 0 ? first()
+              : i == 1 ? second()
+              : NullObjectRef<ThermoType>();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/coefficientPhaseChange.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "coefficientPhaseChange.H"
+#include "fvcGrad.H"
+#include "multicomponentThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(coefficientPhaseChange, 0);
+    addToRunTimeSelectionTable(fvModel, coefficientPhaseChange, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fv::coefficientPhaseChange::readCoeffs()
+{
+    C_.read(coeffs());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::coefficientPhaseChange::coefficientPhaseChange
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    singleComponentPhaseChangeBase
+    (
+        name,
+        modelType,
+        mesh,
+        dict,
+        Pair<bool>(false, false)
+    ),
+    C_("C", dimMass/dimArea/dimTime, NaN)
+{
+    readCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::fv::coefficientPhaseChange::mDot() const
+{
+    const volScalarField& alpha1 =
+        mesh().lookupObject<volScalarField>(alphaNames().first());
+
+    tmp<volScalarField::Internal> tmDot =
+        C_*alpha1()*mag(fvc::grad(alpha1))()();
+
+    if (specieis().first() != -1)
+    {
+        tmDot.ref() *= specieThermos().first().Y()[specieis().first()];
+    }
+
+    return tmDot;
+}
+
+
+void Foam::fv::coefficientPhaseChange::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn
+) const
+{
+    const label i = index(alphaNames(), eqn.psi().name());
+
+    if (i != -1)
+    {
+        const volScalarField& alpha1 =
+            mesh().lookupObject<volScalarField>(alphaNames().first());
+
+        volScalarField::Internal mDotByAlpha1(C_*mag(fvc::grad(alpha1)));
+
+        if (specieis().first() != -1)
+        {
+            mDotByAlpha1 *= specieThermos().first().Y()[specieis().first()];
+        }
+
+        if (i == 0)
+        {
+            eqn -= fvm::Sp(mDotByAlpha1, eqn.psi());
+        }
+        else
+        {
+            eqn +=
+                mDotByAlpha1*alpha1
+              - correction(fvm::Sp(mDotByAlpha1, eqn.psi()));
+        }
+    }
+    else
+    {
+        phaseChangeBase::addSup(alpha, rho, eqn);
+    }
+}
+
+
+void Foam::fv::coefficientPhaseChange::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volScalarField& Yi,
+    fvMatrix<scalar>& eqn
+) const
+{
+    const label i = index(alphaNames(), eqn.psi().name());
+
+    if (i == 0 && specieis().first() != -1 && Yi.member() == specie())
+    {
+        eqn -= fvm::Sp(C_*alpha()*mag(fvc::grad(alpha))()(), Yi);
+    }
+    else
+    {
+        singleComponentPhaseChangeBase::addSup(alpha, rho, Yi, eqn);
+    }
+}
+
+
+bool Foam::fv::coefficientPhaseChange::read(const dictionary& dict)
+{
+    if (singleComponentPhaseChangeBase::read(dict))
+    {
+        readCoeffs();
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/coefficientPhaseChange.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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::coefficientPhaseChange
+
+Description
+    This simple model generates a phase change between two phases calculated
+    from the following expression:
+
+        \f[
+            \dot{m}/V = C \alpha \grad \alpha
+        \f]
+
+    \vartable
+        \dot{m}/V | mass transfer rate per unit volume
+        C         | coefficient
+        \alpha    | volume fraction of the source phase
+    \endvartable
+
+Usage
+    Example usage:
+    \verbatim
+    coefficientPhaseChange
+    {
+        type            coefficientPhaseChange;
+
+        phases          (liquid vapour);
+
+        C               [kg/m^2/s] 0.1;
+    }
+    \endverbatim
+
+SourceFiles
+    coefficientPhaseChange.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef coefficientPhaseChange_H
+#define coefficientPhaseChange_H
+
+#include "singleComponentPhaseChangeBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class coefficientPhaseChange Declaration
+\*---------------------------------------------------------------------------*/
+
+class coefficientPhaseChange
+:
+    public singleComponentPhaseChangeBase
+{
+private:
+
+    // Private Data
+
+        //- Phase change coefficient. Multiplies the estimated interfacial
+        //  area density return the phase change rate. Units of [kg/m^2/s].
+        dimensionedScalar C_;
+
+
+    // Private Member Functions
+
+        //- Non-virtual read
+        void readCoeffs();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("coefficientPhaseChange");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        coefficientPhaseChange
+        (
+            const word& name,
+            const word& modelType,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        // Sources
+
+            //- Return the phase change rate
+            virtual tmp<DimensionedField<scalar, volMesh>> mDot() const;
+
+            //- Override the compressible continuity equation to add
+            //  linearisation w.r.t alpha
+            void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn
+            ) const;
+
+            //- Override the compressible mass fraction equation to add
+            //  linearisation w.r.t the mass fraction
+            void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& Yi,
+                fvMatrix<scalar>& eqn
+            ) const;
+
+
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/multicomponentPhaseChangeBase.C

@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "multicomponentPhaseChangeBase.H"
+#include "basicThermo.H"
+#include "multicomponentThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(multicomponentPhaseChangeBase, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fv::multicomponentPhaseChangeBase::readCoeffs()
+{
+    if (species_ != coeffs().lookup<wordList>("species"))
+    {
+        FatalIOErrorInFunction(coeffs())
+            << "Cannot change the species of a " << typeName << " model "
+            << "at run time" << exit(FatalIOError);
+    }
+
+    energySemiImplicit_ =
+        coeffs().lookup<bool>("energySemiImplicit");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::multicomponentPhaseChangeBase::multicomponentPhaseChangeBase
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    phaseChangeBase(name, modelType, mesh, dict, {true, true}),
+    species_(coeffs().lookup<wordList>("species")),
+    energySemiImplicit_(false)
+{
+    readCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::fv::multicomponentPhaseChangeBase::mDot() const
+{
+    tmp<volScalarField::Internal> tmDot =
+        volScalarField::Internal::New
+        (
+            "mDot",
+            mesh(),
+            dimensionedScalar(dimDensity/dimTime, Zero)
+        );
+
+    forAll(species(), mDoti)
+    {
+        tmDot.ref() += mDot(mDoti);
+    }
+
+    return tmDot;
+}
+
+
+void Foam::fv::multicomponentPhaseChangeBase::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volScalarField& heOrYi,
+    fvMatrix<scalar>& eqn
+) const
+{
+    const label i = index(phaseNames(), alpha.group());
+    const label s = sign(phaseNames(), alpha.group());
+
+    // Energy equation
+    if (index(heNames(), heOrYi.name()) != -1)
+    {
+        const volScalarField& p = this->p();
+        const volScalarField Tchange(vifToVf(this->Tchange()));
+
+        forAll(species(), mDoti)
+        {
+            const label speciei =
+                specieThermos()[i].species()[species()[mDoti]];
+
+            const volScalarField::Internal mDotIn(s*mDot(mDoti));
+
+            // Absolute enthalpies at the interface
+            Pair<tmp<volScalarField::Internal>> has;
+            for (label j = 0; j < 2; ++ j)
+            {
+                has[j] = vfToVif(specieThermos()[j].hai(speciei, p, Tchange));
+            }
+
+            // Direct transfer of energy due to mass transfer
+            if (energySemiImplicit_)
+            {
+                eqn += -fvm::SuSp(-mDotIn, heOrYi) + mDotIn*(has[i]() - heOrYi);
+            }
+            else
+            {
+                eqn += mDotIn*has[i]();
+            }
+
+            // Latent heat of phase change
+            eqn +=
+                (i == 0 ? Lfraction() - 1 : Lfraction())
+               *mDotIn*(has.second() - has.first());
+        }
+    }
+    // Mass fraction equation
+    else if
+    (
+        specieThermos().valid()[i]
+     && specieThermos()[i].containsSpecie(heOrYi.member())
+    )
+    {
+        // A transferring specie
+        if (species().found(heOrYi.member()))
+        {
+            eqn += s*mDot(species()[heOrYi.member()]);
+        }
+        // A non-transferring specie. Do nothing.
+        else
+        {}
+    }
+    // Something else. Fall back.
+    else
+    {
+        phaseChangeBase::addSup(alpha, rho, eqn);
+    }
+}
+
+
+bool Foam::fv::multicomponentPhaseChangeBase::read(const dictionary& dict)
+{
+    if (phaseChangeBase::read(dict))
+    {
+        readCoeffs();
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/multicomponentPhaseChangeBase.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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::multicomponentPhaseChangeBase
+
+Description
+    Base class for phase change models in which multiple components may change
+    phase. This can only be applied between multicomponent phases.
+
+SourceFiles
+    multicomponentPhaseChangeBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multicomponentPhaseChangeBase_H
+#define multicomponentPhaseChangeBase_H
+
+#include "phaseChangeBase.H"
+#include "hashedWordList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class multicomponentPhaseChangeBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class multicomponentPhaseChangeBase
+:
+    public phaseChangeBase
+{
+private:
+
+    // Private Data
+
+        //- The names of the changing species
+        const hashedWordList species_;
+
+        //- Whether or not to linearise the energy source
+        bool energySemiImplicit_;
+
+
+    // Private Member Functions
+
+        //- Non-virtual read
+        void readCoeffs();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multicomponentPhaseChangeBase");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        multicomponentPhaseChangeBase
+        (
+            const word& name,
+            const word& modelType,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the names of the changing species
+            inline const hashedWordList& species() const;
+
+
+        // Sources
+
+            //- Return the mass transfer rates of each specie
+            virtual tmp<DimensionedField<scalar, volMesh>> mDot
+            (
+                const label mDoti
+            ) const = 0;
+
+            //- Return the total mass transfer rate
+            virtual tmp<DimensionedField<scalar, volMesh>> mDot() const;
+
+            //- Override the energy equation to add the phase change heat, or
+            //  the species equations to add the relevant mass sources
+            void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& heOrYi,
+                fvMatrix<scalar>& eqn
+            ) const;
+
+
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "multicomponentPhaseChangeBaseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/multicomponentPhaseChangeBaseI.H

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "multicomponentPhaseChangeBase.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::hashedWordList&
+Foam::fv::multicomponentPhaseChangeBase::species() const
+{
+    return species_;
+}
+
+
+// ************************************************************************* //
+

src/fvModels/derived/phaseChange/phaseChangeBase.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "phaseChangeBase.H"
+#include "fluidThermo.H"
+#include "multicomponentThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(phaseChangeBase, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+const Foam::volScalarField& Foam::fv::phaseChangeBase::p() const
+{
+    for (label i = 0; i < 2; ++ i)
+    {
+        if (isA<fluidThermo>(thermos()[i]))
+        {
+            return refCast<const fluidThermo>(thermos()[i]).p();
+        }
+    }
+
+    return mesh().lookupObject<volScalarField>("p");
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::fv::phaseChangeBase::vifToVf
+(
+    const tmp<volScalarField::Internal>& tvif
+)
+{
+    tmp<volScalarField> tvf =
+        volScalarField::New
+        (
+            tvif().name(),
+            tvif().mesh(),
+            tvif().dimensions(),
+            extrapolatedCalculatedFvPatchField<scalar>::typeName
+        );
+
+    tvf->ref() = tvif();
+    tvf->correctBoundaryConditions();
+
+    tvif.clear();
+
+    return tvf;
+}
+
+
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::fv::phaseChangeBase::vfToVif
+(
+    const tmp<volScalarField>& tvf
+)
+{
+    tmp<volScalarField::Internal> tvif(tvf.ptr());
+
+    tvf.clear();
+
+    return tvif;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::phaseChangeBase::phaseChangeBase
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const Pair<bool> specieThermosRequired
+)
+:
+    massTransferBase(name, modelType, mesh, dict),
+    thermos_(mesh, phaseNames()),
+    specieThermos_(thermos_),
+    heNames_(thermos_.first().he().name(), thermos_.second().he().name())
+{
+    forAll(specieThermos_.valid(), i)
+    {
+        if (!specieThermos_.valid()[i] && specieThermosRequired[i])
+        {
+            FatalErrorInFunction
+                << "Model " << name << " of type " << modelType
+                << " requires a multicomponent thermo for phase "
+                << phaseNames()[i] << exit(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::fv::phaseChangeBase::Tchange() const
+{
+    const volScalarField::Internal mDot(this->mDot());
+
+    return pos0(mDot)*thermos().first().T() + neg(mDot)*thermos().second().T();
+}
+
+
+Foam::tmp<Foam::volScalarField::Internal>
+Foam::fv::phaseChangeBase::Lfraction() const
+{
+    const volScalarField& kappa1 = thermos().first().kappa();
+    const volScalarField& kappa2 = thermos().second().kappa();
+
+    return vfToVif(kappa2/(kappa1 + kappa2));
+}
+
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/phaseChangeBase.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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::phaseChangeBase
+
+Description
+    Base class for phase change models
+
+SourceFiles
+    phaseChangeBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseChangeBase_H
+#define phaseChangeBase_H
+
+#include "massTransferBase.H"
+#include "ThermoRefPair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class basicThermo;
+class multicomponentThermo;
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class phaseChangeBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseChangeBase
+:
+    public massTransferBase
+{
+private:
+
+    // Private Data
+
+        //- The thermo references
+        const ThermoRefPair<basicThermo> thermos_;
+
+        //- The specie thermo references
+        const ThermoRefPair<multicomponentThermo> specieThermos_;
+
+        //- Names of the energy fields
+        const Pair<word> heNames_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Access the pressure field
+        const volScalarField& p() const;
+
+        //- Add a boundary field to the given internal field
+        static tmp<volScalarField> vifToVf
+        (
+            const tmp<DimensionedField<scalar, volMesh>>& tvif
+        );
+
+        //- Remove the boundary field from the given geometric field
+        static tmp<DimensionedField<scalar, volMesh>> vfToVif
+        (
+            const tmp<volScalarField>& tvf
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("phaseChangeBase");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        phaseChangeBase
+        (
+            const word& name,
+            const word& modelType,
+            const fvMesh& mesh,
+            const dictionary& dict,
+            const Pair<bool> specieThermosRequired
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the thermo references
+            inline const ThermoRefPair<basicThermo>& thermos() const;
+
+            //- Return the specie thermo references
+            inline const ThermoRefPair<multicomponentThermo>&
+                specieThermos() const;
+
+            //- Return the names of the energy fields
+            inline const Pair<word>& heNames() const;
+
+
+        // Evaluation
+
+            //- Return the temperature at which the phases are considered to be
+            //  changing. By default this is considered to be the temperature
+            //  of the "source" phase (i.e., the phase for which mDot is
+            //  negative), but this can be overridden to account for heat
+            //  transfer modelling or similar.
+            virtual tmp<DimensionedField<scalar, volMesh>> Tchange() const;
+
+            //- Return the fraction of the latent heat that is transferred into
+            //  the second phase. By default this is weighted by the phase
+            //  thermal conductivities, but this can be overridden to account
+            //  for heat transfer modelling or similar.
+            virtual tmp<DimensionedField<scalar, volMesh>> Lfraction() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "phaseChangeBaseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/phaseChangeBaseI.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "phaseChangeBase.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::ThermoRefPair<Foam::basicThermo>&
+Foam::fv::phaseChangeBase::thermos() const
+{
+    return thermos_;
+}
+
+
+inline const Foam::ThermoRefPair<Foam::multicomponentThermo>&
+Foam::fv::phaseChangeBase::specieThermos() const
+{
+    return specieThermos_;
+}
+
+
+inline const Foam::Pair<Foam::word>&
+Foam::fv::phaseChangeBase::heNames() const
+{
+    return heNames_;
+}
+
+
+// ************************************************************************* //
+

src/fvModels/derived/phaseChange/singleComponentPhaseChangeBase.C

@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "singleComponentPhaseChangeBase.H"
+#include "basicThermo.H"
+#include "multicomponentThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(singleComponentPhaseChangeBase, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fv::singleComponentPhaseChangeBase::readCoeffs()
+{
+    if
+    (
+        (specieThermos().valid().first() || specieThermos().valid().second())
+     && specie_ != coeffs().lookup<word>("specie")
+    )
+    {
+        FatalIOErrorInFunction(coeffs())
+            << "Cannot change the specie of a " << typeName << " model "
+            << "at run time" << exit(FatalIOError);
+    }
+
+    energySemiImplicit_ =
+        coeffs().lookup<bool>("energySemiImplicit");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::singleComponentPhaseChangeBase::singleComponentPhaseChangeBase
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const Pair<bool> specieThermosRequired
+)
+:
+    phaseChangeBase(name, modelType, mesh, dict, specieThermosRequired),
+    specie_
+    (
+        specieThermos().valid().first() || specieThermos().valid().second()
+      ? coeffs().lookup<word>("specie")
+      : word::null
+    ),
+    specieis_
+    (
+        specieThermos().valid().first()
+      ? specieThermos().first().species()[specie_]
+      : -1,
+        specieThermos().valid().second()
+      ? specieThermos().second().species()[specie_]
+      : -1
+    ),
+    energySemiImplicit_(false)
+{
+    readCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::fv::singleComponentPhaseChangeBase::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    const volScalarField& heOrYi,
+    fvMatrix<scalar>& eqn
+) const
+{
+    const label i = index(phaseNames(), alpha.group());
+    const label s = sign(phaseNames(), alpha.group());
+
+    // Energy equation
+    if (index(heNames(), heOrYi.name()) != -1)
+    {
+        const volScalarField& p = this->p();
+        const volScalarField Tchange(vifToVf(this->Tchange()));
+
+        const volScalarField::Internal mDotIn(s*mDot());
+
+        // Absolute enthalpies at the interface
+        Pair<tmp<volScalarField::Internal>> has;
+        for (label j = 0; j < 2; ++ j)
+        {
+            has[j] =
+                specieThermos().valid()[j]
+              ? vfToVif(specieThermos()[j].hai(specieis_[j], p, Tchange))
+              : vfToVif(thermos()[j].ha(p, Tchange));
+        }
+
+        // Direct transfer of energy due to mass transfer
+        if (energySemiImplicit_)
+        {
+            eqn += -fvm::SuSp(-mDotIn, heOrYi) + mDotIn*(has[i]() - heOrYi);
+        }
+        else
+        {
+            eqn += mDotIn*has[i]();
+        }
+
+        // Latent heat of phase change
+        eqn +=
+            (i == 0 ? Lfraction() - 1 : Lfraction())
+           *mDotIn*(has.second() - has.first());
+    }
+    // Mass fraction equation
+    else if
+    (
+        specieThermos().valid()[i]
+     && specieThermos()[i].containsSpecie(heOrYi.member())
+    )
+    {
+        // The transferring specie
+        if (heOrYi.member() == specie())
+        {
+            eqn += s*mDot();
+        }
+        // A non-transferring specie. Do nothing.
+        else
+        {}
+    }
+    // Something else. Fall back.
+    else
+    {
+        massTransferBase::addSup(alpha, rho, heOrYi, eqn);
+    }
+}
+
+
+bool Foam::fv::singleComponentPhaseChangeBase::read(const dictionary& dict)
+{
+    if (phaseChangeBase::read(dict))
+    {
+        readCoeffs();
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/singleComponentPhaseChangeBase.H

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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::singleComponentPhaseChangeBase
+
+Description
+    Base class for phase change models in which only a single component changes
+    phase. Can be applied to any combination of pure and multicomponent phases.
+    If either phase is multicomponent, then a single specie must be identified
+    as the one that changes phase.
+
+SourceFiles
+    singleComponentPhaseChangeBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef singleComponentPhaseChangeBase_H
+#define singleComponentPhaseChangeBase_H
+
+#include "phaseChangeBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class singleComponentPhaseChangeBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class singleComponentPhaseChangeBase
+:
+    public phaseChangeBase
+{
+private:
+
+    // Private Data
+
+        //- The name of the changing specie, or word::null if neither thermo
+        //  is multicomponent
+        const word specie_;
+
+        //- The indices of the changing species in the phases, or -1 if the
+        //  phases are not multicomponent
+        const Pair<label> specieis_;
+
+        //- Whether or not to linearise the energy source
+        bool energySemiImplicit_;
+
+
+    // Private Member Functions
+
+        //- Non-virtual read
+        void readCoeffs();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("singleComponentPhaseChangeBase");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        singleComponentPhaseChangeBase
+        (
+            const word& name,
+            const word& modelType,
+            const fvMesh& mesh,
+            const dictionary& dict,
+            const Pair<bool> specieThermosRequired
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the name of the changing specie
+            inline const word& specie() const;
+
+            //- Return the indices of the changing species in the phases
+            inline const Pair<label>& specieis() const;
+
+
+        // Sources
+
+            //- Override the energy equation to add the phase change heat, or
+            //  the species equations to add the relevant mass sources
+            void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                const volScalarField& heOrYi,
+                fvMatrix<scalar>& eqn
+            ) const;
+
+
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "singleComponentPhaseChangeBaseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvModels/derived/phaseChange/singleComponentPhaseChangeBaseI.H

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021-2023 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 "singleComponentPhaseChangeBase.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::word&
+Foam::fv::singleComponentPhaseChangeBase::specie() const
+{
+    return specie_;
+}
+
+
+inline const Foam::Pair<Foam::label>&
+Foam::fv::singleComponentPhaseChangeBase::specieis() const
+{
+    return specieis_;
+}
+
+
+// ************************************************************************* //
+