ENH: Adding features for phase change solvers

1) Adding interfaceHeight FO
2) Adding interfaceHeatResistance mass transfer model to
   interCondensatingEvaporatingFoam with spread source approach
3) Reworking framework for icoReactingMultiphaseInterFoam

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H

@@ -1,7 +1,7 @@
 {
     radiation->correct();
     rhoCp =  rho*fluid.Cp();
-
+    
     const surfaceScalarField rhoCpPhi(fvc::interpolate(fluid.Cp())*rhoPhi);
 
     const volScalarField kappaEff
@@ -31,7 +31,7 @@
     fvOptions.correct(T);
 
     fluid.correct();
-
+    
     Info<< "min/max(T) = "
         << min(T).value() << ", " << max(T).value() << endl;
 }

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/UEqn.H

@@ -2,7 +2,7 @@ fvVectorMatrix UEqn
 (
     fvm::ddt(rho, U)
   + fvm::div(rhoPhi, U)
-    //- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
+//- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
   + turbulence->divDevRhoReff(U)
   ==
     fvOptions(rho, U)

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.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.
@@ -98,13 +98,13 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        fluid.correctMassSources(T);
         fluid.solve();
         rho = fluid.rho();
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
             #include "UEqn.H"
             #include "YEqns.H"
             #include "TEqn.H"

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C

@@ -50,6 +50,44 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Lee
 template<class Thermo, class OtherThermo>
 Foam::tmp<Foam::volScalarField>
 Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Kexp
+(
+    const volScalarField& refValue
+)
+{
+    {
+        const volScalarField from
+        (
+            min(max(this->pair().from(), scalar(0)), scalar(1))
+        );
+        
+        const volScalarField coeff
+        (
+            C_*from*this->pair().from().rho()*pos(from - alphaMin_)
+           *(refValue - Tactivate_)
+           /Tactivate_
+        );
+        
+        if (sign(C_.value()) > 0)
+        {
+            return
+            (
+                coeff*pos(refValue - Tactivate_)
+            );
+        }
+        else
+        {
+            return
+            (
+               coeff*pos(Tactivate_ - refValue)
+            );
+        }
+    }
+}
+
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSp
 (
     label variable,
     const volScalarField& refValue
@@ -61,31 +99,68 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Kexp
         (
             min(max(this->pair().from(), scalar(0)), scalar(1))
         );
+        
+        const volScalarField coeff
+        (
+            C_*from*this->pair().from().rho()*pos(from - alphaMin_)
+            /Tactivate_
+        );
 
         if (sign(C_.value()) > 0)
         {
             return
             (
-                C_
-              * from
-              * this->pair().from().rho()
-              * (refValue.oldTime() - Tactivate_)
-              * pos(from - alphaMin_)
-              * pos(refValue.oldTime() - Tactivate_)/Tactivate_
+                coeff*pos(refValue - Tactivate_)
             );
         }
         else
         {
             return
             (
-               -C_
-              * from
-              * this->pair().from().rho()
-              * pos(from - alphaMin_)
-              * (Tactivate_ - refValue.oldTime())
-              * pos(Tactivate_ - refValue.oldTime())/Tactivate_
+                coeff*pos(Tactivate_ - refValue)
             );
+        }
+    }
+    else
+    {
+        return tmp<volScalarField> ();
+    }
+}
 
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSu
+(
+    label variable,
+    const volScalarField& refValue
+)
+{
+    if (this->modelVariable_ == variable)
+    {
+        volScalarField from
+        (
+            min(max(this->pair().from(), scalar(0)), scalar(1))
+        );
+        
+        const volScalarField coeff
+        (
+            C_*from*this->pair().from().rho()*pos(from - alphaMin_)
+        );
+
+        if (sign(C_.value()) > 0)
+        {
+            return
+            (
+                -coeff*pos(refValue - Tactivate_)
+            );
+        }
+        else
+        {
+            return
+            (
+                coeff*pos(Tactivate_ - refValue)
+            );
         }
     }
     else
@@ -103,4 +178,12 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Tactivate() const
 }
 
 
+template<class Thermo, class OtherThermo>
+bool 
+Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::includeDivU()
+{
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H

@@ -146,12 +146,29 @@ public:
         //- Explicit mass transfer coefficient
         virtual tmp<volScalarField> Kexp
         (
-            label variable,
+            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 T transition between phases
         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();
 };
 
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C

@@ -90,4 +90,10 @@ const Foam::word Foam::interfaceCompositionModel::variable() const
 }
 
 
+bool Foam::interfaceCompositionModel::includeDivU()
+{
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H

@@ -67,7 +67,8 @@ public:
         {
             T,  /* temperature based */
             P,  /* pressure based */
-            Y   /* mass fraction based */
+            Y,   /* mass fraction based */
+            alpha /* alpha source */
         };
 
         static const Enum<modelVariable> modelVariableNames;
@@ -170,15 +171,33 @@ public:
             const volScalarField& Tf
         ) const = 0;
 
-        //- Explicit mass transfer coefficient
+        //- Explicit full mass transfer
         virtual tmp<volScalarField> Kexp
+        (
+            const volScalarField& field
+        ) = 0;
+        
+        //- Implicit mass transfer
+        virtual tmp<volScalarField> KSp
         (
             label modelVariable,
             const volScalarField& field
         ) = 0;
+        
+        //- Explicit mass transfer
+        virtual tmp<volScalarField> KSu
+        (
+            label modelVariable,
+            const volScalarField& field
+        ) = 0;
+        
 
         //- Reference value
         virtual const dimensionedScalar& Tactivate() const = 0;
+        
+        //- Adds and substract alpha*div(U) as a source term
+        //  for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
+        virtual bool includeDivU();
 
         //- Returns the variable on which the model is based
         const word variable() const;

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C

@@ -88,9 +88,8 @@ Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
 template<class Thermo, class OtherThermo>
 Foam::tmp<Foam::volScalarField>
 Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
-::Kexp(label variable, const volScalarField& field)
+::Kexp(const volScalarField& field)
 {
-    if (this->modelVariable_ == variable)
     {
         const volScalarField& to = this->pair().to();
 
@@ -246,10 +245,29 @@ Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>
 
         return massFluxEvap*areaDensity*Nl*from;
     }
-    else
-    {
-        return tmp<volScalarField> ();
-    }
+}
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>::KSu
+(
+    label variable,
+    const volScalarField& refValue
+)
+{
+    return tmp<volScalarField> ();
+}
+
+
+template<class Thermo, class OtherThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::meltingEvaporationModels::kineticGasEvaporation<Thermo, OtherThermo>::KSp
+(
+    label variable,
+    const volScalarField& refValue
+)
+{
+    return tmp<volScalarField> ();
 }
 
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H

@@ -185,7 +185,20 @@ public:
         //- Explicit mass transfer coefficient
         virtual tmp<volScalarField> Kexp
         (
-            label variable,
+            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
         );
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H

@@ -33,47 +33,9 @@
         (
             fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
+          + fluid.volTransfer(p_rgh)
         );
-
-        forAllConstIters(fluid.totalPhasePairs(), iter)
-        {
-            const phasePair& pair = iter()();
-
-            const phaseModel& phase1 = pair.phase1();
-            const phaseModel& phase2 = pair.phase2();
-
-            const phasePairKey key12
-            (
-                phase1.name(),
-                phase2.name(),
-                true
-            );
-
-            // Mass transfer from phase2 to phase1
-            tmp<volScalarField> tdmdt12(fluid.dmdt(key12));
-            const volScalarField& dmdt12 = tdmdt12();
-
-            const phasePairKey key21
-            (
-                phase2.name(),
-                phase1.name(),
-                true
-            );
-
-            // Mass transfer from phase1 to phase2
-            tmp<volScalarField> tdmdt21(fluid.dmdt(key21));
-            const volScalarField&  dmdt21 = tdmdt21();
-
-            const volScalarField dmdtNet(dmdt21 - dmdt12);
-
-            p_rghEqn +=
-                dmdtNet*
-                (
-                  - fluid.coeffs(phase1.name())
-                  + fluid.coeffs(phase2.name())
-                );
-        }
-
+          
         p_rghEqn.setReference(pRefCell, pRefValue);
 
         p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-202 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,7 @@ License
 #include "fvcDiv.H"
 #include "fvmSup.H"
 #include "fvMatrix.H"
+#include "volFields.H"
 #include "fundamentalConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -152,9 +153,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::dmdt
             (
                 "dmdt",
                 this->mesh().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
+                this->mesh()
             ),
             this->mesh(),
             dimensionedScalar(dimDensity/dimTime, Zero)
@@ -213,15 +212,45 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
                         (
                             "tdmdtYki",
                             this->mesh().time().timeName(),
-                            this->mesh(),
-                            IOobject::NO_READ,
-                            IOobject::NO_WRITE
+                            this->mesh()
                         ),
                         this->mesh(),
                         dimensionedScalar(dimDensity/dimTime, Zero)
                     )
                 );
                 volScalarField& dmdtNetki = tdmdtNetki.ref();
+                
+                tmp<volScalarField> tSp
+                (
+                    new volScalarField
+                    (
+                        IOobject
+                        (
+                            "Sp",
+                            this->mesh().time().timeName(),
+                            this->mesh()
+                        ),
+                        this->mesh(),
+                        dimensionedScalar(dimDensity/dimTime/dimTemperature, Zero)
+                    )
+                );
+                volScalarField& Sp = tSp.ref();
+                
+                tmp<volScalarField> tSu
+                (
+                    new volScalarField
+                    (
+                        IOobject
+                        (
+                            "Su",
+                            this->mesh().time().timeName(),
+                            this->mesh()
+                        ),
+                        this->mesh(),
+                        dimensionedScalar(dimDensity/dimTime, Zero)
+                    )
+                );
+                volScalarField& Su = tSu.ref();
 
 
                 if (massTransferModels_.found(keyik))
@@ -229,14 +258,28 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
                     autoPtr<interfaceCompositionModel>& interfacePtr =
                         massTransferModels_[keyik];
 
-                    // Explicit temperature mass transfer rate
-                    tmp<volScalarField> Kexp =
-                        interfacePtr->Kexp(interfaceCompositionModel::T, T);
-
-                    if (Kexp.valid())
+                    dmdtNetki -= *dmdt_[keyik];
+                    
+                    tmp<volScalarField> KSp =
+                        interfacePtr->KSp(interfaceCompositionModel::T, T);
+                    
+                    if (KSp.valid())
                     {
-                        dmdtNetki -= Kexp.ref();
-                        *dmdt_[keyik] = Kexp.ref();
+                        Sp -= KSp.ref();
+                    }
+                    
+                    tmp<volScalarField> KSu =
+                        interfacePtr->KSu(interfaceCompositionModel::T, T);
+                    
+                    if (KSu.valid())
+                    {
+                        Su -= KSu.ref();
+                    }
+                    
+                    // If linearization is not provided used full explicit
+                    if (!KSp.valid() && !KSu.valid())
+                    {
+                        Su -= *dmdt_[keyik];
                     }
                 }
 
@@ -246,30 +289,36 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
                     autoPtr<interfaceCompositionModel>& interfacePtr =
                         massTransferModels_[keyki];
 
-                    // Explicit temperature mass transfer rate
-                    const tmp<volScalarField> Kexp =
-                        interfacePtr->Kexp(interfaceCompositionModel::T, T);
+                    dmdtNetki += *dmdt_[keyki];
 
-                    if (Kexp.valid())
+                    
+                    tmp<volScalarField> KSp =
+                        interfacePtr->KSp(interfaceCompositionModel::T, T);
+                    
+                    if (KSp.valid())
                     {
-                        dmdtNetki += Kexp.ref();
-                        *dmdt_[keyki] = Kexp.ref();
+                        Sp += KSp.ref();
+                    }
+                    
+                    tmp<volScalarField> KSu =
+                        interfacePtr->KSu(interfaceCompositionModel::T, T);
+                    
+                    if (KSu.valid())
+                    {
+                        Su += KSu.ref();
+                    }
+                    
+                    // If linearization is not provided used full explicit
+                    if (!KSp.valid() && !KSu.valid())
+                    {
+                        Su += *dmdt_[keyki];
                     }
-
                 }
 
-                word keyikName(phasei.name() + phasek.name());
-                word keykiName(phasek.name() + phasei.name());
+                tmp<volScalarField> L = calculateL(dmdtNetki, keyik, keyki, T);
 
-                eqn -=
-                    (
-                        dmdtNetki
-                        *(
-                            calculateL(dmdtNetki, keyik, keyki, T)
-                            - (phasek.Cp() - phasei.Cp())
-                            * constant::standard::Tstd
-                        )
-                    );
+                //eqn -= dmdtNetki*L;
+                eqn -= fvm::Sp(Sp*L.ref(), T) + Su*L.ref();
             }
         }
     }
