waxSolventEvaporation.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) 2017-2018 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 "waxSolventEvaporation.H"
28 #include "thermoSingleLayer.H"
29 #include "zeroField.H"
30 
31 #include "fvmDdt.H"
32 #include "fvmDiv.H"
33 #include "fvcDiv.H"
34 #include "fvmSup.H"
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 namespace regionModels
41 {
42 namespace surfaceFilmModels
43 {
44 
45 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 
47 defineTypeNameAndDebug(waxSolventEvaporation, 0);
48 
50 (
51  phaseChangeModel,
52  waxSolventEvaporation,
53  dictionary
54 );
55 
56 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
57 
59 (
60  const scalar Re,
61  const scalar Sc
62 ) const
63 {
64  if (Re < 5.0E+05)
65  {
66  return 0.664*sqrt(Re)*cbrt(Sc);
67  }
68  else
69  {
70  return 0.037*pow(Re, 0.8)*cbrt(Sc);
71  }
72 }
73 
74 
75 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
76 
77 waxSolventEvaporation::waxSolventEvaporation
78 (
80  const dictionary& dict
81 )
82 :
83  phaseChangeModel(typeName, film, dict),
84  Wwax_
85  (
86  IOobject
87  (
88  typeName + ":Wwax",
89  film.regionMesh().time().constant(),
90  film.regionMesh()
91  ),
92  readScalar(coeffDict_.lookup("Wwax"))
93  ),
94  Wsolvent_
95  (
96  IOobject
97  (
98  typeName + ":Wsolvent",
99  film.regionMesh().time().constant(),
100  film.regionMesh()
101  ),
102  readScalar(coeffDict_.lookup("Wsolvent"))
103  ),
104  Ysolvent0_
105  (
106  IOobject
107  (
108  typeName + ":Ysolvent0",
109  film.regionMesh().time().constant(),
110  film.regionMesh(),
113  )
114  ),
115  Ysolvent_
116  (
117  IOobject
118  (
119  typeName + ":Ysolvent",
120  film.regionMesh().time().timeName(),
121  film.regionMesh(),
124  ),
125  film.regionMesh()
126  ),
127  deltaMin_(readScalar(coeffDict_.lookup("deltaMin"))),
128  L_(readScalar(coeffDict_.lookup("L"))),
129  TbFactor_(coeffDict_.lookupOrDefault<scalar>("TbFactor", 1.1)),
130  YInfZero_(coeffDict_.lookupOrDefault<Switch>("YInfZero", false)),
131  activityCoeff_
132  (
133  Function1<scalar>::New("activityCoeff", coeffDict_)
134  )
135 {}
136 
137 
138 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
139 
141 {}
142 
143 
144 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
145 
146 template<class YInfType>
148 (
149  const scalar dt,
150  scalarField& availableMass,
151  scalarField& dMass,
152  scalarField& dEnergy,
153  YInfType YInf
154 )
155 {
156  const thermoSingleLayer& film = filmType<thermoSingleLayer>();
157 
158  const volScalarField& delta = film.delta();
159  const volScalarField& deltaRho = film.deltaRho();
160  const surfaceScalarField& phi = film.phi();
161 
162  // Set local thermo properties
163  const SLGThermo& thermo = film.thermo();
164  const filmThermoModel& filmThermo = film.filmThermo();
165  const label vapId = thermo.carrierId(filmThermo.name());
166 
167  // Retrieve fields from film model
168  const scalarField& pInf = film.pPrimary();
169  const scalarField& T = film.T();
170  const scalarField& hs = film.hs();
171  const scalarField& rho = film.rho();
172  const scalarField& rhoInf = film.rhoPrimary();
173  const scalarField& muInf = film.muPrimary();
174  const scalarField& magSf = film.magSf();
175  const vectorField dU(film.UPrimary() - film.Us());
176  const scalarField limMass
177  (
178  max(scalar(0), availableMass - deltaMin_*rho*magSf)
179  );
180 
181  // Molecular weight of vapour [kg/kmol]
182  const scalar Wvap = thermo.carrier().W(vapId);
183 
184  const scalar Wwax = Wwax_.value();
185  const scalar Wsolvent = Wsolvent_.value();
186 
187  volScalarField::Internal evapRateCoeff
188  (
189  IOobject
190  (
191  typeName + ":evapRateCoeff",
192  film.regionMesh().time().timeName(),
193  film.regionMesh(),
196  false
197  ),
198  film.regionMesh(),
200  );
201 
202  volScalarField::Internal evapRateInf
203  (
204  IOobject
205  (
206  typeName + ":evapRateInf",
207  film.regionMesh().time().timeName(),
208  film.regionMesh(),
211  false
212  ),
213  film.regionMesh(),
214  dimensionedScalar("zero", dimDensity*dimVelocity, 0)
215  );
216 
217  bool filmPresent = false;
218 
219  forAll(dMass, celli)
220  {
221  if (delta[celli] > deltaMin_)
222  {
223  filmPresent = true;
224 
225  const scalar Ysolvent = Ysolvent_[celli];
226 
227  // Molefraction of solvent in liquid film
228  const scalar Xsolvent
229  (
230  Ysolvent*Wsolvent/((1 - Ysolvent)*Wwax + Ysolvent*Wsolvent)
231  );
232 
233  // Primary region density [kg/m3]
234  const scalar rhoInfc = rhoInf[celli];
235 
236  // Cell pressure [Pa]
237  const scalar pc = pInf[celli];
238 
239  // Calculate the boiling temperature
240  const scalar Tb = filmThermo.Tb(pc);
241 
242  // Local temperature - impose lower limit of 200 K for stability
243  const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[celli]));
244 
245  const scalar pPartialCoeff
246  (
247  filmThermo.pv(pc, Tloc)*activityCoeff_->value(Xsolvent)
248  );
249 
250  scalar XsCoeff = pPartialCoeff/pc;
251 
252  // Vapour phase mole fraction of solvent at interface
253  scalar Xs = XsCoeff*Xsolvent;
254 
255  if (Xs > 1)
256  {
258  << "Solvent vapour pressure > ambient pressure"
259  << endl;
260 
261  XsCoeff /= Xs;
262  Xs = 1;
263  }
264 
265  // Vapour phase mass fraction of solvent at the interface
266  const scalar YsCoeff
267  (
268  XsCoeff/(XsCoeff*Xsolvent*Wsolvent + (1 - Xs)*Wvap)
269  );
270 
271  // Primary region viscosity [Pa.s]
272  const scalar muInfc = muInf[celli];
273 
274  // Reynolds number
275  const scalar Re = rhoInfc*mag(dU[celli])*L_/muInfc;
276 
277  // Vapour diffusivity [m2/s]
278  const scalar Dab = filmThermo.D(pc, Tloc);
279 
280  // Schmidt number
281  const scalar Sc = muInfc/(rhoInfc*(Dab + rootVSmall));
282 
283  // Sherwood number
284  const scalar Sh = this->Sh(Re, Sc);
285 
286  // Mass transfer coefficient [m/s]
287  evapRateCoeff[celli] = rhoInfc*Sh*Dab/(L_ + rootVSmall);
288 
289  // Solvent mass transfer
290  const scalar dm
291  (
292  max
293  (
294  dt*magSf[celli]
295  *evapRateCoeff[celli]*(YsCoeff*Ysolvent - YInf[celli]),
296  0
297  )
298  );
299 
300  if (dm > limMass[celli])
301  {
302  evapRateCoeff[celli] *= limMass[celli]/dm;
303  }
304 
305  evapRateInf[celli] = evapRateCoeff[celli]*YInf[celli];
306  evapRateCoeff[celli] *= YsCoeff;
307 
308  // hVap[celli] = filmThermo.hl(pc, Tloc);
309  }
310  }
311 
312  const dimensionedScalar deltaRho0Bydt
313  (
314  "deltaRho0",
315  deltaRho.dimensions()/dimTime,
316  rootVSmall/dt
317  );
318 
319  volScalarField::Internal impingementRate
320  (
321  max
322  (
323  -film.rhoSp()(),
324  dimensionedScalar("zero", film.rhoSp().dimensions(), 0)
325  )
326  );
327 
328  if (filmPresent)
329  {
330  // Solve for the solvent mass fraction
331  fvScalarMatrix YsolventEqn
332  (
333  fvm::ddt(deltaRho, Ysolvent_)
334  + fvm::div(phi, Ysolvent_)
335  ==
336  deltaRho0Bydt*Ysolvent_()
337 
338  + evapRateInf
339 
340  // Include the effect of the impinging droplets
341  // added with Ysolvent = Ysolvent0
342  + impingementRate*Ysolvent0_
343 
344  - fvm::Sp
345  (
346  deltaRho0Bydt
347  + evapRateCoeff
348  + film.rhoSp()()
349  + impingementRate,
350  Ysolvent_
351  )
352  );
353 
354  YsolventEqn.relax();
355  YsolventEqn.solve();
356 
357  Ysolvent_.min(1);
358  Ysolvent_.max(0);
359 
360  scalarField dm
361  (
362  dt*magSf*rhoInf*(evapRateCoeff*Ysolvent_ + evapRateInf)
363  );
364 
365  dMass += dm;
366 
367  // Heat is assumed to be removed by heat-transfer to the wall
368  // so the energy remains unchanged by the phase-change.
369  dEnergy += dm*hs;
370 
371  // Latent heat [J/kg]
372  // dEnergy += dm*(hs[celli] + hVap);
373  }
374 }
375 
376 
378 (
379  const scalar dt,
380  scalarField& availableMass,
381  scalarField& dMass,
382  scalarField& dEnergy
383 )
384 {
385  if (YInfZero_)
386  {
387  correctModel(dt, availableMass, dMass, dEnergy, zeroField());
388  }
389  else
390  {
391  const thermoSingleLayer& film = filmType<thermoSingleLayer>();
392  const label vapId = film.thermo().carrierId(film.filmThermo().name());
393  const scalarField& YInf = film.YPrimary()[vapId];
394 
395  correctModel(dt, availableMass, dMass, dEnergy, YInf);
396  }
397 }
398 
399 
400 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
401 
402 } // End namespace surfaceFilmModels
403 } // End namespace regionModels
404 } // End namespace Foam
405 
406 // ************************************************************************* //
scalar delta
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
virtual const word & name() const =0
Return the specie name.
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
virtual const volScalarField & deltaRho() const
Return the film thickness*density (helper field) [kg/m3].
surfaceScalarField & phi
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 > &)
autoPtr< Function1< scalar > > activityCoeff_
Activity coefficient as a function of solvent mole fraction.
uniformDimensionedScalarField Wwax_
Molecular weight of wax [kg/kmol].
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
virtual const volVectorField & Us() const
Return the film surface velocity [m/s].
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
const volScalarField & rhoPrimary() const
Density [kg/m3].
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
rhoReactionThermo & thermo
Definition: createFields.H:28
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
const scalar deltaMin_
Minimum film height for model to be active.
addToRunTimeSelectionTable(surfaceFilmRegionModel, kinematicSingleLayer, mesh)
Macros for easy insertion into run-time selection tables.
Base class for surface film phase change models.
const volScalarField & muPrimary() const
Viscosity [Pa.s].
const filmThermoModel & filmThermo() const
Film thermo.
const dimensionSet & dimensions() const
Return dimensions.
const volVectorField & UPrimary() const
Velocity [m/s].
A class representing the concept of a field of 0 used to avoid unnecessary manipulations for objects ...
Definition: zeroField.H:50
dimensionedScalar cbrt(const dimensionedScalar &ds)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calculate the matrix for the first temporal derivative.
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
const Type & value() const
Return const reference to value.
virtual const surfaceScalarField & phi() const
Return the film flux [kg.m/s].
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:61
const PtrList< volScalarField > & YPrimary() const
Specie mass fractions [0-1].
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
void min(const dimensioned< Type > &)
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 successful.
Definition: doubleScalar.H:68
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].
Calculate the divergence of the given field.
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
uniformDimensionedScalarField Wsolvent_
Molecular weight of liquid [kg/kmol].
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:523
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 > &)
const dimensionSet dimDensity
const basicSpecieMixture & carrier() const
Return reference to the gaseous components.
Definition: SLGThermo.C:108
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Calculate the matrix for the divergence of the given field and flux.
label carrierId(const word &cmptName, bool allowNotFound=false) const
Index of carrier component.
Definition: SLGThermo.C:148
virtual const volScalarField & magSf() const
Return the face area magnitudes / [m2].
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
void max(const dimensioned< Type > &)
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
virtual const volScalarField & hs() const
Return the film sensible enthalpy [J/kg].
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar Sh(const scalar Re, const scalar Sc) const
Return Sherwood number as a function of Reynolds and Schmidt numbers.
virtual const volScalarField & rho() const
Return the film density [kg/m3].
virtual scalar W(const label speciei) const =0
Molecular weight of the given specie [kg/kmol].
uniformDimensionedScalarField Ysolvent0_
Initial solvent mass-fraction.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
defineTypeNameAndDebug(kinematicSingleLayer, 0)
const SLGThermo & thermo() const
Return const reference to the SLGThermo object.
void correctModel(const scalar dt, scalarField &availableMass, scalarField &dMass, scalarField &dEnergy, YInfType YInf)
Calculate the matrix for implicit and explicit sources.
virtual scalar pv(const scalar p, const scalar T) const =0
Return vapour pressure [Pa].
Thermodynamic form of single-cell layer surface film model.
Namespace for OpenFOAM.
virtual const volScalarField & T() const
Return the film mean temperature [K].
const dimensionSet dimVelocity