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 | Copyright (C) 2014 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  Pc_(pg.Pc_),
63  Vc_(pg.Vc_),
64  Zc_(pg.Zc_),
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  Istream& is
85 )
86 {
88 }
89 
90 
91 template<class Specie>
94 (
95  const dictionary& dict
96 )
97 {
99  (
101  );
102 }
103 
104 
105 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
106 
107 template<class Specie>
108 inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
109 (
110  scalar p,
111  scalar T
112 ) const
113 {
114  scalar z = Z(p, T);
115  return p/(z*this->R()*T);
116 }
117 
118 
119 template<class Specie>
120 inline Foam::scalar Foam::PengRobinsonGas<Specie>::s
121 (
122  scalar p,
123  scalar T
124 ) const
125 {
126  //***HGW This is the expression for a perfect gas
127  // Need to add the entropy defect for Peng-Robinson
128  return -RR*log(p/Pstd);
129 }
130 
131 
132 template<class Specie>
133 inline Foam::scalar Foam::PengRobinsonGas<Specie>::psi
134 (
135  scalar p,
136  scalar T
137 ) const
138 {
139  scalar z = Z(p, T);
140  return 1.0/(z*this->R()*T);
141 }
142 
143 
144 template<class Specie>
145 inline Foam::scalar Foam::PengRobinsonGas<Specie>::Z
146 (
147  scalar p,
148  scalar T
149 ) const
150 {
151  scalar a = 0.45724*sqr(this->R())*sqr(Tc_)/Pc_;
152  scalar b = 0.07780*this->R()*Tc_/Pc_;
153  scalar Tr = T/Tc_;
154 
155  scalar alpha =
156  sqr
157  (
158  1.0
159  + (0.37464 + 1.54226*omega_- 0.26992*sqr(omega_))
160  * (1.0 - sqrt(Tr))
161  );
162 
163  scalar B = b*p/(this->R()*T);
164  scalar A = a*alpha*p/sqr(this->R()*T);
165 
166  scalar a2 = B - 1.0;
167  scalar a1 = A - 2.0*B - 3.0*sqr(B);
168  scalar a0 = -A*B + sqr(B) + pow3(B);
169 
170  scalar Q = (3.0*a1 - a2*a2)/9.0;
171  scalar Rl = (9.0*a2*a1 - 27.0*a0 - 2.0*a2*a2*a2)/54;
172 
173  scalar Q3 = Q*Q*Q;
174  scalar D = Q3 + Rl*Rl;
175 
176  scalar root = -1.0;
177 
178  if (D <= 0.0)
179  {
180  scalar th = ::acos(Rl/sqrt(-Q3));
181  scalar qm = 2.0*sqrt(-Q);
182  scalar r1 = qm*cos(th/3.0) - a2/3.0;
183  scalar r2 = qm*cos((th + 2.0*constant::mathematical::pi)/3.0) - a2/3.0;
184  scalar r3 = qm*cos((th + 4.0*constant::mathematical::pi)/3.0) - a2/3.0;
185 
186  root = max(r1, max(r2, r3));
187  }
188  else
189  {
190  // one root is real
191  scalar D05 = sqrt(D);
192  scalar S = pow(Rl + D05, 1.0/3.0);
193  scalar Tl = 0;
194  if (D05 > Rl)
195  {
196  Tl = -pow(mag(Rl - D05), 1.0/3.0);
197  }
198  else
199  {
200  Tl = pow(Rl - D05, 1.0/3.0);
201  }
202 
203  root = S + Tl - a2/3.0;
204  }
205 
206  return root;
207 }
208 
209 
210 template<class Specie>
211 inline Foam::scalar Foam::PengRobinsonGas<Specie>::cpMcv
212 (
213  scalar p,
214  scalar T
215 ) const
216 {
217  return RR*Z(p, T);
218 }
219 
220 
221 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
222 
223 template<class Specie>
224 inline void Foam::PengRobinsonGas<Specie>::operator+=
225 (
226  const PengRobinsonGas<Specie>& pg
227 )
228 {
229  scalar molr1 = this->nMoles();
230  Specie::operator+=(pg);
231 
232  molr1 /= this->nMoles();
233  scalar molr2 = pg.nMoles()/this->nMoles();
234 
235  Tc_ = molr1*Tc_ + molr2*pg.Tc_;
236  Vc_ = molr1*Vc_ + molr2*pg.Vc_;
237  Zc_ = molr1*Zc_ + molr2*pg.Zc_;
238  Pc_ = RR*Zc_*Tc_/Vc_;
239  omega_ = molr1*omega_ + molr2*pg.omega_;
240 }
241 
242 
243 template<class Specie>
244 inline void Foam::PengRobinsonGas<Specie>::operator-=
245 (
246  const PengRobinsonGas<Specie>& pg
247 )
248 {
249  scalar molr1 = this->nMoles();
250 
251  Specie::operator-=(pg);
252 
253  molr1 /= this->nMoles();
254  scalar molr2 = pg.nMoles()/this->nMoles();
255 
256  Tc_ = molr1*Tc_ - molr2*pg.Tc_;
257  Vc_ = molr1*Vc_ - molr2*pg.Vc_;
258  Zc_ = molr1*Zc_ - molr2*pg.Zc_;
259  Pc_ = RR*Zc_*Tc_/Vc_;
260  omega_ = molr1*omega_ - molr2*pg.omega_;
261 }
262 
263 
264 template<class Specie>
265 inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
266 {
267  Specie::operator*=(s);
268 }
269 
270 
271 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
272 
273 
274 template<class Specie>
275 Foam::PengRobinsonGas<Specie> Foam::operator+
276 (
277  const PengRobinsonGas<Specie>& pg1,
278  const PengRobinsonGas<Specie>& pg2
279 )
280 {
281  scalar nMoles = pg1.nMoles() + pg2.nMoles();
282  scalar molr1 = pg1.nMoles()/nMoles;
283  scalar molr2 = pg2.nMoles()/nMoles;
284 
285  const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
286  const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
287  const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
288 
290  (
291  static_cast<const Specie&>(pg1)
292  + static_cast<const Specie&>(pg2),
293  Tc,
294  Vc,
295  Zc,
296  RR*Zc*Tc/Vc,
297  molr1*pg1.omega_ + molr2*pg2.omega_
298  );
299 }
300 
301 
302 template<class Specie>
303 Foam::PengRobinsonGas<Specie> Foam::operator-
304 (
305  const PengRobinsonGas<Specie>& pg1,
306  const PengRobinsonGas<Specie>& pg2
307 )
308 {
309  scalar nMoles = pg1.nMoles() + pg2.nMoles();
310  scalar molr1 = pg1.nMoles()/nMoles;
311  scalar molr2 = pg2.nMoles()/nMoles;
312 
313  const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
314  const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
315  const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
316 
318  (
319  static_cast<const Specie&>(pg1)
320  - static_cast<const Specie&>(pg2),
321  Tc,
322  Vc,
323  Zc,
324  RR*Zc*Tc/Vc,
325  molr1*pg1.omega_ - molr2*pg2.omega_
326  );
327 }
328 
329 
330 template<class Specie>
331 Foam::PengRobinsonGas<Specie> Foam::operator*
332 (
333  const scalar s,
334  const PengRobinsonGas<Specie>& pg
335 )
336 {
338  (
339  s*static_cast<const Specie&>(pg),
340  pg.Tc_,
341  pg.Vc_,
342  pg.Zc_,
343  pg.Pc_,
344  pg.omega_
345  );
346 }
347 
348 
349 template<class Specie>
350 Foam::PengRobinsonGas<Specie> Foam::operator==
351 (
352  const PengRobinsonGas<Specie>& pg1,
353  const PengRobinsonGas<Specie>& pg2
354 )
355 {
356  return pg2 - pg1;
357 }
358 
359 
360 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar pow3(const dimensionedScalar &ds)
PengRobinsonGas(const Specie &sp, const scalar &Tc, const scalar &Vc, const scalar &Zc, const scalar &Pc, const scalar &omega)
Construct from components.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
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 ))
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
autoPtr< PengRobinsonGas > clone() const
Construct and return a clone.
dimensioned< scalar > mag(const dimensioned< Type > &)
#define R(A, B, C, D, E, F, K, M)
static autoPtr< PengRobinsonGas > New(Istream &is)
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
A class for handling words, derived from string.
Definition: word.H:59
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const volScalarField & T
Definition: createFields.H:25
const dimensionedScalar Pstd
Standard pressure.
scalar s(const scalar p, const scalar T) const
Return entropy [J/(kmol K)].
dictionary dict
dimensionedScalar log(const dimensionedScalar &ds)
const dimensionedScalar a0
Bohr radius: default SI units: [m].
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
scalar Z(scalar p, scalar T) const
Return compression factor [-].
PengRobinsonGas gas equation of state.
const scalar RR
Universal gas constant (default in [J/(kmol K)])
scalar cpMcv(scalar p, scalar T) const
Return (cp - cv) [J/(kmol K].
scalar rho(scalar p, scalar T) const
Return density [kg/m^3].
void operator*=(const scalar)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
This function object calculates and outputs the second invariant of the velocity gradient tensor [1/s...
Definition: Q.H:66
scalar psi(scalar p, scalar T) const
Return compressibility rho/p [s^2/m^2].
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:117