@@ -277,6 +326,401 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
 }
 
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
+(
+    const volScalarField& p
+)
+{
+    tmp<fvScalarMatrix> tEqnPtr
+    (
+        new fvScalarMatrix(p, dimVolume/dimTime)
+    );
+    
+    fvScalarMatrix& eqn = tEqnPtr.ref();
+    
+    tmp<volScalarField> tSp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Sp",
+                this->mesh().time().timeName(),
+                this->mesh()
+            ),
+            this->mesh(),
+            dimensionedScalar(dimless/dimTime/dimPressure, Zero)
+        )
+    );
+    volScalarField& Sp = tSp.ref();
+    
+    tmp<volScalarField> tSu
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Su",
+                this->mesh().time().timeName(),
+                this->mesh()
+            ),
+            this->mesh(),
+            dimensionedScalar(dimless/dimTime, Zero)
+        )
+    );
+    volScalarField& Su = tSu.ref();
+    
+    forAllConstIters(this->totalPhasePairs(), iter)
+    {
+        const phasePair& pair = iter()();
+
+        const phaseModel& phase1 = pair.phase1();
+        const phaseModel& phase2 = pair.phase2();
+
+        const phasePairKey key12
+        (
+            phase1.name(),
+            phase2.name(),
+            true
+        );
+        
+        if (massTransferModels_.found(key12))
+        {
+            autoPtr<interfaceCompositionModel>& interfacePtr =
+                massTransferModels_[key12];
+           
+            tmp<volScalarField> KSp =
+                interfacePtr->KSp(interfaceCompositionModel::P, p);
+            
+            if (KSp.valid())
+            {
+                Sp += KSp.ref();
+            }
+            
+            tmp<volScalarField> KSu =
+                interfacePtr->KSu(interfaceCompositionModel::P, p);
+            
+            if (KSu.valid())
+            {
+                Su += KSu.ref();
+            }
+            
+            // If linearization is not provided used full explicit
+            if (!KSp.valid() && !KSu.valid())
+            {
+                Su -= 
+                    *dmdt_[key12] 
+                    *(
+                        - this->coeffs(phase1.name())
+                        + this->coeffs(phase2.name())
+                    );
+            }
+        }
+        
+        const phasePairKey key21
+        (
+            phase2.name(),
+            phase1.name(),
+            true
+        );
+        
+        if (massTransferModels_.found(key21))
+        {
+            autoPtr<interfaceCompositionModel>& interfacePtr =
+                massTransferModels_[key21];
+           
+            tmp<volScalarField> KSp =
+                interfacePtr->KSp(interfaceCompositionModel::P, p);
+            
+            if (KSp.valid())
+            {
+                Sp += KSp.ref();
+            }
+            
+            tmp<volScalarField> KSu =
+                interfacePtr->KSu(interfaceCompositionModel::P, p);
+            
+            if (KSu.valid())
+            {
+                Su += KSu.ref();
+            }
+            
+            // If linearization is not provided used full explicit
+            if (!KSp.valid() && !KSu.valid())
+            {
+                Su += 
+                    *dmdt_[key21]
+                    *(
+                        - this->coeffs(phase1.name())
+                        + this->coeffs(phase2.name())
+                    );
+            }
+        }
+
+    }
+    
+    eqn += fvm::Sp(Sp, p) + Su;
+    return tEqnPtr;
+}
+
+
+template<class BasePhaseSystem>
+void Foam::MassTransferPhaseSystem<BasePhaseSystem>::correctMassSources
+(
+    const volScalarField& T
+)
+{
+    forAllConstIters(this->phaseModels_, iteri)
+    {
+        const phaseModel& phasei = iteri()();
+
+        auto iterk = iteri;
+
+        for (++iterk; iterk != this->phaseModels_.end(); ++iterk)
+        {
+            if (iteri()().name() != iterk()().name())
+            {
+                const phaseModel& phasek = iterk()();
+
+                // Phase i to phase k
+                const phasePairKey keyik(phasei.name(), phasek.name(), true);
+
+                // Phase k to phase i
+                const phasePairKey keyki(phasek.name(), phasei.name(), true);
+                
+                if (massTransferModels_.found(keyik))
+                {
+                    autoPtr<interfaceCompositionModel>& interfacePtr =
+                        massTransferModels_[keyik];
+                        
+                    tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
+                    
+                    *dmdt_[keyik] = Kexp.ref();
+                    
+                }
+                
+                if (massTransferModels_.found(keyki))
+                {
+                    autoPtr<interfaceCompositionModel>& interfacePtr =
+                        massTransferModels_[keyki];
+
+                    // Explicit temperature mass transfer rate
+                    const tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
+                    
+                    *dmdt_[keyki] = Kexp.ref();
+                }
+            }
+        }
+    }
+}
+
+
+template<class BasePhaseSystem>
+void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
+(
+    SuSpTable& Su,
+    SuSpTable& Sp
+)
+{
+    // This term adds and substract alpha*div(U) as a source term
+    // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
+    bool includeDivU(true);
+    
+    forAllConstIters(this->totalPhasePairs(), iter)
+    {
+        const phasePair& pair = iter()();
+
+        const phaseModel& phase1 = pair.phase1();
+        const phaseModel& phase2 = pair.phase2();
+
+        const volScalarField& alpha1 = pair.phase1();
+        const volScalarField& alpha2 = pair.phase2();
+
+        tmp<volScalarField> tCoeffs1 = this->coeffs(phase1.name());
+        const volScalarField&  coeffs1 = tCoeffs1();
+
+        tmp<volScalarField> tCoeffs2 = this->coeffs(phase2.name());
+        const volScalarField&  coeffs2 = tCoeffs2();
+
+        // Phase 1 to phase 2
+        const phasePairKey key12
+        (
+            phase1.name(),
+            phase2.name(),
+            true
+        );
+        
+        tmp<volScalarField> tdmdt12(this->dmdt(key12));
+        volScalarField& dmdt12 = tdmdt12.ref();
+        
+        if (massTransferModels_.found(key12))
+        {
+            autoPtr<interfaceCompositionModel>& interfacePtr =
+                massTransferModels_[key12];
+             
+            tmp<volScalarField> KSu =
+                interfacePtr->KSu(interfaceCompositionModel::alpha, phase1);
+            
+            if (KSu.valid())
+            {
+                dmdt12 = KSu.ref();
+            }
+            
+            includeDivU = interfacePtr->includeDivU();
+        }
+       
+
+        // Phase 2 to phase 1
+        const phasePairKey key21
+        (
+            phase2.name(),
+            phase1.name(),
+            true
+        );
+
+        tmp<volScalarField> tdmdt21(this->dmdt(key21));
+        volScalarField& dmdt21 = tdmdt21.ref();
+        
+        if (massTransferModels_.found(key21))
+        {
+            autoPtr<interfaceCompositionModel>& interfacePtr =
+                massTransferModels_[key21];
+            
+            tmp<volScalarField> KSu =
+                interfacePtr->KSu(interfaceCompositionModel::alpha, phase2);
+            
+            if (KSu.valid())
+            {
+                dmdt21 = KSu.ref();
+            }
+            
+            includeDivU = interfacePtr->includeDivU();
+        }
+
+        volScalarField::Internal& SpPhase1 = Sp[phase1.name()];
+
+        volScalarField::Internal& SuPhase1 = Su[phase1.name()];
+
+        volScalarField::Internal& SpPhase2 = Sp[phase2.name()];
+
+        volScalarField::Internal& SuPhase2 = Su[phase2.name()];
+
+        const volScalarField dmdtNet(dmdt21 - dmdt12);
+
+        const volScalarField coeffs12(coeffs1 - coeffs2);
+
+        const surfaceScalarField& phi = this->phi();
+        
+        if (includeDivU)
+        {
+            SuPhase1 +=
+                fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1));
+                
+            SuPhase2 +=
+                fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1));
+        }
+        
+        // NOTE: dmdtNet is distributed in terms =
+        //  Source for phase 1 =
+        //      dmdtNet/rho1
+        //    - alpha1*dmdtNet(1/rho1 - 1/rho2)
+
+        forAll(dmdtNet, celli)
+        {
+            scalar dmdt21 = dmdtNet[celli];
+            scalar coeffs12Cell = coeffs12[celli];
+
+            scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0);
+
+            // exp.
+            SuPhase1[celli] += coeffs1[celli]*dmdt21;
+
+            if (includeDivU)
+            {
+                if (dmdt21 > 0)
+                {
+                    if (coeffs12Cell > 0)
+                    {
+                        // imp
+                        SpPhase1[celli] -= dmdt21*coeffs12Cell;
+                    }
+                    else if (coeffs12Cell < 0)
+                    {
+                        // exp
+                        SuPhase1[celli] -=
+                            dmdt21*coeffs12Cell*alpha1Limited;
+                    }
+                }
+                else if (dmdt21 < 0)
+                {
+                    if (coeffs12Cell > 0)
+                    {
+                        // exp
+                        SuPhase1[celli] -=
+                            dmdt21*coeffs12Cell*alpha1Limited;
+                    }
+                    else if (coeffs12Cell < 0)
+                    {
+                        // imp
+                        SpPhase1[celli] -= dmdt21*coeffs12Cell;
+                    }
+                }
+            }
+        }
+
+        forAll(dmdtNet, celli)
+        {
+            scalar dmdt12 = -dmdtNet[celli];
+            scalar coeffs21Cell = -coeffs12[celli];
+
+            scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0);
+
+            // exp
+            SuPhase2[celli] += coeffs2[celli]*dmdt12;
+
+            if (includeDivU)
+            {
+                if (dmdt12 > 0)
+                {
+                    if (coeffs21Cell > 0)
+                    {
+                        // imp
+                        SpPhase2[celli] -= dmdt12*coeffs21Cell;
+                    }
+                    else if (coeffs21Cell < 0)
+                    {
+                        // exp
+                        SuPhase2[celli] -=
+                            dmdt12*coeffs21Cell*alpha2Limited;
+                    }
+                }
+                else if (dmdt12 < 0)
+                {
+                    if (coeffs21Cell > 0)
+                    {
+                        // exp
+                        SuPhase2[celli] -=
+                            coeffs21Cell*dmdt12*alpha2Limited;
+                    }
+                    else if (coeffs21Cell < 0)
+                    {
+                        // imp
+                        SpPhase2[celli] -= dmdt12*coeffs21Cell;
+                    }
+                }
+            }
+
+        }
+
+        // Update ddtAlphaMax
+        this->ddtAlphaMax_ =  
+            max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
+    }
+}
+
+        
 template<class BasePhaseSystem>
 void Foam::MassTransferPhaseSystem<BasePhaseSystem>::massSpeciesTransfer
 (

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.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.
@@ -59,7 +59,7 @@ public:
 
     // Public typedefs
 
-     typedef
+    typedef
             HashTable
             <
                 autoPtr<interfaceCompositionModel>,
@@ -67,6 +67,9 @@ public:
                 phasePairKey::hash
             >
             massTransferModelTable;
+            
+            
+    typedef HashTable<volScalarField::Internal> SuSpTable;
 
 protected:
 
@@ -123,9 +126,17 @@ public:
 
     // Mass transfer functions
 
-        //- Return the heat transfer matrix and fill dmdt for phases
+        //- Return the heat transfer matrix
         virtual  tmp<fvScalarMatrix> heatTransfer(const volScalarField& T);
-
+        
+        //- Return the volumetric rate transfer matrix
+        virtual  tmp<fvScalarMatrix> volTransfer(const volScalarField& p);
+        
+        //- Correct/calculates mass sources dmdt for phases
+        virtual void correctMassSources(const volScalarField& T);
+        
+        //- Calculate mass transfer for alpha's
+        virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp);
 
         //- Calculate mass transfer for species
         virtual void massSpeciesTransfer

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C

@@ -90,9 +90,7 @@ MultiComponentPhaseModel
                 (
                     IOobject::groupName("X" + species_[i], phaseName),
                     fluid.mesh().time().timeName(),
-                    fluid.mesh(),
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
+                    fluid.mesh()
                 ),
                 Y()[i]
             )

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -120,166 +120,60 @@ Foam::multiphaseSystem::multiphaseSystem
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::multiphaseSystem::calculateSuSp()
-{
-    forAllConstIters(totalPhasePairs_, iter)
-    {
-        const phasePair& pair = iter()();
-
-        const phaseModel& phase1 = pair.phase1();
-        const phaseModel& phase2 = pair.phase2();
-
-        const volScalarField& alpha1 = pair.phase1();
-        const volScalarField& alpha2 = pair.phase2();
-
-        tmp<volScalarField> tCoeffs1 = this->coeffs(phase1.name());
-        const volScalarField&  coeffs1 = tCoeffs1();
-
-        tmp<volScalarField> tCoeffs2 = this->coeffs(phase2.name());
-        const volScalarField&  coeffs2 = tCoeffs2();
-
-        // Phase 1 to phase 2
-        const phasePairKey key12
-        (
-            phase1.name(),
-            phase2.name(),
-            true
-        );
-
-
-        tmp<volScalarField> tdmdt12(this->dmdt(key12));
-        const volScalarField& dmdt12 = tdmdt12();
-
-        // Phase 2 to phase 1
-        const phasePairKey key21
-        (
-            phase2.name(),
-            phase1.name(),
-            true
-        );
-
-        tmp<volScalarField> tdmdt21(this->dmdt(key21));
-        const volScalarField& dmdt21 = tdmdt21();
-
-        volScalarField::Internal& SpPhase1 = Sp_[phase1.name()];
-
-        volScalarField::Internal& SuPhase1 = Su_[phase1.name()];
-
-        volScalarField::Internal& SpPhase2 = Sp_[phase2.name()];
-
-        volScalarField::Internal& SuPhase2 = Su_[phase2.name()];
-
-        const volScalarField dmdtNet(dmdt21 - dmdt12);
-
-        const volScalarField coeffs12(coeffs1 - coeffs2);
-
-        // NOTE: dmdtNet is distributed in terms =
-        //  Source for phase 1 =
-        //      dmdtNet/rho1
-        //    - alpha1*dmdtNet(1/rho1 - 1/rho2)
-
-        forAll(dmdtNet, celli)
-        {
-            scalar dmdt21 = dmdtNet[celli];
-            scalar coeffs12Cell = coeffs12[celli];
-
-            scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0);
-
-            // exp.
-            SuPhase1[celli] += coeffs1[celli]*dmdt21;
-
-            if (dmdt21 > 0)
-            {
-                if (coeffs12Cell > 0)
-                {
-                    // imp
-                    SpPhase1[celli] -= dmdt21*coeffs12Cell;
-                }
-                else if (coeffs12Cell < 0)
-                {
-                    // exp
-                    SuPhase1[celli] -=
-                        dmdt21*coeffs12Cell*alpha1Limited;
-                }
-            }
-            else if (dmdt21 < 0)
-            {
-                if (coeffs12Cell > 0)
-                {
-                    // exp
-                    SuPhase1[celli] -=
-                        dmdt21*coeffs12Cell*alpha1Limited;
-                }
-                else if (coeffs12Cell < 0)
-                {
-                    // imp
-                    SpPhase1[celli] -= dmdt21*coeffs12Cell;
-                }
-            }
-        }
-
-        forAll(dmdtNet, celli)
-        {
-            scalar dmdt12 = -dmdtNet[celli];
-            scalar coeffs21Cell = -coeffs12[celli];
-
-            scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0);
-
-            // exp
-            SuPhase2[celli] += coeffs2[celli]*dmdt12;
-
-            if (dmdt12 > 0)
-            {
-                if (coeffs21Cell > 0)
-                {
-                    // imp
-                    SpPhase2[celli] -= dmdt12*coeffs21Cell;
-                }
-                else if (coeffs21Cell < 0)
-                {
-                    // exp
-                    SuPhase2[celli] -=
-                        dmdt12*coeffs21Cell*alpha2Limited;
-                }
-            }
-            else if (dmdt12 < 0)
-            {
-                if (coeffs21Cell > 0)
-                {
-                    // exp
-                    SuPhase2[celli] -=
-                        coeffs21Cell*dmdt12*alpha2Limited;
-                }
-                else if (coeffs21Cell < 0)
-                {
-                    // imp
-                    SpPhase2[celli] -= dmdt12*coeffs21Cell;
-                }
-            }
-        }
-
-        // Update ddtAlphaMax
-        ddtAlphaMax_ =
-            max
-            (
-                ddtAlphaMax_.value(),
-                max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()))
-            );
-    }
+{   
+    this->alphaTransfer(Su_, Sp_);
 }
 
 
 void Foam::multiphaseSystem::solve()
 {
-    const fvMesh& mesh = this->mesh();
-
-    const dictionary& alphaControls = mesh.solverDict("alpha");
+    const dictionary& alphaControls = mesh_.solverDict("alpha");
     label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
+    
+    volScalarField& alpha = phases_.first();
+    
+    if (nAlphaSubCycles > 1)
+    {
+        surfaceScalarField rhoPhiSum
+        (
+            IOobject
+            (
+                "rhoPhiSum",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar(rhoPhi_.dimensions(), Zero)
+        );
+
+        dimensionedScalar totalDeltaT = mesh_.time().deltaT();
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            solveAlphas();
+            rhoPhiSum += (mesh_.time().deltaT()/totalDeltaT)*rhoPhi_;
+        }
+
+        rhoPhi_ = rhoPhiSum;
+    }
+    else
+    {
+        solveAlphas();
+    }
+    
+}
+
+void Foam::multiphaseSystem::solveAlphas()
+{
+    
+    mesh_.solverDict("alpha").readEntry("cAlphas", cAlphas_);
+    const dictionary& alphaControls = mesh_.solverDict("alpha");
     label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
-    mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
-
-    // Reset ddtAlphaMax
-    ddtAlphaMax_ = dimensionedScalar(dimless, Zero);
-
+    
     PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
 
     const surfaceScalarField& phi = this->phi();
@@ -389,12 +283,11 @@ void Foam::multiphaseSystem::solve()
             Sp_[phase.name()] = dimensionedScalar("Sp", dimless/dimTime, Zero);
 
             // Add alpha*div(U)
-            const volScalarField& alpha = phase;
-            Su_[phase.name()] +=
-                fvc::div(phi)*min(max(alpha, scalar(0)), scalar(1));
+            //const volScalarField& alpha = phase;
+            //Su_[phase.name()] +=
+            //    fvc::div(phi)*min(max(alpha, scalar(0)), scalar(1));
         }
 
