ENH: surfaceFilmModels::standardPhaseChange: Add support for treating YInf as zero

by setting the optional switch YInfZero to true.

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "standardPhaseChange.H"
 #include "addToRunTimeSelectionTable.H"
 #include "thermoSingleLayer.H"
+#include "zeroField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,7 +78,8 @@ standardPhaseChange::standardPhaseChange
     phaseChangeModel(typeName, film, dict),
     deltaMin_(readScalar(coeffDict_.lookup("deltaMin"))),
     L_(readScalar(coeffDict_.lookup("L"))),
-    TbFactor_(coeffDict_.lookupOrDefault<scalar>("TbFactor", 1.1))
+    TbFactor_(coeffDict_.lookupOrDefault<scalar>("TbFactor", 1.1)),
+    YInfZero_(coeffDict_.lookupOrDefault<Switch>("YInfZero", false))
 {}
 
 
@@ -89,26 +91,28 @@ standardPhaseChange::~standardPhaseChange()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class YInfType>
 void standardPhaseChange::correctModel
 (
     const scalar dt,
     scalarField& availableMass,
     scalarField& dMass,
-    scalarField& dEnergy
+    scalarField& dEnergy,
+    YInfType YInf
 )
 {
     const thermoSingleLayer& film = filmType<thermoSingleLayer>();
 
-    // set local thermo properties
+    // Set local thermo properties
     const SLGThermo& thermo = film.thermo();
     const filmThermoModel& filmThermo = film.filmThermo();
     const label vapId = thermo.carrierId(filmThermo.name());
 
-    // retrieve fields from film model
+    // Retrieve fields from film model
     const scalarField& delta = film.delta();
-    const scalarField& YInf = film.YPrimary()[vapId];
     const scalarField& pInf = film.pPrimary();
     const scalarField& T = film.T();
+    const scalarField& hs = film.hs();
     const scalarField& rho = film.rho();
     const scalarField& rhoInf = film.rhoPrimary();
     const scalarField& muInf = film.muPrimary();
@@ -116,36 +120,44 @@ void standardPhaseChange::correctModel
     const vectorField dU(film.UPrimary() - film.Us());
     const scalarField limMass
     (
-        max(scalar(0.0), availableMass - deltaMin_*rho*magSf)
+        max(scalar(0), availableMass - deltaMin_*rho*magSf)
     );
 
+    // Molecular weight of vapour [kg/kmol]
+    const scalar Wvap = thermo.carrier().W(vapId);
+
+    // Molecular weight of liquid [kg/kmol]
+    const scalar Wliq = filmThermo.W();
+
     forAll(dMass, celli)
     {
+        scalar dm = 0;
+
         if (delta[celli] > deltaMin_)
         {
-            // cell pressure [Pa]
+            // Cell pressure [Pa]
             const scalar pc = pInf[celli];
 
-            // calculate the boiling temperature
+            // Calculate the boiling temperature
             const scalar Tb = filmThermo.Tb(pc);
 
-            // local temperature - impose lower limit of 200 K for stability
+            // Local temperature - impose lower limit of 200 K for stability
             const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[celli]));
 
-            // saturation pressure [Pa]
+            // Saturation pressure [Pa]
             const scalar pSat = filmThermo.pv(pc, Tloc);
 
-            // latent heat [J/kg]
+            // Latent heat [J/kg]
             const scalar hVap = filmThermo.hl(pc, Tloc);
 
-            // calculate mass transfer
+            // Calculate mass transfer
             if (pSat >= 0.95*pc)
             {
-                // boiling
+                // Boiling
                 const scalar Cp = filmThermo.Cp(pc, Tloc);
                 const scalar Tcorr = max(0.0, T[celli] - Tb);
                 const scalar qCorr = limMass[celli]*Cp*(Tcorr);
-                dMass[celli] = qCorr/hVap;
+                dm = qCorr/hVap;
             }
             else
             {
@@ -158,39 +170,59 @@ void standardPhaseChange::correctModel
                 // Reynolds number
                 const scalar Re = rhoInfc*mag(dU[celli])*L_/muInfc;
 
-                // molecular weight of vapour [kg/kmol]
+                // Vapour mass fraction at interface
-                const scalar Wvap = thermo.carrier().W(vapId);
-
-                // molecular weight of liquid [kg/kmol]
-                const scalar Wliq = filmThermo.W();
-
-                // vapour mass fraction at interface
                 const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
 
-                // vapour diffusivity [m2/s]
+                // Vapour diffusivity [m2/s]
                 const scalar Dab = filmThermo.D(pc, Tloc);
 
                 // Schmidt number
-                const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));
+                const scalar Sc = muInfc/(rhoInfc*(Dab + rootVSmall));
 
                 // Sherwood number
                 const scalar Sh = this->Sh(Re, Sc);
 
-                // mass transfer coefficient [m/s]
+                // Mass transfer coefficient [m/s]
-                const scalar hm = Sh*Dab/(L_ + ROOTVSMALL);
+                const scalar hm = Sh*Dab/(L_ + rootVSmall);
 
-                // add mass contribution to source
+                // Add mass contribution to source
-                dMass[celli] =
+                dm = dt*magSf[celli]*rhoInfc*hm*(Ys - YInf[celli])/(1.0 - Ys);
-                    dt*magSf[celli]*rhoInfc*hm*(Ys - YInf[celli])/(1.0 - Ys);
             }
 
-            dMass[celli] = min(limMass[celli], max(0.0, dMass[celli]));
+            dMass[celli] += min(limMass[celli], max(dm, 0));
-            dEnergy[celli] = dMass[celli]*hVap;
+
+            // Heat is assumed to be removed by heat-transfer to the wall
+            // so the energy remains unchanged by the phase-change.
+            dEnergy[celli] += dm*hs[celli];
+            // dEnergy[celli] += dm*(hs[celli] + hVap);
         }
     }
 }
 
 
+void standardPhaseChange::correctModel
+(
+    const scalar dt,
+    scalarField& availableMass,
+    scalarField& dMass,
+    scalarField& dEnergy
+)
+{
+    if (YInfZero_)
+    {
+        correctModel(dt, availableMass, dMass, dEnergy, zeroField());
+    }
+    else
+    {
+        const thermoSingleLayer& film = filmType<thermoSingleLayer>();
+        const label vapId = film.thermo().carrierId(film.filmThermo().name());
+        const scalarField& YInf = film.YPrimary()[vapId];
+
+        correctModel(dt, availableMass, dMass, dEnergy, YInf);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace surfaceFilmModels

src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H

@@ -79,12 +79,25 @@ protected:
         //  Used to set max limit on temperature to Tb*TbFactor
         const scalar TbFactor_;
 
+        //- Switch to treat YInf as zero
+        Switch YInfZero_;
+
 
     // Protected member functions
 
         //- Return Sherwood number as a function of Reynolds and Schmidt numbers
         scalar Sh(const scalar Re, const scalar Sc) const;
 
+        template<class YInfType>
+        void correctModel
+        (
+            const scalar dt,
+            scalarField& availableMass,
+            scalarField& dMass,
+            scalarField& dEnergy,
+            YInfType YInf
+        );
+
 
 public: