LiquidEvapFuchsKnudsen.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2020-2021 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "LiquidEvapFuchsKnudsen.H"
29 
30 
31 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
32 
33 template<class CloudType>
35 (
36  const scalar massliq,
37  const scalar masssol,
38  scalar& Xliq,
39  scalar& Xsol
40 ) const
41 {
42  const scalar Yliq = massliq/(massliq + masssol);
43  const scalar Ysol = 1 - Yliq;
44  Xliq = Yliq/liquids_.properties()[liqToLiqMap_].W();
45  Xsol = Ysol/this->owner().thermo().solids().properties()[solToSolMap_].W();
46  Xliq /= (Xliq + Xsol);
47  Xsol = 1 - Xliq;
48 }
49 
50 
51 template<class CloudType>
53 (
54  const label celli
55 ) const
56 {
57  scalarField Xc(this->owner().thermo().carrier().Y().size());
58 
59  forAll(Xc, i)
60  {
61  Xc[i] =
62  this->owner().thermo().carrier().Y()[i][celli]
63  /this->owner().thermo().carrier().W(i);
64  }
65 
66  return Xc/sum(Xc);
67 }
68 
69 
70 template<class CloudType>
72 (
73  const scalar Re,
74  const scalar Sc
75 ) const
76 {
77  return cbrt(1 + Sc*Re)*max(1, pow(Re, 0.077));
78 }
79 
80 
81 template<class CloudType>
83 (
84  const scalar Xliq,
85  const scalar Xsol
86 ) const
87 {
88  switch (method_)
89  {
90  case pUNIFAC:
91  {
93  << "Activity coefficient UNIFAC is not implemented " << nl
94  << abort(FatalError);
95  break;
96  }
97  case pHoff:
98  {
99  const scalar ic = this->coeffDict().getScalar("ic");
100  return inv((1 + ic*Xsol/(Xliq + ROOTVSMALL)));
101  break;
102  }
103  }
104  return -1;
105 }
106 
107 
108 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
109 
110 template<class CloudType>
112 (
113  const dictionary& dict,
114  CloudType& owner
115 )
116 :
117  PhaseChangeModel<CloudType>(dict, owner, typeName),
118  method_(pHoff),
119  gamma_(this->coeffDict().getScalar("gamma")),
120  alpha_(this->coeffDict().getScalar("alpham")),
121  liquids_(owner.thermo().liquids()),
122  solution_(this->coeffDict().lookup("solution")),
123  liqToCarrierMap_(-1),
124  liqToLiqMap_(-1),
125  solToSolMap_(-1)
126 {
127  if (solution_.size() > 2)
128  {
130  << "Solution is not well defined. It should be (liquid solid)"
131  << nl << exit(FatalError);
132  }
133  else
134  {
135  Info<< "Participating liquid-solid species:" << endl;
136 
137  Info<< " " << solution_[0] << endl;
138  liqToCarrierMap_ =
139  owner.composition().carrierId(solution_[0]);
140 
141  // Determine mapping between model active liquids and global liquids
142  const label idLiquid = owner.composition().idLiquid();
143  liqToLiqMap_ =
144  owner.composition().localId(idLiquid, solution_[0]);
145 
146  // Mapping for the solid
147  const label idSolid = owner.composition().idSolid();
148 
149  solToSolMap_ =
150  owner.composition().localId(idSolid, solution_[1]);
151 
152  const word activityCoefficienType
153  (
154  this->coeffDict().getWord("activityCoefficient")
155  );
156 
157  if (activityCoefficienType == "Hoff")
158  {
159  method_ = pHoff;
160  }
161  else if (activityCoefficienType == "UNIFAC")
162  {
163  method_ = pUNIFAC;
164  }
165  else
166  {
168  << "activityCoefficient must be either 'Hoff' or 'UNIFAC'"
169  << nl << exit(FatalError);
170  }
171  }
172 }
173 
174 
175 template<class CloudType>
177 (
179 )
180 :
182  method_(pcm.method_),
183  gamma_(pcm.gamma_),
184  alpha_(pcm.alpha_),
185  liquids_(pcm.owner().thermo().liquids()),
186  solution_(pcm.solution_),
187  liqToCarrierMap_(pcm.liqToCarrierMap_),
188  liqToLiqMap_(pcm.liqToLiqMap_),
189  solToSolMap_(pcm.solToSolMap_)
190 {}
191 
192 
193 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
194 
195 template<class CloudType>
197 (
198  const scalar dt,
199  const label celli,
200  const scalar Re,
201  const scalar Pr,
202  const scalar d,
203  const scalar nu,
204  const scalar rho,
205  const scalar T,
206  const scalar Ts,
207  const scalar pc,
208  const scalar Tc,
209  const scalarField& X,
210  const scalarField& solMass,
211  const scalarField& liqMass,
212  scalarField& dMassPC
213 ) const
214 {
215  const scalar rhog = this->owner().thermo().thermo().rho()()[celli];
216 
217  const label gid = liqToCarrierMap_;
218  const label lid = liqToLiqMap_;
219  const label sid = solToSolMap_;
220 
221  const scalar W = liquids_.properties()[lid].W();
222 
223  const scalar YeInf = this->owner().thermo().carrier().Y()[gid][celli];
224 
225  const scalar sigma = liquids_.properties()[lid].sigma(pc, Ts);
226 
227  // Kelvin effect
228  const scalar Ke = exp(4*sigma*W/(RR*rho*d*T));
229 
230  // vapour diffusivity [m2/s]
231  const scalar Dab = liquids_.properties()[lid].D(pc, Ts);
232 
233  // saturation pressure for species i [pa]
234  const scalar pSat = liquids_.properties()[lid].pv(pc, T);
235 
236  scalar Xliq(0), Xsol(0);
237  calcXcSolution(liqMass[lid], solMass[sid], Xliq, Xsol);
238 
239  // Activity Coefficient (gammaE*Xe)
240  const scalar gamma = activityCoeff(Xliq, Xsol);
241 
242  // water concentration at surface
243  const scalar Rliq = RR/W;
244  const scalar YeSurf = max(gamma*Ke*pSat/(Rliq*T*rhog), 0);
245 
246  const scalar Kn = 2*gamma_/d;
247  const scalar Cm =
248  (1+Kn)/(1+ (4/(3*alpha_) + 0.377)*Kn + sqr(Kn)*4/(3*alpha_));
249 
250  // Schmidt number
251  const scalar Sc = nu/(Dab + ROOTVSMALL);
252 
253  // Sherwood number
254  const scalar Sherwood = Sh(Re, Sc);
255 
256  // mass flux density [kg/m2/s]
257  const scalar Ni = (rhog*Sherwood*Dab*Cm/d)*log((1 - YeInf)/(1 - YeSurf));
258 
259  // mass transfer [kg]
260  const scalar As = Foam::constant::mathematical::pi*d*d;
261  dMassPC[lid] += Ni*As*dt;
262 }
263 
264 
265 template<class CloudType>
267 (
268  const label idc,
269  const label idl,
270  const scalar p,
271  const scalar T
272 ) const
273 {
274  scalar dh = 0;
275 
276  typedef PhaseChangeModel<CloudType> parent;
277  switch (parent::enthalpyTransfer_)
278  {
279  case (parent::etLatentHeat):
280  {
281  dh = liquids_.properties()[idl].hl(p, T);
282  break;
283  }
284  case (parent::etEnthalpyDifference):
285  {
286  scalar hc = this->owner().composition().carrier().Ha(idc, p, T);
287  scalar hp = liquids_.properties()[idl].h(p, T);
288 
289  dh = hc - hp;
290  break;
291  }
292  default:
293  {
295  << "Unknown enthalpyTransfer type" << abort(FatalError);
296  }
297  }
298 
299  return dh;
300 }
301 
302 
303 template<class CloudType>
305 (
306  const scalarField& X
307 ) const
308 {
309  return Zero;
310 }
311 
312 
313 template<class CloudType>
315 (
316  const scalar p,
317  const scalarField& X
318 ) const
319 {
320  // If liquid is present calculates pvInter
321  if (sum(X) > SMALL)
322  {
323  return liquids_.pvInvert(p, X);
324  }
325  else
326  {
327  return GREAT;
328  }
329 }
330 
331 
332 // ************************************************************************* //
Foam::constant::thermodynamic::RR
const scalar RR
Universal gas constant: default in [J/(kmol K)].
Definition: thermodynamicConstants.C:46
Foam::PhaseChangeModel
Templated phase change model class.
Definition: ReactingCloud.H:61
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::LiquidEvapFuchsKnudsen
Liquid evaporation/condensation model for solution of liquid and solid.
Definition: LiquidEvapFuchsKnudsen.H:66
Foam::LiquidEvapFuchsKnudsen::method_
activityCoeffMethodType method_
Method used.
Definition: LiquidEvapFuchsKnudsen.H:87
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::LiquidEvapFuchsKnudsen::solToSolMap_
label solToSolMap_
Mapping between local and global solid species.
Definition: LiquidEvapFuchsKnudsen.H:108
Foam::LiquidEvapFuchsKnudsen::calculate
virtual void calculate(const scalar dt, const label celli, const scalar Re, const scalar Pr, const scalar d, const scalar nu, const scalar rho, const scalar T, const scalar Ts, const scalar pc, const scalar Tc, const scalarField &X, const scalarField &Xsol, const scalarField &liqMass, scalarField &dMassPC) const
Update model.
Definition: LiquidEvapFuchsKnudsen.C:197
thermo
psiReactionThermo & thermo
Definition: createFields.H:28
Foam::LiquidEvapFuchsKnudsen::LiquidEvapFuchsKnudsen
LiquidEvapFuchsKnudsen(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
Definition: LiquidEvapFuchsKnudsen.C:112
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
rho
rho
Definition: readInitialConditions.H:88
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::LiquidEvapFuchsKnudsen::gamma_
scalar gamma_
Mean gas free path.
Definition: LiquidEvapFuchsKnudsen.H:90
Foam::LiquidEvapFuchsKnudsen::calcXcSolution
void calcXcSolution(const scalar massliq, const scalar masssol, scalar &Xliq, scalar &Xsol) const
Calculate volumetric fractions of components in the solution.
Definition: LiquidEvapFuchsKnudsen.C:35
nu
volScalarField & nu
Definition: readMechanicalProperties.H:176
Foam::LiquidEvapFuchsKnudsen::dh
virtual scalar dh(const label idc, const label idl, const scalar p, const scalar T) const
Return the enthalpy per unit mass.
Definition: LiquidEvapFuchsKnudsen.C:267
Foam::LiquidEvapFuchsKnudsen::solution_
List< word > solution_
List of active liquid names i.e (liquidName solidName)
Definition: LiquidEvapFuchsKnudsen.H:99
Foam::Field< scalar >
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:73
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Pr
dimensionedScalar Pr("Pr", dimless, laminarTransport)
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::log
dimensionedScalar log(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:262
Foam::LiquidEvapFuchsKnudsen::calcXc
tmp< scalarField > calcXc(const label celli) const
Calculate the carrier phase component volume fractions at celli.
Definition: LiquidEvapFuchsKnudsen.C:53
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
T
const volScalarField & T
Definition: createFieldRefs.H:2
Y
PtrList< volScalarField > & Y
Definition: createFieldRefs.H:7
rhog
dimensionedScalar rhog("rhog", dimDensity, transportProperties)
Foam::PhaseChangeModel::Sh
scalar Sh() const
Sherwood number.
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::LiquidEvapFuchsKnudsen::activityCoeff
scalar activityCoeff(const scalar Xliq, const scalar Ysol) const
Return activity coefficient.
Definition: LiquidEvapFuchsKnudsen.C:83
Foam::LiquidEvapFuchsKnudsen::liqToLiqMap_
label liqToLiqMap_
Mapping between local and global liquid species.
Definition: LiquidEvapFuchsKnudsen.H:105
Foam::LiquidEvapFuchsKnudsen::TMax
virtual scalar TMax(const scalar p, const scalarField &X) const
Return maximum/limiting temperature.
Definition: LiquidEvapFuchsKnudsen.C:315
Foam::LiquidEvapFuchsKnudsen::liqToCarrierMap_
label liqToCarrierMap_
Mapping between liquid and carrier species.
Definition: LiquidEvapFuchsKnudsen.H:102
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::LiquidEvapFuchsKnudsen::Tvap
virtual scalar Tvap(const scalarField &X) const
Return vapourisation temperature.
Definition: LiquidEvapFuchsKnudsen.C:305
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
gamma
const scalar gamma
Definition: EEqn.H:9
Foam::Re
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
LiquidEvapFuchsKnudsen.H
Foam::LiquidEvapFuchsKnudsen::alpha_
scalar alpha_
The mass thermal accomodation.
Definition: LiquidEvapFuchsKnudsen.H:93
sigma
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155