standardPhaseChange.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 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"
28 #include "thermoSingleLayer.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace regionModels
35 {
36 namespace surfaceFilmModels
37 {
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 defineTypeNameAndDebug(standardPhaseChange, 0);
42 
44 (
45  phaseChangeModel,
46  standardPhaseChange,
47  dictionary
48 );
49 
50 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
51 
53 (
54  const scalar Re,
55  const scalar Sc
56 ) const
57 {
58  if (Re < 5.0E+05)
59  {
60  return 0.664*sqrt(Re)*cbrt(Sc);
61  }
62  else
63  {
64  return 0.037*pow(Re, 0.8)*cbrt(Sc);
65  }
66 }
67 
68 
69 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
70 
71 standardPhaseChange::standardPhaseChange
72 (
73  surfaceFilmModel& owner,
74  const dictionary& dict
75 )
76 :
77  phaseChangeModel(typeName, owner, dict),
78  deltaMin_(readScalar(coeffDict_.lookup("deltaMin"))),
79  L_(readScalar(coeffDict_.lookup("L"))),
80  TbFactor_(coeffDict_.lookupOrDefault<scalar>("TbFactor", 1.1))
81 {}
82 
83 
84 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
85 
87 {}
88 
89 
90 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
91 
93 (
94  const scalar dt,
95  scalarField& availableMass,
96  scalarField& dMass,
97  scalarField& dEnergy
98 )
99 {
100  const thermoSingleLayer& film = filmType<thermoSingleLayer>();
101 
102  // set local thermo properties
103  const SLGThermo& thermo = film.thermo();
104  const filmThermoModel& filmThermo = film.filmThermo();
105  const label vapId = thermo.carrierId(filmThermo.name());
106 
107  // retrieve fields from film model
108  const scalarField& delta = film.delta();
109  const scalarField& YInf = film.YPrimary()[vapId];
110  const scalarField& pInf = film.pPrimary();
111  const scalarField& T = film.T();
112  const scalarField& rho = film.rho();
113  const scalarField& rhoInf = film.rhoPrimary();
114  const scalarField& muInf = film.muPrimary();
115  const scalarField& magSf = film.magSf();
116  const vectorField dU(film.UPrimary() - film.Us());
117  const scalarField limMass
118  (
119  max(scalar(0.0), availableMass - deltaMin_*rho*magSf)
120  );
121 
122  forAll(dMass, cellI)
123  {
124  if (delta[cellI] > deltaMin_)
125  {
126  // cell pressure [Pa]
127  const scalar pc = pInf[cellI];
128 
129  // calculate the boiling temperature
130  const scalar Tb = filmThermo.Tb(pc);
131 
132  // local temperature - impose lower limit of 200 K for stability
133  const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[cellI]));
134 
135  // saturation pressure [Pa]
136  const scalar pSat = filmThermo.pv(pc, Tloc);
137 
138  // latent heat [J/kg]
139  const scalar hVap = filmThermo.hl(pc, Tloc);
140 
141  // calculate mass transfer
142  if (pSat >= 0.95*pc)
143  {
144  // boiling
145  const scalar Cp = filmThermo.Cp(pc, Tloc);
146  const scalar Tcorr = max(0.0, T[cellI] - Tb);
147  const scalar qCorr = limMass[cellI]*Cp*(Tcorr);
148  dMass[cellI] = qCorr/hVap;
149  }
150  else
151  {
152  // Primary region density [kg/m3]
153  const scalar rhoInfc = rhoInf[cellI];
154 
155  // Primary region viscosity [Pa.s]
156  const scalar muInfc = muInf[cellI];
157 
158  // Reynolds number
159  const scalar Re = rhoInfc*mag(dU[cellI])*L_/muInfc;
160 
161  // molecular weight of vapour [kg/kmol]
162  const scalar Wvap = thermo.carrier().W(vapId);
163 
164  // molecular weight of liquid [kg/kmol]
165  const scalar Wliq = filmThermo.W();
166 
167  // vapour mass fraction at interface
168  const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
169 
170  // vapour diffusivity [m2/s]
171  const scalar Dab = filmThermo.D(pc, Tloc);
172 
173  // Schmidt number
174  const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));
175 
176  // Sherwood number
177  const scalar Sh = this->Sh(Re, Sc);
178 
179  // mass transfer coefficient [m/s]
180  const scalar hm = Sh*Dab/(L_ + ROOTVSMALL);
181 
182  // add mass contribution to source
183  dMass[cellI] =
184  dt*magSf[cellI]*rhoInfc*hm*(Ys - YInf[cellI])/(1.0 - Ys);
185  }
186 
187  dMass[cellI] = min(limMass[cellI], max(0.0, dMass[cellI]));
188  dEnergy[cellI] = dMass[cellI]*hVap;
189  }
190  }
191 }
192 
193 
194 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
195 
196 } // End namespace surfaceFilmModels
197 } // End namespace regionModels
198 } // End namespace Foam
199 
200 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual scalar pv(const scalar p, const scalar T) const =0
Return vapour pressure [Pa].
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual scalar Cp(const scalar p, const scalar T) const =0
Return specific heat capacity [J/kg/K].
virtual scalar D(const scalar p, const scalar T) const =0
Return diffusivity [m2/s].
const volScalarField & rhoPrimary() const
Density / [kg/m3].
const scalar TbFactor_
Boiling temperature factor / [].
const dimensionedScalar & pSat
Definition: createFields.H:41
const SLGThermo & thermo() const
Return const reference to the SLGThermo object.
defineTypeNameAndDebug(kinematicSingleLayer, 0)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
virtual const volScalarField & T() const
Return the film mean temperature [K].
virtual const volScalarField & rho() const
Return the film density [kg/m3].
const volScalarField & delta() const
Return const access to the film thickness / [m].
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedScalar cbrt(const dimensionedScalar &ds)
addToRunTimeSelectionTable(surfaceFilmModel, kinematicSingleLayer, mesh)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const PtrList< volScalarField > & YPrimary() const
Specie mass fractions / [0-1].
const volVectorField & UPrimary() const
Velocity / [m/s].
Namespace for OpenFOAM.
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: SLGThermo.H:62
virtual scalar W() const =0
Return molecular weight [kg/kmol].
const volScalarField & muPrimary() const
Viscosity / [Pa.s].
virtual scalar Tb(const scalar p) const =0
Return boiling temperature [K].
const filmThermoModel & filmThermo() const
Film thermo.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual const word & name() const =0
Return the specie name.
virtual const volScalarField & magSf() const
Return the face area magnitudes / [m2].
scalar delta
Macros for easy insertion into run-time selection tables.
scalar Sh(const scalar Re, const scalar Sc) const
Return Sherwood number as a function of Reynolds and Schmidt numbers.
virtual scalar hl(const scalar p, const scalar T) const =0
Return latent heat [J/kg].
const volScalarField & pPrimary() const
Pressure / [Pa].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual const volVectorField & Us() const
Return the film surface velocity [m/s].
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
virtual void correctModel(const scalar dt, scalarField &availableMass, scalarField &dMass, scalarField &dEnergy)
Correct.
const scalar deltaMin_
Minimum film height for model to be active.