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-2017 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& film,
74  const dictionary& dict
75 )
76 :
77  phaseChangeModel(typeName, film, 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& hs = film.hs();
113  const scalarField& rho = film.rho();
114  const scalarField& rhoInf = film.rhoPrimary();
115  const scalarField& muInf = film.muPrimary();
116  const scalarField& magSf = film.magSf();
117  const vectorField dU(film.UPrimary() - film.Us());
118  const scalarField limMass
119  (
120  max(scalar(0), availableMass - deltaMin_*rho*magSf)
121  );
122 
123  forAll(dMass, celli)
124  {
125  scalar dm = 0;
126 
127  if (delta[celli] > deltaMin_)
128  {
129  // cell pressure [Pa]
130  const scalar pc = pInf[celli];
131 
132  // calculate the boiling temperature
133  const scalar Tb = filmThermo.Tb(pc);
134 
135  // local temperature - impose lower limit of 200 K for stability
136  const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[celli]));
137 
138  // saturation pressure [Pa]
139  const scalar pSat = filmThermo.pv(pc, Tloc);
140 
141  // latent heat [J/kg]
142  const scalar hVap = filmThermo.hl(pc, Tloc);
143 
144  // calculate mass transfer
145  if (pSat >= 0.95*pc)
146  {
147  // boiling
148  const scalar Cp = filmThermo.Cp(pc, Tloc);
149  const scalar Tcorr = max(0.0, T[celli] - Tb);
150  const scalar qCorr = limMass[celli]*Cp*(Tcorr);
151  dm = qCorr/hVap;
152  }
153  else
154  {
155  // Primary region density [kg/m3]
156  const scalar rhoInfc = rhoInf[celli];
157 
158  // Primary region viscosity [Pa.s]
159  const scalar muInfc = muInf[celli];
160 
161  // Reynolds number
162  const scalar Re = rhoInfc*mag(dU[celli])*L_/muInfc;
163 
164  // molecular weight of vapour [kg/kmol]
165  const scalar Wvap = thermo.carrier().W(vapId);
166 
167  // molecular weight of liquid [kg/kmol]
168  const scalar Wliq = filmThermo.W();
169 
170  // vapour mass fraction at interface
171  const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
172 
173  // vapour diffusivity [m2/s]
174  const scalar Dab = filmThermo.D(pc, Tloc);
175 
176  // Schmidt number
177  const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));
178 
179  // Sherwood number
180  const scalar Sh = this->Sh(Re, Sc);
181 
182  // mass transfer coefficient [m/s]
183  const scalar hm = Sh*Dab/(L_ + ROOTVSMALL);
184 
185  // add mass contribution to source
186  dm = dt*magSf[celli]*rhoInfc*hm*(Ys - YInf[celli])/(1.0 - Ys);
187  }
188 
189  dMass[celli] += min(limMass[celli], max(dm, 0));
190 
191  // Heat is assumed to be removed by heat-transfer to the wall
192  // so the energy remains unchanged by the phase-change.
193  dEnergy[celli] += dm*hs[celli];
194  // dEnergy[celli] += dm*(hs[celli] + hVap);
195  }
196  }
197 }
198 
199 
200 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
201 
202 } // End namespace surfaceFilmModels
203 } // End namespace regionModels
204 } // End namespace Foam
205 
206 // ************************************************************************* //
scalar delta
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
virtual const word & name() const =0
Return the specie name.
virtual void correctModel(const scalar dt, scalarField &availableMass, scalarField &dMass, scalarField &dEnergy)
Correct.
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual scalar Cp(const scalar p, const scalar T) const =0
Return specific heat capacity [J/kg/K].
virtual scalar hl(const scalar p, const scalar T) const =0
Return latent heat [J/kg].
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual const volVectorField & Us() const
Return the film surface velocity [m/s].
const volScalarField & rhoPrimary() const
Density [kg/m3].
Macros for easy insertion into run-time selection tables.
const volScalarField & muPrimary() const
Viscosity [Pa.s].
const filmThermoModel & filmThermo() const
Film thermo.
addToRunTimeSelectionTable(surfaceFilmModel, kinematicSingleLayer, mesh)
const volVectorField & UPrimary() const
Velocity [m/s].
scalar Sh(const scalar Re, const scalar Sc) const
Return Sherwood number as a function of Reynolds and Schmidt numbers.
const scalar TbFactor_
Boiling temperature factor / [].
dimensionedScalar cbrt(const dimensionedScalar &ds)
const PtrList< volScalarField > & YPrimary() const
Specie mass fractions [0-1].
Thermo package for (S)olids (L)iquids and (G)ases Takes reference to thermo package, and provides:
Definition: SLGThermo.H:62
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual scalar D(const scalar p, const scalar T) const =0
Return diffusivity [m2/s].
virtual scalar Tb(const scalar p) const =0
Return boiling temperature [K].
const volScalarField & pPrimary() const
Pressure [Pa].
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const volScalarField & delta() const
Return const access to the film thickness [m].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual scalar W() const =0
Return molecular weight [kg/kmol].
virtual const volScalarField & magSf() const
Return the face area magnitudes / [m2].
virtual const volScalarField & hs() const
Return the film sensible enthalpy [J/kg].
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual const volScalarField & rho() const
Return the film density [kg/m3].
const scalar deltaMin_
Minimum film height for model to be active.
defineTypeNameAndDebug(kinematicSingleLayer, 0)
const SLGThermo & thermo() const
Return const reference to the SLGThermo object.
virtual scalar pv(const scalar p, const scalar T) const =0
Return vapour pressure [Pa].
Namespace for OpenFOAM.
virtual const volScalarField & T() const
Return the film mean temperature [K].