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-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 "waxSolventEvaporation.H"
27 #include "thermoSingleLayer.H"
28 #include "liquidThermo.H"
29 #include "basicSpecieMixture.H"
30 
31 #include "fvmDdt.H"
32 #include "fvmDiv.H"
33 #include "fvcDiv.H"
34 #include "fvmSup.H"
35 
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 namespace regionModels
43 {
44 namespace surfaceFilmModels
45 {
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 defineTypeNameAndDebug(waxSolventEvaporation, 0);
50 
52 (
53  phaseChangeModel,
54  waxSolventEvaporation,
55  dictionary
56 );
57 
58 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
59 
61 (
62  const scalar Re,
63  const scalar Sc
64 ) const
65 {
66  if (Re < 5.0E+05)
67  {
68  return 0.664*sqrt(Re)*cbrt(Sc);
69  }
70  else
71  {
72  return 0.037*pow(Re, 0.8)*cbrt(Sc);
73  }
74 }
75 
76 
77 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
78 
80 (
82  const dictionary& dict
83 )
84 :
85  speciePhaseChange(typeName, film, dict),
86  Wwax_
87  (
88  IOobject
89  (
90  IOobject::modelName("Wwax", typeName),
91  film.regionMesh().time().constant(),
92  film.regionMesh()
93  ),
94  coeffDict_.lookup<scalar>("Wwax")
95  ),
96  Wsolvent_
97  (
98  IOobject
99  (
100  IOobject::modelName("Wsolvent", typeName),
101  film.regionMesh().time().constant(),
102  film.regionMesh()
103  ),
104  coeffDict_.lookup<scalar>("Wsolvent")
105  ),
106  Ysolvent0_
107  (
108  IOobject
109  (
110  IOobject::modelName("Ysolvent0", typeName),
111  film.regionMesh().time().constant(),
112  film.regionMesh(),
115  )
116  ),
117  Ysolvent_
118  (
119  IOobject
120  (
121  IOobject::modelName("Ysolvent", typeName),
122  film.regionMesh().time().timeName(),
123  film.regionMesh(),
126  ),
127  film.regionMesh()
128  ),
129  deltaMin_(coeffDict_.lookup<scalar>("deltaMin")),
130  L_(coeffDict_.lookup<scalar>("L")),
131  TbFactor_(coeffDict_.lookupOrDefault<scalar>("TbFactor", 1.1)),
132  YInfZero_(coeffDict_.lookupOrDefault<Switch>("YInfZero", false)),
133  activityCoeff_
134  (
135  Function1<scalar>::New("activityCoeff", coeffDict_)
136  )
137 {}
138 
139 
140 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
141 
143 {}
144 
145 
146 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
147 
148 template<class YInfType>
150 (
151  const scalar dt,
152  scalarField& availableMass,
153  scalarField& dMass,
154  scalarField& dEnergy,
155  YInfType YInf
156 )
157 {
158  const thermoSingleLayer& film = filmType<thermoSingleLayer>();
159 
160  const volScalarField& alpha = film.alpha();
161  const volScalarField& delta = film.delta();
162  const volScalarField& rho = film.rho();
163  const surfaceScalarField& phi = film.phi();
164 
165  // Set local liquidThermo properties
166  const liquidProperties& liquidThermo =
167  refCast<const heRhoThermopureMixtureliquidProperties>(film.thermo())
168  .cellThermoMixture(0).properties();
169 
170  const basicSpecieMixture& primarySpecieThermo =
171  refCast<const basicSpecieMixture>(film.primaryThermo());
172 
173  // Retrieve fields from film model
174  const scalarField& pInf = film.pPrimary();
175  const scalarField& T = film.thermo().T();
176  const scalarField& he = film.thermo().he();
177  const scalarField& rhoInf = film.rhoPrimary();
178  const scalarField& muInf = film.muPrimary();
179  const scalarField& V = film.regionMesh().V();
180  const scalarField& magSf = film.magSf();
181  const scalarField& VbyA = film.VbyA();
182  const vectorField dU(film.UPrimary() - film.Us());
183  const scalarField limMass
184  (
185  max(scalar(0), availableMass - deltaMin_*rho*magSf)
186  );
187 
188  // Molecular weight of vapour [kg/kmol]
189  const scalar Wvap = this->Wvap();
190 
191  const scalar Wwax = Wwax_.value();
192  const scalar Wsolvent = Wsolvent_.value();
193 
194  volScalarField::Internal evapRateCoeff
195  (
196  IOobject
197  (
198  IOobject::modelName("evapRateCoeff", typeName),
199  film.regionMesh().time().timeName(),
200  film.regionMesh(),
203  false
204  ),
205  film.regionMesh(),
207  );
208 
209  volScalarField::Internal evapRateInf
210  (
211  IOobject
212  (
213  IOobject::modelName("evapRateInf", typeName),
214  film.regionMesh().time().timeName(),
215  film.regionMesh(),
218  false
219  ),
220  film.regionMesh(),
221  dimensionedScalar(evapRateCoeff.dimensions(), 0)
222  );
223 
224  // Local surface temperature at which evaporation takes place
225  scalarField Tloc(dMass.size());
226 
227  bool filmPresent = false;
228 
229  forAll(dMass, celli)
230  {
231  if (delta[celli] > deltaMin_)
232  {
233  filmPresent = true;
234 
235  const scalar Ysolvent = Ysolvent_[celli];
236 
237  // Molefraction of solvent in liquid film
238  const scalar Xsolvent
239  (
240  Ysolvent*Wsolvent/((1 - Ysolvent)*Wwax + Ysolvent*Wsolvent)
241  );
242 
243  // Primary region density [kg/m^3]
244  const scalar rhoInfc = rhoInf[celli];
245 
246  // Cell pressure [Pa]
247  const scalar pc = pInf[celli];
248 
249  // Calculate the boiling temperature
250  const scalar Tb = liquidThermo.pvInvert(pc);
251 
252  // Local temperature - impose lower limit of 200 K for stability
253  Tloc[celli] = min(TbFactor_*Tb, max(200.0, T[celli]));
254 
255  const scalar pPartialCoeff
256  (
257  liquidThermo.pv(pc, Tloc[celli])*activityCoeff_->value(Xsolvent)
258  );
259 
260  scalar XsCoeff = pPartialCoeff/pc;
261 
262  // Vapour phase mole fraction of solvent at interface
263  scalar Xs = XsCoeff*Xsolvent;
264 
265  if (Xs > 1)
266  {
268  << "Solvent vapour pressure > ambient pressure"
269  << endl;
270 
271  XsCoeff /= Xs;
272  Xs = 1;
273  }
274 
275  // Vapour phase mass fraction of solvent at the interface
276  const scalar YsCoeff
277  (
278  XsCoeff/(XsCoeff*Xsolvent*Wsolvent + (1 - Xs)*Wvap)
279  );
280 
281  // Primary region viscosity [Pa.s]
282  const scalar muInfc = muInf[celli];
283 
284  // Reynolds number
285  const scalar Re = rhoInfc*mag(dU[celli])*L_/muInfc;
286 
287  // Vapour diffusivity [m^2/s]
288  const scalar Dab = liquidThermo.D(pc, Tloc[celli]);
289 
290  // Schmidt number
291  const scalar Sc = muInfc/(rhoInfc*(Dab + rootVSmall));
292 
293  // Sherwood number
294  const scalar Sh = this->Sh(Re, Sc);
295 
296  // Mass transfer coefficient [kg/m^3 s]
297  evapRateCoeff[celli] =
298  rhoInfc*Sh*Dab/(VbyA[celli]*(L_ + rootVSmall));
299 
300  // Solvent mass transfer
301  const scalar dm
302  (
303  max
304  (
305  dt*V[celli]
306  *evapRateCoeff[celli]*(YsCoeff*Ysolvent - YInf[celli]),
307  0
308  )
309  );
310 
311  if (dm > limMass[celli])
312  {
313  evapRateCoeff[celli] *= limMass[celli]/dm;
314  }
315 
316  evapRateInf[celli] = evapRateCoeff[celli]*YInf[celli];
317  evapRateCoeff[celli] *= YsCoeff;
318  }
319  }
320 
321  reduce(filmPresent, orOp<bool>());
322 
323  const dimensionedScalar rho0Bydt
324  (
325  "rho0Bydt",
327  rootVSmall/dt
328  );
329 
330  volScalarField::Internal impingementRate
331  (
332  max
333  (
334  -film.rhoSp(),
335  dimensionedScalar(film.rhoSp().dimensions(), 0)
336  )
337  );
338 
339  if (filmPresent)
340  {
341  // Solve for the solvent mass fraction
342  fvScalarMatrix YsolventEqn
343  (
344  fvm::ddt(alpha, rho, Ysolvent_) + fvm::div(phi, Ysolvent_)
345  - fvm::Sp(film.continuityErr(), Ysolvent_)
346  ==
347  rho0Bydt*Ysolvent_()
348 
349  + evapRateInf
350 
351  // Include the effect of the impinging droplets
352  // added with Ysolvent = Ysolvent0
353  + impingementRate*Ysolvent0_
354 
355  - fvm::Sp
356  (
357  rho0Bydt
358  + evapRateCoeff
359  + film.rhoSp()
360  + impingementRate,
361  Ysolvent_
362  )
363  );
364 
365  YsolventEqn.relax();
366  YsolventEqn.solve();
367 
368  Ysolvent_.min(1);
369  Ysolvent_.max(0);
370 
371  scalarField dm
372  (
373  dt*V*rhoInf*(evapRateCoeff*Ysolvent_ + evapRateInf)
374  );
375 
376  dMass += dm;
377 
378  // Assume that the vapour transferred to the primary region is
379  // already at temperature Tloc so that all heat required for
380  // the phase-change is provided by the film
381  dEnergy += dm*primarySpecieThermo.Hs(vapId(), pInf, Tloc);
382 
383  // Heat is assumed to be removed by heat-transfer to the wall
384  // so the energy remains unchanged by the phase-change.
385  dEnergy += dm*he;
386  // dEnergy += dm*(h[celli] + hVap);
387  }
388 }
389 
390 
392 (
393  const scalar dt,
394  scalarField& availableMass,
395  scalarField& dMass,
396  scalarField& dEnergy
397 )
398 {
399  if (YInfZero_)
400  {
401  correctModel(dt, availableMass, dMass, dEnergy, zeroField());
402  }
403  else
404  {
405  const thermoSingleLayer& film = filmType<thermoSingleLayer>();
406  const scalarField& YInf = film.YPrimary()[vapId()];
407 
408  correctModel(dt, availableMass, dMass, dEnergy, YInf);
409  }
410 }
411 
412 
413 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
414 
415 } // End namespace surfaceFilmModels
416 } // End namespace regionModels
417 } // End namespace Foam
418 
419 // ************************************************************************* //
scalar delta
Run-time selectable general function of one variable.
Definition: Function1.H:52
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const volScalarField & VbyA() const
Return the cell layer volume/area [m].
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 > &)
virtual scalar pvInvert(scalar p) const
Invert the vapour pressure relationship to retrieve the.
autoPtr< Function1< scalar > > activityCoeff_
Activity coefficient as a function of solvent mole fraction.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
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:251
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/any.
Definition: Switch.H:60
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
const volScalarField & rhoPrimary() const
Density [kg/m^3].
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
Specialisation of basicMixture for a mixture consisting of a number for molecular species...
tmp< volVectorField::Internal > Us() const
Return the film surface velocity [m/s].
const scalar deltaMin_
Minimum film height for model to be active.
volScalarField::Internal & rhoSp()
Mass [kg/m^2/s].
Macros for easy insertion into run-time selection tables.
virtual volScalarField & he()=0
Enthalpy/Internal energy [J/kg].
const volScalarField & muPrimary() const
Viscosity [Pa.s].
const dimensionSet dimTime
const dimensionSet & dimensions() const
Return dimensions.
const volVectorField & UPrimary() const
Velocity [m/s].
virtual Type value(const scalar x) const =0
Return value as a function of scalar x.
A class representing the concept of a field of 0 used to avoid unnecessary manipulations for objects ...
Definition: zeroField.H:50
const dimensionSet dimDensity
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:123
const Type & value() const
Return const reference to value.
const surfaceScalarField & phi() const
Return the film flux [kg m/s].
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:55
The thermophysical properties of a liquid.
const PtrList< volScalarField > & YPrimary() const
Specie mass fractions [0-1].
phi
Definition: correctPhi.H:3
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
virtual scalar pv(scalar p, scalar T) const =0
Vapour pressure [Pa].
const volScalarField::Internal & continuityErr() const
Return the current continuity error.
void min(const dimensioned< Type > &)
const volScalarField & pPrimary() const
Pressure [Pa].
Calculate the divergence of the given field.
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
thermo he()
uniformDimensionedScalarField Wsolvent_
Molecular weight of liquid [kg/kmol].
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:521
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 > &)
virtual scalar D(scalar p, scalar T) const =0
Vapour diffusivity [m^2/s].
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Calculate the matrix for the divergence of the given field and flux.
const volScalarField::Internal & magSf() const
Return the face area magnitudes [m^2].
#define WarningInFunction
Report a warning using Foam::Warning.
virtual const volScalarField & T() const =0
Temperature [K].
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
void max(const dimensioned< Type > &)
virtual scalar Hs(const label speciei, const scalar p, const scalar T) const =0
Sensible enthalpy [J/kg].
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
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.
const volScalarField & rho() const
Return the film density [kg/m^3].
static word modelName(Name name, const word &model)
Return the name of the object within the given model.
uniformDimensionedScalarField Ysolvent0_
Initial solvent mass-fraction.
Specie phase change model abstract base class.
waxSolventEvaporation(surfaceFilmRegionModel &film, const dictionary &dict)
Construct from surface film model.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
const fluidThermo & primaryThermo() const
Return const reference to the primary region thermo object.
defineTypeNameAndDebug(kinematicSingleLayer, 0)
void correctModel(const scalar dt, scalarField &availableMass, scalarField &dMass, scalarField &dEnergy, YInfType YInf)
Calculate the matrix for implicit and explicit sources.
addToRunTimeSelectionTable(surfaceFilmModel, noFilm, mesh)
Thermodynamic form of single-cell layer surface film model.
const volScalarField & alpha() const
Return const access to the film volume fraction [].
Namespace for OpenFOAM.