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-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 "PengRobinsonGas.H"
27 #include "mathematicalConstants.H"
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 template<class Specie>
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 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
52 
53 template<class Specie>
55 (
56  const word& name,
57  const PengRobinsonGas& pg
58 )
59 :
60  Specie(name, pg),
61  Tc_(pg.Tc_),
62  Vc_(pg.Vc_),
63  Zc_(pg.Zc_),
64  Pc_(pg.Pc_),
65  omega_(pg.omega_)
66 {}
67 
68 
69 template<class Specie>
72 {
74  (
75  new PengRobinsonGas<Specie>(*this)
76  );
77 }
78 
79 
80 template<class Specie>
83 (
84  const dictionary& dict
85 )
86 {
88  (
90  );
91 }
92 
93 
94 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
95 
96 template<class Specie>
97 inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
98 (
99  scalar p,
100  scalar T
101 ) const
102 {
103  const scalar Z = this->Z(p, T);
104  return p/(Z*this->R()*T);
105 }
106 
107 
108 template<class Specie>
109 inline Foam::scalar Foam::PengRobinsonGas<Specie>::H(scalar p, scalar T) const
110 {
111  const scalar Pr = p/Pc_;
112  const scalar Tr = T/Tc_;
113  const scalar B = 0.07780*Pr/Tr;
114  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
115  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
116 
117  const scalar Z = this->Z(p, T);
118 
119  return
120  this->R()
121  *Tc_
122  *(
123  Tr*(Z - 1)
124  - 2.078*(1 + kappa)*sqrt(alpha)
125  *log((Z + 2.414*B)/(Z - 0.414*B))
126  );
127 }
128 
129 
130 template<class Specie>
131 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Cp(scalar p, scalar T) const
132 {
133  const scalar Tr = T/Tc_;
134  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
135  const scalar b = 0.07780*RR*Tc_/Pc_;
136  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
137  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
138 
139  const scalar A = a*alpha*p/sqr(RR*T);
140  const scalar B = b*p/(RR*T);
141 
142  const scalar Z = this->Z(p, T);
143 
144  const scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_));
145  const scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_));
146 
147  const scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B);
148  const scalar N = ap*B/(b*RR);
149 
150  const scalar root2 = sqrt(2.0);
151 
152  return
153  (
154  app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B))
155  + RR*sqr(M - N)/(sqr(M) - 2*A*(Z + B))
156  - RR
157  )/this->W();
158 }
159 
160 
161 template<class Specie>
162 inline Foam::scalar Foam::PengRobinsonGas<Specie>::E(scalar p, scalar T) const
163 {
164  const scalar Pr = p/Pc_;
165  const scalar Tr = T/Tc_;
166  const scalar B = 0.07780*Pr/Tr;
167  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
168  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
169 
170  const scalar Z = this->Z(p, T);
171 
172  return
173  this->R()
174  *Tc_
175  *(
176  - 2.078*(1 + kappa)*sqrt(alpha)
177  *log((Z + 2.414*B)/(Z - 0.414*B))
178  );
179 }
180 
181 
182 template<class Specie>
183 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Cv(scalar p, scalar T) const
184 {
185  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
186  const scalar b = 0.07780*RR*Tc_/Pc_;
187  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
188 
189  const scalar B = b*p/(RR*T);
190 
191  const scalar Z = this->Z(p, T);
192 
193  const scalar app = kappa*a*(1 + kappa)/(2*sqrt(pow3(T)*Tc_));
194 
195  const scalar root2 = sqrt(2.0);
196 
197  return
198  (
199  app*(T/(2*root2*b))*log((Z + (root2 + 1)*B)/(Z - (root2 - 1)*B))
200  - RR
201  )/this->W();
202 }
203 
204 
205 template<class Specie>
206 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Sp
207 (
208  scalar p,
209  scalar T
210 ) const
211 {
212  const scalar Pr = p/Pc_;
213  const scalar Tr = T/Tc_;
214  const scalar B = 0.07780*Pr/Tr;
215  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
216 
217  const scalar Z = this->Z(p, T);
218 
219  return
220  this->R()
221  *(
222  - log(p/Pstd)
223  + (
224  log(Z - B)
225  - 2.078*kappa*((1 + kappa)/sqrt(Tr) - kappa)
226  *log((Z + 2.414*B)/(Z - 0.414*B))
227  )
228  );
229 }
230 
231 
232 template<class Specie>
233 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Sv
234 (
235  scalar p,
236  scalar T
237 ) const
238 {
240  return 0;
241 }
242 
243 
244 template<class Specie>
245 inline Foam::scalar Foam::PengRobinsonGas<Specie>::psi
246 (
247  scalar p,
248  scalar T
249 ) const
250 {
251  const scalar Z = this->Z(p, T);
252 
253  return 1.0/(Z*this->R()*T);
254 }
255 
256 
257 template<class Specie>
258 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Z
259 (
260  scalar p,
261  scalar T
262 ) const
263 {
264  const scalar Tr = T/Tc_;
265  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
266  const scalar b = 0.07780*RR*Tc_/Pc_;
267  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
268  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
269 
270  const scalar A = a*alpha*p/sqr(RR*T);
271  const scalar B = b*p/(RR*T);
272 
273  const scalar a2 = B - 1;
274  const scalar a1 = A - 2*B - 3*sqr(B);
275  const scalar a0 = -A*B + sqr(B) + pow3(B);
276 
277  const scalar Q = (3*a1 - a2*a2)/9.0;
278  const scalar Rl = (9*a2*a1 - 27*a0 - 2*a2*a2*a2)/54.0;
279 
280  const scalar Q3 = Q*Q*Q;
281  const scalar D = Q3 + Rl*Rl;
282 
283  scalar root = -1;
284 
285  if (D <= 0)
286  {
287  const scalar th = ::acos(Rl/sqrt(-Q3));
288  const scalar qm = 2*sqrt(-Q);
289  const scalar r1 = qm*cos(th/3.0) - a2/3.0;
290  const scalar r2 =
291  qm*cos((th + 2*constant::mathematical::pi)/3.0) - a2/3.0;
292  const scalar r3 =
293  qm*cos((th + 4*constant::mathematical::pi)/3.0) - a2/3.0;
294 
295  root = max(r1, max(r2, r3));
296  }
297  else
298  {
299  // One root is real
300  const scalar D05 = sqrt(D);
301  const scalar S = pow(Rl + D05, 1.0/3.0);
302  scalar Tl = 0;
303  if (D05 > Rl)
304  {
305  Tl = -pow(mag(Rl - D05), 1.0/3.0);
306  }
307  else
308  {
309  Tl = pow(Rl - D05, 1.0/3.0);
310  }
311 
312  root = S + Tl - a2/3.0;
313  }
314 
315  return root;
316 }
317 
318 
319 template<class Specie>
320 inline Foam::scalar Foam::PengRobinsonGas<Specie>::CpMCv
321 (
322  scalar p,
323  scalar T
324 ) const
325 {
326  const scalar Tr = T/Tc_;
327  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
328  const scalar b = 0.07780*RR*Tc_/Pc_;
329  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
330  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
331 
332  const scalar A = alpha*a*p/sqr(RR*T);
333  const scalar B = b*p/(RR*T);
334 
335  const scalar Z = this->Z(p, T);
336 
337  const scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_));
338  const scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B);
339  const scalar N = ap*B/(b*RR);
340 
341  return this->R()*sqr(M - N)/(sqr(M) - 2*A*(Z + B));
342 }
343 
344 
345 template<class Specie>
346 inline Foam::scalar Foam::PengRobinsonGas<Specie>::alphav
347 (
348  scalar p,
349  scalar T
350 ) const
351 {
353  return 0;
354 }
355 
356 
357 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
358 
359 template<class Specie>
360 inline void Foam::PengRobinsonGas<Specie>::operator+=
361 (
362  const PengRobinsonGas<Specie>& pg
363 )
364 {
365  scalar X1 = this->Y()/this->W();
366  Specie::operator+=(pg);
367 
368  if (mag(this->Y()) > small)
369  {
370  X1 *= this->W()/this->Y();
371  const scalar X2 = this->W()*pg.Y()/(pg.W()*this->Y());
372 
373  Tc_ = X1*Tc_ + X2*pg.Tc_;
374  Vc_ = X1*Vc_ + X2*pg.Vc_;
375  Zc_ = X1*Zc_ + X2*pg.Zc_;
376  Pc_ = RR*Zc_*Tc_/Vc_;
377  omega_ = X1*omega_ + X2*pg.omega_;
378  }
379 }
380 
381 
382 template<class Specie>
383 inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
384 {
385  Specie::operator*=(s);
386 }
387 
388 
389 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
390 
391 
392 template<class Specie>
393 Foam::PengRobinsonGas<Specie> Foam::operator+
394 (
395  const PengRobinsonGas<Specie>& pg1,
396  const PengRobinsonGas<Specie>& pg2
397 )
398 {
399  Specie sp
400  (
401  static_cast<const Specie&>(pg1)
402  + static_cast<const Specie&>(pg2)
403  );
404 
405  if (mag(sp.Y()) < small)
406  {
408  (
409  sp,
410  pg1.Tc_,
411  pg1.Vc_,
412  pg1.Zc_,
413  pg1.Pc_,
414  pg1.omega_
415  );
416  }
417  else
418  {
419  const scalar X1 = sp.W()*pg1.Y()/(pg1.W()*sp.Y());
420  const scalar X2 = sp.W()*pg2.Y()/(pg2.W()*sp.Y());
421 
422  const scalar Tc = X1*pg1.Tc_ + X2*pg2.Tc_;
423  const scalar Vc = X1*pg1.Vc_ + X2*pg2.Vc_;
424  const scalar Zc = X1*pg1.Zc_ + X2*pg2.Zc_;
425 
427  (
428  sp,
429  Tc,
430  Vc,
431  Zc,
432  RR*Zc*Tc/Vc,
433  X1*pg1.omega_ + X2*pg2.omega_
434  );
435  }
436 }
437 
438 
439 template<class Specie>
440 Foam::PengRobinsonGas<Specie> Foam::operator*
441 (
442  const scalar s,
443  const PengRobinsonGas<Specie>& pg
444 )
445 {
447  (
448  s*static_cast<const Specie&>(pg),
449  pg.Tc_,
450  pg.Vc_,
451  pg.Zc_,
452  pg.Pc_,
453  pg.omega_
454  );
455 }
456 
457 
458 template<class Specie>
459 Foam::PengRobinsonGas<Specie> Foam::operator==
460 (
461  const PengRobinsonGas<Specie>& pg1,
462  const PengRobinsonGas<Specie>& pg2
463 )
464 {
466  (
467  static_cast<const Specie&>(pg1) == static_cast<const Specie&>(pg2),
468  NaN,
469  NaN,
470  NaN,
471  NaN,
472  NaN
473  );
474 }
475 
476 
477 // ************************************************************************* //
dictionary dict
dimensionedScalar acos(const dimensionedScalar &ds)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
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
scalar alphav(const scalar p, const scalar T) const
Return volumetric coefficient of thermal expansion [1/T].
dimensionedSymmTensor sqr(const dimensionedVector &dv)
scalar psi(scalar p, scalar T) const
Return compressibility [s^2/m^2].
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
dimensionedScalar sqrt(const dimensionedScalar &ds)
scalar Z(scalar p, scalar T) const
Return compression factor [].
void operator*=(const scalar)
scalar H(const scalar p, const scalar T) const
Return enthalpy contribution [J/kg].
scalar Sp(const scalar p, const scalar T) const
Return entropy contribution to the integral of Cp/T [J/kg/K].
scalar Cv(scalar p, scalar T) const
Return Cv contribution [J/(kg K].
autoPtr< PengRobinsonGas > clone() const
Construct and return a clone.
scalar rho(scalar p, scalar T) const
Return density [kg/m^3].
scalar Cp(scalar p, scalar T) const
Return Cp contribution [J/(kg K].
const dimensionedScalar RR
Universal gas constant: default SI units: [J/kmol/K].
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
dimensionedScalar cos(const dimensionedScalar &ds)
scalar Sv(const scalar p, const scalar T) const
Return entropy contribution to the integral of Cv/T [J/kg/K].
scalar CpMCv(scalar p, scalar T) const
Return (Cp - Cv) [J/(kg K].
A class for handling words, derived from string.
Definition: word.H:59
tmp< fvMatrix< Type > > S(const Pair< tmp< volScalarField::Internal >> &, const GeometricField< Type, fvPatchField, volMesh > &)
static autoPtr< PengRobinsonGas > New(const dictionary &dict)
PengRobinsonGas cubic equation of state for gases.
PengRobinsonGas(const Specie &sp, const scalar &Tc, const scalar &Vc, const scalar &Zc, const scalar &Pc, const scalar &omega)
Construct from components.
const dimensionedScalar Pstd
Standard pressure.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
const volScalarField & T
#define R(A, B, C, D, E, F, K, M)
const scalarList W(::W(thermo))
PtrList< volScalarField > & Y
label N
Definition: createFields.H:9
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
scalar E(const scalar p, const scalar T) const
Return internal energy contribution [J/kg].
const dimensionedScalar a0
Bohr radius: default SI units: [m].
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:353
#define M(I)