diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/InterfaceCompositionModel/InterfaceCompositionModels.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/InterfaceCompositionModel/InterfaceCompositionModels.C
index a459e0ee7182727d46be24550c6a9a94fa3920b3..ff9bcb7d2d66ef7436ea2b17cbd1cd4ef2b9445c 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/InterfaceCompositionModel/InterfaceCompositionModels.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/InterfaceCompositionModel/InterfaceCompositionModels.C
@@ -50,6 +50,7 @@ License
 
 #include "kineticGasEvaporation.H"
 #include "Lee.H"
+#include "interfaceHeatResistance.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -328,6 +329,152 @@ namespace Foam
         );
 
 
+
+        // interfaceHeatResistance model definitions
+
+        // From pure phase (rho const) to phase (rho const)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            constRhoHThermoPhysics,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            constRhoHThermoPhysics
+        );
+
+        // From pure phase (rho const) to phase (Boussinesq)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            constRhoHThermoPhysics,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            BoussinesqFluidEThermoPhysics
+        );
+
+
+        // From pure phase (solidThermo) to phase (Boussinesq)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hConstSolidThermoPhysics,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            BoussinesqFluidEThermoPhysics
+        );
+
+        // From pure phase (solidThermo) to phase (rho const)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hConstSolidThermoPhysics,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            constRhoHThermoPhysics
+        );
+
+        // From pure phase (all-poly solidThermo) to phase (ico-rho)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hPolyTranspPolyIcoSolidThermoPhysics,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            icoPoly8HThermoPhysics
+        );
+
+        // From pure phase (exp-Transp, hPower solidThermo) to phase (ico-rho)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hPowerSolidThermoPhysics,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            icoPoly8HThermoPhysics
+        );
+
+
+        // From pure phase (const rho) to multi phase (incomp ideal gas)
+        makeInterfaceContSpecieMixtureType
+        (
+            interfaceHeatResistance,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            constRhoHThermoPhysics,
+            heRhoThermo,
+            rhoReactionThermo,
+            multiComponentMixture,
+            constIncompressibleGasHThermoPhysics
+        );
+
+
+        // From pure phase (Boussinesq) to phase (solidThermo)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            BoussinesqFluidEThermoPhysics,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hConstSolidThermoPhysics
+        );
+
+        // From pure phase (rho const) to phase (solidThermo)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            constRhoHThermoPhysics,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hConstSolidThermoPhysics
+        );
+
+        //From pure liquid phase (ico-rho) to pure phase (exp-Transp, hPower solidThermo)
+        makeInterfacePureType
+        (
+            interfaceHeatResistance,
+            heRhoThermo,
+            rhoThermo,
+            pureMixture,
+            icoPoly8HThermoPhysics,
+            heSolidThermo,
+            solidThermo,
+            pureMixture,
+            hPowerSolidThermoPhysics
+        );
+
 }
 
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H
index e1f60133733b60e0755f1d0157c645b1cb8e7751..04d6f7ce7158051225e02c21125bb7dbd7738ebd 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,6 +80,8 @@ Usage
         Property    | Description            | Required    | Default value
         Tactivate   | Activation temperature | yes
         C           | Model constant         | yes
+        includeVolChange    | Volumen change  | no          | yes
+        species     | Specie name on the other phase | no   | none
     \endtable
 
 SourceFiles
@@ -143,7 +145,7 @@ public:
 
     // Member Functions
 
-        //- Explicit mass transfer coefficient
+        //- Explicit total mass transfer coefficient
         virtual tmp<volScalarField> Kexp
         (
             const volScalarField& field
@@ -169,6 +171,7 @@ public:
         //- Adds and substract alpha*div(U) as a source term
         //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
         virtual bool includeDivU();
+
 };
 
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C
index 8d54a861079e5b839c62d3c0f1049f456414f831..3f84cf11ceaddddbe28d55ec2cdc28d471ff3e4d 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +44,7 @@ Foam::interfaceCompositionModel::modelVariableNames
     { modelVariable::T, "temperature" },
     { modelVariable::P, "pressure" },
     { modelVariable::Y, "massFraction" },
+    { modelVariable::alpha, "alphaVolumeFraction" },
 };
 
 
@@ -64,6 +65,7 @@ Foam::interfaceCompositionModel::interfaceCompositionModel
             modelVariable::T
         )
     ),
+    includeVolChange_(dict.lookupOrDefault<bool>("includeVolChange", true)),
     pair_(pair),
     speciesName_(dict.lookupOrDefault<word>("species", "none")),
     mesh_(pair_.from().mesh())
@@ -96,4 +98,9 @@ bool Foam::interfaceCompositionModel::includeDivU()
 }
 
 
+bool Foam::interfaceCompositionModel::includeVolChange()
+{
+    return includeVolChange_;
+}
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
index 2260b1db1dffa6156fcbf223b7f6e48ad0b829a7..9962a95add446dbe4249f5c4d164471f52f34400 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,6 +76,9 @@ public:
         //- Enumeration for model variables
         modelVariable modelVariable_;
 
+        //- Add volume change in pEq
+        bool includeVolChange_;
+
 
 protected:
 
@@ -198,6 +201,9 @@ public:
         //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
         virtual bool includeDivU();
 
