diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index 71c5a06724c5829f79184cbdf3de358629c481fc..ff4fb10e9dddaf69fb27da3418c39a9516ba56a8 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -28,7 +28,6 @@ License
 #include "BlendedInterfacialModel.H"
 #include "heatTransferModel.H"
 #include "massTransferModel.H"
-#include "interfaceCompositionModel.H"
 
 #include "HashPtrTable.H"
 
@@ -59,12 +58,6 @@ HeatAndMassTransferPhaseSystem
         massTransferModels_
     );
 
-    this->generatePairsAndSubModels
-    (
-        "interfaceComposition",
-        interfaceCompositionModels_
-    );
-
     forAllConstIter
     (
         phaseSystem::phasePairTable,
@@ -332,263 +325,6 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
 }
 
 
-template<class BasePhaseSystem>
-Foam::autoPtr<Foam::phaseSystem::massTransferTable>
-Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::massTransfer() const
-{
-    // Create a mass transfer matrix for each species of each phase
-    autoPtr<phaseSystem::massTransferTable> eqnsPtr
-    (
-        new phaseSystem::massTransferTable()
-    );
-
-    phaseSystem::massTransferTable& eqns = eqnsPtr();
-
-    forAllConstIter
-    (
-        phaseSystem::phaseModelTable,
-        this->phaseModels_,
-        phaseModelIter
-    )
-    {
-        const phaseModel& phase(phaseModelIter());
-
-        const PtrList<volScalarField>& Yi = phase.Y();
-
-        forAll(Yi, i)
-        {
-            eqns.insert
-            (
-                Yi[i].name(),
-                new fvScalarMatrix(Yi[i], dimMass/dimTime)
-            );
-        }
-    }
-
-    // Reset the interfacial mass flow rates
-    forAllConstIter
-    (
-        phaseSystem::phasePairTable,
-        this->phasePairs_,
-        phasePairIter
-    )
-    {
-        const phasePair& pair(phasePairIter());
-
-        if (pair.ordered())
-        {
-            continue;
-        }
-
-        *dmdt_[pair] =
-            *dmdtExplicit_[pair];
-
-        *dmdtExplicit_[pair] =
-            dimensionedScalar("zero", dimDensity/dimTime, 0);
-    }
-
-    // Sum up the contribution from each interface composition model
-    forAllConstIter
-    (
-        interfaceCompositionModelTable,
-        interfaceCompositionModels_,
-        interfaceCompositionModelIter
-    )
-    {
-        const interfaceCompositionModel& compositionModel
-        (
-            interfaceCompositionModelIter()
-        );
-
-        const phasePair& pair
-        (
-            this->phasePairs_[interfaceCompositionModelIter.key()]
-        );
-        const phaseModel& phase = pair.phase1();
-        const phaseModel& otherPhase = pair.phase2();
-        const phasePairKey key(phase.name(), otherPhase.name());
-
-        const volScalarField& Tf(*Tf_[key]);
-
-        volScalarField& dmdtExplicit(*dmdtExplicit_[key]);
-        volScalarField& dmdt(*dmdt_[key]);
-
-        scalar dmdtSign(Pair<word>::compare(dmdt_.find(key).key(), key));
-
-        const volScalarField K
-        (
-            massTransferModels_[key][phase.name()]->K()
-        );
-
-        forAllConstIter
-        (
-            hashedWordList,
-            compositionModel.species(),
-            memberIter
-        )
-        {
-            const word& member = *memberIter;
-
-            const word name
-            (
-                IOobject::groupName(member, phase.name())
-            );
-
-            const word otherName
-            (
-                IOobject::groupName(member, otherPhase.name())
-            );
-
-            const volScalarField KD
-            (
-                K*compositionModel.D(member)
-            );
-
-            const volScalarField Yf
-            (
-                compositionModel.Yf(member, Tf)
-            );
-
-            // Implicit transport through the phase
-            *eqns[name] +=
-                phase.rho()*KD*Yf
-              - fvm::Sp(phase.rho()*KD, eqns[name]->psi());
-
-            // Sum the mass transfer rate
-            dmdtExplicit += dmdtSign*phase.rho()*KD*Yf;
-            dmdt -= dmdtSign*phase.rho()*KD*eqns[name]->psi();
-
-            // Explicit transport out of the other phase
-            if (eqns.found(otherName))
-            {
-                *eqns[otherName] -=
-                    otherPhase.rho()*KD*compositionModel.dY(member, Tf);
-            }
-        }
-    }
-
-    return eqnsPtr;
-}
-
-
-template<class BasePhaseSystem>
-void Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::correctThermo()
-{
-    BasePhaseSystem::correctThermo();
-
-    // This loop solves for the interface temperatures, Tf, and updates the
-    // interface composition models.
-    //
-    // The rate of heat transfer to the interface must equal the latent heat
-    // consumed at the interface, i.e.:
-    //
-    // H1*(T1 - Tf) + H2*(T2 - Tf) == mDotL
-    //                             == K*rho*(Yfi - Yi)*Li
-    //
-    // Yfi is likely to be a strong non-linear (typically exponential) function
-    // of Tf, so the solution for the temperature is newton-accelerated
-
-    forAllConstIter
-    (
-        phaseSystem::phasePairTable,
-        this->phasePairs_,
-        phasePairIter
-    )
-    {
-        const phasePair& pair(phasePairIter());
-
-        if (pair.ordered())
-        {
-            continue;
-        }
-
-        const phasePairKey key12(pair.first(), pair.second(), true);
-        const phasePairKey key21(pair.second(), pair.first(), true);
-
-        volScalarField H1(heatTransferModels_[pair][pair.first()]->K());
-        volScalarField H2(heatTransferModels_[pair][pair.second()]->K());
-        dimensionedScalar HSmall("small", heatTransferModel::dimK, SMALL);
-
-        volScalarField mDotL
-        (
-            IOobject
-            (
-                "mDotL",
-                this->mesh().time().timeName(),
-                this->mesh()
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0)
-        );
-        volScalarField mDotLPrime
-        (
-            IOobject
-            (
-                "mDotLPrime",
-                this->mesh().time().timeName(),
-                this->mesh()
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", mDotL.dimensions()/dimTemperature, 0)
-        );
-
-        volScalarField& Tf = *Tf_[pair];
-
-        // Add latent heats from forward and backward models
-        if (interfaceCompositionModels_.found(key12))
-        {
-            interfaceCompositionModels_[key12]->addMDotL
-            (
-                massTransferModels_[pair][pair.first()]->K(),
-                Tf,
-                mDotL,
-                mDotLPrime
-            );
-        }
-        if (interfaceCompositionModels_.found(key21))
-        {
-            interfaceCompositionModels_[key21]->addMDotL
-            (
-                massTransferModels_[pair][pair.second()]->K(),
-                Tf,
-                mDotL,
-                mDotLPrime
-            );
-        }
-
-        // Update the interface temperature by applying one step of newton's
-        // method to the interface relation
-        Tf -=
-            (
-                H1*(Tf - pair.phase1().thermo().T())
-              + H2*(Tf - pair.phase2().thermo().T())
-              + mDotL
-            )
-           /(
-                max(H1 + H2 + mDotLPrime, HSmall)
-            );
-
-        // Update the interface compositions
-        if (interfaceCompositionModels_.found(key12))
-        {
-            interfaceCompositionModels_[key12]->update(Tf);
-        }
-        if (interfaceCompositionModels_.found(key21))
-        {
-            interfaceCompositionModels_[key21]->update(Tf);
-        }
-
-        Tf.correctBoundaryConditions();
-
-        Info<< "Tf." << pair.name()
-            << ": min = " << min(Tf.internalField())
-            << ", mean = " << average(Tf.internalField())
-            << ", max = " << max(Tf.internalField())
-            << endl;
-    }
-}
-
-
 template<class BasePhaseSystem>
 bool Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::read()
 {
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
index 27e4f32a77034a47c7977d700907bf2d41bb6cd0..3f20c575f83d4abccd2d1504c1134373eb7c2325 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H
@@ -25,12 +25,8 @@ Class
     Foam::HeatAndMassTransferPhaseSystem
 
 Description
-    Class which models interfacial heat and mass transfer between a number of
-    phases. Mass is transferred to or from a phase according to a composition
-    model at the surface of another phase. Heat is transferred from a phase to
-    an interfacial temperature. The interface temperature is calculated such
-    that the net rate at which the heat is transferred to the interface is
-    equal to the latent heat consumed by the mass transfer.
+    Base class to support interfacial heat and mass transfer between a number
+    of phases.
 
 SourceFiles
     HeatAndMassTransferPhaseSystem.C
@@ -53,7 +49,6 @@ class BlendedInterfacialModel;
 class blendingMethod;
 class heatTransferModel;
 class massTransferModel;
-class interfaceCompositionModel;
 
 /*---------------------------------------------------------------------------*\
                  Class HeatAndMassTransferPhaseSystem Declaration
@@ -88,13 +83,6 @@ protected:
             phasePairKey::hash
         > massTransferModelTable;
 
-        typedef HashTable
-        <
-            autoPtr<interfaceCompositionModel>,
-            phasePairKey,
-            phasePairKey::hash
-        > interfaceCompositionModelTable;
-
 
     // Protected data
 
@@ -117,9 +105,6 @@ protected:
             //- Mass transfer models
             massTransferModelTable massTransferModels_;
 
-            //- Interface composition models
-            interfaceCompositionModelTable interfaceCompositionModels_;
-
 
 public:
 
@@ -151,10 +136,10 @@ public:
 
         //- Return the mass transfer matrices
         virtual autoPtr<phaseSystem::massTransferTable>
-            massTransfer() const;
+            massTransfer() const = 0;
 
         //- Correct the thermodynamics
-        virtual void correctThermo();
+        virtual void correctThermo() = 0;
 
         //- Read base phaseProperties dictionary
         virtual bool read();
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
new file mode 100644
index 0000000000000000000000000000000000000000..47c994334d63bcc48c9be3fd696b4f21fb68647c
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
@@ -0,0 +1,336 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "InterfaceCompositionPhaseChangePhaseSystem.H"
+#include "interfaceCompositionModel.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
+InterfaceCompositionPhaseChangePhaseSystem
+(
+    const fvMesh& mesh
+)
+:
+    HeatAndMassTransferPhaseSystem<BasePhaseSystem>(mesh)
+{
+    this->generatePairsAndSubModels
+    (
+        "interfaceComposition",
+        interfaceCompositionModels_
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
+~InterfaceCompositionPhaseChangePhaseSystem()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::phaseSystem::massTransferTable>
+Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
+massTransfer() const
+{
+    // Create a mass transfer matrix for each species of each phase
+    autoPtr<phaseSystem::massTransferTable> eqnsPtr
+    (
+        new phaseSystem::massTransferTable()
+    );
+
+    phaseSystem::massTransferTable& eqns = eqnsPtr();
+
+    forAllConstIter
+    (
+        phaseSystem::phaseModelTable,
+        this->phaseModels_,
+        phaseModelIter
+    )
+    {
+        const phaseModel& phase(phaseModelIter());
+
+        const PtrList<volScalarField>& Yi = phase.Y();
+
+        forAll(Yi, i)
+        {
+            eqns.insert
+            (
+                Yi[i].name(),
+                new fvScalarMatrix(Yi[i], dimMass/dimTime)
+            );
+        }
+    }
+
+    // Reset the interfacial mass flow rates
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        *this->dmdt_[pair] =
+            *this->dmdtExplicit_[pair];
+
+        *this->dmdtExplicit_[pair] =
+            dimensionedScalar("zero", dimDensity/dimTime, 0);
+    }
+
+    // Sum up the contribution from each interface composition model
+    forAllConstIter
+    (
+        interfaceCompositionModelTable,
+        interfaceCompositionModels_,
+        interfaceCompositionModelIter
+    )
+    {
+        const interfaceCompositionModel& compositionModel
+        (
+            interfaceCompositionModelIter()
+        );
+
+        const phasePair& pair
+        (
+            this->phasePairs_[interfaceCompositionModelIter.key()]
+        );
+        const phaseModel& phase = pair.phase1();
+        const phaseModel& otherPhase = pair.phase2();
+        const phasePairKey key(phase.name(), otherPhase.name());
+
+        const volScalarField& Tf(*this->Tf_[key]);
+
+        volScalarField& dmdtExplicit(*this->dmdtExplicit_[key]);
+        volScalarField& dmdt(*this->dmdt_[key]);
+
+        scalar dmdtSign(Pair<word>::compare(this->dmdt_.find(key).key(), key));
+
+        const volScalarField K
+        (
+            this->massTransferModels_[key][phase.name()]->K()
+        );
+
+        forAllConstIter
+        (
+            hashedWordList,
+            compositionModel.species(),
+            memberIter
+        )
+        {
+            const word& member = *memberIter;
+
+            const word name
+            (
+                IOobject::groupName(member, phase.name())
+            );
+
+            const word otherName
+            (
+                IOobject::groupName(member, otherPhase.name())
+            );
+
+            const volScalarField KD
+            (
+                K*compositionModel.D(member)
+            );
+
+            const volScalarField Yf
+            (
+                compositionModel.Yf(member, Tf)
+            );
+
+            // Implicit transport through the phase
+            *eqns[name] +=
+                phase.rho()*KD*Yf
+              - fvm::Sp(phase.rho()*KD, eqns[name]->psi());
+
+            // Sum the mass transfer rate
+            dmdtExplicit += dmdtSign*phase.rho()*KD*Yf;
+            dmdt -= dmdtSign*phase.rho()*KD*eqns[name]->psi();
+
+            // Explicit transport out of the other phase
+            if (eqns.found(otherName))
+            {
+                *eqns[otherName] -=
+                    otherPhase.rho()*KD*compositionModel.dY(member, Tf);
+            }
+        }
+    }
+
+    return eqnsPtr;
+}
+
+
+template<class BasePhaseSystem>
+void Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
+correctThermo()
+{
+    BasePhaseSystem::correctThermo();
+
+    // This loop solves for the interface temperatures, Tf, and updates the
+    // interface composition models.
+    //
+    // The rate of heat transfer to the interface must equal the latent heat
+    // consumed at the interface, i.e.:
+    //
+    // H1*(T1 - Tf) + H2*(T2 - Tf) == mDotL
+    //                             == K*rho*(Yfi - Yi)*Li
+    //
+    // Yfi is likely to be a strong non-linear (typically exponential) function
+    // of Tf, so the solution for the temperature is newton-accelerated
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const phasePairKey key12(pair.first(), pair.second(), true);
+        const phasePairKey key21(pair.second(), pair.first(), true);
+
+        volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
+        volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
+        dimensionedScalar HSmall("small", heatTransferModel::dimK, SMALL);
+
+        volScalarField mDotL
+        (
+            IOobject
+            (
+                "mDotL",
+                this->mesh().time().timeName(),
+                this->mesh()
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0)
+        );
+        volScalarField mDotLPrime
+        (
+            IOobject
+            (
+                "mDotLPrime",
+                this->mesh().time().timeName(),
+                this->mesh()
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", mDotL.dimensions()/dimTemperature, 0)
+        );
+
+        volScalarField& Tf = *this->Tf_[pair];
+
+        // Add latent heats from forward and backward models
+        if (this->interfaceCompositionModels_.found(key12))
+        {
+            this->interfaceCompositionModels_[key12]->addMDotL
+            (
+                this->massTransferModels_[pair][pair.first()]->K(),
+                Tf,
+                mDotL,
+                mDotLPrime
+            );
+        }
+        if (this->interfaceCompositionModels_.found(key21))
+        {
+            this->interfaceCompositionModels_[key21]->addMDotL
+            (
+                this->massTransferModels_[pair][pair.second()]->K(),
+                Tf,
+                mDotL,
+                mDotLPrime
+            );
+        }
+
+        // Update the interface temperature by applying one step of newton's
+        // method to the interface relation
+        Tf -=
+            (
+                H1*(Tf - pair.phase1().thermo().T())
+              + H2*(Tf - pair.phase2().thermo().T())
+              + mDotL
+            )
+           /(
+                max(H1 + H2 + mDotLPrime, HSmall)
+            );
+
+        Tf.correctBoundaryConditions();
+
+        Info<< "Tf." << pair.name()
+            << ": min = " << min(Tf.internalField())
+            << ", mean = " << average(Tf.internalField())
+            << ", max = " << max(Tf.internalField())
+            << endl;
+
+        // Update the interface compositions
+        if (this->interfaceCompositionModels_.found(key12))
+        {
+            this->interfaceCompositionModels_[key12]->update(Tf);
+        }
+        if (this->interfaceCompositionModels_.found(key21))
+        {
+            this->interfaceCompositionModels_[key21]->update(Tf);
+        }
+    }
+}
+
+
+template<class BasePhaseSystem>
+bool Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::read()
+{
+    if (BasePhaseSystem::read())
+    {
+        bool readOK = true;
+
+        // Models ...
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H
new file mode 100644
index 0000000000000000000000000000000000000000..b87c5a709da68677ed6c29e919831a93f04e59c6
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::InterfaceCompositionPhaseChangePhaseSystem
+
+Description
+    Class to provide interfacial heat and mass transfer between a number of
+    phases according to a interface composition model.
+
+    The interface temperature is calculated such that the net rate at which the
+    heat is transferred to the interface is equal to the latent heat consumed by
+    the mass transfer.
+
+SourceFiles
+    InterfaceCompositionPhaseChangePhaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef InterfaceCompositionPhaseChangePhaseSystem_H
+#define InterfaceCompositionPhaseChangePhaseSystem_H
+
+#include "HeatAndMassTransferPhaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class interfaceCompositionModel;
+
+/*---------------------------------------------------------------------------*\
+                 Class InterfaceCompositionPhaseChangePhaseSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasePhaseSystem>
+class InterfaceCompositionPhaseChangePhaseSystem
+:
+    public HeatAndMassTransferPhaseSystem<BasePhaseSystem>
+{
+protected:
+
+    // Protected typedefs
+
+        typedef HashTable
+        <
+            autoPtr<interfaceCompositionModel>,
+            phasePairKey,
+            phasePairKey::hash
+        > interfaceCompositionModelTable;
+
+
+    // Protected data
+
+        // Sub Models
+
+            //- Interface composition models
+            interfaceCompositionModelTable interfaceCompositionModels_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        InterfaceCompositionPhaseChangePhaseSystem(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~InterfaceCompositionPhaseChangePhaseSystem();
+
+
+    // Member Functions
+
+        //- Return the mass transfer matrices
+        virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
+
+        //- Correct the thermodynamics
+        virtual void correctThermo();
+
+        //- Read base phaseProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "InterfaceCompositionPhaseChangePhaseSystem.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
new file mode 100644
index 0000000000000000000000000000000000000000..bbce1962a12a3e105be2176f723349e63abd14a1
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ThermalPhaseChangePhaseSystem.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
+ThermalPhaseChangePhaseSystem
+(
+    const fvMesh& mesh
+)
+:
+    HeatAndMassTransferPhaseSystem<BasePhaseSystem>(mesh),
+    volatile_(this->lookup("volatile")),
+    saturationModel_(saturationModel::New(this->subDict("saturationModel")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
+~ThermalPhaseChangePhaseSystem()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::phaseSystem::massTransferTable>
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
+{
+    // Create a mass transfer matrix for each species of each phase
+    autoPtr<phaseSystem::massTransferTable> eqnsPtr
+    (
+        new phaseSystem::massTransferTable()
+    );
+
+    phaseSystem::massTransferTable& eqns = eqnsPtr();
+
+    forAllConstIter
+    (
+        phaseSystem::phaseModelTable,
+        this->phaseModels_,
+        phaseModelIter
+    )
+    {
+        const phaseModel& phase(phaseModelIter());
+
+        const PtrList<volScalarField>& Yi = phase.Y();
+
+        forAll(Yi, i)
+        {
+            eqns.insert
+            (
+                Yi[i].name(),
+                new fvScalarMatrix(Yi[i], dimMass/dimTime)
+            );
+        }
+    }
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+        const phaseModel& phase = pair.phase1();
+        const phaseModel& otherPhase = pair.phase2();
+
+        const word name
+        (
+            IOobject::groupName(volatile_, phase.name())
+        );
+
+        const word otherName
+        (
+            IOobject::groupName(volatile_, otherPhase.name())
+        );
+
+        const volScalarField dmdt(this->dmdt(pair));
+        const volScalarField dmdt12(posPart(dmdt));
+        const volScalarField dmdt21(negPart(dmdt));
+
+        *eqns[name] += fvm::Sp(dmdt21, eqns[name]->psi()) - dmdt21;
+        *eqns[otherName] += dmdt12 - fvm::Sp(dmdt12, eqns[otherName]->psi());
+    }
+
+    return eqnsPtr;
+}
+
+
+template<class BasePhaseSystem>
+void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
+{
+    BasePhaseSystem::correctThermo();
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const phaseModel& phase1 = pair.phase1();
+        const phaseModel& phase2 = pair.phase2();
+
+        volScalarField& Tf = *this->Tf_[pair];
+        Tf = saturationModel_->Tsat(phase1.thermo().p());
+
+        Info<< "Tf." << pair.name()
+            << ": min = " << min(Tf.internalField())
+            << ", mean = " << average(Tf.internalField())
+            << ", max = " << max(Tf.internalField())
+            << endl;
+
+        volScalarField& dmdt(*this->dmdt_[pair]);
+
+        volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
+        volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
+
+        const volScalarField& T1(phase1.thermo().T());
+        const volScalarField& T2(phase2.thermo().T());
+
+        const volScalarField& he1(phase1.thermo().he());
+        const volScalarField& he2(phase2.thermo().he());
+
+        volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
+        volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
+
+        dmdt =
+            (H2*(Tf - T1) + H1*(Tf - T2))
+           /min
+            (
+                (pos(dmdt)*he2 + neg(dmdt)*hef2)
+              - (neg(dmdt)*he1 + pos(dmdt)*hef1),
+                0.3*(hef1 - hef2)
+            );
+    }
+}
+
+
+template<class BasePhaseSystem>
+bool Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::read()
+{
+    if (BasePhaseSystem::read())
+    {
+        bool readOK = true;
+
+        // Models ...
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
new file mode 100644
index 0000000000000000000000000000000000000000..9a0a0d109a446a09a79c29dfd62cce912ccb3453
--- /dev/null
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ThermalPhaseChangePhaseSystem
+
+Description
+    Class to provide interfacial heat and mass transfer between a number of
+    phases according the interfacial temperature approximated by the saturation
+    temperature.
+
+    Currently only a single specified specie is considered volatile and changes
+    phase, all other species are considered nonvolatile and do not
+    affect the mass-transfer.
+
+SourceFiles
+    ThermalPhaseChangePhaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ThermalPhaseChangePhaseSystem_H
+#define ThermalPhaseChangePhaseSystem_H
+
+#include "HeatAndMassTransferPhaseSystem.H"
+#include "saturationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class ThermalPhaseChangePhaseSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasePhaseSystem>
+class ThermalPhaseChangePhaseSystem
+:
+    public HeatAndMassTransferPhaseSystem<BasePhaseSystem>
+{
+
+protected:
+
+    // Protected data
+
+        //- Name of the volatile specie
+        word volatile_;
+
+        //- The saturation model used to evaluate Tsat = Tf
+        autoPtr<saturationModel> saturationModel_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        ThermalPhaseChangePhaseSystem(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~ThermalPhaseChangePhaseSystem();
+
+
+    // Member Functions
+
+        //- Return the mass transfer matrices
+        virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
+
+        //- Correct the thermodynamics
+        virtual void correctThermo();
+
+        //- Read base phaseProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ThermalPhaseChangePhaseSystem.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystems.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystems.C
index 6b649afb6410a077c24217622632de4ecdcc7d27..075a578c2330cd3ff6df6764d089ef772ed716e6 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystems.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystems.C
@@ -29,7 +29,8 @@ License
 #include "twoPhaseSystem.H"
 #include "MomentumTransferPhaseSystem.H"
 #include "HeatTransferPhaseSystem.H"
-#include "HeatAndMassTransferPhaseSystem.H"
+#include "InterfaceCompositionPhaseChangePhaseSystem.H"
+#include "ThermalPhaseChangePhaseSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,19 +52,35 @@ namespace Foam
     );
 
     typedef
-        HeatAndMassTransferPhaseSystem
+        InterfaceCompositionPhaseChangePhaseSystem
         <
             MomentumTransferPhaseSystem<twoPhaseSystem>
         >
-        heatMassAndMomentumTransferTwoPhaseSystem;
+        interfaceCompositionPhaseChangeTwoPhaseSystem;
 
     addNamedToRunTimeSelectionTable
     (
         twoPhaseSystem,
-        heatMassAndMomentumTransferTwoPhaseSystem,
+        interfaceCompositionPhaseChangeTwoPhaseSystem,
         dictionary,
-        heatMassAndMomentumTransferTwoPhaseSystem
+        interfaceCompositionPhaseChangeTwoPhaseSystem
+    );
+
+    typedef
+        ThermalPhaseChangePhaseSystem
+        <
+            MomentumTransferPhaseSystem<twoPhaseSystem>
+        >
+        thermalPhaseChangeTwoPhaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        twoPhaseSystem,
+        thermalPhaseChangeTwoPhaseSystem,
+        dictionary,
+        thermalPhaseChangeTwoPhaseSystem
     );
 }
 
+
 // ************************************************************************* //