PengRobinsonGasI.H
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) 2014-2017 OpenFOAM Foundation
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 "PengRobinsonGas.H"
29 #include "mathematicalConstants.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 template<class Specie>
34 inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
35 (
36  const Specie& sp,
37  const scalar& Tc,
38  const scalar& Vc,
39  const scalar& Zc,
40  const scalar& Pc,
41  const scalar& omega
42 )
43 :
44  Specie(sp),
45  Tc_(Tc),
46  Vc_(Vc),
47  Zc_(Zc),
48  Pc_(Pc),
49  omega_(omega)
50 {}
51 
52 
53 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
54 
55 template<class Specie>
56 inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
57 (
58  const word& name,
59  const PengRobinsonGas& pg
60 )
61 :
62  Specie(name, pg),
63  Tc_(pg.Tc_),
64  Vc_(pg.Vc_),
65  Zc_(pg.Zc_),
66  Pc_(pg.Pc_),
67  omega_(pg.omega_)
68 {}
69 
70 
71 template<class Specie>
72 inline Foam::autoPtr<Foam::PengRobinsonGas <Specie>>
73 Foam::PengRobinsonGas<Specie>::clone() const
74 {
75  return autoPtr<PengRobinsonGas<Specie>>::New(*this);
76 }
77 
78 
79 template<class Specie>
80 inline Foam::autoPtr<Foam::PengRobinsonGas<Specie>>
81 Foam::PengRobinsonGas<Specie>::New
82 (
83  const dictionary& dict
84 )
85 {
86  return autoPtr<PengRobinsonGas<Specie>>::New(dict);
87 }
88 
89 
90 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
91 
92 template<class Specie>
93 inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
94 (
95  scalar p,
96  scalar T
97 ) const
98 {
99  const scalar Z = this->Z(p, T);
100  return p/(Z*this->R()*T);
101 }
102 
103 
104 template<class Specie>
105 inline Foam::scalar Foam::PengRobinsonGas<Specie>::H(scalar p, scalar T) const
106 {
107  const scalar Pr = p/Pc_;
108  const scalar Tr = T/Tc_;
109  const scalar B = 0.07780*Pr/Tr;
110  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
111  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
112 
113  const scalar Z = this->Z(p, T);
114 
115  return
116  this->R()
117  *Tc_
118  *(
119  Tr*(Z - 1)
120  - 2.078*(1 + kappa)*sqrt(alpha)
121  *log((Z + 2.414*B)/(Z - 0.414*B))
122  );
123 }
124 
125 
126 template<class Specie>
127 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Cp(scalar p, scalar T) const
128 {
129  const scalar Tr = T/Tc_;
130  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
131  const scalar b = 0.07780*RR*Tc_/Pc_;
132  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
133  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
134 
135  const scalar A = a*alpha*p/sqr(RR*T);
136  const scalar B = b*p/(RR*T);
137 
138  const scalar Z = this->Z(p, T);
139 
140  const scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_));
141  const scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_));
142 
143  const scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B);
144  const scalar N = ap*B/(b*RR);
145 
146  const scalar root2 = sqrt(2.0);
147 
148  return
149  (
150  app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B))
151  + RR*sqr(M - N)/(sqr(M) - 2*A*(Z + B))
152  - RR
153  )/this->W();
154 }
155 
156 
157 template<class Specie>
158 inline Foam::scalar Foam::PengRobinsonGas<Specie>::E(scalar p, scalar T) const
159 {
160  const scalar Pr = p/Pc_;
161  const scalar Tr = T/Tc_;
162  const scalar B = 0.07780*Pr/Tr;
163  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
164  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
165 
166  const scalar Z = this->Z(p, T);
167 
168  return
169  this->R()
170  *Tc_
171  *(
172  - 2.078*(1 + kappa)*sqrt(alpha)
173  *log((Z + 2.414*B)/(Z - 0.414*B))
174  );
175 }
176 
177 
178 template<class Specie>
179 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Cv(scalar p, scalar T) const
180 {
181  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
182  const scalar b = 0.07780*RR*Tc_/Pc_;
183  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
184 
185  const scalar B = b*p/(RR*T);
186 
187  const scalar Z = this->Z(p, T);
188 
189  const scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_));
190 
191  const scalar root2 = sqrt(2.0);
192 
193  return
194  (
195  app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B))
196  - RR
197  )/this->W();
198 }
199 
200 
201 template<class Specie>
202 inline Foam::scalar Foam::PengRobinsonGas<Specie>::S
203 (
204  scalar p,
205  scalar T
206 ) const
207 {
208  const scalar Pr = p/Pc_;
209  const scalar Tr = T/Tc_;
210  const scalar B = 0.07780*Pr/Tr;
211  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
212 
213  const scalar Z = this->Z(p, T);
214 
215  return
216  this->R()
217  *(
218  - log(p/Pstd)
219  + (
220  log(Z - B)
221  - 2.078*kappa*((1 + kappa)/sqrt(Tr) - kappa)
222  *log((Z + 2.414*B)/(Z - 0.414*B))
223  )
224  );
225 }
226 
227 
228 template<class Specie>
229 inline Foam::scalar Foam::PengRobinsonGas<Specie>::psi
230 (
231  scalar p,
232  scalar T
233 ) const
234 {
235  const scalar Z = this->Z(p, T);
236 
237  return 1.0/(Z*this->R()*T);
238 }
239 
240 
241 template<class Specie>
242 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Z
243 (
244  scalar p,
245  scalar T
246 ) const
247 {
248  const scalar Tr = T/Tc_;
249  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
250  const scalar b = 0.07780*RR*Tc_/Pc_;
251  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
252  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
253 
254  const scalar A = a*alpha*p/sqr(RR*T);
255  const scalar B = b*p/(RR*T);
256 
257  const scalar a2 = B - 1;
258  const scalar a1 = A - 2*B - 3*sqr(B);
259  const scalar a0 = -A*B + sqr(B) + pow3(B);
260 
261  const scalar Q = (3*a1 - a2*a2)/9.0;
262  const scalar Rl = (9*a2*a1 - 27*a0 - 2*a2*a2*a2)/54.0;
263 
264  const scalar Q3 = Q*Q*Q;
265  const scalar D = Q3 + Rl*Rl;
266 
267  scalar root = -1;
268 
269  if (D <= 0)
270  {
271  const scalar th = ::acos(Rl/sqrt(-Q3));
272  const scalar qm = 2*sqrt(-Q);
273  const scalar r1 = qm*cos(th/3.0) - a2/3.0;
274  const scalar r2 =
275  qm*cos((th + 2*constant::mathematical::pi)/3.0) - a2/3.0;
276  const scalar r3 =
277  qm*cos((th + 4*constant::mathematical::pi)/3.0) - a2/3.0;
278 
279  root = max(r1, max(r2, r3));
280  }
281  else
282  {
283  // One root is real
284  const scalar D05 = sqrt(D);
285  const scalar S = cbrt(Rl + D05);
286  scalar Tl = 0;
287  if (D05 > Rl)
288  {
289  Tl = -cbrt(mag(Rl - D05));
290  }
291  else
292  {
293  Tl = cbrt(Rl - D05);
294  }
295 
296  root = S + Tl - a2/3.0;
297  }
298 
299  return root;
300 }
301 
302 
303 template<class Specie>
304 inline Foam::scalar Foam::PengRobinsonGas<Specie>::CpMCv
305 (
306  scalar p,
307  scalar T
308 ) const
309 {
310  const scalar Tr = T/Tc_;
311  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
312  const scalar b = 0.07780*RR*Tc_/Pc_;
313  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
314  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
315 
316  const scalar A = alpha*a*p/sqr(RR*T);
317  const scalar B = b*p/(RR*T);
318 
319  const scalar Z = this->Z(p, T);
320 
321  const scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_));
322  const scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B);
323  const scalar N = ap*B/(b*RR);
324 
325  return this->R()*sqr(M - N)/(sqr(M) - 2*A*(Z + B));
326 }
327 
328 
329 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
330 
331 template<class Specie>
332 inline void Foam::PengRobinsonGas<Specie>::operator+=
333 (
334  const PengRobinsonGas<Specie>& pg
335 )
336 {
337  scalar Y1 = this->Y();
338  Specie::operator+=(pg);
339 
340  if (mag(this->Y()) > SMALL)
341  {
342  Y1 /= this->Y();
343  const scalar Y2 = pg.Y()/this->Y();
344 
345  Tc_ = Y1*Tc_ + Y2*pg.Tc_;
346  Vc_ = Y1*Vc_ + Y2*pg.Vc_;
347  Zc_ = Y1*Zc_ + Y2*pg.Zc_;
348  Pc_ = RR*Zc_*Tc_/Vc_;
349  omega_ = Y1*omega_ + Y2*pg.omega_;
350  }
351 }
352 
353 
354 template<class Specie>
355 inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
356 {
357  Specie::operator*=(s);
358 }
359 
360 
361 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
362 
363 
364 template<class Specie>
365 Foam::PengRobinsonGas<Specie> Foam::operator+
366 (
367  const PengRobinsonGas<Specie>& pg1,
368  const PengRobinsonGas<Specie>& pg2
369 )
370 {
371  Specie sp
372  (
373  static_cast<const Specie&>(pg1)
374  + static_cast<const Specie&>(pg2)
375  );
376 
377  if (mag(sp.Y()) < SMALL)
378  {
379  return PengRobinsonGas<Specie>
380  (
381  sp,
382  pg1.Tc_,
383  pg1.Vc_,
384  pg1.Zc_,
385  pg1.Pc_,
386  pg1.omega_
387  );
388  }
389  else
390  {
391  const scalar Y1 = pg1.Y()/sp.Y();
392  const scalar Y2 = pg2.Y()/sp.Y();
393 
394  const scalar Tc = Y1*pg1.Tc_ + Y2*pg2.Tc_;
395  const scalar Vc = Y1*pg1.Vc_ + Y2*pg2.Vc_;
396  const scalar Zc = Y1*pg1.Zc_ + Y2*pg2.Zc_;
397 
398  return PengRobinsonGas<Specie>
399  (
400  sp,
401  Tc,
402  Vc,
403  Zc,
404  RR*Zc*Tc/Vc,
405  Y1*pg1.omega_ + Y2*pg2.omega_
406  );
407  }
408 }
409 
410 
411 template<class Specie>
412 Foam::PengRobinsonGas<Specie> Foam::operator*
413 (
414  const scalar s,
415  const PengRobinsonGas<Specie>& pg
416 )
417 {
418  return PengRobinsonGas<Specie>
419  (
420  s*static_cast<const Specie&>(pg),
421  pg.Tc_,
422  pg.Vc_,
423  pg.Zc_,
424  pg.Pc_,
425  pg.omega_
426  );
427 }
428 
429 
430 template<class Specie>
431 Foam::PengRobinsonGas<Specie> Foam::operator==
432 (
433  const PengRobinsonGas<Specie>& pg1,
434  const PengRobinsonGas<Specie>& pg2
435 )
436 {
437  Specie sp
438  (
439  static_cast<const Specie&>(pg1)
440  == static_cast<const Specie&>(pg2)
441  );
442 
443  const scalar Y1 = pg1.Y()/sp.Y();
444  const scalar Y2 = pg2.Y()/sp.Y();
445 
446  const scalar Tc = Y2*pg2.Tc_ - Y1*pg1.Tc_;
447  const scalar Vc = Y2*pg2.Vc_ - Y1*pg1.Vc_;
448  const scalar Zc = Y2*pg2.Zc_ - Y1*pg1.Zc_;
449 
450  return PengRobinsonGas<Specie>
451  (
452  sp,
453  Tc,
454  Vc,
455  Zc,
456  RR*Zc*Tc/Vc,
457  Y2*pg2.omega_ - Y1*pg1.omega_
458  );
459 }
460 
461 
462 // ************************************************************************* //
Foam::constant::thermodynamic::RR
const scalar RR
Universal gas constant: default in [J/(kmol K)].
Definition: thermodynamicConstants.C:46
Foam::PengRobinsonGas::psi
scalar psi(scalar p, scalar T) const
Return compressibility rho/p [s^2/m^2].
Definition: PengRobinsonGasI.H:230
mathematicalConstants.H
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::constant::standard::Pstd
const dimensionedScalar Pstd
Standard pressure.
Definition: thermodynamicConstants.C:48
Foam::PengRobinsonGas::Z
scalar Z(scalar p, scalar T) const
Return compression factor [].
Definition: PengRobinsonGasI.H:243
Foam::PengRobinsonGas::operator*=
void operator*=(const scalar)
Definition: PengRobinsonGasI.H:355
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:212
Foam::PengRobinsonGas::Cv
scalar Cv(scalar p, scalar T) const
Return Cv departure [J/(kg K].
Definition: PengRobinsonGasI.H:179
B
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
Foam::PengRobinsonGas::rho
scalar rho(scalar p, scalar T) const
Return density [kg/m^3].
Definition: PengRobinsonGasI.H:94
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::PengRobinsonGas::CpMCv
scalar CpMCv(scalar p, scalar T) const
Return (Cp - Cv) [J/(kg K].
Definition: PengRobinsonGasI.H:305
Foam::constant::electromagnetic::kappa
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
Foam::PengRobinsonGas::clone
autoPtr< PengRobinsonGas > clone() const
Construct and return a clone.
Definition: PengRobinsonGasI.H:73
Foam::PengRobinsonGas::H
scalar H(const scalar p, const scalar T) const
Return enthalpy departure [J/kg].
Definition: PengRobinsonGasI.H:105
R
#define R(A, B, C, D, E, F, K, M)
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::pow3
dimensionedScalar pow3(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:89
Foam::PengRobinsonGas::New
static autoPtr< PengRobinsonGas > New(const dictionary &dict)
Definition: PengRobinsonGasI.H:82
Pr
dimensionedScalar Pr("Pr", dimless, laminarTransport)
Foam::PengRobinsonGas
PengRobinsonGas gas equation of state.
Definition: PengRobinsonGas.H:52
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
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::PengRobinsonGas::PengRobinsonGas
PengRobinsonGas(const Specie &sp, const scalar &Tc, const scalar &Vc, const scalar &Zc, const scalar &Pc, const scalar &omega)
Construct from components.
Definition: PengRobinsonGasI.H:35
Foam::PengRobinsonGas::Cp
scalar Cp(scalar p, scalar T) const
Return Cp departure [J/(kg K].
Definition: PengRobinsonGasI.H:127
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
PengRobinsonGas.H
Foam::log
dimensionedScalar log(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:262
Foam::PengRobinsonGas::S
scalar S(const scalar p, const scalar T) const
Return entropy [J/(kg K)].
Definition: PengRobinsonGasI.H:203
T
const volScalarField & T
Definition: createFieldRefs.H:2
Y
PtrList< volScalarField > & Y
Definition: createFieldRefs.H:7
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::constant::atomic::a0
const dimensionedScalar a0
Bohr radius: default SI units: [m].
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:268
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
M
#define M(I)
D
const dimensionedScalar & D
Definition: solveBulkSurfactant.H:4
Foam::PengRobinsonGas::E
scalar E(const scalar p, const scalar T) const
Return internal energy departure [J/kg].
Definition: PengRobinsonGasI.H:158
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
N
const Vector< label > N(dict.get< Vector< label >>("N"))
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265