liquidMixtureProperties.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-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 
27 #include "dictionary.H"
28 #include "specie.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 const Foam::scalar Foam::liquidMixtureProperties::TrMax = 0.999;
33 
34 
35 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
36 
38 (
39  const dictionary& dict
40 )
41 :
42  components_(),
43  properties_()
44 {
45  components_ = dict.toc();
46  properties_.setSize(components_.size());
47 
48  forAll(components_, i)
49  {
50  if (dict.isDict(components_[i]))
51  {
52  properties_.set
53  (
54  i,
55  liquidProperties::New(dict.subDict(components_[i]))
56  );
57  }
58  else
59  {
60  properties_.set
61  (
62  i,
63  liquidProperties::New(components_[i])
64  );
65  }
66  }
67 }
68 
69 
71 (
72  const liquidMixtureProperties& lm
73 )
74 :
75  components_(lm.components_),
76  properties_(lm.properties_.size())
77 {
78  forAll(properties_, i)
79  {
80  properties_.set(i, lm.properties_(i)->clone());
81  }
82 }
83 
84 
85 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
86 
89 (
90  const dictionary& thermophysicalProperties
91 )
92 {
94  (
95  new liquidMixtureProperties(thermophysicalProperties)
96  );
97 }
98 
99 
100 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
101 
102 Foam::scalar Foam::liquidMixtureProperties::Tc(const scalarField& X) const
103 {
104  scalar vTc = 0.0;
105  scalar vc = 0.0;
106 
107  forAll(properties_, i)
108  {
109  scalar x1 = X[i]*properties_[i].Vc();
110  vc += x1;
111  vTc += x1*properties_[i].Tc();
112  }
113 
114  return vTc/vc;
115 }
116 
117 
119 {
120  scalar Tpt = 0.0;
121 
122  forAll(properties_, i)
123  {
124  Tpt += X[i]*properties_[i].Tt();
125  }
126 
127  return Tpt;
128 }
129 
130 
132 (
133  const scalar p,
134  const scalarField& X
135 ) const
136 {
137  // Set upper and lower bounds
138  scalar Thi = Tc(X);
139  scalar Tlo = Tpt(X);
140 
141  // Check for critical and solid phase conditions
142  if (p >= pv(p, Thi, X))
143  {
144  return Thi;
145  }
146  else if (p < pv(p, Tlo, X))
147  {
149  << "Pressure below triple point pressure: "
150  << "p = " << p << " < Pt = " << pv(p, Tlo, X) << nl << endl;
151  return -1;
152  }
153 
154  // Set initial guess
155  scalar T = (Thi + Tlo)*0.5;
156 
157  while ((Thi - Tlo) > 1.0e-4)
158  {
159  if ((pv(p, T, X) - p) <= 0.0)
160  {
161  Tlo = T;
162  }
163  else
164  {
165  Thi = T;
166  }
167 
168  T = (Thi + Tlo)*0.5;
169  }
170 
171  return T;
172 }
173 
174 
176 {
177  scalar Tpc = 0.0;
178 
179  forAll(properties_, i)
180  {
181  Tpc += X[i]*properties_[i].Tc();
182  }
183 
184  return Tpc;
185 }
186 
187 
189 {
190  scalar Vc = 0.0;
191  scalar Zc = 0.0;
192 
193  forAll(properties_, i)
194  {
195  Vc += X[i]*properties_[i].Vc();
196  Zc += X[i]*properties_[i].Zc();
197  }
198 
199  return RR*Zc*Tpc(X)/Vc;
200 }
201 
202 
204 {
205  scalar omega = 0.0;
206 
207  forAll(properties_, i)
208  {
209  omega += X[i]*properties_[i].omega();
210  }
211 
212  return omega;
213 }
214 
215 
217 (
218  const scalar p,
219  const scalar Tg,
220  const scalar Tl,
221  const scalarField& Xg,
222  const scalarField& Xl
223 ) const
224 {
225  scalarField Xs(Xl.size());
226 
227  // Raoult's Law
228  forAll(Xs, i)
229  {
230  scalar Ti = min(TrMax*properties_[i].Tc(), Tl);
231  Xs[i] = properties_[i].pv(p, Ti)*Xl[i]/p;
232  }
233 
234  return Xs;
235 }
236 
237 
238 Foam::scalar Foam::liquidMixtureProperties::W(const scalarField& X) const
239 {
240  scalar W = 0.0;
241 
242  forAll(properties_, i)
243  {
244  W += X[i]*properties_[i].W();
245  }
246 
247  return W;
248 }
249 
250 
252 {
253  scalarField Y(X.size());
254  scalar sumY = 0.0;
255 
256  forAll(Y, i)
257  {
258  Y[i] = X[i]*properties_[i].W();
259  sumY += Y[i];
260  }
261 
262  Y /= sumY;
263 
264  return Y;
265 }
266 
267 
269 {
270  scalarField X(Y.size());
271  scalar sumX = 0.0;
272 
273  forAll(X, i)
274  {
275  X[i] = Y[i]/properties_[i].W();
276  sumX += X[i];
277  }
278 
279  X /= sumX;
280 
281  return X;
282 }
283 
284 
286 (
287  const scalar p,
288  const scalar T,
289  const scalarField& X
290 ) const
291 {
292  scalar sumY = 0.0;
293  scalar v = 0.0;
294 
295  forAll(properties_, i)
296  {
297  if (X[i] > small)
298  {
299  scalar Ti = min(TrMax*properties_[i].Tc(), T);
300  scalar rho = properties_[i].rho(p, Ti);
301 
302  if (rho > small)
303  {
304  scalar Yi = X[i]*properties_[i].W();
305  sumY += Yi;
306  v += Yi/rho;
307  }
308  }
309  }
310 
311  return sumY/v;
312 }
313 
314 
316 (
317  const scalar p,
318  const scalar T,
319  const scalarField& X
320 ) const
321 {
322  scalar sumY = 0.0;
323  scalar pv = 0.0;
324 
325  forAll(properties_, i)
326  {
327  if (X[i] > small)
328  {
329  scalar Yi = X[i]*properties_[i].W();
330  sumY += Yi;
331 
332  scalar Ti = min(TrMax*properties_[i].Tc(), T);
333  pv += Yi*properties_[i].pv(p, Ti);
334  }
335  }
336 
337  return pv/sumY;
338 }
339 
340 
342 (
343  const scalar p,
344  const scalar T,
345  const scalarField& X
346 ) const
347 {
348  scalar sumY = 0.0;
349  scalar hl = 0.0;
350 
351  forAll(properties_, i)
352  {
353  if (X[i] > small)
354  {
355  scalar Yi = X[i]*properties_[i].W();
356  sumY += Yi;
357 
358  scalar Ti = min(TrMax*properties_[i].Tc(), T);
359  hl += Yi*properties_[i].hl(p, Ti);
360  }
361  }
362 
363  return hl/sumY;
364 }
365 
366 
368 (
369  const scalar p,
370  const scalar T,
371  const scalarField& X
372 ) const
373 {
374  scalar sumY = 0.0;
375  scalar Cp = 0.0;
376 
377  forAll(properties_, i)
378  {
379  if (X[i] > small)
380  {
381  scalar Yi = X[i]*properties_[i].W();
382  sumY += Yi;
383 
384  scalar Ti = min(TrMax*properties_[i].Tc(), T);
385  Cp += Yi*properties_[i].Cp(p, Ti);
386  }
387  }
388 
389  return Cp/sumY;
390 }
391 
392 
394 (
395  const scalar p,
396  const scalar T,
397  const scalarField& X
398 ) const
399 {
400  // sigma is based on surface mole fractions
401  // which are estimated from Raoult's Law
402  scalar sigma = 0.0;
403  scalarField Xs(X.size());
404  scalar XsSum = 0.0;
405 
406  forAll(properties_, i)
407  {
408  scalar Ti = min(TrMax*properties_[i].Tc(), T);
409  scalar Pvs = properties_[i].pv(p, Ti);
410 
411  Xs[i] = X[i]*Pvs/p;
412  XsSum += Xs[i];
413  }
414 
415  Xs /= XsSum;
416 
417  forAll(properties_, i)
418  {
419  if (Xs[i] > small)
420  {
421  scalar Ti = min(TrMax*properties_[i].Tc(), T);
422  sigma += Xs[i]*properties_[i].sigma(p, Ti);
423  }
424  }
425 
426  return sigma;
427 }
428 
429 
431 (
432  const scalar p,
433  const scalar T,
434  const scalarField& X
435 ) const
436 {
437  scalar mu = 0.0;
438 
439  forAll(properties_, i)
440  {
441  if (X[i] > small)
442  {
443  scalar Ti = min(TrMax*properties_[i].Tc(), T);
444  mu += X[i]*log(properties_[i].mu(p, Ti));
445  }
446  }
447 
448  return exp(mu);
449 }
450 
451 
453 (
454  const scalar p,
455  const scalar T,
456  const scalarField& X
457 ) const
458 {
459  // Calculate superficial volume fractions phii
460  scalarField phii(X.size());
461  scalar pSum = 0.0;
462 
463  forAll(properties_, i)
464  {
465  scalar Ti = min(TrMax*properties_[i].Tc(), T);
466 
467  scalar Vi = properties_[i].W()/properties_[i].rho(p, Ti);
468  phii[i] = X[i]*Vi;
469  pSum += phii[i];
470  }
471 
472  phii /= pSum;
473 
474  scalar K = 0.0;
475 
476  forAll(properties_, i)
477  {
478  scalar Ti = min(TrMax*properties_[i].Tc(), T);
479 
480  forAll(properties_, j)
481  {
482  scalar Tj = min(TrMax*properties_[j].Tc(), T);
483 
484  scalar Kij =
485  2.0
486  /(
487  1.0/properties_[i].kappa(p, Ti)
488  + 1.0/properties_[j].kappa(p, Tj)
489  );
490  K += phii[i]*phii[j]*Kij;
491  }
492  }
493 
494  return K;
495 }
496 
497 
499 (
500  const scalar p,
501  const scalar T,
502  const scalarField& X
503 ) const
504 {
505  // Blanc's law
506  scalar Dinv = 0.0;
507 
508  forAll(properties_, i)
509  {
510  if (X[i] > small)
511  {
512  scalar Ti = min(TrMax*properties_[i].Tc(), T);
513  Dinv += X[i]/properties_[i].D(p, Ti);
514  }
515  }
516 
517  return 1.0/Dinv;
518 }
519 
520 
521 // ************************************************************************* //
scalar Cp(const scalar p, const scalar T, const scalarField &X) const
Calculate the mixture heat capacity [J/kg/K].
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
dimensionedScalar log(const dimensionedScalar &ds)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
liquidMixtureProperties(const dictionary &dict)
Construct from dictionary.
scalar D(const scalar p, const scalar T, const scalarField &X) const
Vapour diffusivity [m^2/s].
CGAL::Exact_predicates_exact_constructions_kernel K
wordList toc() const
Return the table of contents.
Definition: dictionary.C:1115
scalar mu(const scalar p, const scalar T, const scalarField &X) const
Calculate the mixture viscosity [Pa s].
bool isDict(const word &) const
Check if entry is a sub-dictionary.
Definition: dictionary.C:956
scalar kappa(const scalar p, const scalar T, const scalarField &X) const
Estimate thermal conductivity [W/m/K].
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:1002
scalarField Xs(const scalar p, const scalar Tg, const scalar Tl, const scalarField &Xg, const scalarField &Xl) const
Return the surface molar fractions.
scalar rho(const scalar p, const scalar T, const scalarField &X) const
Calculate the mixture density [kg/m^3].
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
scalar Tpc(const scalarField &X) const
Return pseudocritical temperature according to Kay&#39;s rule.
static autoPtr< liquidMixtureProperties > New(const dictionary &)
Select construct from dictionary.
scalar Ppc(const scalarField &X) const
Return pseudocritical pressure (modified Prausnitz and Gunn)
dimensionedScalar exp(const dimensionedScalar &ds)
scalarField X(const scalarField &Y) const
Returns the mole fractions corresponding to the given mass fractions.
scalar sigma(const scalar p, const scalar T, const scalarField &X) const
Estimate mixture surface tension [N/m].
scalar Tc(const scalarField &X) const
Calculate the critical temperature of mixture.
scalarField Y(const scalarField &X) const
Returns the mass fractions corresponding to the given mole fractions.
scalar Tpt(const scalarField &X) const
Return pseudo triple point temperature (mole averaged formulation)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
static const char nl
Definition: Ostream.H:260
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
scalar hl(const scalar p, const scalar T, const scalarField &X) const
Calculate the mixture latent heat [J/kg].
void setSize(const label)
Reset size of List.
Definition: List.C:281
#define WarningInFunction
Report a warning using Foam::Warning.
scalar pv(const scalar p, const scalar T, const scalarField &X) const
Calculate the mixture vapour pressure [Pa].
scalar W(const scalarField &X) const
Calculate the mean molecular weight [kg/kmol].
volScalarField & p
scalar omega(const scalarField &X) const
Return mixture accentric factor.
scalar pvInvert(const scalar p, const scalarField &X) const
Invert the vapour pressure relationship to retrieve the boiling.