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

 {
     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,37 +289,438 @@ 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();
+            }
+        }
+    }