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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "standardPhaseChange.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "thermoSingleLayer.H"
  29 #include "specie.H"
  30 #include "heatTransferModel.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36 namespace regionModels
  37 {
  38 namespace surfaceFilmModels
  39 {
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 defineTypeNameAndDebug(standardPhaseChange, 0);
  44
  45 addToRunTimeSelectionTable
  46 (
  47     phaseChangeModel,
  48     standardPhaseChange,
  49     dictionary
  50 );
  51
  52 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  53
  54 scalar standardPhaseChange::Sh
  55 (
  56     const scalar Re,
  57     const scalar Sc
  58 ) const
  59 {
  60     if (Re < 5.0E+05)
  61     {
  62         return 0.664*sqrt(Re)*cbrt(Sc);
  63     }
  64     else
  65     {
  66         return 0.037*pow(Re, 0.8)*cbrt(Sc);
  67     }
  68 }
  69
  70
  71 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  72
  73 standardPhaseChange::standardPhaseChange
  74 (
  75     const surfaceFilmModel& owner,
  76     const dictionary& dict
  77 )
  78 :
  79     phaseChangeModel(typeName, owner, dict),
  80     Tb_(readScalar(coeffs_.lookup("Tb"))),
  81     deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
  82     L_(readScalar(coeffs_.lookup("L"))),
  83     TbFactor_(coeffs_.lookupOrDefault<scalar>("TbFactor", 1.1))
  84 {}
  85
  86
  87 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  88
  89 standardPhaseChange::~standardPhaseChange()
  90 {}
  91
  92
  93 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  94
  95 void standardPhaseChange::correctModel
  96 (
  97     const scalar dt,
  98     scalarField& availableMass,
  99     scalarField& dMass,
 100     scalarField& dEnergy
 101 )
 102 {
 103     const thermoSingleLayer& film = refCast<const thermoSingleLayer>(owner_);
 104
 105     // set local thermo properties
 106     const SLGThermo& thermo = film.thermo();
 107     const label liqId = film.liquidId();
 108     const liquidProperties& liq = thermo.liquids().properties()[liqId];
 109     const label vapId = thermo.carrierId(thermo.liquids().components()[liqId]);
 110
 111     // retrieve fields from film model
 112     const scalarField& delta = film.delta();
 113     const scalarField& YInf = film.YPrimary()[vapId];
 114     const scalarField& pInf = film.pPrimary();
 115     const scalarField& T = film.T();
 116     const scalarField& Tw = film.Tw();
 117     const scalarField& rho = film.rho();
 118     const scalarField& TInf = film.TPrimary();
 119     const scalarField& rhoInf = film.rhoPrimary();
 120     const scalarField& muInf = film.muPrimary();
 121     const scalarField& magSf = film.magSf();
 122     const scalarField hInf(film.htcs().h());
 123     const scalarField hFilm(film.htcw().h());
 124     const vectorField dU(film.UPrimary() - film.Us());
 125     const scalarField limMass
 126     (
 127         max(scalar(0.0), availableMass - deltaMin_*rho*magSf)
 128     );
 129
 130     forAll(dMass, cellI)
 131     {
 132         if (delta[cellI] > deltaMin_)
 133         {
 134             // cell pressure [Pa]
 135             const scalar pc = pInf[cellI];
 136
 137             // local temperature - impose lower limit of 200 K for stability
 138             const scalar Tloc = min(TbFactor_*Tb_, max(200.0, T[cellI]));
 139
 140             // saturation pressure [Pa]
 141             const scalar pSat = liq.pv(pc, Tloc);
 142
 143             // latent heat [J/kg]
 144             const scalar hVap = liq.hl(pc, Tloc);
 145
 146             // calculate mass transfer
 147             if (pSat >= 0.95*pc)
 148             {
 149                 // boiling
 150                 const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]);
 151                 const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]);
 152
 153                 const scalar Cp = liq.Cp(pc, Tloc);
 154                 const scalar Tcorr = max(0.0, T[cellI] - Tb_);
 155                 const scalar qCorr = limMass[cellI]*Cp*(Tcorr);
 156                 dMass[cellI] =
 157                     dt*magSf[cellI]/hVap*(qDotInf + qDotFilm)
 158                   + qCorr/hVap;
 159             }
 160             else
 161             {
 162                 // Primary region density [kg/m3]
 163                 const scalar rhoInfc = rhoInf[cellI];
 164
 165                 // Primary region viscosity [Pa.s]
 166                 const scalar muInfc = muInf[cellI];
 167
 168                 // Reynolds number
 169                 const scalar Re = rhoInfc*mag(dU[cellI])*L_/muInfc;
 170
 171                 // molecular weight of vapour [kg/kmol]
 172                 const scalar Wvap = thermo.carrier().W(vapId);
 173
 174                 // molecular weight of liquid [kg/kmol]
 175                 const scalar Wliq = liq.W();
 176
 177                 // vapour mass fraction at interface
 178                 const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
 179
 180                 // vapour diffusivity [m2/s]
 181                 const scalar Dab = liq.D(pc, Tloc);
 182
 183                 // Schmidt number
 184                 const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));
 185
 186                 // Sherwood number
 187                 const scalar Sh = this->Sh(Re, Sc);
 188
 189                 // mass transfer coefficient [m/s]
 190                 const scalar hm = Sh*Dab/(L_ + ROOTVSMALL);
 191
 192                 // add mass contribution to source
 193                 dMass[cellI] =
 194                     dt*magSf[cellI]*rhoInfc*hm*(Ys - YInf[cellI])/(1.0 - Ys);
 195             }
 196
 197             dMass[cellI] = min(limMass[cellI], max(0.0, dMass[cellI]));
 198             dEnergy[cellI] = dMass[cellI]*hVap;
 199         }
 200     }
 201 }
 202
 203
 204 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 205
 206 } // End namespace surfaceFilmModels
 207 } // End namespace regionModels
 208 } // End namespace Foam
 209
 210 // ************************************************************************* //