+        //- Add volume change in pEq
+        bool includeVolChange();
+
         //- Returns the variable on which the model is based
         const word variable() const;
 };
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.C
new file mode 100644
index 0000000000000000000000000000000000000000..797d71ceef410f7ce4d0ce5bc6766ac67e20f83f
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.C
@@ -0,0 +1,343 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 Henning Scheufler
+-------------------------------------------------------------------------------
+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 "interfaceHeatResistance.H"
+#include "constants.H"
+#include "isoCutCell.H"
+#include "volPointInterpolation.H"
+#include "wallPolyPatch.H"
+#include "fvcSmooth.H"
+
+using namespace Foam::constant;
+
+
+// * * * * * * * * * * *  Private Member Functions * * * * * * * * * * * * * //
+
+template<class Thermo, class OtherThermo>
+void Foam::meltingEvaporationModels::
+interfaceHeatResistance<Thermo, OtherThermo>
+::updateInterface(const volScalarField& T)
+{
+    const fvMesh& mesh = this->mesh_;
+
+    const volScalarField& alpha = this->pair().from();
+
+    scalarField ap
+    (
+        volPointInterpolation::New(mesh).interpolate(alpha)
+    );
+
+    isoCutCell cutCell(mesh, ap);
+
+    forAll(interfaceArea_, celli)
+    {
+        label status = cutCell.calcSubCell(celli, isoAlpha_);
+        interfaceArea_[celli] = 0;
+        if (status == 0) // cell is cut
+        {
+            interfaceArea_[celli] =
+                mag(cutCell.isoFaceArea())/mesh.V()[celli];
+        }
+    }
+
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<wallPolyPatch>(pbm[patchi]))
+        {
+            const polyPatch& pp = pbm[patchi];
+            forAll(pp.faceCells(), faceI)
+            {
+                const label pCelli = pp.faceCells()[faceI];
+                bool interface(false);
+                if
+                (
+                     sign(R_.value()) > 0
+                 && (T[pCelli] - Tactivate_.value()) > 0
+                )
+                {
+                    interface = true;
+                }
+
+                if
+                (
+                    sign(R_.value()) < 0
+                 && (T[pCelli] - Tactivate_.value()) < 0
+                )
+                {
+                    interface = true;
+                }
+
+                if
+                (
+                    interface
+                 &&
+                    (
+                        alpha[pCelli] < 2*isoAlpha_
+                     && alpha[pCelli] > 0.5*isoAlpha_
+                    )
+                )
+                {
+                    interfaceArea_[pCelli] =
+                        mag(pp.faceAreas()[faceI])/mesh.V()[pCelli];
+                }
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Thermo, class OtherThermo>
+Foam::meltingEvaporationModels::interfaceHeatResistance<Thermo, OtherThermo>
+::interfaceHeatResistance
+(
+    const dictionary& dict,
+    const phasePair& pair
+)
+:
+    InterfaceCompositionModel<Thermo, OtherThermo>(dict, pair),
+    R_("R", dimPower/dimArea/dimTemperature, dict),
+    Tactivate_("Tactivate", dimTemperature, dict),
+    interfaceArea_
+    (
+        IOobject
+        (
+            "interfaceArea",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimless/dimLength, Zero)
+    ),
+    mDotc_
+    (
+        IOobject
+        (
+            "mDotc",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+    mDotcSpread_
+    (
+        IOobject
+        (
+            "mDotcSpread",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+    htc_
+    (
+        IOobject
+        (
+            "htc",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimMass/dimArea/dimTemperature/dimTime, Zero)
+    ),
+    isoAlpha_(dict.lookupOrDefault<scalar>("isoAlpha", 0.5)),
+    spread_(dict.lookupOrDefault<scalar>("spread", 3))
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::interfaceHeatResistance<Thermo, OtherThermo>
+::Kexp(const volScalarField& T)
+{
+
+    const fvMesh& mesh = this->mesh_;
+
+    updateInterface(T);
+
+    tmp<volScalarField> tdeltaT
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "tdeltaT",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(dimTemperature, Zero)
+        )
+    );
+    volScalarField& deltaT = tdeltaT.ref();
+
+    dimensionedScalar T0("T0", dimTemperature, Zero);
+
+    if (sign(R_.value()) > 0)
+    {
+        deltaT = max(T - Tactivate_, T0);
+    }
+    else
+    {
+        deltaT = max(Tactivate_ - T, T0);
+    }
+
+    word fullSpeciesName = this->transferSpecie();
+    auto tempOpen = fullSpeciesName.find('.');
+    const word speciesName(fullSpeciesName.substr(0, tempOpen));
+
+    tmp<volScalarField> L = this->L(speciesName, T);
+
+    htc_ = R_/L();
+
+    const volScalarField& to = this->pair().to();
+    const volScalarField& from = this->pair().from();
+
+    dimensionedScalar D
+    (
+        "D",
+        dimArea,
+        spread_/sqr(gAverage(this->mesh_.nonOrthDeltaCoeffs()))
+    );
+
+    const dimensionedScalar MdotMin("MdotMin", mDotc_.dimensions(), 1e-3);
+
+    if (max(mDotc_) > MdotMin)
+    {
+        fvc::spreadSource
+        (
+            mDotcSpread_,
+            mDotc_,
+            from,
+            to,
+            D,
+            1e-3
+        );
+    }
+
+    mDotc_ = interfaceArea_*htc_*deltaT;
+
+    return tmp<volScalarField>(new volScalarField(mDotc_));
+}
+
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::interfaceHeatResistance<Thermo, OtherThermo>
+::KSp
+(
+    label variable,
+    const volScalarField& refValue
+)
+{
+    if (this->modelVariable_ == variable)
+    {
+        const volScalarField coeff(htc_*interfaceArea_);
+
+        if (sign(R_.value()) > 0)
+        {
+            return(coeff*pos(refValue - Tactivate_));
+        }
+        else
+        {
+            return(coeff*pos(Tactivate_ - refValue));
+        }
+    }
+    else
+    {
+        return tmp<volScalarField> ();
+    }
+}
+
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::interfaceHeatResistance<Thermo, OtherThermo>
+::KSu
+(
+    label variable,
+    const volScalarField& refValue
+)
+{
+    if (this->modelVariable_ == variable)
+    {
+        const volScalarField coeff(htc_*interfaceArea_*Tactivate_);
+
+        if (sign(R_.value()) > 0)
+        {
+            return(-coeff*pos(refValue - Tactivate_));
+        }
+        else
+        {
+            return(coeff*pos(Tactivate_ - refValue));
+        }
+    }
+    else if (interfaceCompositionModel::P == variable)
+    {
+        return tmp<volScalarField>(new volScalarField(mDotcSpread_));
+    }
+    else
+    {
+        return tmp<volScalarField> ();
+    }
+}
+
+
+template<class Thermo, class OtherThermo>
+const Foam::dimensionedScalar&
+Foam::meltingEvaporationModels::interfaceHeatResistance<Thermo, OtherThermo>
+::Tactivate() const
+{
+    return Tactivate_;
+}
+
+
+template<class Thermo, class OtherThermo>
+bool
+Foam::meltingEvaporationModels::
+interfaceHeatResistance<Thermo, OtherThermo>::includeDivU()
+{
+    return true;
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.H
new file mode 100644
index 0000000000000000000000000000000000000000..3c49e725107e52257b1ef916c7077abc88944399
--- /dev/null
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceHeatResistance/interfaceHeatResistance.H
@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 Henning Scheufler
+-------------------------------------------------------------------------------
+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::meltingEvaporationModels::interfaceHeatResistance
+
+Description
+    Interface Heat Resistance type of condensation/saturation model using
+    spread source distribution following:
+
+    References:
+    \verbatim
+        Hardt, S., Wondra, F. (2008).
+        Evaporation model for interfacial flows based on a continuum-
+        field representation of the source term
+        Journal of Computational Physics 227 (2008), 5871-5895
+    \endverbatim
+
+
+Usage
+
+    Example usage:
+    \verbatim
+    massTransferModel
+    (
+        (liquid to gas)
+        {
+            type                interfaceHeatResistance;
+            R                   2e6;
+            Tactivate           373;
+        }
+    );
+    \endverbatim
+
+    where:
+    \table
+        Property    | Description             | Required    | Default value
+        R           | Heat transfer coefficient | yes
+        includeVolChange    | Volumen change  | no          | yes
+        isoAlpha    | iso-alpha for interface | no          | 0.5
+        Tactivate   | Saturation temperature  | yes
+        species     | Specie name on the other phase | no   | none
+        spread      | Cells to spread the source for pEq | no   | 3
+    \endtable
+
+SourceFiles
+    interfaceHeatResistance.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meltingEvaporationModels_interfaceHeatResistance_H
+#define meltingEvaporationModels_interfaceHeatResistance_H
+
+
+#include "InterfaceCompositionModel.H"
+
+// * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * * * * *//
+
+namespace Foam
+{
+
+class phasePair;
+
+namespace meltingEvaporationModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class interfaceHeatResistance
+\*---------------------------------------------------------------------------*/
+
+template<class Thermo, class OtherThermo>
+class interfaceHeatResistance
+:
+    public InterfaceCompositionModel<Thermo, OtherThermo>
+{
+    // Private data
+
+        //- Heat transfer coefficient [1/s/K]
+        dimensionedScalar R_;
+
+        //- Activation temperature
+        const dimensionedScalar Tactivate_;
+
+        //- Interface area
+        volScalarField interfaceArea_;
+
+        //- Mass source
+        volScalarField mDotc_;
+
+        //- Spread mass source
+        volScalarField mDotcSpread_;
+
+        //- Heat transfer coefficient
+        volScalarField htc_;
+
+        //- Interface Iso-value
+        scalar isoAlpha_;
+
+        //- Spread for mass source
+        scalar spread_;
+
+
+    // Private member
+
+        //- Update interface
+        void updateInterface(const volScalarField& T);
+
+public:
+
+    //- Runtime type information
+    TypeName("interfaceHeatResistance");
+
+
+    // Constructors
+
+        //- Construct from components
+        interfaceHeatResistance
+        (
+            const dictionary& dict,
+            const phasePair& pair
+        );
+
+
+    //- Destructor
+    virtual ~interfaceHeatResistance() = default;
+
+
+    // Member Functions
+
+        //- Explicit total mass transfer coefficient
+        virtual tmp<volScalarField> Kexp
+        (
+            const volScalarField& field
+        );
+
+        //- Implicit mass transfer coefficient
+        virtual tmp<volScalarField> KSp
+        (
+            label modelVariable,
+            const volScalarField& field
+        );
+
+        //- Explicit mass transfer coefficient
+        virtual tmp<volScalarField> KSu
+        (
+            label modelVariable,
+            const volScalarField& field
+        );
+
+        //- Return Tactivate
+        virtual const dimensionedScalar& Tactivate() const;
+
+        //- Adds and substract alpha*div(U) as a source term
+        //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
+        virtual bool includeDivU();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace meltingEvaporationModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "interfaceHeatResistance.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C
index 3c544fdaeaa6923e8f7f79ac80ccab298d7fc1f1..814fa9d8622d3ae40883b6a655177ea8c81b2973 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,16 +27,90 @@ License
 
 #include "kineticGasEvaporation.H"
 #include "constants.H"
-#include "fvcGrad.H"
-#include "fvcSnGrad.H"
-#include "fvcDiv.H"
-#include "surfaceInterpolate.H"
-#include "fvcReconstruct.H"
-#include "fvm.H"
-#include "zeroGradientFvPatchFields.H"
+#include "isoCutCell.H"
+#include "volPointInterpolation.H"
+#include "wallPolyPatch.H"
+#include "fvcSmooth.H"
 
 using namespace Foam::constant;
 
+
+// * * * * * * * * * * *  Private Member Functions * * * * * * * * * * * * * //
+
+template<class Thermo, class OtherThermo>
+void Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
+::updateInterface(const volScalarField& T)
+{
+    const fvMesh& mesh = this->mesh_;
+
+    const volScalarField& alpha = this->pair().from();
+
+    scalarField ap
+    (
+        volPointInterpolation::New(mesh).interpolate(alpha)
+    );
+
+    isoCutCell cutCell(mesh, ap);
+
+    forAll(interfaceArea_, celli)
+    {
+        label status = cutCell.calcSubCell(celli, isoAlpha_);
+        interfaceArea_[celli] = 0;
+        if (status == 0) // cell is cut
+        {
+            interfaceArea_[celli] =
+                mag(cutCell.isoFaceArea())/mesh.V()[celli];
+        }
+    }
+
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<wallPolyPatch>(pbm[patchi]))
+        {
+            const polyPatch& pp = pbm[patchi];
+            forAll(pp.faceCells(), faceI)
+            {
+                const label pCelli = pp.faceCells()[faceI];
+                bool interface(false);
+                if
+                (
+                     sign(C_.value()) > 0
+                 && (T[pCelli] - Tactivate_.value()) > 0
+                )
+                {
+                    interface = true;
+                }
+
+                if
+                (
+                    sign(C_.value()) < 0
+                 && (T[pCelli] - Tactivate_.value()) < 0
+                )
+                {
+                    interface = true;
+                }
+
+                if
+                (
+                    interface
+                 &&
+                    (
+                        alpha[pCelli] < 2*isoAlpha_
+                     && alpha[pCelli] > 0.5*isoAlpha_
+                    )
+                )
+                {
+                    interfaceArea_[pCelli] =
+                        mag(pp.faceAreas()[faceI])/mesh.V()[pCelli];
+                }
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Thermo, class OtherThermo>
@@ -51,28 +125,59 @@ Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
     C_("C", dimless, dict),
     Tactivate_("Tactivate", dimTemperature, dict),
     Mv_("Mv", dimMass/dimMoles, -1, dict),
-    alphaMax_(dict.lookupOrDefault<scalar>("alphaMax", 1.0)),
-    alphaMin_(dict.lookupOrDefault<scalar>("alphaMin", 0.5)),
-    alphaRestMax_(dict.lookupOrDefault<scalar>("alphaRestMax", 0.01))
+    interfaceArea_
+    (
+        IOobject
+        (
+            "interfaceArea",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimless/dimLength, Zero)
+    ),
+    htc_
+    (
+        IOobject
+        (
+            "htc",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimMass/dimArea/dimTemperature/dimTime, Zero)
+    ),
+    mDotc_
+    (
+        IOobject
+        (
+            "mDotc",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+    isoAlpha_(dict.lookupOrDefault<scalar>("isoAlpha", 0.5))
 {
-    if (this->transferSpecie() != "none")
-    {
-        word fullSpeciesName = this->transferSpecie();
-        auto tempOpen = fullSpeciesName.find('.');
-        const word speciesName(fullSpeciesName.substr(0, tempOpen));
+    word speciesName = IOobject::member(this->transferSpecie());
 
-        // Get the "to" thermo
-        const typename OtherThermo::thermoType& toThermo =
-            this->getLocalThermo
-            (
-                speciesName,
-                this->toThermo_
-            );
-
-         // Convert from g/mol to Kg/mol
-        Mv_.value() = toThermo.W()*1e-3;
-    }
+    // Get the "to" thermo
+    const typename OtherThermo::thermoType& toThermo =
+        this->getLocalThermo
+        (
+            speciesName,
+            this->toThermo_
+        );
 
+    // Convert from g/mol to Kg/mol
+    Mv_.value() = toThermo.W()*1e-3;
 
     if (Mv_.value() == -1)
     {
@@ -88,186 +193,132 @@ Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
 template<class Thermo, class OtherThermo>
 Foam::tmp<Foam::volScalarField>
 Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
-::Kexp(const volScalarField& field)
+::Kexp(const volScalarField& T)
 {
-    {
-        const volScalarField& to = this->pair().to();
-
-        const volScalarField& from = this->pair().from();
 
-        const fvMesh& mesh = this->mesh_;
+    const fvMesh& mesh = this->mesh_;
 
-        const volScalarField& T =
-            mesh.lookupObject<volScalarField>("T").oldTime();
-
-        const dimensionedScalar HerztKnudsConst
+    const dimensionedScalar HerztKnudsConst
+    (
+        sqrt
         (
-            sqrt
-            (
-                Mv_
-               /2.0
-               /constant::physicoChemical::R
-               /mathematical::pi
-               /pow3(Tactivate_)
-            )
-        );
-
-        word fullSpeciesName = this->transferSpecie();
-        auto tempOpen = fullSpeciesName.find('.');
-        const word speciesName(fullSpeciesName.substr(0, tempOpen));
+            2.0*mathematical::pi
+          * pow3(Tactivate_)
+          * constant::physicoChemical::R/Mv_
+        )
+    );
 
-        tmp<volScalarField> L = this->L(speciesName, field);
+    word speciesName = IOobject::member(this->transferSpecie());
+    tmp<volScalarField> L = this->L(speciesName, T);
 
-        const volVectorField gradFrom(fvc::grad(from));
-        const volVectorField gradTo(fvc::grad(to));
+    updateInterface(T);
 
-        const volScalarField areaDensity("areaDensity", mag(gradFrom));
-
-        const volScalarField gradAlphaf(gradFrom & gradTo);
-
-        volScalarField Tmask("Tmask", from*0.0);
-
-        forAll(Tmask, celli)
-        {
-            if (gradAlphaf[celli] < 0)
-            {
-                if (from[celli] > alphaMin_ && from[celli] < alphaMax_)
-                {
-                    {
-                        scalar alphaRes = 1.0 - from[celli] - to[celli];
-                        if (alphaRes < alphaRestMax_)
-                        {
-                            Tmask[celli] = 1.0;
-                        }
-                    }
-                }
-            }
-        }
-
-        tmp<volScalarField> tRhom
+    tmp<volScalarField> tRhov
+    (
+        new volScalarField
         (
-            new volScalarField
+            IOobject
             (
-                IOobject
-                (
-                    "trhom",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh,
-                dimensionedScalar(dimDensity, Zero)
-            )
-        );
-        volScalarField& rhom = tRhom.ref();
-
-        tmp<volScalarField> tTdelta
+                "tRhov",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(dimDensity, Zero)
+        )
+    );
+    volScalarField& rhov = tRhov.ref();
+
+    tmp<volScalarField> tdeltaT
+    (
+        new volScalarField
         (
-            new volScalarField
+            IOobject
             (
-                IOobject
-                (
-                    "trhom",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh,
-                dimensionedScalar(dimTemperature, Zero)
-            )
-        );
-        volScalarField& tDelta = tTdelta.ref();
-
-        if (sign(C_.value()) > 0)
-        {
-            rhom =
-                this->pair().to().rho()*this->pair().from().rho()
-              / (this->pair().from().rho() - this->pair().to().rho());
-
-            tDelta = max
-            (
-                (T*Tmask - Tactivate_),
-                dimensionedScalar("T0", dimTemperature, Zero)
-            );
-        }
-        else
-        {
-            rhom =
-                this->pair().to().rho()*this->pair().from().rho()
-              / (this->pair().to().rho() - this->pair().from().rho());
-
-            tDelta = max
-            (
-                Tmask*(Tactivate_ - T),
-                dimensionedScalar("T0", dimTemperature, Zero)
-            );
-        }
-
-        volScalarField massFluxEvap
-        (
-            "massFluxEvap",
-            2*mag(C_)/(2 - mag(C_))
-          * HerztKnudsConst
-          * L()
-          * rhom
-          * tDelta
-        );
-
-        // 'from' phase normalization
-        // WIP: Normalization could be convinient for cases where the area were
-        // the source term is calculated is uniform
-        const dimensionedScalar Nl
-        (
-            gSum((areaDensity*mesh.V())())
-           /(
-               gSum
-               (
-                   ((areaDensity*from)*mesh.V())()
-               )
-             + dimensionedScalar("SMALL", dimless, VSMALL)
-            )
-        );
+                "tdeltaT",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(dimTemperature, Zero)
+        )
+    );
+    volScalarField& deltaT = tdeltaT.ref();
+
+    dimensionedScalar T0("T0", dimTemperature, Zero);
+
+    if (sign(C_.value()) > 0)
+    {
+        rhov = this->pair().to().rho();
+        deltaT = max(T - Tactivate_, T0);
+    }
+    else
+    {
+        rhov = this->pair().from().rho();
+        deltaT = max(Tactivate_ - T, T0);
+    }
 
+    htc_ = 2*mag(C_)/(2-mag(C_))*(L()*rhov/HerztKnudsConst);
 
-        if (mesh.time().outputTime() && debug)
-        {
-            areaDensity.write();
-            Tmask.write();
-            volScalarField mKGasDot
-            (
-                "mKGasDot",
-                massFluxEvap*areaDensity*Nl*from
-            );
-            mKGasDot.write();
-        }
+    mDotc_ = htc_*deltaT*interfaceArea_;
 
-        return massFluxEvap*areaDensity*Nl*from;
-    }
+    return tmp<volScalarField>(new volScalarField(mDotc_));
 }
 
+
 template<class Thermo, class OtherThermo>
 Foam::tmp<Foam::volScalarField>
-Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>::KSu
+Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>::KSp
 (
     label variable,
     const volScalarField& refValue
 )
 {
-    return tmp<volScalarField> ();
+    if (this->modelVariable_ == variable)
+    {
+        const volScalarField coeff(htc_*interfaceArea_);
+
+        if (sign(C_.value()) > 0)
+        {
+            return(coeff*pos(refValue - Tactivate_));
+        }
+        else
+        {
+            return(coeff*pos(Tactivate_ - refValue));
+        }
+    }
+    else
+    {
+        return tmp<volScalarField> ();
+    }
 }
 
 
 template<class Thermo, class OtherThermo>
 Foam::tmp<Foam::volScalarField>
-Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>::KSp
+Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>::KSu
 (
     label variable,
     const volScalarField& refValue
 )
 {
-    return tmp<volScalarField> ();
+    if (this->modelVariable_ == variable)
+    {
+        const volScalarField coeff(htc_*interfaceArea_*Tactivate_);
+
+        if (sign(C_.value()) > 0)
+        {
+            return(-coeff*pos(refValue - Tactivate_));
+        }
+        else
+        {
+            return(coeff*pos(Tactivate_ - refValue));
+        }
+    }
+    else
+    {
+        return tmp<volScalarField> ();
+    }
 }
 
 
@@ -280,4 +331,12 @@ Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
 }
 
 
+template<class Thermo, class OtherThermo>
+bool
+Foam::meltingEvaporationModels::
+kineticGasEvaporation<Thermo, OtherThermo>::includeDivU()
+{
+    return true;
+}
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H
index 417c51f2269687cf94dad33ce5bf41d83f19aa3d..26f293fe07fe850fccb5a72d07d469bf1c0dbd2a 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,9 +66,8 @@ Description
     \f[
         Flux =
            2 \frac{C}{2 - C}
-           \sqrt{\frac{M}{2 \pi R T_{activate}}}
-           L (\rho_{v}*\rho_{l}/(\rho_{l} - \rho_{v}))
-           (T - T_{activate})/T_{activate}
+           \sqrt{\frac{M}{2 \pi R {T_activate}^3}} L (\rho_{v})
+           (T - T_{activate})
     \f]
 
     This assumes liquid and vapour are in equilibrium, then the accomodation
@@ -91,8 +90,7 @@ Usage
             type                kineticGasEvaporation;
             species             vapour.gas;
             C                   0.1;
-            alphaMin            0.0;
-            alphaMax            0.2;
+            isoAlpha            0.1;
             Tactivate           373;
         }
     );
@@ -101,12 +99,12 @@ Usage
     where:
     \table
         Property    | Description             | Required    | Default value
-        C           | Accomodation coefficient (C > 0 for evaporation, C < 0 for
+        C           | Coefficient (C > 0 for evaporation, C < 0 for
             condensation) | yes
-        alphaMin    | Minimum value of alpha  | no          | 0.0
-        alphaMax    | Maximum values of alpha | no          | 0.5
+        includeVolChange    | Volumen change  | no          | yes
+        isoAlpha    | iso-alpha for interface | no          | 0.5
         Tactivate   | Saturation temperature  | yes
-        species     | Specie name on the other phase | yes
+        species     | Specie name on the other phase | no   | none
     \endtable
 
 SourceFiles
@@ -150,14 +148,23 @@ class kineticGasEvaporation
         //- Molar weight of the vapour in the continous phase
         dimensionedScalar Mv_;
 
-        //- 'To' phase maximum value for the mass transfer
-        scalar alphaMax_;
+        //- Interface area
+        volScalarField interfaceArea_;
 
-        //- 'To' phase minumum value for the mass transfer
-        scalar alphaMin_;
+        //- Heat transfer coefficient
+        volScalarField htc_;
 
-        //- Alpha maximum for the rest of phases
-        scalar alphaRestMax_;
+        //- Mass source
+        volScalarField mDotc_;
+
+        //- Interface Iso-value
+        scalar isoAlpha_;
+
+
+    // Private member
+
+        //- Update interface
+        void updateInterface(const volScalarField& T);
 
 
 public:
@@ -182,7 +189,7 @@ public:
 
     // Member Functions
 
-        //- Explicit mass transfer coefficient
+        //- Explicit total mass transfer coefficient
         virtual tmp<volScalarField> Kexp
         (
             const volScalarField& field
@@ -204,6 +211,10 @@ public:
 
         //- Return Tactivate
         virtual const dimensionedScalar& Tactivate() const;
+
+        //- Adds and substract alpha*div(U) as a source term
+        //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
+        virtual bool includeDivU();
 };
 
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
index 8128eb068a40027d6394c2ef07d18ae97ca53f69..ad860631ce800e85a8dce6d8553f81ffb301ce82 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H
@@ -33,9 +33,13 @@
         (
             fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
-          + fluid.volTransfer(p_rgh)
         );
 
+        if (fluid.includeVolChange())
+        {
+            p_rghEqn += fluid.volTransfer(p_rgh);
+        }
+
         p_rghEqn.setReference(pRefCell, pRefValue);
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
index 99777692bf1747cd5c08552a311be458536acff5..0abc85f561b2b78e063a095d2d35ca64bd891cfa 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-202 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -396,7 +396,12 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
 
             if (KSp.valid())
             {
-                Sp += KSp.ref();
+                Sp -=
+                    KSp.ref()
+                   *(
+                        - this->coeffs(phase1.name())
+                        + this->coeffs(phase2.name())
+                    );
             }
 
             tmp<volScalarField> KSu =
@@ -404,7 +409,12 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
 
             if (KSu.valid())
             {
-                Su += KSu.ref();
+                Su -=
+                    KSu.ref()
+                   *(
+                        - this->coeffs(phase1.name())
+                        + this->coeffs(phase2.name())
+                    );
             }
 
             // If linearization is not provided used full explicit
@@ -436,7 +446,12 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
 
             if (KSp.valid())
             {
-                Sp += KSp.ref();
+                Sp +=
+                    KSp.ref()
+                   *(
+                        - this->coeffs(phase1.name())
+                        + this->coeffs(phase2.name())
+                    );
             }
 
             tmp<volScalarField> KSu =
@@ -444,7 +459,12 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
 
             if (KSu.valid())
             {
-                Su += KSu.ref();
+                Su +=
+                    KSu.ref()
+                   *(
+                        - this->coeffs(phase1.name())
+                        + this->coeffs(phase2.name())
+                    );
             }
 
             // If linearization is not provided used full explicit
@@ -744,4 +764,18 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::massSpeciesTransfer
 }
 
 
+template<class BasePhaseSystem>
+bool Foam::MassTransferPhaseSystem<BasePhaseSystem>::includeVolChange()
+{
+    bool includeVolChange(true);
+    forAllIters(massTransferModels_, iter)
+    {
+        if (!iter()->includeVolChange())
+        {
+            includeVolChange = false;
+        }
+    }
+    return includeVolChange;
+}
+
 // ************************************************************************* //
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H
index 030d02e99290c339588aca3e97dfdf66ccd45298..0b2aa5b6f35d688e43e1d6f9e6ef6217e7302089 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H
@@ -127,12 +127,15 @@ public:
     // Mass transfer functions
 
         //- Return the heat transfer matrix
+        // NOTE: Call KSu and KSp with T as variable,if not provided uses dmdt.
         virtual  tmp<fvScalarMatrix> heatTransfer(const volScalarField& T);
 
         //- Return the volumetric rate transfer matrix
+        // NOTE: Call KSu and KSp with p as variable,if not provided uses dmdt.
         virtual  tmp<fvScalarMatrix> volTransfer(const volScalarField& p);
 
         //- Correct/calculates mass sources dmdt for phases
+        // NOTE: Call the kexp() for all the mass transfer models.
         virtual void correctMassSources(const volScalarField& T);
 
         //- Calculate mass transfer for alpha's
@@ -147,6 +150,8 @@ public:
             const word speciesName
         );
 
+        //- Add volume change in pEq
+        virtual bool includeVolChange();
 };
 
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H
index a5394939baf6b36efcc85e574357ae7239862a7d..3d12b3308dbc81fe7b75052dbe509bf3e8058200 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.H
@@ -529,6 +529,9 @@ public:
             const word speciesName
         ) = 0;
 
+        //- Add volume change in pEq
+        virtual bool includeVolChange() = 0;
+
 
     // Solve phases and correct models
 
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
index b3442ed9b83cbf8ea9094dafa346529ebb000532..8ca6ef800af05583f0079bfae5f6d68988891f11 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C
@@ -110,6 +110,34 @@ interfaceHeatResistance
         dimensionedScalar(dimDensity/dimTime, Zero)
     ),
 
+    mDotcSpread_
+    (
+        IOobject
+        (
+            "mDotcSpread",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+
+    mDoteSpread_
+    (
+        IOobject
+        (
+            "mDoteSpread",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+
     spread_
     (
         optionalSubDict(type() + "Coeffs").get<scalar>("spread")
@@ -259,6 +287,47 @@ correct()
         mDote_[celli] = min(max(mDote_[celli], maxCond), maxEvap);
         mDotc_[celli] = min(max(mDotc_[celli], maxCond), maxEvap);
     }
+
+    // Calculate the spread sources
+
+    dimensionedScalar D
+    (
+        "D",
+        dimArea,
+        spread_/sqr(gAverage(mesh_.nonOrthDeltaCoeffs()))
+    );
+
+
+    const volScalarField& alpha1 = mixture_.alpha1();
+    const volScalarField& alpha2 = mixture_.alpha2();
+
+    const dimensionedScalar MDotMin("MdotMin", mDotc_.dimensions(), 1e-3);
+
+    if (max(mDotc_) > MDotMin)
+    {
+        fvc::spreadSource
+        (
+            mDotcSpread_,
+            mDotc_,
+            alpha1,
+            alpha2,
+            D,
+            1e-3
+        );
+    }
+
+    if (max(mDote_) > MDotMin)
+    {
+        fvc::spreadSource
+        (
+            mDoteSpread_,
+            mDote_,
+            alpha1,
+            alpha2,
+            D,
+            1e-3
+        );
+    }
 }
 
 
@@ -325,63 +394,12 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
 vDot() const
 {
 
-    dimensionedScalar D
-    (
-        "D",
-        dimArea,
-        spread_/sqr(gAverage(mesh_.nonOrthDeltaCoeffs()))
-    );
-
-
-    const volScalarField& alpha1 = mixture_.alpha1();
-    const volScalarField& alpha2 = mixture_.alpha2();
-
-    const dimensionedScalar MDotMin("MdotMin", mDotc_.dimensions(), 1e-3);
-
-    Pair<tmp<volScalarField>> mDotSpread
-    (
-        tmp<volScalarField>(mDotc_*0.0),
-        tmp<volScalarField>(mDote_*0.0)
-    );
-
-    if (max(mDotc_) > MDotMin)
-    {
-        fvc::spreadSource
-        (
-            mDotSpread[0].ref(),
-            mDotc_,
-            alpha1,
-            alpha2,
-            D,
-            1e-3
-        );
-    }
-
-    if (max(mDote_) > MDotMin)
-    {
-        fvc::spreadSource
-        (
-            mDotSpread[1].ref(),
-            mDote_,
-            alpha1,
-            alpha2,
-            D,
-            1e-3
-        );
-    }
-
     dimensionedScalar pCoeff(1.0/mixture_.rho1() - 1.0/mixture_.rho2());
 
-    if (mesh_.time().outputTime())
-    {
-        volScalarField mDotS("mDotSpread", mDotSpread[1].ref());
-        mDotS.write();
-    }
-
     return Pair<tmp<volScalarField>>
     (
-        pCoeff*mDotSpread[0],
-       -pCoeff*mDotSpread[1]
+        pCoeff*mDotcSpread_,
+       -pCoeff*mDoteSpread_
     );
 }
 
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H
index 67c0aa5676d9d2f3b363de4b9e868d1afd58a264..c29cedb8201842a6b95f62b30edaaa6c1056b183 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H
@@ -78,6 +78,12 @@ class interfaceHeatResistance
         //- Mass evaporation source
         volScalarField mDote_;
 
+        //- Spread condensation mass source
+        volScalarField mDotcSpread_;
+
+        //- Spread evaporation mass source
+        volScalarField mDoteSpread_;
+
         //- Spread for mass source
         scalar spread_;
 
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/T b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/T
new file mode 100644
index 0000000000000000000000000000000000000000..6f96fca3087046ed16daad7a73169060cc1c2751
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/T
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 371;
+
+boundaryField
+{
+    left
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+    right
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+    bottom
+    {
+        type            externalWallHeatFluxTemperature;
+        mode            flux;
+        q               uniform  18.5e3;
+        kappaMethod     fluidThermo;
+        value           $internalField;
+    }
+    top
+    {
+        type            inletOutlet;
+        value           uniform 371;
+        inletValue      uniform 371;
+    }
+    frontAndBack
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/U b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..0ac96679a9a2f3c37be5fa98bd33e7891553dc4a
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    left
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    right
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    bottom
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    top
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+    frontAndBack
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/air.gas b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/air.gas
new file mode 100644
index 0000000000000000000000000000000000000000..0baf7b0d99e792ed076ab433a29186aee4bdd502
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/air.gas
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      air.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    left
+    {
+        type            zeroGradient;
+    }
+    right
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/alpha.gas b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/alpha.gas
new file mode 100644
index 0000000000000000000000000000000000000000..308f1b3ccd65ef7c2e0204b2cde9419a0a41cfb6
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/alpha.gas
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+     left
+    {
+        type            zeroGradient;
+    }
+    right
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/alpha.liquid b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/alpha.liquid
new file mode 100644
index 0000000000000000000000000000000000000000..660caff82bf63322290385d3ff150502f4d47ba0
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/alpha.liquid
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.liquid;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    left
+    {
+        type            zeroGradient;
+    }
+    right
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/p_rgh b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..7603b7ee489c82912846d3bc2935a64a30c8f38a
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/p_rgh
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [ 1 -1 -2 0 0 0 0 ];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    left
+    {
+        type            fixedFluxPressure;
+        value           uniform 1e5;
+    }
+    right
+    {
+        type            fixedFluxPressure;
+        value           uniform 1e5;
+    }
+    bottom
+    {
+        type            fixedFluxPressure;
+        value           uniform 1e5;
+    }
+    top
+    {
+        type            prghTotalPressure;
+        p0              $internalField;
+        value           uniform 100000;
+    }
+    frontAndBack
+    {
+        type            fixedFluxPressure;
+        value           uniform 1e5;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/vapour.gas b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/vapour.gas
new file mode 100644
index 0000000000000000000000000000000000000000..15613b0907b17af7095b4075e5565557415ef3e5
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/0.orig/vapour.gas
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      vapour.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    left
+    {
+        type            zeroGradient;
+    }
+    right
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/Allclean b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..fb1f3847301c377e02e12439ba58cbf303af3ef9
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/Allclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/Allrun b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..e2f1c31545c467e93adcf482c783e8e908c4ce9c
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/Allrun
@@ -0,0 +1,14 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+runApplication blockMesh
+runApplication setFields
+runApplication decomposePar
+runParallel $(getApplication)
+runApplication reconstructPar
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/README b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/README
new file mode 100644
index 0000000000000000000000000000000000000000..7043c87bd6ccda684491a7b2635aad3e49754657
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/README
@@ -0,0 +1,7 @@
+kineticGasEvaporation model test case
+
+Still evaporating pool with heat flux from bottom.
+
+The theoretical area-averaged mass flux evaporating is 1.2e-4 Kg/s.
+This averaged value is reached approximately at 110 secs 
+
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/g b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/g
new file mode 100755
index 0000000000000000000000000000000000000000..2dbd4439740a48b20a2813a31aa2af62ac525c60
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/g
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/phaseProperties b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/phaseProperties
new file mode 100755
index 0000000000000000000000000000000000000000..3d06c5bf0165a3a7b815f0dab97590b6120084fb
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/phaseProperties
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      phaseProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+type    massTransferMultiphaseSystem;
+
+phases  (gas liquid);
+
+liquid
+{
+    type            pureMovingPhaseModel;
+}
+
+gas
+{
+    type            multiComponentMovingPhaseModel;
+}
+
+surfaceTension
+(
+    (gas and liquid)
+    {
+        type            constant;
+        sigma           0.00;
+    }
+);
+
+
+massTransferModel
+(
+    (liquid to gas)
+    {
+        type                kineticGasEvaporation;
+        species             vapour.gas;
+        C                   0.1;
+        isoAlpha            0.1;
+        Tactivate           372;
+   }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/thermophysicalProperties.gas b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/thermophysicalProperties.gas
new file mode 100755
index 0000000000000000000000000000000000000000..27c24984eb4ff3179384e45b6a4d86ebc4af42ea
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/thermophysicalProperties.gas
@@ -0,0 +1,83 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         multiComponentMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState incompressiblePerfectGas;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+// ************************************************************************* //
+
+species
+(
+    air
+    vapour
+);
+
+inertSpecie air;
+
+vapour
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18.9;
+    }
+    equationOfState
+    {
+        pRef         1e5;
+    }
+    thermodynamics
+    {
+        Hf          0;
+        Cp          2030;
+    }
+    transport
+    {
+        mu          0.9e-05;
+        Pr          0.7;
+    }
+}
+
+air
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    equationOfState
+    {
+        pRef         1e5;
+    }
+    thermodynamics
+    {
+        Hf          0;
+        Cp          900;
+    }
+    transport
+    {
+        mu          0.9e-05;
+        Pr          0.7;
+    }
+}
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/thermophysicalProperties.liquid b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/thermophysicalProperties.liquid
new file mode 100755
index 0000000000000000000000000000000000000000..1c3f30ed7138e2e42a64b18e1f0d43be1e4b7833
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/thermophysicalProperties.liquid
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.liquid;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       18.9;
+    }
+    equationOfState
+    {
+        rho             958.4;
+    }
+    thermodynamics
+    {
+        Cp              4216;
+        Hf              2.45e6;
+    }
+    transport
+    {
+        Pr              6.62;
+        mu              959e-6;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/turbulenceProperties b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/turbulenceProperties
new file mode 100755
index 0000000000000000000000000000000000000000..e35639d9fc92bd8d62ed792ef33fab2811f2a1ed
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/blockMeshDict b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/blockMeshDict
new file mode 100755
index 0000000000000000000000000000000000000000..58fc138313f970898b43ed13da469a21996f9e41
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/blockMeshDict
@@ -0,0 +1,89 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.25;
+
+vertices
+(
+    (0 0 0)
+    (0.5 0 0)
+    (0.5 0.5 0)
+    (0 0.5 0)
+    (0 0 0.5)
+    (0.5 0 0.5)
+    (0.5 0.5 0.5)
+    (0 0.5 0.5)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (50 50 50) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    left
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    right
+    {
+        type wall;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    bottom
+    {
+        type wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    top
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    frontAndBack
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/controlDict b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/controlDict
new file mode 100755
index 0000000000000000000000000000000000000000..98b0ee2ffc7e1b9d6ebf51cb52dd6a90ddab9071
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/controlDict
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoReactingMultiphaseInterFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          1e-3;
+
+writeControl    adjustableRunTime;
+
+writeInterval   5;
+
+purgeWrite      4;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+compression     off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxDeltaT       1e-1;
+
+maxCo           3;
+maxAlphaCo      2;
+maxAlphaDdt     1;
+
+functions
+{
+
+    mass
+    {
+        type            volFieldValue;
+        libs            (fieldFunctionObjects);
+
+        writeControl    timeStep;
+        writeInterval   10;
+        writeFields     false;
+        log             true;
+
+        operation       volIntegrate;
+
+        fields
+        (
+            dmdt.liquidToGas
+        );
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/decomposeParDict b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/decomposeParDict
new file mode 100755
index 0000000000000000000000000000000000000000..087c88f826f9b943192bf07b36853e87951c37e3
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          simple;
+
+simpleCoeffs
+{
+    n               (2 2 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (1 1 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/fvSchemes b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/fvSchemes
new file mode 100755
index 0000000000000000000000000000000000000000..76863fc98b4cd86cff77dd8a5c923093f7ee2709
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/fvSchemes
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+
+    div(rhoPhi,U)           Gauss upwind;
+
+    "div\(phi,alpha.*\)"    Gauss vanLeer;
+    "div\(phir,alpha.*\)"   Gauss linear;
+
+    "div\(Yiphir,alpha.*\)" Gauss linear;
+    "div\(phi,.*\.gas.*\)"  Gauss vanLeer;
+
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+
+    div(phi,T)              Gauss upwind;
+
+    div((alpha.gas*U))      Gauss linear;
+    div((alpha.liquid*U))   Gauss linear;
+
+    div((p*U))              Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         orthogonal;
+}
+
+fluxRequired
+{
+    default         no;
+    p_rgh           ;
+    alpha.gas       ;
+    alpha.liquid    ;
+    Xvapour.gas     ;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/fvSolution b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/fvSolution
new file mode 100755
index 0000000000000000000000000000000000000000..d57847a803c67b02d2cf2cf8d2b00cb63901f48e
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/fvSolution
@@ -0,0 +1,114 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+
+    "alpha.*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.0;
+
+        cAlphas          ((liquid and gas) 1);
+
+        nAlphaCorr      1;
+        nAlphaSubCycles 1;
+
+        // Compressiion factor for species in each alpha phase
+        // NOTE: It should be similar to cAlpha
+        cYi             1;
+    }
+
+    p_rgh
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-09;
+        relTol          0.02;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "U.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    "Yi.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-09;
+        relTol          0;
+        residualAlpha   1e-6;
+    }
+
+    "T.*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+
+    "mDotSmear.*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DIC;
+        tolerance       1e-08;
+        relTol          0.0;
+    }
+
+    "Xvapour.gas.*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-08;
+        relTol          0.0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   no;
+    nOuterCorrectors    1;
+    nCorrectors         2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+    }
+    equations
+    {
+        ".*"     1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/setFieldsDict b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/setFieldsDict
new file mode 100755
index 0000000000000000000000000000000000000000..48a4d62c0d91005b59fc29a57fc13ae322e407f9
--- /dev/null
+++ b/tutorials/multiphase/icoReactingMultiPhaseInterFoam/poolEvaporation/system/setFieldsDict
@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.gas 1
+    volScalarFieldValue alpha.liquid 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (-1 -1 -1) (1 0.012 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha.liquid 1
+            volScalarFieldValue alpha.gas 0
+        );
+    }
+);
+// ************************************************************************* //