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-2023 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().WiValue(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
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();
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 typename CloudType::parcelType& p,
135  const typename CloudType::parcelType::trackingData& td,
136  const scalar dt,
137  const scalar Re,
138  const scalar Pr,
139  const scalar d,
140  const scalar nu,
141  const scalar T,
142  const scalar Ts,
143  const scalar pc,
144  const scalar Tc,
145  const scalarField& X,
146  scalarField& dMassPC
147 ) const
148 {
149  // immediately evaporate mass that has reached critical condition
150  if ((liquids_.Tc(X) - T) < small)
151  {
152  if (debug)
153  {
155  << "Parcel reached critical conditions: "
156  << "evaporating all available mass" << endl;
157  }
158 
159  forAll(activeLiquids_, i)
160  {
161  const label lid = liqToLiqMap_[i];
162  dMassPC[lid] = great;
163  }
164 
165  return;
166  }
167 
168  // droplet surface pressure assumed to surface vapour pressure
169  scalar ps = liquids_.pv(pc, Ts, X);
170 
171  // vapour density at droplet surface [kg/m^3]
172  scalar rhos = ps*liquids_.W(X)/(RR*Ts);
173 
174  // construct carrier phase species volume fractions for cell, celli
175  const scalarField XcMix(calcXc(p.cell()));
176 
177  // carrier thermo properties
178  scalar hsc = 0.0;
179  scalar hc = 0.0;
180  scalar Cpc = 0.0;
181  scalar kappac = 0.0;
182  forAll(this->owner().composition().carrier().Y(), i)
183  {
184  scalar Yc = this->owner().composition().carrier().Y()[i][p.cell()];
185  hc += Yc*this->owner().composition().carrier().hai(i, pc, Tc);
186  hsc += Yc*this->owner().composition().carrier().hai(i, ps, Ts);
187  Cpc += Yc*this->owner().composition().carrier().Cpi(i, ps, Ts);
188  kappac += Yc*this->owner().composition().carrier().kappai(i, ps, Ts);
189  }
190 
191  // calculate mass transfer of each specie in liquid
192  forAll(activeLiquids_, i)
193  {
194  const label gid = liqToCarrierMap_[i];
195  const label lid = liqToLiqMap_[i];
196 
197  // boiling temperature at cell pressure for liquid species lid [K]
198  const scalar TBoil = liquids_.properties()[lid].pvInvert(pc);
199 
200  // limit droplet temperature to boiling/critical temperature
201  const scalar Td = min(T, 0.999*TBoil);
202 
203  // saturation pressure for liquid species lid [Pa]
204  const scalar pSat = liquids_.properties()[lid].pv(pc, Td);
205 
206  // carrier phase concentration
207  const scalar Xc = XcMix[gid];
208 
209 
210  if (Xc*pc > pSat)
211  {
212  // saturated vapour - no phase change
213  }
214  else
215  {
216  // vapour diffusivity [m^2/s]
217  const scalar Dab = liquids_.properties()[lid].D(ps, Ts);
218 
219  // Schmidt number
220  const scalar Sc = nu/(Dab + rootVSmall);
221 
222  // Sherwood number
223  const scalar Sh = this->Sh(Re, Sc);
224 
225 
226  if (pSat > 0.999*pc)
227  {
228  // boiling
229 
230  const scalar deltaT = max(T - TBoil, 0.5);
231 
232  // vapour heat of formation
233  const scalar hv = liquids_.properties()[lid].hl(pc, Td);
234 
235  // empirical heat transfer coefficient W/m2/K
236  scalar alphaS = 0.0;
237  if (deltaT < 5.0)
238  {
239  alphaS = 760.0*pow(deltaT, 0.26);
240  }
241  else if (deltaT < 25.0)
242  {
243  alphaS = 27.0*pow(deltaT, 2.33);
244  }
245  else
246  {
247  alphaS = 13800.0*pow(deltaT, 0.39);
248  }
249 
250  // flash-boil vaporisation rate
251  const scalar Gf = alphaS*deltaT*pi*sqr(d)/hv;
252 
253  // model constants
254  // NOTE: using Sherwood number instead of Nusselt number
255  const scalar A = (hc - hsc)/hv;
256  const scalar B = pi*kappac/Cpc*d*Sh;
257 
258  scalar G = 0.0;
259  if (A > 0.0)
260  {
261  // heat transfer from the surroundings contributes
262  // to the vaporisation process
263  scalar Gr = 1e-5;
264 
265  for (label i=0; i<50; i++)
266  {
267  scalar GrDash = Gr;
268 
269  G = B/(1.0 + Gr)*log(1.0 + A*(1.0 + Gr));
270  Gr = Gf/G;
271 
272  if (mag(Gr - GrDash)/GrDash < 1e-3)
273  {
274  break;
275  }
276  }
277  }
278 
279  dMassPC[lid] += (G + Gf)*dt;
280  }
281  else
282  {
283  // evaporation
284 
285  // surface molar fraction - Raoult's Law
286  const scalar Xs = X[lid]*pSat/pc;
287 
288  // molar ratio
289  const scalar Xr = (Xs - Xc)/max(small, 1.0 - Xs);
290 
291  if (Xr > 0)
292  {
293  // mass transfer [kg]
294  dMassPC[lid] += pi*d*Sh*Dab*rhos*log(1.0 + Xr)*dt;
295  }
296  }
297  }
298  }
299 }
300 
301 
302 template<class CloudType>
304 (
305  const label idc,
306  const label idl,
307  const scalar p,
308  const scalar T
309 ) const
310 {
311  scalar dh = 0;
312 
313  scalar TDash = T;
314  if (liquids_.properties()[idl].pv(p, T) >= 0.999*p)
315  {
316  TDash = liquids_.properties()[idl].pvInvert(p);
317  }
318 
319  typedef PhaseChangeModel<CloudType> parent;
320  switch (parent::enthalpyTransfer_)
321  {
322  case (parent::etLatentHeat):
323  {
324  dh = liquids_.properties()[idl].hl(p, TDash);
325  break;
326  }
327  case (parent::etEnthalpyDifference):
328  {
329  scalar hc =
330  this->owner().composition().carrier().hai(idc, p, TDash);
331  scalar hp = liquids_.properties()[idl].ha(p, TDash);
332 
333  dh = hc - hp;
334  break;
335  }
336  default:
337  {
339  << "Unknown enthalpyTransfer type" << abort(FatalError);
340  }
341  }
342 
343  return dh;
344 }
345 
346 
347 template<class CloudType>
349 (
350  const scalarField& X
351 ) const
352 {
353  return liquids_.Tpt(X);
354 }
355 
356 
357 template<class CloudType>
359 (
360  const scalar p,
361  const scalarField& X
362 ) const
363 {
364  return liquids_.pvInvert(p, X);
365 }
366 
367 
368 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const CloudType & owner() const
Return const access to the owner cloud.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:225
Liquid evaporation model.
List< word > activeLiquids_
List of active liquid names.
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.
List< label > liqToLiqMap_
Mapping between local and global liquid species.
List< label > liqToCarrierMap_
Mapping between liquid and carrier species.
virtual scalar Tvap(const scalarField &X) const
Return vapourisation temperature.
virtual void calculate(const typename CloudType::parcelType &p, const typename CloudType::parcelType::trackingData &td, const scalar dt, 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.
tmp< scalarField > calcXc(const label celli) const
Calculate the carrier phase component volume fractions at celli.
virtual scalar TMax(const scalar p, const scalarField &X) const
Return maximum/limiting temperature.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Templated phase change model class.
scalar Sh() const
Sherwood number.
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
const dictionary & properties() const
Return const access to the properties dictionary.
Definition: subModelBase.C:134
A class for managing temporary objects.
Definition: tmp.H:55
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
#define WarningInFunction
Report a warning using Foam::Warning.
mathematical constants.
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
const dimensionedScalar G
Newtonian constant of gravitation.
const doubleScalar e
Definition: doubleScalar.H:106
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
spatialTransform Xr(const vector &a, const scalar omega)
Rotational spatial transformation tensor about axis a by omega radians.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
dimensionedScalar log(const dimensionedScalar &ds)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
error FatalError
dimensionedScalar cbrt(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:266
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
scalarField Re(const UList< complex > &cf)
Definition: complexFields.C:97
dictionary dict
volScalarField & p
PtrList< volScalarField > & Y
fluidMulticomponentThermo & thermo
Definition: createFields.H:31