diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
index 128b529e11e0710eb6904539a64239d40eaf6f98..85c6d473a7a2c86f024d41bd06fdbf2c940bb421 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+++ b/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]
- const scalar Wvap = thermo.carrier().W(vapId);
-
- // molecular weight of liquid [kg/kmol]
- const scalar Wliq = filmThermo.W();
-
- // vapour mass fraction at interface
+ // 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]
- const scalar hm = Sh*Dab/(L_ + ROOTVSMALL);
+ // Mass transfer coefficient [m/s]
+ const scalar hm = Sh*Dab/(L_ + rootVSmall);
- // add mass contribution to source
- dMass[celli] =
- dt*magSf[celli]*rhoInfc*hm*(Ys - YInf[celli])/(1.0 - Ys);
+ // Add mass contribution to source
+ dm = dt*magSf[celli]*rhoInfc*hm*(Ys - YInf[celli])/(1.0 - Ys);
}
- dMass[celli] = min(limMass[celli], max(0.0, dMass[celli]));
- dEnergy[celli] = dMass[celli]*hVap;
+ dMass[celli] += min(limMass[celli], max(dm, 0));
+
+ // 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
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
index 55007cc6c6aa6a8f82e4970c9217cef52ec1b5e6..96d15ed899b92da7b5430d505b646b9a36a29d23 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
+++ b/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: