My Project
Loading...
Searching...
No Matches
PengRobinson.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_PENG_ROBINSON_HPP
29#define OPM_PENG_ROBINSON_HPP
30
32
36
38
39#include <cassert>
40#include <csignal>
41#include <string>
42
43namespace Opm {
44
58template <class Scalar, bool UseLegacy=true>
60{
62 static const Scalar R;
63
65 { }
66
67public:
68 static void init(Scalar /*aMin*/, Scalar /*aMax*/, unsigned /*na*/,
69 Scalar /*bMin*/, Scalar /*bMax*/, unsigned /*nb*/)
70 {
71/*
72 // resize the tabulation for the critical points
73 criticalTemperature_.resize(aMin, aMax, na, bMin, bMax, nb);
74 criticalPressure_.resize(aMin, aMax, na, bMin, bMax, nb);
75 criticalMolarVolume_.resize(aMin, aMax, na, bMin, bMax, nb);
76
77 Scalar VmCrit, pCrit, TCrit;
78 for (unsigned i = 0; i < na; ++i) {
79 Scalar a = criticalTemperature_.iToX(i);
80 assert(std::abs(criticalTemperature_.xToI(criticalTemperature_.iToX(i)) - i) < 1e-10);
81
82 for (unsigned j = 0; j < nb; ++j) {
83 Scalar b = criticalTemperature_.jToY(j);
84 assert(std::abs(criticalTemperature_.yToJ(criticalTemperature_.jToY(j)) - j) < 1e-10);
85
86 findCriticalPoint_(TCrit, pCrit, VmCrit, a, b);
87
88 criticalTemperature_.setSamplePoint(i, j, TCrit);
89 criticalPressure_.setSamplePoint(i, j, pCrit);
90 criticalMolarVolume_.setSamplePoint(i, j, VmCrit);
91 }
92 }
93 */
94 }
95
104 template <class Evaluation, class Params>
105 static Evaluation computeVaporPressure(const Params& params, const Evaluation& T)
106 {
107 typedef typename Params::Component Component;
110
111 // initial guess of the vapor pressure
112 Evaluation Vm[3];
113 const Scalar eps = Component::criticalPressure()*1e-10;
114
115 // use the Ambrose-Walton method to get an initial guess of
116 // the vapor pressure
117 Evaluation pVap = ambroseWalton_(params, T);
118
119 // Newton-Raphson method
120 for (unsigned i = 0; i < 5; ++i) {
121 // calculate the molar densities
122 assert(molarVolumes(Vm, params, T, pVap) == 3);
123
124 const Evaluation& f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]);
125 Evaluation df_dp =
126 fugacityDifference_(params, T, pVap + eps, Vm[0], Vm[2])
127 -
128 fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]);
129 df_dp /= 2*eps;
130
131 const Evaluation& delta = f/df_dp;
132 pVap = pVap - delta;
133
134 if (std::abs(scalarValue(delta/pVap)) < 1e-10)
135 break;
136 }
137
138 return pVap;
139 }
140
145 template <class FluidState, class Params>
146 static
147 typename FluidState::Scalar
148 computeMolarVolume(const FluidState& fs,
149 Params& params,
150 unsigned phaseIdx,
151 bool isGasPhase)
152 {
153 Valgrind::CheckDefined(fs.temperature(phaseIdx));
154 Valgrind::CheckDefined(fs.pressure(phaseIdx));
155
156 typedef typename FluidState::Scalar Evaluation;
157
158 Evaluation Vm = 0;
159 Valgrind::SetUndefined(Vm);
160
161 const Evaluation& T = fs.temperature(phaseIdx);
162 const Evaluation& p = fs.pressure(phaseIdx);
163
164 const Evaluation& a = params.a(phaseIdx); // "attractive factor"
165 const Evaluation& b = params.b(phaseIdx); // "co-volume"
166
167 if (!std::isfinite(scalarValue(a))
168 || std::abs(scalarValue(a)) < 1e-30)
169 return std::numeric_limits<Scalar>::quiet_NaN();
170 if (!std::isfinite(scalarValue(b)) || b <= 0)
171 return std::numeric_limits<Scalar>::quiet_NaN();
172
173 const Evaluation& RT= Constants<Scalar>::R*T;
174 const Evaluation& Astar = a*p/(RT*RT);
175 const Evaluation& Bstar = b*p/RT;
176
177 const Evaluation& a1 = 1.0;
178 const Evaluation& a2 = - (1 - Bstar);
179 const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
180 const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
181
182 // ignore the first two results if the smallest
183 // compressibility factor is <= 0.0. (this means that if we
184 // would get negative molar volumes for the liquid phase, we
185 // consider the liquid phase non-existant.)
186 Evaluation Z[3] = {0.0,0.0,0.0};
187 Valgrind::CheckDefined(a1);
188 Valgrind::CheckDefined(a2);
189 Valgrind::CheckDefined(a3);
190 Valgrind::CheckDefined(a4);
191
192 int numSol = cubicRoots(Z, a1, a2, a3, a4);
193 if (numSol == 3) {
194 // the EOS has three intersections with the pressure,
195 // i.e. the molar volume of gas is the largest one and the
196 // molar volume of liquid is the smallest one
197 if (isGasPhase)
198 Vm = max(1e-7, Z[2]*RT/p);
199 else
200 Vm = max(1e-7, Z[0]*RT/p);
201 }
202 else if (numSol == 1) {
203 // the EOS only has one intersection with the pressure,
204 // for the other phase, we take the extremum of the EOS
205 // with the largest distance from the intersection.
206 Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
207 Vm = VmCubic;
208
209 if (UseLegacy) {
210 // find the extrema (if they are present)
211 Evaluation Vmin, Vmax, pmin, pmax;
212 if (findExtrema_(Vmin, Vmax, pmin, pmax, a, b, T)) {
213 if (isGasPhase)
214 Vm = std::max(Vmax, VmCubic);
215 else {
216 if (Vmin > 0)
217 Vm = std::min(Vmin, VmCubic);
218 else
219 Vm = VmCubic;
220 }
221 } else {
222 // the EOS does not exhibit any physically meaningful
223 // extrema, and the fluid is critical...
224 Vm = VmCubic;
225 handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
226 }
227 }
228 }
229
230 Valgrind::CheckDefined(Vm);
231 assert(std::isfinite(scalarValue(Vm)));
232 assert(Vm > 0);
233 return Vm;
234 }
235
246 template <class Evaluation, class Params>
247 static Evaluation computeFugacityCoeffient(const Params& params)
248 {
249 const Evaluation& T = params.temperature();
250 const Evaluation& p = params.pressure();
251 const Evaluation& Vm = params.molarVolume();
252
253 const Evaluation& RT = Constants<Scalar>::R*T;
254 const Evaluation& Z = p*Vm/RT;
255 const Evaluation& Bstar = p*params.b() / RT;
256
257 const Evaluation& tmp =
258 (Vm + params.b()*(1 + std::sqrt(2))) /
259 (Vm + params.b()*(1 - std::sqrt(2)));
260 const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
261 const Evaluation& fugCoeff =
262 exp(Z - 1) / (Z - Bstar) *
263 pow(tmp, expo);
264
265 return fugCoeff;
266 }
267
278 template <class Evaluation, class Params>
279 static Evaluation computeFugacity(const Params& params)
280 { return params.pressure()*computeFugacityCoeff(params); }
281
282protected:
283 template <class FluidState, class Params, class Evaluation = typename FluidState::Scalar>
284 static void handleCriticalFluid_(Evaluation& Vm,
285 const FluidState& /*fs*/,
286 const Params& params,
287 unsigned phaseIdx,
288 bool isGasPhase)
289 {
290 Evaluation Tcrit, pcrit, Vcrit;
291 findCriticalPoint_(Tcrit,
292 pcrit,
293 Vcrit,
294 Evaluation(params.a(phaseIdx)),
295 Evaluation(params.b(phaseIdx)) );
296
297
298 //Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx));
299
300 if (isGasPhase)
301 Vm = max(Vm, Vcrit);
302 else
303 Vm = min(Vm, Vcrit);
304 }
305
306 template <class Evaluation>
307 static void findCriticalPoint_(Evaluation& Tcrit,
308 Evaluation& pcrit,
309 Evaluation& Vcrit,
310 const Evaluation& a,
311 const Evaluation& b)
312 {
313 Evaluation minVm(0);
314 Evaluation maxVm(1e30);
315
316 Evaluation minP(0);
317 Evaluation maxP(1e30);
318
319 // first, we need to find an isotherm where the EOS exhibits
320 // a maximum and a minimum
321 Evaluation Tmin = 250; // [K]
322 for (unsigned i = 0; i < 30; ++i) {
323 bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
324 if (hasExtrema)
325 break;
326 Tmin /= 2;
327 };
328
329 Evaluation T = Tmin;
330
331 // Newton's method: Start at minimum temperature and minimize
332 // the "gap" between the extrema of the EOS
333 unsigned iMax = 100;
334 for (unsigned i = 0; i < iMax; ++i) {
335 // calculate function we would like to minimize
336 Evaluation f = maxVm - minVm;
337
338 // check if we're converged
339 if (f < 1e-10 || (i == iMax - 1 && f < 1e-8)) {
340 Tcrit = T;
341 pcrit = (minP + maxP)/2;
342 Vcrit = (maxVm + minVm)/2;
343 return;
344 }
345
346 // backward differences. Forward differences are not
347 // robust, since the isotherm could be critical if an
348 // epsilon was added to the temperature. (this is case
349 // rarely happens, though)
350 const Scalar eps = - 1e-11;
351 [[maybe_unused]] bool found = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
352 assert(found);
353 assert(std::isfinite(scalarValue(maxVm)));
354 Evaluation fStar = maxVm - minVm;
355
356 // derivative of the difference between the maximum's
357 // molar volume and the minimum's molar volume regarding
358 // temperature
359 Evaluation fPrime = (fStar - f)/eps;
360 if (std::abs(scalarValue(fPrime)) < 1e-40) {
361 Tcrit = T;
362 pcrit = (minP + maxP)/2;
363 Vcrit = (maxVm + minVm)/2;
364 return;
365 }
366
367 // update value for the current iteration
368 Evaluation delta = f/fPrime;
369 assert(std::isfinite(scalarValue(delta)));
370 if (delta > 0)
371 delta = -10;
372
373 // line search (we have to make sure that both extrema
374 // still exist after the update)
375 for (unsigned j = 0; ; ++j) {
376 if (j >= 20) {
377 if (f < 1e-8) {
378 Tcrit = T;
379 pcrit = (minP + maxP)/2;
380 Vcrit = (maxVm + minVm)/2;
381 return;
382 }
383
384 const std::string msg =
385 "Could not determine the critical point for a="
386 + std::to_string(getValue(a))
387 + ", b=" + std::to_string(getValue(b));
388 throw NumericalProblem(msg);
389 }
390
391 if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
392 // if the isotherm for T - delta exhibits two
393 // extrema the update is finished
394 T -= delta;
395 break;
396 }
397 else
398 delta /= 2;
399 };
400 }
401
402 assert(false);
403 }
404
405 // find the two molar volumes where the EOS exhibits extrema and
406 // which are larger than the covolume of the phase
407 template <class Evaluation>
408 static bool findExtrema_(Evaluation& Vmin,
409 Evaluation& Vmax,
410 Evaluation& /*pMin*/,
411 Evaluation& /*pMax*/,
412 const Evaluation& a,
413 const Evaluation& b,
414 const Evaluation& T)
415 {
416 Scalar u = 2;
417 Scalar w = -1;
418
419 const Evaluation& RT = Constants<Scalar>::R*T;
420 // calculate coefficients of the 4th order polynominal in
421 // monomial basis
422 const Evaluation& a1 = RT;
423 const Evaluation& a2 = 2*RT*u*b - 2*a;
424 const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
425 const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
426 const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
427
428 assert(std::isfinite(scalarValue(a1)));
429 assert(std::isfinite(scalarValue(a2)));
430 assert(std::isfinite(scalarValue(a3)));
431 assert(std::isfinite(scalarValue(a4)));
432 assert(std::isfinite(scalarValue(a5)));
433
434 // Newton method to find first root
435
436 // if the values which we got on Vmin and Vmax are usefull, we
437 // will reuse them as initial value, else we will start 10%
438 // above the covolume
439 Evaluation V = b*1.1;
440 Evaluation delta = 1.0;
441 for (unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
442 const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
443 const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
444
445 if (std::abs(scalarValue(fPrime)) < 1e-20) {
446 // give up if the derivative is zero
447 return false;
448 }
449
450
451 delta = f/fPrime;
452 V -= delta;
453
454 if (i > 200) {
455 // give up after 200 iterations...
456 return false;
457 }
458 }
459 assert(std::isfinite(scalarValue(V)));
460
461 // polynomial division
462 Evaluation b1 = a1;
463 Evaluation b2 = a2 + V*b1;
464 Evaluation b3 = a3 + V*b2;
465 Evaluation b4 = a4 + V*b3;
466
467 // invert resulting cubic polynomial analytically
468 std::array<Evaluation,4> allV;
469 allV[0] = V;
470 const decltype(allV.size()) numSol = 1 + invertCubicPolynomial<Evaluation>(allV.data() + 1, b1, b2, b3, b4);
471
472 assert (numSol <= allV.size());
473
474 // sort all roots of the derivative
475 std::sort(allV.begin(), allV.begin() + numSol);
476
477 // check whether the result is physical
478 if (numSol != 4 || allV[numSol - 2] < b) {
479 // the second largest extremum is smaller than the phase's
480 // covolume which is physically impossible
481 return false;
482 }
483
484
485 // it seems that everything is okay...
486 Vmin = allV[numSol - 2];
487 Vmax = allV[numSol - 1];
488 return true;
489 }
490
503 template <class Evaluation, class Params>
504 static Evaluation ambroseWalton_(const Params& /*params*/, const Evaluation& T)
505 {
506 typedef typename Params::Component Component;
507
508 const Evaluation& Tr = T / Component::criticalTemperature();
509 const Evaluation& tau = 1 - Tr;
510 const Evaluation& omega = Component::acentricFactor();
511
512 const Evaluation& f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
513 const Evaluation& f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
514 const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
515
516 return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
517 }
518
529 template <class Evaluation, class Params>
530 static Evaluation fugacityDifference_(const Params& params,
531 const Evaluation& T,
532 const Evaluation& p,
533 const Evaluation& VmLiquid,
534 const Evaluation& VmGas)
535 { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
536
537
538};
539
540} // namespace Opm
541
542#endif
Provides the OPM specific exception classes.
Relations valid for an ideal gas.
Provides free functions to invert polynomials of degree 1, 2 and 3.
This is a fluid state which allows to set the fluid temperatures and takes all other quantities from ...
Implements a scalar function that depends on two variables and which is sampled on an uniform X-Y gri...
Abstract base class of a pure chemical species.
Definition Component.hpp:44
static Scalar acentricFactor()
Returns the acentric factor of the component.
Definition Component.hpp:112
static Scalar criticalPressure()
Returns the critical pressure in of the component.
Definition Component.hpp:105
static Scalar criticalTemperature()
Returns the critical temperature in of the component.
Definition Component.hpp:99
A central place for various physical constants occuring in some equations.
Definition Constants.hpp:41
static const Scalar R
The ideal gas constant [J/(mol K)].
Definition Constants.hpp:45
Implements the Peng-Robinson equation of state for liquids and gases.
Definition PengRobinson.hpp:60
static Evaluation computeVaporPressure(const Params &params, const Evaluation &T)
Predicts the vapor pressure for the temperature given in setTP().
Definition PengRobinson.hpp:105
static Evaluation fugacityDifference_(const Params &params, const Evaluation &T, const Evaluation &p, const Evaluation &VmLiquid, const Evaluation &VmGas)
Returns the difference between the liquid and the gas phase fugacities in [bar].
Definition PengRobinson.hpp:530
static FluidState::Scalar computeMolarVolume(const FluidState &fs, Params &params, unsigned phaseIdx, bool isGasPhase)
Computes molar volumes where the Peng-Robinson EOS is true.
Definition PengRobinson.hpp:148
static Evaluation computeFugacity(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition PengRobinson.hpp:279
static Evaluation computeFugacityCoeffient(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition PengRobinson.hpp:247
static Evaluation ambroseWalton_(const Params &, const Evaluation &T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition PengRobinson.hpp:504
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition PolynomialUtils.hpp:331