LiquidEvaporationBoil.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "LiquidEvaporationBoil.H"
27 #include "specie.H"
28 #include "mathematicalConstants.H"
29 
30 using namespace Foam::constant::mathematical;
31 
32 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
33 
34 template<class CloudType>
36 (
37  const label celli
38 ) const
39 {
40  scalarField Xc(this->owner().composition().carrier().Y().size());
41 
42  forAll(Xc, i)
43  {
44  Xc[i] =
45  this->owner().composition().carrier().Y()[i][celli]
46  /this->owner().composition().carrier().Wi(i);
47  }
48 
49  return Xc/sum(Xc);
50 }
51 
52 
53 template<class CloudType>
55 (
56  const scalar Re,
57  const scalar Sc
58 ) const
59 {
60  return 2.0 + 0.6*Foam::sqrt(Re)*cbrt(Sc);
61 }
62 
63 
64 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
65 
66 template<class CloudType>
68 (
69  const dictionary& dict,
70  CloudType& owner
71 )
72 :
73  PhaseChangeModel<CloudType>(dict, owner, typeName),
74  liquids_(owner.thermo().liquids()),
75  activeLiquids_(this->coeffDict().lookup("activeLiquids")),
76  liqToCarrierMap_(activeLiquids_.size(), -1),
77  liqToLiqMap_(activeLiquids_.size(), -1)
78 {
79  if (activeLiquids_.size() == 0)
80  {
82  << "Evaporation model selected, but no active liquids defined"
83  << nl << endl;
84  }
85  else
86  {
87  Info<< "Participating liquid species:" << endl;
88 
89  // Determine mapping between liquid and carrier phase species
90  forAll(activeLiquids_, i)
91  {
92  Info<< " " << activeLiquids_[i] << endl;
93  liqToCarrierMap_[i] =
94  owner.composition().carrierId(activeLiquids_[i]);
95  }
96 
97  // Determine mapping between model active liquids and global liquids
98  const label idLiquid = owner.composition().idLiquid();
99  forAll(activeLiquids_, i)
100  {
101  liqToLiqMap_[i] =
102  owner.composition().localId(idLiquid, activeLiquids_[i]);
103  }
104  }
105 }
106 
107 
108 template<class CloudType>
110 (
112 )
113 :
115  liquids_(pcm.owner().thermo().liquids()),
116  activeLiquids_(pcm.activeLiquids_),
117  liqToCarrierMap_(pcm.liqToCarrierMap_),
118  liqToLiqMap_(pcm.liqToLiqMap_)
119 {}
120 
121 
122 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
123 
124 template<class CloudType>
126 {}
127 
128 
129 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
130 
131 template<class CloudType>
133 (
134  const scalar dt,
135  const label celli,
136  const scalar Re,
137  const scalar Pr,
138  const scalar d,
139  const scalar nu,
140  const scalar T,
141  const scalar Ts,
142  const scalar pc,
143  const scalar Tc,
144  const scalarField& X,
145  scalarField& dMassPC
146 ) const
147 {
148  // immediately evaporate mass that has reached critical condition
149  if ((liquids_.Tc(X) - T) < small)
150  {
151  if (debug)
152  {
154  << "Parcel reached critical conditions: "
155  << "evaporating all available mass" << endl;
156  }
157 
158  forAll(activeLiquids_, i)
159  {
160  const label lid = liqToLiqMap_[i];
161  dMassPC[lid] = great;
162  }
163 
164  return;
165  }
166 
167  // droplet surface pressure assumed to surface vapour pressure
168  scalar ps = liquids_.pv(pc, Ts, X);
169 
170  // vapour density at droplet surface [kg/m^3]
171  scalar rhos = ps*liquids_.W(X)/(RR*Ts);
172 
173  // construct carrier phase species volume fractions for cell, celli
174  const scalarField XcMix(calcXc(celli));
175 
176  // carrier thermo properties
177  scalar Hsc = 0.0;
178  scalar Hc = 0.0;
179  scalar Cpc = 0.0;
180  scalar kappac = 0.0;
181  forAll(this->owner().composition().carrier().Y(), i)
182  {
183  scalar Yc = this->owner().composition().carrier().Y()[i][celli];
184  Hc += Yc*this->owner().composition().carrier().Ha(i, pc, Tc);
185  Hsc += Yc*this->owner().composition().carrier().Ha(i, ps, Ts);
186  Cpc += Yc*this->owner().composition().carrier().Cp(i, ps, Ts);
187  kappac += Yc*this->owner().composition().carrier().kappa(i, ps, Ts);
188  }
189 
190  // calculate mass transfer of each specie in liquid
191  forAll(activeLiquids_, i)
192  {
193  const label gid = liqToCarrierMap_[i];
194  const label lid = liqToLiqMap_[i];
195 
196  // boiling temperature at cell pressure for liquid species lid [K]
197  const scalar TBoil = liquids_.properties()[lid].pvInvert(pc);
198 
199  // limit droplet temperature to boiling/critical temperature
200  const scalar Td = min(T, 0.999*TBoil);
201 
202  // saturation pressure for liquid species lid [Pa]
203  const scalar pSat = liquids_.properties()[lid].pv(pc, Td);
204 
205  // carrier phase concentration
206  const scalar Xc = XcMix[gid];
207 
208 
209  if (Xc*pc > pSat)
210  {
211  // saturated vapour - no phase change
212  }
213  else
214  {
215  // vapour diffusivity [m^2/s]
216  const scalar Dab = liquids_.properties()[lid].D(ps, Ts);
217 
218  // Schmidt number
219  const scalar Sc = nu/(Dab + rootVSmall);
220 
221  // Sherwood number
222  const scalar Sh = this->Sh(Re, Sc);
223 
224 
225  if (pSat > 0.999*pc)
226  {
227  // boiling
228 
229  const scalar deltaT = max(T - TBoil, 0.5);
230 
231  // vapour heat of formation
232  const scalar hv = liquids_.properties()[lid].hl(pc, Td);
233 
234  // empirical heat transfer coefficient W/m2/K
235  scalar alphaS = 0.0;
236  if (deltaT < 5.0)
237  {
238  alphaS = 760.0*pow(deltaT, 0.26);
239  }
240  else if (deltaT < 25.0)
241  {
242  alphaS = 27.0*pow(deltaT, 2.33);
243  }
244  else
245  {
246  alphaS = 13800.0*pow(deltaT, 0.39);
247  }
248 
249  // flash-boil vaporisation rate
250  const scalar Gf = alphaS*deltaT*pi*sqr(d)/hv;
251 
252  // model constants
253  // NOTE: using Sherwood number instead of Nusselt number
254  const scalar A = (Hc - Hsc)/hv;
255  const scalar B = pi*kappac/Cpc*d*Sh;
256 
257  scalar G = 0.0;
258  if (A > 0.0)
259  {
260  // heat transfer from the surroundings contributes
261  // to the vaporisation process
262  scalar Gr = 1e-5;
263 
264  for (label i=0; i<50; i++)
265  {
266  scalar GrDash = Gr;
267 
268  G = B/(1.0 + Gr)*log(1.0 + A*(1.0 + Gr));
269  Gr = Gf/G;
270 
271  if (mag(Gr - GrDash)/GrDash < 1e-3)
272  {
273  break;
274  }
275  }
276  }
277 
278  dMassPC[lid] += (G + Gf)*dt;
279  }
280  else
281  {
282  // evaporation
283 
284  // surface molar fraction - Raoult's Law
285  const scalar Xs = X[lid]*pSat/pc;
286 
287  // molar ratio
288  const scalar Xr = (Xs - Xc)/max(small, 1.0 - Xs);
289 
290  if (Xr > 0)
291  {
292  // mass transfer [kg]
293  dMassPC[lid] += pi*d*Sh*Dab*rhos*log(1.0 + Xr)*dt;
294  }
295  }
296  }
297  }
298 }
299 
300 
301 template<class CloudType>
303 (
304  const label idc,
305  const label idl,
306  const scalar p,
307  const scalar T
308 ) const
309 {
310  scalar dh = 0;
311 
312  scalar TDash = T;
313  if (liquids_.properties()[idl].pv(p, T) >= 0.999*p)
314  {
315  TDash = liquids_.properties()[idl].pvInvert(p);
316  }
317 
318  typedef PhaseChangeModel<CloudType> parent;
319  switch (parent::enthalpyTransfer_)
320  {
321  case (parent::etLatentHeat):
322  {
323  dh = liquids_.properties()[idl].hl(p, TDash);
324  break;
325  }
326  case (parent::etEnthalpyDifference):
327  {
328  scalar hc = this->owner().composition().carrier().Ha(idc, p, TDash);
329  scalar hp = liquids_.properties()[idl].Ha(p, TDash);
330 
331  dh = hc - hp;
332  break;
333  }
334  default:
335  {
337  << "Unknown enthalpyTransfer type" << abort(FatalError);
338  }
339  }
340 
341  return dh;
342 }
343 
344 
345 template<class CloudType>
347 (
348  const scalarField& X
349 ) const
350 {
351  return liquids_.Tpt(X);
352 }
353 
354 
355 template<class CloudType>
357 (
358  const scalar p,
359  const scalarField& X
360 ) const
361 {
362  return liquids_.pvInvert(p, X);
363 }
364 
365 
366 // ************************************************************************* //
const dictionary & properties() const
Return const access to the properties dictionary.
Definition: subModelBase.C:134
tmp< scalarField > calcXc(const label celli) const
Calculate the carrier phase component volume fractions at celli.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
List< word > activeLiquids_
List of active liquid names.
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
Templated phase change model class.
Definition: ReactingCloud.H:57
dimensionedScalar log(const dimensionedScalar &ds)
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
basicSpecieMixture & composition
dimensionedSymmTensor sqr(const dimensionedVector &dv)
List< label > liqToCarrierMap_
Mapping between liquid and carrier species.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const dimensionedScalar G
Newtonian constant of gravitation.
virtual scalar Tvap(const scalarField &X) const
Return vapourisation temperature.
List< label > liqToLiqMap_
Mapping between local and global liquid species.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const CloudType & owner() const
Return const access to the owner cloud.
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
virtual scalar dh(const label idc, const label idl, const scalar p, const scalar T) const
Return the enthalpy per unit mass.
LiquidEvaporationBoil(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
mathematical constants.
Liquid evaporation model.
dimensionedScalar cbrt(const dimensionedScalar &ds)
spatialTransform Xr(const vector &a, const scalar omega)
Rotational spatial transformation tensor about axis a by omega radians.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const char nl
Definition: Ostream.H:260
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual scalar TMax(const scalar p, const scalarField &X) const
Return maximum/limiting temperature.
#define WarningInFunction
Report a warning using Foam::Warning.
virtual void calculate(const scalar dt, const label celli, const scalar Re, const scalar Pr, const scalar d, const scalar nu, const scalar T, const scalar Ts, const scalar pc, const scalar Tc, const scalarField &X, scalarField &dMassPC) const
Update model.
PtrList< volScalarField > & Y
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:53
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
virtual ~LiquidEvaporationBoil()
Destructor.
scalar Sh() const
Sherwood number.