-
         // Fill Su and Sp
         calculateSuSp();
 
@@ -454,12 +347,12 @@ void Foam::multiphaseSystem::solve()
             //phiAlpha += upwind<scalar>(mesh_, phi).flux(alpha1);
             fvScalarMatrix alpha1Eqn
             (
-                fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+                fv::EulerDdtScheme<scalar>(mesh_).fvmDdt(alpha1)
               + fv::gaussConvectionScheme<scalar>
                 (
-                    mesh,
+                    mesh_,
                     phi,
-                    upwind<scalar>(mesh, phi)
+                    upwind<scalar>(mesh_, phi)
                 ).fvmDiv(phi, alpha1)
               ==
                  Su + fvm::Sp(Sp, alpha1)
@@ -468,60 +361,21 @@ void Foam::multiphaseSystem::solve()
             alpha1Eqn.solve();
 
             phiAlpha += alpha1Eqn.flux();
+            
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                alpha1,
+                phi,
+                phiAlpha,
+                Sp,
+                Su,
+                oneField(),
+                zeroField()
+            );
 
-            if (nAlphaSubCycles > 1)
-            {
-                for
-                (
-                    subCycle<volScalarField> alphaSubCycle
-                    (
-                        alpha1,
-                        nAlphaSubCycles
-                    );
-                    !(++alphaSubCycle).end();
-                )
-                {
-                    MULES::explicitSolve
-                    (
-                        geometricOneField(),
-                        alpha1,
-                        phi,
-                        phiAlpha,
-                        (alphaSubCycle.index()*Sp)(),
-                        (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
-                        oneField(),
-                        zeroField()
-                    );
-
-                    if (alphaSubCycle.index() == 1)
-                    {
-                        phase.alphaPhi() = phiAlpha;
-                    }
-                    else
-                    {
-                        phase.alphaPhi() += phiAlpha;
-                    }
-                }
-
-                phase.alphaPhi() /= nAlphaSubCycles;
-            }
-            else
-            {
-                MULES::explicitSolve
-                (
-                    geometricOneField(),
-                    alpha1,
-                    phi,
-                    phiAlpha,
-                    Sp,
-                    Su,
-                    oneField(),
-                    zeroField()
-                );
-
-                phase.alphaPhi() = phiAlpha;
-            }
-
+            phase.alphaPhi() = phiAlpha;
+            
             ++phasei;
         }
 
@@ -549,7 +403,6 @@ void Foam::multiphaseSystem::solve()
 
                 // Update rhoPhi
                 rhoPhi_ += fvc::interpolate(phase.rho()) * phase.alphaPhi();
-
             }
 
             Info<< "Phase-sum volume fraction, min, max = "
@@ -566,7 +419,7 @@ void Foam::multiphaseSystem::solve()
                 alpha += alpha*sumCorr;
 
                 Info<< alpha.name() << " volume fraction = "
-                    << alpha.weightedAverage(mesh.V()).value()
+                    << alpha.weightedAverage(mesh_.V()).value()
                     << "  Min(alpha) = " << min(alpha).value()
                     << "  Max(alpha) = " << max(alpha).value()
                     << endl;

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.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.
@@ -89,6 +89,9 @@ protected:
 
         //- Calculate Sp and Su
         void calculateSuSp();
+        
+        //- Solve alphas
+        void solveAlphas();
 
 
 public:

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem_backUp.C

