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