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-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 
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>::S
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>::psi
234 (
235  scalar p,
236  scalar T
237 ) const
238 {
239  const scalar Z = this->Z(p, T);
240 
241  return 1.0/(Z*this->R()*T);
242 }
243 
244 
245 template<class Specie>
246 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Z
247 (
248  scalar p,
249  scalar T
250 ) const
251 {
252  const scalar Tr = T/Tc_;
253  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
254  const scalar b = 0.07780*RR*Tc_/Pc_;
255  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
256  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
257 
258  const scalar A = a*alpha*p/sqr(RR*T);
259  const scalar B = b*p/(RR*T);
260 
261  const scalar a2 = B - 1;
262  const scalar a1 = A - 2*B - 3*sqr(B);
263  const scalar a0 = -A*B + sqr(B) + pow3(B);
264 
265  const scalar Q = (3*a1 - a2*a2)/9.0;
266  const scalar Rl = (9*a2*a1 - 27*a0 - 2*a2*a2*a2)/54.0;
267 
268  const scalar Q3 = Q*Q*Q;
269  const scalar D = Q3 + Rl*Rl;
270 
271  scalar root = -1;
272 
273  if (D <= 0)
274  {
275  const scalar th = ::acos(Rl/sqrt(-Q3));
276  const scalar qm = 2*sqrt(-Q);
277  const scalar r1 = qm*cos(th/3.0) - a2/3.0;
278  const scalar r2 =
279  qm*cos((th + 2*constant::mathematical::pi)/3.0) - a2/3.0;
280  const scalar r3 =
281  qm*cos((th + 4*constant::mathematical::pi)/3.0) - a2/3.0;
282 
283  root = max(r1, max(r2, r3));
284  }
285  else
286  {
287  // One root is real
288  const scalar D05 = sqrt(D);
289  const scalar S = pow(Rl + D05, 1.0/3.0);
290  scalar Tl = 0;
291  if (D05 > Rl)
292  {
293  Tl = -pow(mag(Rl - D05), 1.0/3.0);
294  }
295  else
296  {
297  Tl = pow(Rl - D05, 1.0/3.0);
298  }
299 
300  root = S + Tl - a2/3.0;
301  }
302 
303  return root;
304 }
305 
306 
307 template<class Specie>
308 inline Foam::scalar Foam::PengRobinsonGas<Specie>::CpMCv
309 (
310  scalar p,
311  scalar T
312 ) const
313 {
314  const scalar Tr = T/Tc_;
315  const scalar a = 0.45724*sqr(RR*Tc_)/Pc_;
316  const scalar b = 0.07780*RR*Tc_/Pc_;
317  const scalar kappa = 0.37464 + 1.54226*omega_ - 0.26992*sqr(omega_);
318  const scalar alpha = sqr(1 + kappa*(1 - sqrt(Tr)));
319 
320  const scalar A = alpha*a*p/sqr(RR*T);
321  const scalar B = b*p/(RR*T);
322 
323  const scalar Z = this->Z(p, T);
324 
325  const scalar ap = kappa*a*(kappa/Tc_ - (1 + kappa)/sqrt(T*Tc_));
326  const scalar M = (sqr(Z) + 2*B*Z - sqr(B))/(Z - B);
327  const scalar N = ap*B/(b*RR);
328 
329  return this->R()*sqr(M - N)/(sqr(M) - 2*A*(Z + B));
330 }
331 
332 
333 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
334 
335 template<class Specie>
336 inline void Foam::PengRobinsonGas<Specie>::operator+=
337 (
338  const PengRobinsonGas<Specie>& pg
339 )
340 {
341  scalar Y1 = this->Y();
342  Specie::operator+=(pg);
343 
344  if (mag(this->Y()) > small)
345  {
346  Y1 /= this->Y();
347  const scalar Y2 = pg.Y()/this->Y();
348 
349  Tc_ = Y1*Tc_ + Y2*pg.Tc_;
350  Vc_ = Y1*Vc_ + Y2*pg.Vc_;
351  Zc_ = Y1*Zc_ + Y2*pg.Zc_;
352  Pc_ = RR*Zc_*Tc_/Vc_;
353  omega_ = Y1*omega_ + Y2*pg.omega_;
354  }
355 }
356 
357 
358 template<class Specie>
359 inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
360 {
361  Specie::operator*=(s);
362 }
363 
364 
365 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
366 
367 
368 template<class Specie>
369 Foam::PengRobinsonGas<Specie> Foam::operator+
370 (
371  const PengRobinsonGas<Specie>& pg1,
372  const PengRobinsonGas<Specie>& pg2
373 )
374 {
375  Specie sp
376  (
377  static_cast<const Specie&>(pg1)
378  + static_cast<const Specie&>(pg2)
379  );
380 
381  if (mag(sp.Y()) < small)
382  {
384  (
385  sp,
386  pg1.Tc_,
387  pg1.Vc_,
388  pg1.Zc_,
389  pg1.Pc_,
390  pg1.omega_
391  );
392  }
393  else
394  {
395  const scalar Y1 = pg1.Y()/sp.Y();
396  const scalar Y2 = pg2.Y()/sp.Y();
397 
398  const scalar Tc = Y1*pg1.Tc_ + Y2*pg2.Tc_;
399  const scalar Vc = Y1*pg1.Vc_ + Y2*pg2.Vc_;
400  const scalar Zc = Y1*pg1.Zc_ + Y2*pg2.Zc_;
401 
403  (
404  sp,
405  Tc,
406  Vc,
407  Zc,
408  RR*Zc*Tc/Vc,
409  Y1*pg1.omega_ + Y2*pg2.omega_
410  );
411  }
412 }
413 
414 
415 template<class Specie>
416 Foam::PengRobinsonGas<Specie> Foam::operator*
417 (
418  const scalar s,
419  const PengRobinsonGas<Specie>& pg
420 )
421 {
423  (
424  s*static_cast<const Specie&>(pg),
425  pg.Tc_,
426  pg.Vc_,
427  pg.Zc_,
428  pg.Pc_,
429  pg.omega_
430  );
431 }
432 
433 
434 template<class Specie>
435 Foam::PengRobinsonGas<Specie> Foam::operator==
436 (
437  const PengRobinsonGas<Specie>& pg1,
438  const PengRobinsonGas<Specie>& pg2
439 )
440 {
441  Specie sp
442  (
443  static_cast<const Specie&>(pg1)
444  == static_cast<const Specie&>(pg2)
445  );
446 
447  const scalar Y1 = pg1.Y()/sp.Y();
448  const scalar Y2 = pg2.Y()/sp.Y();
449 
450  const scalar Tc = Y2*pg2.Tc_ - Y1*pg1.Tc_;
451  const scalar Vc = Y2*pg2.Vc_ - Y1*pg1.Vc_;
452  const scalar Zc = Y2*pg2.Zc_ - Y1*pg1.Zc_;
453 
455  (
456  sp,
457  Tc,
458  Vc,
459  Zc,
460  RR*Zc*Tc/Vc,
461  Y2*pg2.omega_ - Y1*pg1.omega_
462  );
463 }
464 
465 
466 // ************************************************************************* //
dictionary dict
dimensionedScalar acos(const dimensionedScalar &ds)
dimensionedScalar log(const dimensionedScalar &ds)
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
scalar psi(scalar p, scalar T) const
Return compressibility rho/p [s^2/m^2].
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 departure [J/kg].
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
scalar Cv(scalar p, scalar T) const
Return Cv departure [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 departure [J/(kg K].
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 CpMCv(scalar p, scalar T) const
Return (Cp - Cv) [J/(kg K].
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
A class for handling words, derived from string.
Definition: word.H:59
static autoPtr< PengRobinsonGas > New(const dictionary &dict)
PengRobinsonGas gas equation of state.
const dimensionedScalar Pstd
Standard pressure.
PengRobinsonGas(const Specie &sp, const scalar &Tc, const scalar &Vc, const scalar &Zc, const scalar &Pc, const scalar &omega)
Construct from components.
scalar S(const scalar p, const scalar T) const
Return entropy [J/kg/K].
const volScalarField & T
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
#define R(A, B, C, D, E, F, K, M)
const scalarList W(::W(thermo))
PtrList< volScalarField > & Y
const dimensionedScalar a0
Bohr radius: default SI units: [m].
label N
Definition: createFields.H:22
dimensioned< scalar > mag(const dimensioned< Type > &)
const scalar RR
Universal gas constant (default in [J/kmol/K])
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
scalar E(const scalar p, const scalar T) const
Return internal energy departure [J/kg].
#define M(I)