@@ -0,0 +1,642 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017-2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "multiphaseSystem.H"
+
+#include "fixedValueFvsPatchFields.H"
+#include "Time.H"
+#include "subCycle.H"
+#include "fvcMeshPhi.H"
+
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
+#include "fvcSnGrad.H"
+#include "fvcDiv.H"
+#include "fvcDdt.H"
+#include "fvcFlux.H"
+#include "fvmDdt.H"
+#include "fvcAverage.H"
+#include "fvMatrix.H"
+#include "fvmSup.H"
+#include "CMULES.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(multiphaseSystem, 0);
+    defineRunTimeSelectionTable(multiphaseSystem, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseSystem::multiphaseSystem
+(
+     const fvMesh& mesh
+)
+:
+    phaseSystem(mesh),
+    cAlphas_(),
+    ddtAlphaMax_(0.0),
+    limitedPhiAlphas_(phaseModels_.size()),
+    Su_(phaseModels_.size()),
+    Sp_(phaseModels_.size())
+{
+    label phasei = 0;
+    phases_.setSize(phaseModels_.size());
+    forAllIters(phaseModels_, iter)
+    {
+        phaseModel& pm = iter()();
+        phases_.set(phasei++, &pm);
+    }
+
+    mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
+
+    // Initiate Su and Sp
+    forAllConstIters(phaseModels_, iter)
+    {
+        const phaseModel& pm = iter()();
+
+        Su_.insert
+        (
+            pm.name(),
+            volScalarField::Internal
+            (
+                IOobject
+                (
+                    "Su" + pm.name(),
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                dimensionedScalar(dimless/dimTime, Zero)
+            )
+        );
+
+        Sp_.insert
+        (
+            pm.name(),
+            volScalarField::Internal
+            (
+                IOobject
+                (
+                    "Sp" + pm.name(),
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                dimensionedScalar(dimless/dimTime, Zero)
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::multiphaseSystem::calculateSuSp()
+{
+    forAllConstIters(totalPhasePairs_, iter)
+    {
+        const phasePair& pair = iter()();
+
+        const phaseModel& phase1 = pair.phase1();
+        const phaseModel& phase2 = pair.phase2();
+
+        const volScalarField& alpha1 = pair.phase1();
+        const volScalarField& alpha2 = pair.phase2();
+
+        tmp<volScalarField> tCoeffs1 = this->coeffs(phase1.name());
+        const volScalarField&  coeffs1 = tCoeffs1();
+
+        tmp<volScalarField> tCoeffs2 = this->coeffs(phase2.name());
+        const volScalarField&  coeffs2 = tCoeffs2();
+
+        // Phase 1 to phase 2
+        const phasePairKey key12
+        (
+            phase1.name(),
+            phase2.name(),
+            true
+        );
+
+
+        tmp<volScalarField> tdmdt12(this->dmdt(key12));
+        const volScalarField& dmdt12 = tdmdt12();
+
+        // Phase 2 to phase 1
+        const phasePairKey key21
+        (
+            phase2.name(),
+            phase1.name(),
+            true
+        );
+
+        tmp<volScalarField> tdmdt21(this->dmdt(key21));
+        const volScalarField& dmdt21 = tdmdt21();
+
+        volScalarField::Internal& SpPhase1 = Sp_[phase1.name()];
+
+        volScalarField::Internal& SuPhase1 = Su_[phase1.name()];
+
+        volScalarField::Internal& SpPhase2 = Sp_[phase2.name()];
+
+        volScalarField::Internal& SuPhase2 = Su_[phase2.name()];
+
+        const volScalarField dmdtNet(dmdt21 - dmdt12);
+
+        const volScalarField coeffs12(coeffs1 - coeffs2);
+
+        // NOTE: dmdtNet is distributed in terms =
+        //  Source for phase 1 =
+        //      dmdtNet/rho1
+        //    - alpha1*dmdtNet(1/rho1 - 1/rho2)
+
+        forAll(dmdtNet, celli)
+        {
+            scalar dmdt21 = dmdtNet[celli];
+            scalar coeffs12Cell = coeffs12[celli];
+
+            scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0);
+
+            // exp.
+            SuPhase1[celli] += coeffs1[celli]*dmdt21;
+
+            if (dmdt21 > 0)
+            {
+                if (coeffs12Cell > 0)
+                {
+                    // imp
+                    SpPhase1[celli] -= dmdt21*coeffs12Cell;
+                }
+                else if (coeffs12Cell < 0)
+                {
+                    // exp
+                    SuPhase1[celli] -=
+                        dmdt21*coeffs12Cell*alpha1Limited;
+                }
+            }
+            else if (dmdt21 < 0)
+            {
+                if (coeffs12Cell > 0)
+                {
+                    // exp
+                    SuPhase1[celli] -=
+                        dmdt21*coeffs12Cell*alpha1Limited;
+                }
+                else if (coeffs12Cell < 0)
+                {
+                    // imp
+                    SpPhase1[celli] -= dmdt21*coeffs12Cell;
+                }
+            }
+        }
+
+        forAll(dmdtNet, celli)
+        {
+            scalar dmdt12 = -dmdtNet[celli];
+            scalar coeffs21Cell = -coeffs12[celli];
+
+            scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0);
+
+            // exp
+            SuPhase2[celli] += coeffs2[celli]*dmdt12;
+
+            if (dmdt12 > 0)
+            {
+                if (coeffs21Cell > 0)
+                {
+                    // imp
+                    SpPhase2[celli] -= dmdt12*coeffs21Cell;
+                }
+                else if (coeffs21Cell < 0)
+                {
+                    // exp
+                    SuPhase2[celli] -=
+                        dmdt12*coeffs21Cell*alpha2Limited;
+                }
+            }
+            else if (dmdt12 < 0)
+            {
+                if (coeffs21Cell > 0)
+                {
+                    // exp
+                    SuPhase2[celli] -=
+                        coeffs21Cell*dmdt12*alpha2Limited;
+                }
+                else if (coeffs21Cell < 0)
+                {
+                    // imp
+                    SpPhase2[celli] -= dmdt12*coeffs21Cell;
+                }
+            }
+
+        }
+
+        // Update ddtAlphaMax
+        ddtAlphaMax_ =  
+            max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
+    }
+}
+
+
+void Foam::multiphaseSystem::solve()
+{
+    const fvMesh& mesh = this->mesh();
+
+    const dictionary& alphaControls = mesh.solverDict("alpha");
+    label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
+    label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
+    mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
+
+    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
+
+    const surfaceScalarField& phi = this->phi();
+
+    surfaceScalarField phic(mag((phi)/mesh_.magSf()));
+
+    // Do not compress interface at non-coupled boundary faces
+    // (inlets, outlets etc.)
+    surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef();
+    forAll(phic.boundaryField(), patchi)
+    {
+        fvsPatchScalarField& phicp = phicBf[patchi];
+
+        if (!phicp.coupled())
+        {
+            phicp == 0;
+        }
+    }
+
+    for (int acorr=0; acorr<nAlphaCorr; acorr++)
+    {
+        label phasei = 0;
+        for (phaseModel& phase1 : phases_)
+        {
+            const volScalarField& alpha1 = phase1;
+
+            phiAlphaCorrs.set
+            (
+                phasei,
+                new surfaceScalarField
+                (
+                    "phi" + alpha1.name() + "Corr",
+                    fvc::flux
+                    (
+                        phi,
+                        alpha1,
+                        "div(phi," + alpha1.name() + ')'
+                    )
+                )
+            );
+
+            surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
+
+            for (phaseModel& phase2 : phases_)
+            {
+                const volScalarField& alpha2 = phase2;
+
+                if (&phase2 == &phase1) continue;
+
+                const phasePairKey key12(phase1.name(), phase2.name());
+
+                if (!cAlphas_.found(key12))
+                {
+                    FatalErrorInFunction
+                        << "Phase compression factor (cAlpha) not found for : "
+                        << key12
+                        << exit(FatalError);
+                }
+                scalar cAlpha = cAlphas_.find(key12)();
+
+                phic = min(cAlpha*phic, max(phic));
+
+                surfaceScalarField phir(phic*nHatf(alpha1, alpha2));
+
+                word phirScheme
+                (
+                    "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
+                );
+
+                phiAlphaCorr += fvc::flux
+                (
+                   -fvc::flux(-phir, alpha2, phirScheme),
+                    alpha1,
+                    phirScheme
+                );
+            }
+
+            // Ensure that the flux at inflow BCs is preserved
+            forAll(phiAlphaCorr.boundaryField(), patchi)
+            {
+                fvsPatchScalarField& phiAlphaCorrp =
+                    phiAlphaCorr.boundaryFieldRef()[patchi];
+
+                if (!phiAlphaCorrp.coupled())
+                {
+                    const scalarField& phi1p = phi.boundaryField()[patchi];
+                    const scalarField& alpha1p =
+                        alpha1.boundaryField()[patchi];
+
+                    forAll(phiAlphaCorrp, facei)
+                    {
+                        if (phi1p[facei] < 0)
+                        {
+                            phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
+                        }
+                    }
+                }
+            }
+
+            ++phasei;
+        }
+
+        // Set Su and Sp to zero
+        for (const phaseModel& phase : phases_)
+        {
+            Su_[phase.name()] = dimensionedScalar("Su", dimless/dimTime, Zero);
+            Sp_[phase.name()] = dimensionedScalar("Sp", dimless/dimTime, Zero);
+
+            // Add alpha*div(U)
+            const volScalarField& alpha = phase;
+            Sp_[phase.name()] +=
+                fvc::div(phi);//*min(max(alpha, scalar(0)), scalar(1));
+        }
+
+
+        // Fill Su and Sp
+        calculateSuSp();
+
+        // Limit phiAlphaCorr on each phase
+        phasei = 0;
+        for (phaseModel& phase : phases_)
+        {
+            volScalarField& alpha1 = phase;
+
+            surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
+
+            volScalarField::Internal& Su = Su_[phase.name()];
+            volScalarField::Internal& Sp = Sp_[phase.name()];
+
+            MULES::limit
+            (
+                1.0/mesh_.time().deltaT().value(),
+                geometricOneField(),
+                alpha1,
+                phi,
+                phiAlphaCorr,
+                Sp,
+                Su,
+                oneField(),
+                zeroField(),
+                true
+            );
+            ++phasei;
+        }
+
+        MULES::limitSum(phiAlphaCorrs);
+
+        volScalarField sumAlpha
+        (
+            IOobject
+            (
+                "sumAlpha",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar(dimless, Zero)
+        );
+
+        phasei = 0;
+        for (phaseModel& phase : phases_)
+        {
+            volScalarField& alpha1 = phase;
+
+            const volScalarField::Internal& Su = Su_[phase.name()];
+
+            const volScalarField::Internal& Sp = Sp_[phase.name()];
+
+            surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
+
+            // Add a bounded upwind U-mean flux
+            phiAlpha += upwind<scalar>(mesh_, phi).flux(alpha1);
+//             fvScalarMatrix alpha1Eqn
+//             (
+//                 fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+//               + fv::gaussConvectionScheme<scalar>
+//                 (
+//                     mesh,
+//                     phi,
+//                     upwind<scalar>(mesh, phi)
+//                 ).fvmDiv(phi, alpha1)
+//               ==
+//                  Su + fvm::Sp(Sp, alpha1)
+//             );
+// 
+//             alpha1Eqn.solve();
+
+            //phiAlpha += alpha1Eqn.flux();
+            
+            if (nAlphaSubCycles > 1)
+            {
+                for
+                (
+                    subCycle<volScalarField> alphaSubCycle
+                    (
+                        alpha1,
+                        nAlphaSubCycles
+                    );
+                    !(++alphaSubCycle).end();
+                )
+                {
+                    MULES::explicitSolve
+                    (
+                        geometricOneField(),
+                        alpha1,
+                        phi,
+                        phiAlpha,
+                        (alphaSubCycle.index()*Sp)(),
+                        (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
+                        oneField(),
+                        zeroField()
+                    );
+
+                    if (alphaSubCycle.index() == 1)
+                    {
+                        phase.alphaPhi() = phiAlpha;
+                    }
+                    else
+                    {
+                        phase.alphaPhi() += phiAlpha;
+                    }
+                }
+
+                phase.alphaPhi() /= nAlphaSubCycles;
+            }
+            else
+            {
+                MULES::explicitSolve
+                (
+                    geometricOneField(),
+                    alpha1,
+                    phi,
+                    phiAlpha,
+                    Sp,
+                    Su,
+                    oneField(),
+                    zeroField()
+                );
+
+                phase.alphaPhi() = phiAlpha;
+            }
+
+            ++phasei;
+        }
+
+        if (acorr == nAlphaCorr - 1)
+        {
+            volScalarField sumAlpha
+            (
+                IOobject
+                (
+                    "sumAlpha",
+                    mesh_.time().timeName(),
+                    mesh_
+                ),
+                mesh_,
+                dimensionedScalar(dimless, Zero)
+            );
+
+            // Reset rhoPhi
+            rhoPhi_ = dimensionedScalar("rhoPhi", dimMass/dimTime, Zero);
+
+            for (phaseModel& phase : phases_)
+            {
+                volScalarField& alpha1 = phase;
+                sumAlpha += alpha1;
+
+                // Update rhoPhi
+                rhoPhi_ += fvc::interpolate(phase.rho()) * phase.alphaPhi();
+            }
+
+            Info<< "Phase-sum volume fraction, min, max = "
+                << sumAlpha.weightedAverage(mesh_.V()).value()
+                << ' ' << min(sumAlpha).value()
+                << ' ' << max(sumAlpha).value()
+                << endl;
+
+            volScalarField sumCorr(1.0 - sumAlpha);
+
+            for (phaseModel& phase : phases_)
+            {
+                volScalarField& alpha = phase;
+                //alpha += alpha*sumCorr;
+
+                Info<< alpha.name() << " volume fraction = "
+                    << alpha.weightedAverage(mesh.V()).value()
+                    << "  Min(alpha) = " << min(alpha).value()
+                    << "  Max(alpha) = " << max(alpha).value()
+                    << endl;
+            }
+        }
+    }
+}
+
+
+const Foam::UPtrList<Foam::phaseModel>& Foam::multiphaseSystem::phases() const
+{
+    return phases_;
+}
+
+
+Foam::UPtrList<Foam::phaseModel>& Foam::multiphaseSystem::phases()
+{
+    return phases_;
+}
+
+
+const Foam::phaseModel& Foam::multiphaseSystem::phase(const label i) const
+{
+    return phases_[i];
+}
+
+
+Foam::phaseModel& Foam::multiphaseSystem::phase(const label i)
+{
+    return phases_[i];
+}
+
+
+Foam::dimensionedScalar Foam::multiphaseSystem::ddtAlphaMax() const
+{
+    return ddtAlphaMax_;
+}
+
+
+Foam::scalar Foam::multiphaseSystem::maxDiffNo() const
+{
+    auto iter = phaseModels_.cbegin();
+
+    scalar maxVal = max(iter()->diffNo()).value();
+
+    for (++iter; iter != phaseModels_.cend(); ++iter)
+    {
+        maxVal = max(maxVal, max(iter()->diffNo()).value());
+    }
+
+    return maxVal * mesh_.time().deltaT().value();
+}
+
+
+const Foam::multiphaseSystem::compressionFluxTable&
+Foam::multiphaseSystem::limitedPhiAlphas() const
+{
+    return limitedPhiAlphas_;
+}
+
+
+Foam::multiphaseSystem::SuSpTable& Foam::multiphaseSystem::Su()
+{
+    return Su_;
+}
+
+
+Foam::multiphaseSystem::SuSpTable& Foam::multiphaseSystem::Sp()
+{
+    return Sp_;
+}
+
+
+bool Foam::multiphaseSystem::read()
+{
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,6 +55,12 @@ const Foam::word Foam::phaseSystem::phasePropertiesName("phaseProperties");
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+void Foam::phaseSystem::calcMu()
+{
+    mu_ = mu()();
+}
+
+
 Foam::phaseSystem::phaseModelTable
 Foam::phaseSystem::generatePhaseModels(const wordList& phaseNames) const
 {
@@ -202,6 +208,18 @@ Foam::phaseSystem::phaseSystem
 :
     basicThermo(mesh, word::null, phasePropertiesName),
     mesh_(mesh),
+    mu_
+    (
+        IOobject
+        (
+            "mu",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar(dimViscosity*dimDensity, Zero),
+        calculatedFvPatchScalarField::typeName
+    ),
     phaseNames_(lookup("phases")),
     phi_
     (
@@ -222,9 +240,7 @@ Foam::phaseSystem::phaseSystem
         (
             "rhoPhi",
             mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
+            mesh_
         ),
         mesh_,
         dimensionedScalar(dimMass/dimTime, Zero)
@@ -265,6 +281,9 @@ Foam::phaseSystem::phaseSystem
 
     // Total phase pair
     generatePairsTable();
+    
+    // Update mu_
+    calcMu();
 }
 
 
@@ -884,6 +903,8 @@ void Foam::phaseSystem::correct()
     {
         iter()->correct();
     }
+    
+    calcMu();
 }
 
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.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.
@@ -78,6 +78,9 @@ public:
 
 
         typedef HashTable<autoPtr<phaseModel>> phaseModelTable;
+        
+
+        typedef HashTable<volScalarField::Internal> SuSpTable;
 
 
 protected:
@@ -113,6 +116,9 @@ protected:
 
         //- Reference to the mesh
         const fvMesh& mesh_;
+        
+        //- Dynamic viscocity
+        volScalarField mu_;
 
         //- Phase names
         wordList phaseNames_;
@@ -148,6 +154,9 @@ protected:
     // Protected member functions
 
 
+        //- Calculate and return the laminar viscosity
+        void calcMu();
+
         //- Generate the phases
         HashTable<autoPtr<phaseModel>> generatePhaseModels
         (
@@ -501,9 +510,15 @@ public:
         (
             const volScalarField& T
         ) = 0;
-
-        //- Calculate mass transfer
-        //virtual void massTransfer(const volScalarField& T) = 0;
+        
+        //- Return the volumetric rate transfer matrix
+        virtual tmp<fvScalarMatrix> volTransfer
+        (
+            const volScalarField& p
+        ) = 0;
+        
+        //- Calculate mass transfer for alpha's
+        virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp) = 0;
 
         //- Calculate mass transfer for species
         virtual void massSpeciesTransfer
@@ -522,6 +537,9 @@ public:
 
         //- Correct the mixture thermos
         virtual void correct();
+        
+        //- Correct mass sources
+        virtual void correctMassSources(const volScalarField& T) = 0;
 
         //- Return the name of the thermo physics
         virtual word thermoName() const
@@ -531,7 +549,7 @@ public:
         }
 
         //- Correct the turbulence
-        // (NOTE: Each phase could help its own turbulence)
+        // (NOTE: Each phase could have its own turbulence)
         virtual void correctTurbulence();
 
         //- Read base phaseProperties dictionary

applications/solvers/multiphase/interCondensatingEvaporatingFoam/Make/options

@@ -1,13 +1,18 @@
-interPhaseChangePath = $(FOAM_SOLVERS)/multiphase/interPhaseChangeFoam
+interPhaseChangeFoam = $(FOAM_SOLVERS)/multiphase/interPhaseChangeFoam
+interFoam = $(FOAM_SOLVERS)/multiphase/interFoam
+VoF = $(FOAM_SOLVERS)/multiphase/VoF
 
 EXE_INC = \
     -I. \
+    -I$(interPhaseChangeFoam) \
+    -I$(interFoam) \
+    -I$(VoF) \
     -ItemperaturePhaseChangeTwoPhaseMixtures/lnInclude \
-    -I$(interPhaseChangePath) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude\
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
@@ -21,6 +26,7 @@ EXE_LIBS = \
     -lfvOptions \
     -lmeshTools \
     -lsampling \
+    -ldynamicFvMesh \
     -lphaseTemperatureChangeTwoPhaseMixtures \
     -ltwoPhaseMixture \
     -linterfaceProperties \
@@ -29,3 +35,5 @@ EXE_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
     -lfluidThermophysicalModels
+
+    

applications/solvers/multiphase/interCondensatingEvaporatingFoam/TEqn.H

@@ -1,18 +1,19 @@
 {
     tmp<volScalarField> tcp(thermo->Cp());
     const volScalarField& cp = tcp();
+    
+    const dimensionedScalar Cp1 = thermo->Cp1();
+    const dimensionedScalar Cp2 = thermo->Cp2();
 
     rhoCp = rho*cp;
 
     kappaEff = thermo->kappa() + rho*cp*turbulence->nut()/Prt;
-
-    const surfaceScalarField rhoCpPhi(fvc::interpolate(cp)*rhoPhi);
-
-    Pair<tmp<volScalarField>> vDotAlphal = mixture->mDot();
-
-    const volScalarField& vDotcAlphal = vDotAlphal[0]();
-    const volScalarField& vDotvAlphal = vDotAlphal[1]();
-    const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
+    
+    const surfaceScalarField rhoCpPhi
+    (
+        "rhoCpPhi",
+        rhoPhi*(Cp1 - Cp2) + phi*rho2*Cp2
+    );
 
     fvScalarMatrix TEqn
     (
@@ -20,7 +21,7 @@
        + fvm::div(rhoCpPhi, T)
        - fvm::Sp(fvc::ddt(rhoCp) + fvc::div(rhoCpPhi), T)
        - fvm::laplacian(kappaEff, T)
-       + thermo->hc()*vDotvmcAlphal
+       + mixture->TSource()
     );
 
 

applications/solvers/multiphase/interCondensatingEvaporatingFoam/alphaCourantNo.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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/>.
+
+Global
+    alphaCourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    runTime.controlDict().get<scalar>("maxAlphaCo")
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        interface.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H

@@ -61,6 +61,20 @@ volScalarField rho
 );
 rho.oldTime();
 
+// Mass flux
+surfaceScalarField rhoPhi
+(
+    IOobject
+    (
+        "rhoPhi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::interpolate(rho)*phi
+);
+
 // Construct interface from alpha1 distribution
 interfaceProperties interface
 (
@@ -106,6 +120,24 @@ if (p_rgh.needReference())
 mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
+// MULES compressed flux is registered in case scalarTransport FO needs it.
+surfaceScalarField alphaPhiUn
+(
+    IOobject
+    (
+        "alphaPhiUn",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(phi.dimensions(), Zero)
+);
+
+#include "createMRF.H"
+#include "createFvOptions.H"
+
 // Turbulent Prandtl number
 dimensionedScalar Prt("Prt", dimless, thermo->transportPropertiesDict());
 

applications/solvers/multiphase/interCondensatingEvaporatingFoam/interCondensatingEvaporatingFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,7 +43,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
 #include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixtureEThermo.H"
@@ -52,7 +56,8 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fixedFluxPressureFvPatchScalarField.H"
+#include "CorrectPhi.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,11 +77,14 @@ int main(int argc, char *argv[])
     #include "addCheckCaseOptions.H"
     #include "setRootCaseLists.H"
     #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
     #include "createFields.H"
-    #include "createFvOptions.H"
-    #include "createTimeControls.H"
+    #include "createAlphaFluxes.H"
+    #include "initCorrectPhi.H"
+    #include "createUfIfPresent.H"
+    
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -90,9 +98,19 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "createTimeControls.H"
-        #include "CourantNo.H"
-        #include "setDeltaT.H"
+        #include "readDyMControls.H"
+        
+        // Store divU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+
+        volScalarField divU("divU", fvc::div(fvc::absolute(phi, U)));
+        
+        {
+            #include "CourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
 
         ++runTime;
 
@@ -101,28 +119,53 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            #include "alphaControls.H"
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                mesh.update();
 
-            surfaceScalarField rhoPhi
-            (
-                IOobject
-                (
-                    "rhoPhi",
-                    runTime.timeName(),
-                    mesh
-                ),
-                mesh,
-                dimensionedScalar(dimMass/dimTime, Zero)
-            );
+                if (mesh.changing())
+                {
+                    // Do not apply previous time-step mesh compression flux
+                    // if the mesh topology changed
+                    if (mesh.topoChanging())
+                    {
+                        talphaPhi1Corr0.clear();
+                    }
 
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    MRF.update();
+
+                    if (correctPhi)
+                    {
+                        // Calculate absolute flux
+                        // from the mapped surface velocity
+                        phi = mesh.Sf() & Uf();
+
+                        #include "correctPhi.H"
+
+                        // Make the flux relative to the mesh motion
+                        fvc::makeRelative(phi, U);
+
+                        interface.correct();
+                    }
+
+                    if (checkMeshCourantNo)
+                    {
+                        #include "meshCourantNo.H"
+                    }
+                }
+            }
+            
             mixture->correct();
 
+            #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 
-            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
             #include "UEqn.H"
             #include "TEqn.H"
-
+            
             // --- Pressure corrector loop
             while (pimple.correct())
             {

applications/solvers/multiphase/interCondensatingEvaporatingFoam/pEqn.H

@@ -33,7 +33,8 @@
         (
             fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
-          - (vDotc - vDotv)
+          ==
+            vDotv - vDotc
         );
 
         p_rghEqn.setReference(pRefCell, pRefValue);

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/Make/files

@@ -3,5 +3,6 @@ temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C
 thermoIncompressibleTwoPhaseMixture/thermoIncompressibleTwoPhaseMixture.C
 twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C
 constant/constant.C
+interfaceHeatResistance/interfaceHeatResistance.C
 
 LIB = $(FOAM_LIBBIN)/libphaseTemperatureChangeTwoPhaseMixtures

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvcGrad.H"
 #include "twoPhaseMixtureEThermo.H"
+#include "fvmSup.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -57,11 +58,15 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::constant
     temperaturePhaseChangeTwoPhaseMixture(mixture, mesh),
     coeffC_
     (
-        "coeffC", dimless/dimTime/dimTemperature, subDict(type() + "Coeffs")
+        "coeffC", 
+        dimless/dimTime/dimTemperature, 
+        optionalSubDict(type() + "Coeffs")
     ),
     coeffE_
     (
-        "coeffE", dimless/dimTime/dimTemperature, subDict(type() + "Coeffs")
+        "coeffE", 
+        dimless/dimTime/dimTemperature, 
+        optionalSubDict(type() + "Coeffs")
     )
 {}
 
@@ -85,8 +90,8 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotAlphal() const
 
     return Pair<tmp<volScalarField>>
     (
-        coeffC_*mixture_.rho2()*max(TSat - T.oldTime(), T0),
-       -coeffE_*mixture_.rho1()*max(T.oldTime() - TSat, T0)
+        coeffC_*mixture_.rho2()*max(TSat - T, T0),
+       -coeffE_*mixture_.rho1()*max(T - TSat, T0)
     );
 }
 
@@ -116,11 +121,20 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
     const dimensionedScalar& TSat = thermo.TSat();
 
     const dimensionedScalar T0(dimTemperature, Zero);
+    
+    if (mesh_.time().outputTime())
+    {
+        volScalarField mDot
+        (
+            "mDot", coeffE_*mixture_.rho1()*limitedAlpha1*max(T - TSat, T0)
+        );
+        mDot.write();
+    }
 
     return Pair<tmp<volScalarField>>
     (
-        coeffC_*mixture_.rho2()*limitedAlpha2*max(TSat - T.oldTime(), T0),
-        coeffE_*mixture_.rho1()*limitedAlpha1*max(T.oldTime() - TSat, T0)
+        coeffC_*mixture_.rho2()*limitedAlpha2*max(TSat - T, T0),
+       -coeffE_*mixture_.rho1()*limitedAlpha1*max(T - TSat, T0)
     );
 }
 
@@ -150,12 +164,66 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotDeltaT() const
 
     return Pair<tmp<volScalarField>>
     (
-        coeffC_*mixture_.rho2()*limitedAlpha2*pos(TSat - T.oldTime()),
-        coeffE_*mixture_.rho1()*limitedAlpha1*pos(T.oldTime() - TSat)
+        coeffC_*mixture_.rho2()*limitedAlpha2*pos(TSat - T),
+        coeffE_*mixture_.rho1()*limitedAlpha1*pos(T - TSat)
     );
 }
 
 
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
+{
+
+    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
+
+    tmp<fvScalarMatrix> tTSource
+    (
+        new fvScalarMatrix
+        (
+            T,
+            dimEnergy/dimTime
+        )
+    );
+
+    fvScalarMatrix& TSource = tTSource.ref();
+
+    const twoPhaseMixtureEThermo& thermo =
+        refCast<const twoPhaseMixtureEThermo>
+        (
+            mesh_.lookupObject<basicThermo>(basicThermo::dictName)
+        );
+
+    const dimensionedScalar& TSat = thermo.TSat();
+    
+    dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
+    
+    volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    volScalarField limitedAlpha2
+    (
+        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
+    );
+   
+    const volScalarField Vcoeff
+    (
+        coeffE_*mixture_.rho1()*limitedAlpha1*L
+    );
+    const volScalarField Ccoeff
+    (
+        coeffC_*mixture_.rho2()*limitedAlpha2*L
+    );
+
+    TSource = 
+        fvm::Sp(Vcoeff, T) - Vcoeff*TSat
+      - fvm::Sp(Ccoeff, T) + Ccoeff*TSat;
+
+    return tTSource;
+}
+
+
 void Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::correct()
 {
 }

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.H

@@ -98,6 +98,9 @@ public:
         //  coefficient to multiply (Tsat - T) for the condensation rate
         //  and a coefficient to multiply  (T - Tsat) for the vaporisation rate
         virtual Pair<tmp<volScalarField>> mDotDeltaT() const;
+        
+        //- Source for T equarion
+        virtual tmp<fvScalarMatrix> TSource() const;
 
         //- Correct the constant phaseChange model
         virtual void correct();

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C

@@ -0,0 +1,409 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "addToRunTimeSelectionTable.H"
+#include "twoPhaseMixtureEThermo.H"
+#include "isoCutCell.H"
+#include "volPointInterpolation.H"
+#include "calculatedFvPatchFields.H"
+#include "wallPolyPatch.H"
+#include "fvcSmooth.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace temperaturePhaseChangeTwoPhaseMixtures
+{
+    defineTypeNameAndDebug(interfaceHeatResistance, 0);
+    addToRunTimeSelectionTable
+    (
+        temperaturePhaseChangeTwoPhaseMixture,
+        interfaceHeatResistance,
+        components
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+interfaceHeatResistance
+(
+    const thermoIncompressibleTwoPhaseMixture& mixture,
+    const fvMesh& mesh
+)
+:
+    temperaturePhaseChangeTwoPhaseMixture(mixture, mesh),
+    R_
+    (
+        "R", 
+        dimPower/dimArea/dimTemperature, optionalSubDict(type() + "Coeffs")
+    ),
+    
+    interfaceArea_
+    (
+        IOobject
+        (
+            "interfaceArea",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimless/dimLength, Zero)
+    ),
+    
+    mDotc_
+    (
+        IOobject
+        (
+            "mDotc",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+    
+    mDote_
+    (
+        IOobject
+        (
+            "mDote",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimDensity/dimTime, Zero)
+    ),
+    
+    spread_
+    (
+        optionalSubDict(type() + "Coeffs").get<scalar>("spread")
+    )
+{
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+vDotAlphal() const
+{
+    dimensionedScalar alphalCoeff(1.0/mixture_.rho1());
+
+    return Pair<tmp<volScalarField>>
+    (
+        (alphalCoeff*mDotc_)/(mixture_.alpha2() + SMALL),
+       -(alphalCoeff*mDote_)/(mixture_.alpha1() + SMALL)
+    );
+}
+
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+mDotAlphal() const
+{
+    volScalarField limitedAlpha1
+    (
+        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
+    );
+
+    volScalarField limitedAlpha2
+    (
+        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
+    );
+    
+    return Pair<tmp<volScalarField>>
+    (
+        (mDotc_/(limitedAlpha2 + SMALL)),
+       -(mDote_/(limitedAlpha1 + SMALL))
+    );
+}
+
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+mDot() const
+{
+    
+    return Pair<tmp<volScalarField>>
+    (
+        tmp<volScalarField>(mDotc_), 
+        tmp<volScalarField>(mDote_)
+    );
+}
+
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+mDotDeltaT() const
+{
+   const twoPhaseMixtureEThermo& thermo =
+        refCast<const twoPhaseMixtureEThermo>
+        (
+            mesh_.lookupObject<basicThermo>(basicThermo::dictName)
+        );
+        
+    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
+
+    const dimensionedScalar& TSat = thermo.TSat();
+
+    Pair<tmp<volScalarField>> mDotce(mDot());
+    
+    return Pair<tmp<volScalarField>>
+    (
+        mDotc_*pos(TSat - T.oldTime())/(TSat - T.oldTime()),
+       -mDote_*pos(T.oldTime() - TSat)/(T.oldTime() - TSat)
+    );
+
+    
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+TSource() const
+{
+
+    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
+
+    tmp<fvScalarMatrix> tTSource
+    (
+        new fvScalarMatrix
+        (
+            T,
+            dimEnergy/dimTime
+        )
+    );
+
+    fvScalarMatrix& TSource = tTSource.ref();
+
+    const twoPhaseMixtureEThermo& thermo =
+        refCast<const twoPhaseMixtureEThermo>
+        (
+            mesh_.lookupObject<basicThermo>(basicThermo::dictName)
+        );
+
+    const dimensionedScalar& TSat = thermo.TSat();
+   
+    // interface heat resistance
+    volScalarField IHRcoeff = interfaceArea_*R_;
+
+    TSource = fvm::Sp(IHRcoeff, T) - IHRcoeff*TSat;
+
+    return tTSource;
+}
+
+
+void Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+correct()
+{
+   
+    // Update Interface
+    updateInterface();
+    
+    // Update mDotc_ and mDote_
+    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
+
+    const twoPhaseMixtureEThermo& thermo =
+        refCast<const twoPhaseMixtureEThermo>
+        (
+            mesh_.lookupObject<basicThermo>(basicThermo::dictName)
+        );
+
+    const dimensionedScalar& TSat = thermo.TSat();
+    const dimensionedScalar T0(dimTemperature, Zero);
+
+    dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
+   
+    // interface heat resistance
+    mDotc_ = interfaceArea_*R_*max(TSat - T, T0)/L;
+    mDote_ = interfaceArea_*R_*(T - TSat)/L;
+       
+    forAll(mDotc_, celli)
+    {
+        scalar rhobyDt = mixture_.rho1().value()/mesh_.time().deltaTValue();
+        scalar maxEvap = mixture_.alpha1()[celli]*rhobyDt; // positive
+        scalar maxCond = -mixture_.alpha2()[celli]*rhobyDt; // negative 
+        mDote_[celli] = min(max(mDote_[celli], maxCond), maxEvap);
+        mDotc_[celli] = min(max(mDotc_[celli], maxCond), maxEvap);
+    }
+}
+
+
+void Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+updateInterface()
+{
+    const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
+    const twoPhaseMixtureEThermo& thermo =
+    refCast<const twoPhaseMixtureEThermo>
+    (
+        mesh_.lookupObject<basicThermo>(basicThermo::dictName)
+    );
+
+    const dimensionedScalar& TSat = thermo.TSat();
+   
+    // interface heat resistance
+    // Interpolating alpha1 cell centre values to mesh points (vertices)
+    scalarField ap
+    (
+        volPointInterpolation::New(mesh_).interpolate(mixture_.alpha1())
+    );
+
+    isoCutCell cutCell(mesh_, ap);
+
+    forAll(interfaceArea_, celli)
+    {
+        label status = cutCell.calcSubCell(celli, 0.5);
+        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(),i)
+            {
+                const label pCelli = pp.faceCells()[i];
+ 
+                if
+                (
+                    (TSat.value() - T[pCelli]) > 0 
+                  && mixture_.alpha1()[pCelli] < 0.9
+                )
+                {
+                    interfaceArea_[pCelli] = 
+                        mag(pp.faceAreas()[i])/mesh_.V()[pCelli];
+                }
+            }
+        }
+    }
+}
+
+Foam::Pair<Foam::tmp<Foam::volScalarField>>
+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]
+    );
+}
+
+
+
+bool Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
+read()
+{
+    if (temperaturePhaseChangeTwoPhaseMixture::read())
+    {
+        optionalSubDict(type() + "Coeffs").readEntry("R", R_);
+        optionalSubDict(type() + "Coeffs").readEntry("spread", spread_);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::phaseChangeTwoPhaseMixtures::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
+
+
+
+SourceFiles
+    interfaceHeatResistance.C
+
+\*--------------------------------------------------------------------*/
+
+#ifndef interfaceHeatResistance_H
+#define interfaceHeatResistance_H
+
+#include "temperaturePhaseChangeTwoPhaseMixture.H"
+
+// * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace temperaturePhaseChangeTwoPhaseMixtures
+{
+
+/*--------------------------------------------------------------------*\
+                    Class interfaceHeatResistance
+\*--------------------------------------------------------------------*/
+
+class interfaceHeatResistance
+:
+    public temperaturePhaseChangeTwoPhaseMixture
+{
+    // Private data
+
+        //- Heat transfer coefficient [1/s/K]
+        dimensionedScalar R_;
+
+        //- Interface area
+        volScalarField interfaceArea_;
+        
+        //- Mass condensation source
+        volScalarField mDotc_;
+        
+        //- Mass evaporation source
+        volScalarField mDote_;
+              
+        //- Spread for mass source
+        scalar spread_;
+
+        
+    // Private member functions
+
+        //- Update interface area
+        void updateInterface();
+        
+
+public:
+
+    //- Runtime type information
+    TypeName("interfaceHeatResistance");
+
+
+    // Constructors
+
+        //- Construct from components
+        interfaceHeatResistance
+        (
+            const thermoIncompressibleTwoPhaseMixture& mixture,
+            const fvMesh& mesh
+        );
+
+    //- Destructor
+    virtual ~interfaceHeatResistance()
+    {}
+
+
+    // Member Functions
+
+        //- Return the mass condensation and vaporisation rates as a
+        //  coefficient to multiply (1 - alphal) for the condensation rate
+        //  and a coefficient to multiply  alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField>> mDotAlphal() const;
+
+        //- Return the mass condensation and vaporisation rates as coefficients
+        virtual Pair<tmp<volScalarField>> mDot() const;
+
+        //- Return the mass condensation and vaporisation rates as a
+        //  coefficient to multiply (Tsat - T) for the condensation rate
+        //  and a coefficient to multiply  (T - Tsat) for the vaporisation rate
+        virtual Pair<tmp<volScalarField>> mDotDeltaT() const;
+        
+        //- Source for T equarion
+        virtual tmp<fvScalarMatrix> TSource() const;
+        
+        //- Volumetric source for alpha (used by alphaEq)
+        virtual Pair<tmp<volScalarField>> vDotAlphal() const;
+        
+        //- Return the volumetric condensation and vaporisation rates as
+        //  coefficients (used by p_rghEq)
+        virtual Pair<tmp<volScalarField>> vDot() const;
+        
+        //- Correct the interfaceHeatResistance phaseChange model
+        virtual void correct();
+
+        //- Read the transportProperties dictionary and update
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace temperaturePhaseChangeTwoPhaseMixtures
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -68,7 +68,7 @@ temperaturePhaseChangeTwoPhaseMixture
 
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixture::vDotAlphal() const
-{
+{ 
     volScalarField alphalCoeff
     (
         1.0/mixture_.rho1() - mixture_.alpha1()

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixtures/temperaturePhaseChangeTwoPhaseMixture.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "volFields.H"
 #include "dimensionedScalar.H"
+#include "fvMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -140,15 +141,18 @@ public:
         //  coefficient to multiply (Tsat - T) for the condensation rate
         //  and a coefficient to multiply  (T - Tsat) for the vaporisation rate
         virtual Pair<tmp<volScalarField>> mDotDeltaT() const = 0;
+        
+        //- Source for T equarion
+        virtual tmp<fvScalarMatrix> TSource() const = 0;
 
         //- Return the volumetric condensation and vaporisation rates as a
         //  coefficient to multiply (1 - alphal) for the condensation rate
         //  and a coefficient to multiply  alphal for the vaporisation rate
-        Pair<tmp<volScalarField>> vDotAlphal() const;
+        virtual Pair<tmp<volScalarField>> vDotAlphal() const;
 
         //- Return the volumetric condensation and vaporisation rates as
         //  coefficients
-        Pair<tmp<volScalarField>> vDot() const;
+        virtual Pair<tmp<volScalarField>> vDot() const;
 
         //- Correct the phaseChange model
         virtual void correct() = 0;

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/thermoIncompressibleTwoPhaseMixture/thermoIncompressibleTwoPhaseMixture.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/thermoIncompressibleTwoPhaseMixture/thermoIncompressibleTwoPhaseMixture.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,44 +40,6 @@ namespace Foam
     defineTypeNameAndDebug(twoPhaseMixtureEThermo, 0);
 }
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-void Foam::twoPhaseMixtureEThermo::eBoundaryCorrection(volScalarField& h)
-{
-    volScalarField::Boundary& hbf = h.boundaryFieldRef();
-
-    forAll(hbf, patchi)
-    {
-        if (isA<gradientEnergyFvPatchScalarField>(hbf[patchi]))
-        {
-            refCast<gradientEnergyFvPatchScalarField>(hbf[patchi]).gradient()
-                = hbf[patchi].fvPatchField::snGrad();
-        }
-        else if (isA<mixedEnergyFvPatchScalarField>(hbf[patchi]))
-        {
-            refCast<mixedEnergyFvPatchScalarField>(hbf[patchi]).refGrad()
-                = hbf[patchi].fvPatchField::snGrad();
-        }
-    }
-}
-
-
-void Foam::twoPhaseMixtureEThermo::init()
-{
-    const volScalarField alpha1Rho1(alpha1()*rho1());
-    const volScalarField alpha2Rho2(alpha2()*rho2());
-
-    e_ =
-        (
-            (T_ - TSat_)*(alpha1Rho1*Cv1() + alpha2Rho2*Cv2())
-          + (alpha1Rho1*Hf1() + alpha2Rho2*Hf2())
-        )
-       /(alpha1Rho1 + alpha2Rho2);
-
-    e_.correctBoundaryConditions();
-}
-
-
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
 Foam::twoPhaseMixtureEThermo::twoPhaseMixtureEThermo
@@ -89,32 +51,8 @@ Foam::twoPhaseMixtureEThermo::twoPhaseMixtureEThermo
     basicThermo(U.mesh(), word::null),
     thermoIncompressibleTwoPhaseMixture(U, phi),
 
-    e_
-    (
-        volScalarField
-        (
-            IOobject
-            (
-                "e",
-                U.mesh().time().timeName(),
-                U.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            U.mesh(),
-            dimensionedScalar(dimEnergy/dimMass, Zero),
-            heBoundaryTypes()
-        )
-    ),
-
-    TSat_("TSat", dimTemperature, static_cast<const basicThermo&>(*this)),
-
-    pDivU_(basicThermo::lookupOrDefault<Switch>("pDivU", true))
-
-{
-    // Initialise e
-    init();
-}
+    TSat_("TSat", dimTemperature, static_cast<const basicThermo&>(*this))
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -122,19 +60,6 @@ Foam::twoPhaseMixtureEThermo::twoPhaseMixtureEThermo
 void Foam::twoPhaseMixtureEThermo::correct()
 {
     incompressibleTwoPhaseMixture::correct();
-
-    const volScalarField alpha1Rho1(alpha1()*rho1());
-    const volScalarField alpha2Rho2(alpha2()*rho2());
-
-    T_ =
-        (
-            (e_*(alpha1Rho1 + alpha2Rho2))
-         -  (alpha1Rho1*Hf1() + alpha2Rho2*Hf2())
-        )
-       /(alpha1Rho1*Cv1() + alpha2Rho2*Cv2())
-       + TSat_;
-
-    T().correctBoundaryConditions();
 }
 
 
@@ -151,15 +76,8 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::he
     const volScalarField& T
 ) const
 {
-    const volScalarField alpha1Rho1(alpha1()*rho1());
-    const volScalarField alpha2Rho2(alpha2()*rho2());
-
-    return
-    (
-        (T - TSat_)*(alpha1Rho1*Cv1() + alpha2Rho2*Cv2())
-        + (alpha1Rho1*Hf1() + alpha2Rho2*Hf2())
-    )
-    / (alpha1Rho1 + alpha2Rho2);
+    NotImplemented;
+    return nullptr;
 }
 
 
@@ -170,31 +88,8 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::he
     const labelList& cells
 ) const
 {
-    tmp<scalarField> the(new scalarField(T.size()));
-    scalarField& he = the.ref();
-
-    const volScalarField alpha1Rho1(alpha1()*rho1());
-    const volScalarField alpha2Rho2(alpha2()*rho2());
-
-    forAll(T, i)
-    {
-        const label celli = cells[i];
-        he[i] =
-            (
-                (T[i] - TSat_.value())
-               *(
-                   alpha1Rho1[celli]*Cv1().value()
-                 + alpha2Rho2[celli]*Cv2().value()
-                )
-              + (
-                    alpha1Rho1[celli]*Hf1().value()
-                  + alpha2Rho2[celli]*Hf2().value()
-                )
-            )
-            / (alpha1Rho1[celli] + alpha2Rho2[celli]);
-    }
-
-    return the;
+    NotImplemented;
+    return nullptr;
 }
 
 
@@ -205,26 +100,8 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::he
     const label patchi
 ) const
 {
-    const scalarField& alpha1p = alpha1().boundaryField()[patchi];
-    const scalarField& alpha2p = alpha2().boundaryField()[patchi];
-
-    const scalarField& Tp = T_.boundaryField()[patchi];
-
-    return
-    (
-        (
-            (Tp - TSat_.value())
-           *(
-               alpha1p*rho1().value()*Cv1().value()
-             + alpha2p*rho2().value()*Cv2().value()
-            )
-          + (
-               alpha1p*rho1().value()*Hf1().value()
-             + alpha2p*rho2().value()*Hf2().value()
-            )
-        )
-        / (alpha1p*rho1().value() + alpha2p*rho2().value())
-    );
+    NotImplemented;
+    return nullptr;
 }
 
 
@@ -252,7 +129,7 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::THE
 (
     const scalarField& h,
     const scalarField& p,
-    const scalarField& T0,      // starting temperature
+    const scalarField& T0,
     const labelList& cells
 ) const
 {
@@ -265,7 +142,7 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::THE
 (
     const scalarField& h,
     const scalarField& p,
-    const scalarField& T0,      // starting temperature
+    const scalarField& T0,
     const label patchi
 ) const
 {
@@ -588,7 +465,6 @@ bool Foam::twoPhaseMixtureEThermo::read()
 {
     if (basicThermo::read() && thermoIncompressibleTwoPhaseMixture::read())
     {
-        basicThermo::readIfPresent("pDivU", pDivU_);
         basicThermo::readEntry("TSat", TSat_);
         return true;
     }

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,10 +27,6 @@ Class
     Foam::twoPhaseMixtureEThermo
 
 Description
-    Two phases thermo Internal energy mixture Defined as:
-    e1 = Cv1(T - Tsat) + Hv1
-    e2 = Cv2(T - Tsat) + Hv2
-    e = (alpha1*rho1*e1 + alpha2*rho2*e2) / (alpha1*rho1 + alpha2*rho2)
 
 SourceFiles
     twoPhaseMixtureEThermo.C
@@ -64,25 +60,9 @@ protected:
 
     // Protected Data
 
-        //- Internal energy field [J]
-        volScalarField e_;
-
         //- Saturation Temperature
         dimensionedScalar TSat_;
 
-        //- Should the p*div(U) term be included in the energy equation
-        Switch pDivU_;
-
-
-    // Protected Member Functions
-
-        //- Correct boundary for e
-        void eBoundaryCorrection (volScalarField& h);
-
-        //- Initialize
-        void init();
-
-
 public:
 
     TypeName("twoPhaseMixtureEThermo");
@@ -105,13 +85,15 @@ public:
         //- Return access to the internal energy field [J/Kg]
         virtual volScalarField& he()
         {
-            return e_;
+            NotImplemented;
+            return p();
         }
 
         //- Return access to the internal energy field [J/Kg]
         virtual const volScalarField& he() const
         {
-            return e_;
+            NotImplemented;
+            return p();
         }
 
         //- Enthalpy/Internal energy
@@ -306,12 +288,6 @@ public:
             {
                 return *this;
             }
-
-            //- Should the dpdt term be included in the enthalpy equation
-            Switch pDivU() const
-            {
-                return pDivU_;
-            }
 };
 
 

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -114,6 +114,9 @@
     }
 
     rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
+    
+    // Cache alphaPhi
+    alphaPhi10 = talphaPhi();
 
     Info<< "Liquid phase volume fraction = "
         << alpha1.weightedAverage(mesh.V()).value()

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -14,7 +14,17 @@
     if (nAlphaSubCycles > 1)
     {
         dimensionedScalar totalDeltaT = runTime.deltaT();
-        surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi);
+        surfaceScalarField rhoPhiSum
+        (
+            IOobject
+            (
+                "rhoPhiSum",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(rhoPhi.dimensions(), Zero)
+        );
 
         for
         (

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -109,3 +109,19 @@ mesh.setFluxRequired(p_rgh.name());
 mesh.setFluxRequired(alpha1.name());
 
 #include "createFvOptions.H"
+
+IOobject alphaPhi10Header
+(
+    IOobject::groupName("alphaPhi0", alpha1.group()),
+    runTime.timeName(),
+    mesh,
+    IOobject::NO_READ,
+    IOobject::NO_WRITE
+);
+
+// MULES flux from previous time-step
+surfaceScalarField alphaPhi10
+(
+    alphaPhi10Header,
+    phi*fvc::interpolate(alpha1)
+);

src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C

@@ -6,6 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 Henning Scheufler
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +32,11 @@ License
 #include "FaceCellWave.H"
 #include "smoothData.H"
 #include "sweepData.H"
+#include "fvMatrices.H"
+#include "fvcVolumeIntegrate.H"
+#include "fvmLaplacian.H"
+#include "fvmSup.H"
+#include "zeroGradientFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -314,4 +321,100 @@ void Foam::fvc::sweep
 }
 
 
+void Foam::fvc::spreadSource
+(
+    volScalarField& mDotOut,
+    const volScalarField& mDotIn,
+    const volScalarField& alpha1,
+    const volScalarField& alpha2,
+    const dimensionedScalar& D,
+    const scalar cutoff
+)
+{
+    const fvMesh& mesh = alpha1.mesh();
+   
+    volScalarField mDotSmear
+    (
+        IOobject
+        (
+            "mDotSmear",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(mDotOut.dimensions(), Zero),
+        zeroGradientFvPatchField<scalar>::typeName
+    );
+
+    //- Smearing of source term field
+    fvScalarMatrix mSourceEqn
+    (
+        fvm::Sp(scalar(1), mDotSmear) 
+      - fvm::laplacian(D, mDotSmear) 
+      == 
+        mDotIn
+    );
+
+    mSourceEqn.solve();
+
+    // Cut cells with cutoff < alpha1 < 1-cutoff and rescale remaining 
+    // source term field
+    
+    dimensionedScalar intvDotLiquid("intvDotLiquid", dimMass/dimTime, 0.0);
+    dimensionedScalar intvDotVapor ("intvDotVapor", dimMass/dimTime, 0.0);
+    
+    const scalarField& Vol = mesh.V();
+    
+    forAll(mesh.C(), celli)
+    {
+        if (alpha1[celli] < cutoff)
+        {
+            intvDotVapor.value() += 
+                alpha2[celli]*mDotSmear[celli]*Vol[celli];
+        }
+        else if (alpha1[celli] > 1.0 - cutoff)
+        {
+            intvDotLiquid.value() += 
+                alpha1[celli]*mDotSmear[celli]*Vol[celli];
+        }
+    }
+    
+    reduce(intvDotVapor.value(), sumOp<scalar>());
+    reduce(intvDotLiquid.value(), sumOp<scalar>());
+
+    //- Calculate Nl and Nv
+    dimensionedScalar Nl ("Nl", dimless, Zero);
+    dimensionedScalar Nv ("Nv", dimless, Zero); 
+    
+    const dimensionedScalar intmSource0(fvc::domainIntegrate(mDotIn));
+
+    if (intvDotVapor.value() > VSMALL)
+    {
+        Nv = intmSource0/intvDotVapor;
+    }
+    if (intvDotLiquid.value() > VSMALL)
+    {
+        Nl = intmSource0/intvDotLiquid;
+    }
+
+    //- Set source terms in cells with alpha1 < cutoff or alpha1 > 1-cutoff
+    forAll(mesh.C(), celli)
+    {
+        if (alpha1[celli] < cutoff)
+        {
+            mDotOut[celli] = Nv.value()*(1 - alpha1[celli])*mDotSmear[celli];
+        }
+        else if (alpha1[celli] > 1.0 - cutoff)
+        {
+            //mDotOut[celli] = 0;
+            mDotOut[celli] = -Nl.value()*alpha1[celli]*mDotSmear[celli];
+        }
+        else
+        {
+            mDotOut[celli] = 0;
+        }
+    }
+}
 // ************************************************************************* //

src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H

@@ -6,6 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 Henning Scheufler
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +44,12 @@ Description
            gradient of alpha is large (where the difference between the values
            in neighbouring cells is larger than alphaDiff) away from that
            starting point of the sweep.
+           
+           
+    spreadSource: spread a source field (mDotIn) for two phase multiphase using
+                  a laplacian operator and diffussivity D.
+                  The spread source (mDotOut) is distributed from alpha1 < cutoff 
+                  to alpha1 > 1 - cutoff, and it is zero across the interface
 
 SourceFiles
     fvcSmooth.C
@@ -52,6 +60,7 @@ SourceFiles
 #define fvcSmooth_H
 
 #include "volFieldsFwd.H"
+#include "dimensionedScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,6 +91,16 @@ namespace fvc
         const label nLayers,
         const scalar alphaDiff = 0.2
     );
+    
+    void spreadSource
+    (
+        volScalarField& mDotOut,
+        const volScalarField& mDotIn,
+        const volScalarField& alpha1,
+        const volScalarField& alpha2,
+        const dimensionedScalar& D,
+        const scalar cutoff
+    );
 }
 }
 

src/functionObjects/field/Make/files

@@ -115,4 +115,6 @@ zeroGradient/zeroGradient.C
 
 stabilityBlendingFactor/stabilityBlendingFactor.C
 
+interfaceHeight/interfaceHeight.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects

src/functionObjects/field/interfaceHeight/interfaceHeight.C

@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017-2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "interfaceHeight.H"
+#include "fvMesh.H"
+#include "interpolation.H"
+#include "IOmanip.H"
+#include "meshSearch.H"
+#include "midPointAndFaceSet.H"
+#include "Time.H"
+#include "uniformDimensionedFields.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(interfaceHeight, 0);
+    addToRunTimeSelectionTable(functionObject, interfaceHeight, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::interfaceHeight::writePositions()
+{
+    const uniformDimensionedVectorField& g =
+        mesh_.time().lookupObject<uniformDimensionedVectorField>("g");
+    vector gHat = vector::zero;
+    
+    if (mag(direction_) > 0.0)
+    {
+        gHat = direction_/mag(direction_);
+    }
+    else
+    {
+        gHat = g.value()/mag(g.value());
+    }
+
+    const volScalarField& alpha =
+        mesh_.lookupObject<volScalarField>(alphaName_);
+
+    autoPtr<interpolation<scalar>>
+        interpolator
+        (
+            interpolation<scalar>::New(interpolationScheme_, alpha)
+        );
+
+    if (Pstream::master())
+    {
+        file(fileID::heightFile) << mesh_.time().timeName() << tab;
+        file(fileID::positionFile) << mesh_.time().timeName() << tab;
+    }
+
+    forAll(locations_, li)
+    {
+        // Create a set along a ray projected in the direction of gravity
+        const midPointAndFaceSet set
+        (
+            "",
+            mesh_,
+            meshSearch(mesh_),
+            "xyz",
+            locations_[li] + gHat*mesh_.bounds().mag(),
+            locations_[li] - gHat*mesh_.bounds().mag()
+        );
+
+        // Find the height of the location above the boundary
+        scalar hLB = set.size() ? - gHat & (locations_[li] - set[0]) : - GREAT;
+        reduce(hLB, maxOp<scalar>());
+
+        // Calculate the integrals of length and length*alpha along the sampling
+        // line. The latter is equal to the equivalent length with alpha equal
+        // to one.
+        scalar sumLength = 0, sumLengthAlpha = 0;
+        for(label si = 0; si < set.size() - 1; ++ si)
+        {
+            if (set.segments()[si] != set.segments()[si+1])
+            {
+                continue;
+            }
+
+            const vector& p0 = set[si], p1 = set[si+1];
+            const label c0 = set.cells()[si], c1 = set.cells()[si+1];
+            const label f0 = set.faces()[si], f1 = set.faces()[si+1];
+            const scalar a0 = interpolator->interpolate(p0, c0, f0);
+            const scalar a1 = interpolator->interpolate(p1, c1, f1);
+
+            const scalar l = - gHat & (p1 - p0);
+            sumLength += l;
+            sumLengthAlpha += l*(a0 + a1)/2;
+        }
+
+        reduce(sumLength, sumOp<scalar>());
+        reduce(sumLengthAlpha, sumOp<scalar>());
+
+        // Write out
+        if (Pstream::master())
+        {
+            // Interface heights above the boundary and location
+            const scalar hIB =
+                liquid_ ? sumLengthAlpha : sumLength - sumLengthAlpha;
+            const scalar hIL = hIB - hLB;
+
+            // Position of the interface
+            const point p = locations_[li] - gHat*hIL;
+
+            const Foam::Omanip<int> w = valueWidth(1);
+
+            file(fileID::heightFile) << w << hIB << w << hIL;
+            file(fileID::positionFile) << '(' << w << p.x() << w << p.y()
+                << valueWidth() << p.z() << ") ";
+        }
+    }
+
+    if (Pstream::master())
+    {
+        file(fileID::heightFile).endl();
+        file(fileID::positionFile).endl();
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::interfaceHeight::writeFileHeader(const fileID i)
+{
+    forAll(locations_, li)
+    {
+        writeHeaderValue
+        (
+            file(i),
+            "Location " + Foam::name(li),
+            locations_[li]
+        );
+    }
+
+    switch (fileID(i))
+    {
+        case fileID::heightFile:
+            writeHeaderValue
+            (
+                file(fileID::heightFile),
+                "hB",
+                "Interface height above the boundary"
+            );
+            writeHeaderValue
+            (
+                file(fileID::heightFile),
+                "hL",
+                "Interface height above the location"
+            );
+            break;
+        case fileID::positionFile:
+            writeHeaderValue(file(i), "p", "Interface position");
+            break;
+    }
+
+    const Foam::Omanip<int> w = valueWidth(1);
+
+    writeCommented(file(i), "Location");
+    forAll(locations_, li)
+    {
+        switch (fileID(i))
+        {
+            case fileID::heightFile:
+                file(i) << w << li << w << ' ';
+                break;
+            case fileID::positionFile:
+                file(i) << w << li << w << ' ' << w << ' ' << "  ";
+                break;
+        }
+    }
+    file(i).endl();
+
+    writeCommented(file(i), "Time");
+    forAll(locations_, li)
+    {
+        switch (fileID(i))
+        {
+            case fileID::heightFile:
+                file(i) << w << "hB" << w << "hL";
+                break;
+            case fileID::positionFile:
+                file(i) << w << "p" << w << ' ' << w << ' ' << "  ";
+                break;
+        }
+    }
+    file(i).endl();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::interfaceHeight::interfaceHeight
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    logFiles(obr_, name),
+    alphaName_("alpha"),
+    liquid_(true),
+    locations_(),
+    interpolationScheme_("cellPoint"),
+    direction_(vector::zero)
+{
+    read(dict);
+    resetNames({"height", "position"});
+    
+    writeFileHeader(fileID::heightFile);
+    writeFileHeader(fileID::positionFile);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::interfaceHeight::~interfaceHeight()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::interfaceHeight::read(const dictionary& dict)
+{
+    
+    dict.readIfPresent("alpha", alphaName_);
+    dict.readIfPresent("liquid", liquid_);
+    dict.lookup("locations") >> locations_;
+    dict.readIfPresent("interpolationScheme", interpolationScheme_);
+    dict.readIfPresent("direction", direction_);
+
+    return true;
+}
+
+
+bool Foam::functionObjects::interfaceHeight::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::interfaceHeight::end()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::interfaceHeight::write()
+{
+    logFiles::write();
+
+    writePositions();
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/interfaceHeight/interfaceHeight.H

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017-2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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::functionObjects::interfaceHeight
+
+Description
+    This function object reports the height of the interface above a set of
+    locations. For each location, it writes the vertical distance of the
+    interface above both the location and the lowest boundary. It also writes
+    the point on the interface from which these heights are computed. It uses
+    an integral approach, so if there are multiple interfaces above or below a
+    location then this method will generate average values.
+
+    Example of function object specification:
+    \verbatim
+    interfaceHeight1
+    {
+        type           interfaceHeight;
+        libs           ("libfieldFunctionObjects.so");
+        alpha          alpha.water;
+        locations      ((0 0 0) (10 0 0) (20 0 0));
+    }
+    \endverbatim
+
+Usage
+    \table
+        Property     | Description               | Required | Default value
+        type         | type name                 | yes      |
+        alpha        | name of the alpha field   | no       | alpha
+        locations    | list of locations to report the height at | yes |
+        liquid       | is the alpha field that of the liquid | no | true
+        direction    | direction of interface | no | g
+    \endtable
+
+SourceFiles
+    interfaceHeight.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interfaceHeight_H
+#define interfaceHeight_H
+
+#include "fvMeshFunctionObject.H"
+#include "logFiles.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class interfaceHeight Declaration
+\*---------------------------------------------------------------------------*/
+
+class interfaceHeight
+:
+    public fvMeshFunctionObject,
+    public logFiles
+{
+    // Private Data
+
+        //- Name of the alpha field
+        word alphaName_;
+
+        //- Is the alpha field that of the liquid under the wave?
+        bool liquid_;
+
+        //- List of locations to report the height at
+        List<point> locations_;
+
+        //- Interpolation scheme
+        word interpolationScheme_;
+        
+        
+        //- Direction of interface motion 
+        vector direction_;
+
+
+    // Private Member Functions
+
+        //- Output positions
+        void writePositions();
+
+
+    // Private Enumerations
+
+        //- File enumeration
+        enum class fileID
+        {
+            heightFile = 0,
+            positionFile = 1
+        };
+
+
+        Ostream& file(const fileID fid)
+        {
+            return logFiles::files(label(fid));
+        }
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Output file header information
+        virtual void writeFileHeader(const fileID i);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("interfaceHeight");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        interfaceHeight
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~interfaceHeight();
+
+
+    // Member Functions
+
+        //- Read
+        virtual bool read(const dictionary&);
+
+        //- Execute
+        virtual bool execute();
+
+        //- Execute at the final time-loop
+        virtual bool end();
+
+        //- Write
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/T

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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 373.15;
+
+boundaryField
+{
+    sideWalls
+    {
+        type            zeroGradient;
+    }
+
+    faceWall
+    {
+        type            fixedValue;
+        value           uniform 378.15;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 373.15;
+        value           uniform 373.15;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/U

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    sideWalls
+    {
+        type            slip;
+        value           uniform (0 0 0);
+    }
+
+    faceWall
+    {
+        type            noSlip;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.gas

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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
+{
+    ".*"
+    {
+        type            zeroGradient;
+    }
+    
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/alpha.liquid

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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 0;
+
+boundaryField
+{
+    ".*"
+    {
+        type            zeroGradient;
+    }
+    
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/p

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+  sideWalls
+  {
+      type            fixedFluxPressure;
+  }
+
+  faceWall
+  {
+      type            fixedFluxPressure;
+  }
+
+  outlet
+  {
+      //type            fixedValue;
+      type            totalPressure;
+      rho             rho;
+      p0              uniform 1e5;
+      value           uniform 1e5;
+  }
+
+  frontAndBack
+  {
+      type            empty;
+  }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/0.orig/p_rgh

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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
+{
+  sideWalls
+  {
+      type            fixedFluxPressure;
+  }
+
+  faceWall
+  {
+      type            fixedFluxPressure;
+  }
+
+  outlet
+  {
+      //type            fixedValue;
+      type            totalPressure;
+      rho             rho;
+      p0              uniform 1e5;
+      value           uniform 1e5;
+  }
+
+  frontAndBack
+  {
+      type            empty;
+  }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/Allclean

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+rm -rf 0
+rm positionClean.dat
+rm OF_vs_Exact.eps
+
+#------------------------------------------------------------------------------

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/Allrun

@@ -0,0 +1,19 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+restore0Dir
+cp -r 0 1.36
+cp system/setAlphaFieldDict.liquid system/setAlphaFieldDict
+runApplication setAlphaField
+rm log.setAlphaField
+
+cp system/setAlphaFieldDict.gas system/setAlphaFieldDict
+runApplication setAlphaField
+
+#runApplication $(getApplication)
+
+#------------------------------------------------------------------------------

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      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 0 0);
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/constant/phaseProperties

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      phaseChangeProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+type    massTransferMultiphaseSystem;
+
+phases  (liquid gas);
+
+liquid
+{
+    type            pureMovingPhaseModel;
+}
+
+gas
+{
+    type            pureMovingPhaseModel;
+}
+
+surfaceTension
+(
+    (gas and liquid)
+    {
+        type            constant;
+        sigma           0.0;
+    }
+);
+
+massTransferModel
+(
+    (liquid to gas)
+    {
+        type            Lee;
+        C               1700;
+        Tactivate       373.25;
+    }
+);
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/constant/thermophysicalProperties.gas

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+// ************************************************************************* //
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18.9;
+    }
+    equationOfState
+    {
+        rho         0.581;
+    }
+    thermodynamics
+    {
+        Hf          0;//-1.338e7; //[J/kg]
+        Cp          2030;
+    }
+    transport
+    {
+        mu          0.9e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/constant/thermophysicalProperties.liquid

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+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.26e6;//-1.5833e7;//deltaHv 2.45e6; //[J/Kg]
+    }
+    transport
+    {
+        mu              959e-6;
+        Pr              6.62;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/constant/transportProperties

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (liquid vapour); // FC-72
+
+liquid
+{
+    transportModel  Newtonian;
+    nu              1e-6;
+    rho             958.4;
+
+
+    Cp              4216;  // irrelevant
+    Cv              4216;  // irrelevant
+    kappa           0.671; // irrelevant
+    hf              0;
+}
+
+vapour
+{
+    transportModel  Newtonian;
+    nu              1e-5; // irrelevant
+    rho             0.581;
+
+
+    Cp              2030;      // FC72 vapour
+    Cv              2030;      // Cv = Cp - R/w
+    kappa           0.025;       // FC72 vapour // 0.01;
+    hf              2260.0e3;
+}
+
+sigma 0;
+
+Prt             0.7;
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/constant/turbulenceProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/plotGraph

@@ -0,0 +1,25 @@
+createGraphs()
+{
+    OF=$1
+    EXPT=$2
+
+    gnuplot<<EOF
+    set terminal postscript default
+    set output "OF_vs_Exact.eps"
+    set xlabel "t [sec]"
+    set ylabel "x [mm]"
+    set grid
+
+
+    plot \
+        "$EXPT" u 1:2 title "Exact", \
+        "$OF" u 1:2 title "OpenFOAM" with line lt -1 lw 1
+EOF
+}
+
+sed -e 's/[()]//g' "postProcessing/interfaceHeight1/1.36/position.dat" > "positionClean.dat"
+
+OF="positionClean.dat"
+EXPT="data.dat"
+
+createGraphs $OF $EXPT

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/blockMeshDict

@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    (0 0 -0.0)
+    (10 0 -0.0)
+    (10 0.2 -0.0)
+    (0 0.2 -0.0)
+    (0 0 0.02)
+    (10 0 0.02)
+    (10 0.2 0.02)
+    (0 0.2 0.02)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (400 2 1)
+    simpleGrading (1 1 1)
+);
+
+patches
+(
+   /* wall sideWalls
+    (
+        (1 5 4 0)
+        (3 7 6 2)
+    )*/
+    wall faceWall
+    (
+      (0 4 7 3)
+    )
+
+    patch outlet
+    (
+      (2 6 5 1)
+    )
+    empty frontAndBack
+    (
+        (0 3 2 1)
+        (4 5 6 7)
+        (1 5 4 0)
+        (3 7 6 2)
+
+    )
+);
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/controlDict

@@ -0,0 +1,68 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoReactingMultiphaseInterFoam;
+
+startFrom       startTime;
+
+startTime       1.36;
+
+stopAt          endTime;
+
+endTime         40;
+
+deltaT          1e-5;
+
+writeControl    adjustableRunTime;
+
+writeInterval   4;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  12;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   10;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.01;
+maxAlphaCo      0.01;
+maxAlphaDdt     0.01;
+
+functions
+{
+    interfaceHeight1
+    {
+        type           interfaceHeight;
+        libs           ("libfieldFunctionObjects.so");
+        alpha          alpha.liquid;
+        locations      ((0 0 0));
+        direction       (1 0 0);
+        writeControl    timeStep;
+        writeInterval   3;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/fvSchemes

@@ -0,0 +1,66 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(rhoPhi,U)           Gauss linearUpwind grad(U);
+    "div\(phi,alpha.*\)"    Gauss vanLeer;
+    "div\(phir,alpha.*\)"   Gauss linear;
+
+    "div\(phi,.*\.gas.*\)"  Gauss vanLeer;
+
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+
+    div(phi,T)              Gauss vanLeer;
+
+    div(rhoPhi,epsilon)     Gauss upwind;
+    div(rhoPhi,k)           Gauss upwind;
+}
+
+laplacianSchemes
+{
+    default        Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+
+fluxRequired
+{
+    default         no;
+    p_rgh;
+    "alpha.*";
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/fvSolution

@@ -0,0 +1,111 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 1;
+        cAlphas          ((liquid and gas) 0.1);
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    rho
+    {
+        solver           diagonal;
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-7;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        solver           PCG;
+        preconditioner   DIC;
+        tolerance        1e-9;
+        relTol           0.001;
+        smoother         DIC;
+    }
+
+    mDotSmearFinal
+    {
+        solver           PCG;
+        tolerance        1e-6;
+        preconditioner   DIC;
+        relTol           0.00;
+        smoother         DIC;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance        1e-9;
+        relTol           0;
+        minIter 10;
+    }
+
+    pcorrFinal
+    {
+        $p_rgh;
+        tolerance        1e-9;
+        relTol           0;
+    }
+
+
+    "(U|h|T.*|k|epsilon|R)"
+    {
+        solver           smoothSolver;
+        smoother         symGaussSeidel;
+        tolerance        1e-8;
+        relTol           0.01;
+    }
+
+    "(U|h|T.*|k|epsilon|R)Final"
+    {
+        $U;
+        tolerance        1e-8;
+        relTol           0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   yes;
+    nCorrectors         4;
+    nOuterCorrectors    1;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"           1;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/setAlphaFieldDict

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system/fluid";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+field "alpha.gas";
+type plane;
+origin (0.503e-3 0 0);
+direction (1 0 0);
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/setAlphaFieldDict.gas

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system/fluid";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+field "alpha.gas";
+type plane;
+origin (0.503e-3 0 0);
+direction (1 0 0);
+
+// ************************************************************************* //

tutorials/verificationAndValidation/icoReactingMultiphaseInterFoam/stefanProblem/system/setAlphaFieldDict.liquid

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system/fluid";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+field "alpha.liquid";
+type plane;
+origin (0.503e-3 0 0);
+direction (-1 0 0);
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/0.orig/T

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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 373.15;
+
+boundaryField
+{
+    sideWalls
+    {
+        type            zeroGradient;
+    }
+
+    faceWall
+    {
+        type            fixedValue;
+        value           uniform 378.15;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/0.orig/U

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    sideWalls
+    {
+        type            slip;
+        value           uniform (0 0 0);
+    }
+
+    faceWall
+    {
+        type            noSlip;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/0.orig/alpha.liquid

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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 0;
+
+boundaryField
+{
+    ".*"
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/0.orig/p

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+  sideWalls
+  {
+      type            fixedFluxPressure;
+  }
+
+  faceWall
+  {
+      type            fixedFluxPressure;
+  }
+
+  outlet
+  {
+      //type            fixedValue;
+      type            totalPressure;
+      rho             rho;
+      p0              uniform 1e5;
+      value           uniform 1e5;
+  }
+
+  frontAndBack
+  {
+      type            empty;
+  }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/0.orig/p_rgh

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      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
+{
+  sideWalls
+  {
+      type            fixedFluxPressure;
+  }
+
+  faceWall
+  {
+      type            fixedFluxPressure;
+  }
+
+  outlet
+  {
+      //type            fixedValue;
+      type            totalPressure;
+      rho             rho;
+      p0              uniform 1e5;
+      value           uniform 1e5;
+  }
+
+  frontAndBack
+  {
+      type            empty;
+  }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/Allclean

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+rm -rf 0
+rm positionClean.dat
+rm OF_vs_Exact.eps
+
+#------------------------------------------------------------------------------

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/Allrun

@@ -0,0 +1,14 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+restore0Dir
+cp -r 0 1.36
+runApplication setAlphaField
+
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      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 0 0);
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/constant/phaseChangeProperties

@@ -0,0 +1,33 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      phaseChangeProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phaseChangeTwoPhaseModel constant;//interfaceHeatResistance;
+
+
+R               1e6;
+maxAlphaRate    1;
+spread          3;
+
+
+coeffC          0;
+coeffE          500;
+
+// interfacePhaseChangeCoeffs
+// {
+//     R_               1e6;
+// }
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/constant/thermophysicalProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+TSat             TSat       [0 0 0 1 0]    373.15;   // saturation temperature
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/constant/transportProperties

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1906                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (liquid vapour); // FC-72
+
+liquid
+{
+    transportModel  Newtonian;
+    nu              1e-6;
+    rho             958.4;
+
+
+    Cp              4216;  // irrelevant
+    Cv              4216;  // irrelevant
+    kappa           0.671; // irrelevant
+    hf              0;
+}
+
+vapour
+{
+    transportModel  Newtonian;
+    nu              1e-5; // irrelevant
+    rho             0.581;
+
+
+    Cp              2030;      // FC72 vapour
+    Cv              2030;      // Cv = Cp - R/w
+    kappa           0.025;       // FC72 vapour // 0.01;
+    hf              2260.0e3;
+}
+
+sigma 0;
+
+Prt             0.7;
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/constant/turbulenceProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/plotGraph

@@ -0,0 +1,25 @@
+createGraphs()
+{
+    OF=$1
+    EXPT=$2
+
+    gnuplot<<EOF
+    set terminal postscript default
+    set output "OF_vs_Exact.eps"
+    set xlabel "t [sec]"
+    set ylabel "x [mm]"
+    set grid
+
+
+    plot \
+        "$EXPT" u 1:2 title "Exact", \
+        "$OF" u 1:2 title "OpenFOAM" with line lt -1 lw 1
+EOF
+}
+
+sed -e 's/[()]//g' "postProcessing/interfaceHeight1/1.36/position.dat" > "positionClean.dat"
+
+OF="positionClean.dat"
+EXPT="data.dat"
+
+createGraphs $OF $EXPT

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/blockMeshDict

@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    (0 0 -0.0)
+    (10 0 -0.0)
+    (10 0.2 -0.0)
+    (0 0.2 -0.0)
+    (0 0 0.02)
+    (10 0 0.02)
+    (10 0.2 0.02)
+    (0 0.2 0.02)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (400 2 1)
+    simpleGrading (1 1 1)
+);
+
+patches
+(
+   /* wall sideWalls
+    (
+        (1 5 4 0)
+        (3 7 6 2)
+    )*/
+    wall faceWall
+    (
+      (0 4 7 3)
+    )
+
+    patch outlet
+    (
+      (2 6 5 1)
+    )
+    empty frontAndBack
+    (
+        (0 3 2 1)
+        (4 5 6 7)
+        (1 5 4 0)
+        (3 7 6 2)
+
+    )
+);
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/controlDict

@@ -0,0 +1,68 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     interCondensatingEvaporatingFoam;
+
+startFrom       startTime;
+
+startTime       1.36;
+
+stopAt          endTime;
+
+endTime         5;
+
+deltaT          1e-5;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.5;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  12;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   10;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.01;
+maxAlphaCo      0.01;
+maxDeltaT       0.01;
+
+functions
+{
+    interfaceHeight1
+    {
+        type           interfaceHeight;
+        libs           ("libfieldFunctionObjects.so");
+        alpha          alpha.liquid;
+        locations      ((0 0 0));
+        direction       (1 0 0);
+        writeControl    timeStep;
+        writeInterval   3;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/fvSchemes

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss linear;
+    div(rhoPhi,U)      Gauss vanLeerV;
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+    div(rhoCpPhi,T) Gauss vanLeer;
+    
+    div((interpolate(cp)*rhoPhi),T) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default        Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+fluxRequired
+{
+    p_rgh;
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/fvSolution

@@ -0,0 +1,118 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    alpha.liquid
+    {
+        nAlphaCorr      2;
+        nAlphaSubCycles 2;
+        cAlpha          1;
+
+        MULESCorr       no;
+        nLimiterIter    5;
+
+        //solver          smoothSolver;
+        //smoother        symGaussSeidel;
+        //tolerance       1e-8;
+        //relTol          0;
+    }
+
+    rho
+    {
+        solver           diagonal;
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-7;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        //solver           GAMG;
+        tolerance        1e-9;
+        solver          PCG;
+        preconditioner  DIC;
+        relTol           0.001;
+        smoother         DIC;
+    }
+
+    mDotSmearFinal
+    {
+        solver           PCG;
+        tolerance        1e-6;
+        preconditioner   DIC;
+        relTol           0.00;
+        smoother         DIC;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance        1e-9;
+        relTol           0;
+        minIter 10;
+    }
+
+    pcorrFinal
+    {
+        $p_rgh;
+        tolerance        1e-9;
+        relTol           0;
+    }
+
+
+    "(U|h|T.*|k|epsilon|R)"
+    {
+        solver           smoothSolver; //PBiCGStab;
+        smoother        symGaussSeidel;
+        //preconditioner   DILU;
+        tolerance        1e-7;
+        relTol           0.;
+        minIter          15;
+        maxIter          50;
+    }
+
+    "(U|h|T.*|k|epsilon|R)Final"
+    {
+        $U;
+        tolerance        1e-7;
+        relTol           0;
+        maxIter          50;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   yes;
+    nCorrectors         4;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"           1;
+    }
+}
+
+// ************************************************************************* //

tutorials/verificationAndValidation/interCondensatingEvaporatingFoam/stefanProblem/system/setAlphaFieldDict

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system/fluid";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+field "alpha.liquid";
+type plane;
+origin (0.503e-3 0 0);
+direction (-1 0 0);
+
+// ************************************************************************* //