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