PengRobinsonGasI.H
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) 2014-2023 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 "PengRobinsonGas.H"
27 #include "mathematicalConstants.H"
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
31 template<class Specie>
32 inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
33 (
34  const Specie& sp,
35  const scalar& Tc,
36  const scalar& Vc,
37  const scalar& Zc,
38  const scalar& Pc,
39  const scalar& omega
40 )
41 :
42  Specie(sp),
43  Tc_(Tc),
44  Vc_(Vc),
45  Zc_(Zc),
46  Pc_(Pc),
47  omega_(omega)
48 {}
49 
50 
51 template<class Specie>
52 inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
53 (
54  const word& name,
55  const PengRobinsonGas& pg
56 )
57 :
58  Specie(name, pg),
59  Tc_(pg.Tc_),
60  Vc_(pg.Vc_),
61  Zc_(pg.Zc_),
62  Pc_(pg.Pc_),
63  omega_(pg.omega_)
64 {}
65 
66 
67 template<class Specie>
68 inline Foam::autoPtr<Foam::PengRobinsonGas <Specie>>
69 Foam::PengRobinsonGas<Specie>::clone() const
70 {
71  return autoPtr<PengRobinsonGas<Specie>>
72  (
73  new PengRobinsonGas<Specie>(*this)
74  );
75 }
76 
77 
78 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
79 
80 template<class Specie>
81 inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
82 (
83  scalar p,
84  scalar T
85 ) const
86 {
87  const scalar Z = this->Z(p, T);
88  return p/(Z*this->R()*T);
89 }
90 
91 
92 template<class Specie>
93 inline Foam::scalar Foam::PengRobinsonGas<Specie>::h(scalar p, scalar T) const
94 {
95  const scalar Pr = p/Pc_;
96  const scalar Tr = T/Tc_;
97  const scalar B = 0.07780*Pr/Tr;
98  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
99  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
100 
101  const scalar Z = this->Z(p, T);
102 
103  return
104  this->R()
105  *Tc_
106  *(
107  Tr*(Z - 1)
108  - 2.078*(1 + kappa)*sqrt(alpha)
109  *log((Z + 2.414*B)/(Z - 0.414*B))
110  );
111 }
112 
113 
114 template<class Specie>
115 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Cp(scalar p, scalar T) const
116 {
117  const scalar Tr = T/Tc_;
118  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
119  const scalar b = 0.07780*RR*Tc_/Pc_;
120  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
121  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
122 
123  const scalar A = a*alpha*p/sqr(RR*T);
124  const scalar B = b*p/(RR*T);
125 
126  const scalar Z = this->Z(p, T);
127 
128  const scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_));
129  const scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_));
130 
131  const scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B);
132  const scalar N = ap*B/(b*RR);
133 
134  const scalar root2 = sqrt(2.0);
135 
136  return
137  (
138  app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B))
139  + RR*sqr(M - N)/(sqr(M) - 2*A*(Z + B))
140  - RR
141  )/this->W();
142 }
143 
144 
145 template<class Specie>
146 inline Foam::scalar Foam::PengRobinsonGas<Specie>::e(scalar p, scalar T) const
147 {
148  const scalar Pr = p/Pc_;
149  const scalar Tr = T/Tc_;
150  const scalar B = 0.07780*Pr/Tr;
151  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
152  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
153 
154  const scalar Z = this->Z(p, T);
155 
156  return
157  this->R()
158  *Tc_
159  *(
160  - 2.078*(1 + kappa)*sqrt(alpha)
161  *log((Z + 2.414*B)/(Z - 0.414*B))
162  );
163 }
164 
165 
166 template<class Specie>
167 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Cv(scalar p, scalar T) const
168 {
169  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
170  const scalar b = 0.07780*RR*Tc_/Pc_;
171  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
172 
173  const scalar B = b*p/(RR*T);
174 
175  const scalar Z = this->Z(p, T);
176 
177  const scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_));
178 
179  const scalar root2 = sqrt(2.0);
180 
181  return
182  (
183  app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B))
184  - RR
185  )/this->W();
186 }
187 
188 
189 template<class Specie>
190 inline Foam::scalar Foam::PengRobinsonGas<Specie>::sp
191 (
192  scalar p,
193  scalar T
194 ) const
195 {
196  const scalar Pr = p/Pc_;
197  const scalar Tr = T/Tc_;
198  const scalar B = 0.07780*Pr/Tr;
199  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
200 
201  const scalar Z = this->Z(p, T);
202 
203  return
204  this->R()
205  *(
206  - log(p/Pstd)
207  + (
208  log(Z - B)
209  - 2.078*kappa*((1 + kappa)/sqrt(Tr) - kappa)
210  *log((Z + 2.414*B)/(Z - 0.414*B))
211  )
212  );
213 }
214 
215 
216 template<class Specie>
217 inline Foam::scalar Foam::PengRobinsonGas<Specie>::sv
218 (
219  scalar p,
220  scalar T
221 ) const
222 {
223  NotImplemented;
224  return 0;
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 = pow(Rl + D05, 1.0/3.0);
286  scalar Tl = 0;
287  if (D05 > Rl)
288  {
289  Tl = -pow(mag(Rl - D05), 1.0/3.0);
290  }
291  else
292  {
293  Tl = pow(Rl - D05, 1.0/3.0);
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 template<class Specie>
330 inline Foam::scalar Foam::PengRobinsonGas<Specie>::alphav
331 (
332  scalar p,
333  scalar T
334 ) const
335 {
336  NotImplemented;
337  return 0;
338 }
339 
340 
341 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
342 
343 template<class Specie>
344 inline void Foam::PengRobinsonGas<Specie>::operator+=
345 (
346  const PengRobinsonGas<Specie>& pg
347 )
348 {
349  scalar X1 = this->Y()/this->W();
350  Specie::operator+=(pg);
351 
352  if (mag(this->Y()) > small)
353  {
354  X1 *= this->W()/this->Y();
355  const scalar X2 = this->W()*pg.Y()/(pg.W()*this->Y());
356 
357  Tc_ = X1*Tc_ + X2*pg.Tc_;
358  Vc_ = X1*Vc_ + X2*pg.Vc_;
359  Zc_ = X1*Zc_ + X2*pg.Zc_;
360  Pc_ = RR*Zc_*Tc_/Vc_;
361  omega_ = X1*omega_ + X2*pg.omega_;
362  }
363 }
364 
365 
366 template<class Specie>
367 inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
368 {
369  Specie::operator*=(s);
370 }
371 
372 
373 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
374 
375 
376 template<class Specie>
377 Foam::PengRobinsonGas<Specie> Foam::operator+
378 (
379  const PengRobinsonGas<Specie>& pg1,
380  const PengRobinsonGas<Specie>& pg2
381 )
382 {
383  Specie sp
384  (
385  static_cast<const Specie&>(pg1)
386  + static_cast<const Specie&>(pg2)
387  );
388 
389  if (mag(sp.Y()) < small)
390  {
391  return PengRobinsonGas<Specie>
392  (
393  sp,
394  pg1.Tc_,
395  pg1.Vc_,
396  pg1.Zc_,
397  pg1.Pc_,
398  pg1.omega_
399  );
400  }
401  else
402  {
403  const scalar X1 = sp.W()*pg1.Y()/(pg1.W()*sp.Y());
404  const scalar X2 = sp.W()*pg2.Y()/(pg2.W()*sp.Y());
405 
406  const scalar Tc = X1*pg1.Tc_ + X2*pg2.Tc_;
407  const scalar Vc = X1*pg1.Vc_ + X2*pg2.Vc_;
408  const scalar Zc = X1*pg1.Zc_ + X2*pg2.Zc_;
409 
410  return PengRobinsonGas<Specie>
411  (
412  sp,
413  Tc,
414  Vc,
415  Zc,
416  RR*Zc*Tc/Vc,
417  X1*pg1.omega_ + X2*pg2.omega_
418  );
419  }
420 }
421 
422 
423 template<class Specie>
424 Foam::PengRobinsonGas<Specie> Foam::operator*
425 (
426  const scalar s,
427  const PengRobinsonGas<Specie>& pg
428 )
429 {
430  return PengRobinsonGas<Specie>
431  (
432  s*static_cast<const Specie&>(pg),
433  pg.Tc_,
434  pg.Vc_,
435  pg.Zc_,
436  pg.Pc_,
437  pg.omega_
438  );
439 }
440 
441 
442 template<class Specie>
443 Foam::PengRobinsonGas<Specie> Foam::operator==
444 (
445  const PengRobinsonGas<Specie>& pg1,
446  const PengRobinsonGas<Specie>& pg2
447 )
448 {
449  return PengRobinsonGas<Specie>
450  (
451  static_cast<const Specie&>(pg1) == static_cast<const Specie&>(pg2),
452  NaN,
453  NaN,
454  NaN,
455  NaN,
456  NaN
457  );
458 }
459 
460 
461 // ************************************************************************* //
A
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
B
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
D
static const Foam::dimensionedScalar D("D", Foam::dimTemperature, 257.14)
M
#define M(I)
Foam::PengRobinsonGas
PengRobinsonGas cubic equation of state for gases.
Definition: PengRobinsonGas.H:127
Foam::PengRobinsonGas::Cv
scalar Cv(scalar p, scalar T) const
Return Cv contribution [J/(kg K].
Definition: PengRobinsonGasI.H:167
Foam::PengRobinsonGas::psi
scalar psi(scalar p, scalar T) const
Return compressibility [s^2/m^2].
Definition: PengRobinsonGasI.H:230
Foam::PengRobinsonGas::alphav
scalar alphav(const scalar p, const scalar T) const
Return volumetric coefficient of thermal expansion [1/T].
Definition: PengRobinsonGasI.H:331
Foam::PengRobinsonGas::e
scalar e(const scalar p, const scalar T) const
Return internal energy contribution [J/kg].
Definition: PengRobinsonGasI.H:146
Foam::PengRobinsonGas::rho
scalar rho(scalar p, scalar T) const
Return density [kg/m^3].
Definition: PengRobinsonGasI.H:82
Foam::PengRobinsonGas::CpMCv
scalar CpMCv(scalar p, scalar T) const
Return (Cp - Cv) [J/(kg K].
Definition: PengRobinsonGasI.H:305
Foam::PengRobinsonGas::h
scalar h(const scalar p, const scalar T) const
Return enthalpy contribution [J/kg].
Definition: PengRobinsonGasI.H:93
Foam::PengRobinsonGas::Cp
scalar Cp(scalar p, scalar T) const
Return Cp contribution [J/(kg K].
Definition: PengRobinsonGasI.H:115
Foam::PengRobinsonGas::clone
autoPtr< PengRobinsonGas > clone() const
Construct and return a clone.
Definition: PengRobinsonGasI.H:69
Foam::PengRobinsonGas::sv
scalar sv(const scalar p, const scalar T) const
Return entropy contribution to the integral of Cv/T [J/kg/K].
Definition: PengRobinsonGasI.H:218
Foam::PengRobinsonGas::sp
scalar sp(const scalar p, const scalar T) const
Return entropy contribution to the integral of Cp/T [J/kg/K].
Definition: PengRobinsonGasI.H:191
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:367
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:33
Foam::autoPtr
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
Foam::omega
Convenience class to handle the input of constant rotational speed. Reads an omega entry with default...
Definition: omega.H:54
Foam::word
A class for handling words, derived from string.
Definition: word.H:62
NotImplemented
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:381
s
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
b
volScalarField & b
Definition: createFields.H:25
alpha
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
Foam::constant::atomic::a0
const dimensionedScalar a0
Bohr radius: default SI units: [m].
Foam::constant::electromagnetic::kappa
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
Foam::constant::physicoChemical::RR
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
Foam::constant::standard::Pstd
const dimensionedScalar Pstd
Standard pressure.
Foam::fvm::S
tmp< fvMatrix< Type > > S(const Pair< tmp< volScalarField::Internal >> &, const VolField< Type > &)
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:49
Foam::pow3
dimensionedScalar pow3(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:87
Foam::name
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
Foam::log
dimensionedScalar log(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:275
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:73
Foam::W
scalarList W(const fluidMulticomponentThermo &thermo)
Definition: thermoTypeFunctions.H:29
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:142
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
Foam::max
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:104
Foam::R
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:55
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:278
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:281
p
volScalarField & p
Definition: createFieldRefs.H:4
Y
PtrList< volScalarField > & Y
Definition: createFieldRefs.H:3