28 #ifndef OPM_PENG_ROBINSON_HPP
29 #define OPM_PENG_ROBINSON_HPP
54 template <
class Scalar>
58 static const Scalar R;
64 static void init(Scalar , Scalar ,
unsigned ,
65 Scalar , Scalar ,
unsigned )
100 template <
class Evaluation,
class Params>
103 typedef typename Params::Component
Component;
116 for (
unsigned i = 0; i < 5; ++i) {
118 assert(molarVolumes(Vm, params, T, pVap) == 3);
127 const Evaluation& delta = f/df_dp;
130 if (std::abs(scalarValue(delta/pVap)) < 1e-10)
141 template <
class Flu
idState,
class Params>
143 typename FluidState::Scalar
149 Valgrind::CheckDefined(fs.temperature(phaseIdx));
150 Valgrind::CheckDefined(fs.pressure(phaseIdx));
152 typedef typename FluidState::Scalar Evaluation;
155 Valgrind::SetUndefined(Vm);
157 const Evaluation& T = fs.temperature(phaseIdx);
158 const Evaluation& p = fs.pressure(phaseIdx);
160 const Evaluation& a = params.a(phaseIdx);
161 const Evaluation& b = params.b(phaseIdx);
163 if (!std::isfinite(scalarValue(a))
164 || std::abs(scalarValue(a)) < 1e-30)
165 return std::numeric_limits<Scalar>::quiet_NaN();
166 if (!std::isfinite(scalarValue(b)) || b <= 0)
167 return std::numeric_limits<Scalar>::quiet_NaN();
170 const Evaluation& Astar = a*p/(RT*RT);
171 const Evaluation& Bstar = b*p/RT;
173 const Evaluation& a1 = 1.0;
174 const Evaluation& a2 = - (1 - Bstar);
175 const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
176 const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
182 Evaluation Z[3] = {0.0,0.0,0.0};
183 Valgrind::CheckDefined(a1);
184 Valgrind::CheckDefined(a2);
185 Valgrind::CheckDefined(a3);
186 Valgrind::CheckDefined(a4);
194 Vm = max(1e-7, Z[2]*RT/p);
196 Vm = max(1e-7, Z[0]*RT/p);
198 else if (numSol == 1) {
202 Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
206 Evaluation Vmin, Vmax, pmin, pmax;
207 if (findExtrema_(Vmin, Vmax,
212 Vm = std::max(Vmax, VmCubic);
215 Vm = std::min(Vmin, VmCubic);
224 handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
228 Valgrind::CheckDefined(Vm);
229 assert(std::isfinite(scalarValue(Vm)));
244 template <
class Evaluation,
class Params>
247 const Evaluation& T = params.temperature();
248 const Evaluation& p = params.pressure();
249 const Evaluation& Vm = params.molarVolume();
252 const Evaluation& Z = p*Vm/RT;
253 const Evaluation& Bstar = p*params.b() / RT;
255 const Evaluation& tmp =
256 (Vm + params.b()*(1 + std::sqrt(2))) /
257 (Vm + params.b()*(1 - std::sqrt(2)));
258 const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
259 const Evaluation& fugCoeff =
260 exp(Z - 1) / (Z - Bstar) *
276 template <
class Evaluation,
class Params>
278 {
return params.pressure()*computeFugacityCoeff(params); }
281 template <
class Flu
idState,
class Params,
class Evaluation =
typename Flu
idState::Scalar>
282 static void handleCriticalFluid_(Evaluation& Vm,
284 const Params& params,
288 Evaluation Tcrit, pcrit, Vcrit;
289 findCriticalPoint_(Tcrit,
304 template <
class Evaluation>
305 static void findCriticalPoint_(Evaluation& Tcrit,
312 Evaluation maxVm(1e30);
315 Evaluation maxP(1e30);
319 Evaluation Tmin = 250;
320 for (
unsigned i = 0; i < 30; ++i) {
321 bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
332 for (
unsigned i = 0; i < iMax; ++i) {
334 Evaluation f = maxVm - minVm;
337 if (f < 1e-10 || (i == iMax - 1 && f < 1e-8)) {
339 pcrit = (minP + maxP)/2;
340 Vcrit = (maxVm + minVm)/2;
348 const Scalar eps = - 1e-11;
349 assert(findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps));
350 assert(std::isfinite(scalarValue(maxVm)));
351 Evaluation fStar = maxVm - minVm;
356 Evaluation fPrime = (fStar - f)/eps;
357 if (std::abs(scalarValue(fPrime)) < 1e-40) {
359 pcrit = (minP + maxP)/2;
360 Vcrit = (maxVm + minVm)/2;
365 Evaluation delta = f/fPrime;
366 assert(std::isfinite(scalarValue(delta)));
372 for (
unsigned j = 0; ; ++j) {
376 pcrit = (minP + maxP)/2;
377 Vcrit = (maxVm + minVm)/2;
381 std::ostringstream oss;
382 oss <<
"Could not determine the critical point for a=" << a <<
", b=" << b;
383 throw NumericalIssue(oss.str());
386 if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
402 template <
class Evaluation>
403 static bool findExtrema_(Evaluation& Vmin,
417 const Evaluation& a1 = RT;
418 const Evaluation& a2 = 2*RT*u*b - 2*a;
419 const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
420 const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
421 const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
423 assert(std::isfinite(scalarValue(a1)));
424 assert(std::isfinite(scalarValue(a2)));
425 assert(std::isfinite(scalarValue(a3)));
426 assert(std::isfinite(scalarValue(a4)));
427 assert(std::isfinite(scalarValue(a5)));
434 Evaluation V = b*1.1;
435 Evaluation delta = 1.0;
436 for (
unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
437 const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
438 const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
440 if (std::abs(scalarValue(fPrime)) < 1e-20) {
454 assert(std::isfinite(scalarValue(V)));
458 Evaluation b2 = a2 + V*b1;
459 Evaluation b3 = a3 + V*b2;
460 Evaluation b4 = a4 + V*b3;
465 int numSol = 1 + invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
468 std::sort(allV + 0, allV + numSol);
471 if (numSol != 4 || allV[numSol - 2] < b) {
479 Vmin = allV[numSol - 2];
480 Vmax = allV[numSol - 1];
496 template <
class Evaluation,
class Params>
499 typedef typename Params::Component
Component;
502 const Evaluation& tau = 1 - Tr;
505 const Evaluation& f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
506 const Evaluation& f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
507 const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
522 template <
class Evaluation,
class Params>
526 const Evaluation& VmLiquid,
527 const Evaluation& VmGas)
528 {
return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
Relations valid for an ideal gas.
Provides free functions to invert polynomials of degree 1, 2 and 3.
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:331
This is a fluid state which allows to set the fluid temperatures and takes all other quantities from ...
Abstract base class of a pure chemical species.
Definition: Component.hpp:42
static Scalar acentricFactor()
Returns the acentric factor of the component.
Definition: Component.hpp:110
static Scalar criticalPressure()
Returns the critical pressure in of the component.
Definition: Component.hpp:103
static Scalar criticalTemperature()
Returns the critical temperature in of the component.
Definition: Component.hpp:97
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:56
static Evaluation computeFugacityCoeffient(const Params ¶ms)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:245
static Evaluation fugacityDifference_(const Params ¶ms, 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:523
static Evaluation computeVaporPressure(const Params ¶ms, const Evaluation &T)
Predicts the vapor pressure for the temperature given in setTP().
Definition: PengRobinson.hpp:101
static Evaluation ambroseWalton_(const Params &, const Evaluation &T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition: PengRobinson.hpp:497
static FluidState::Scalar computeMolarVolume(const FluidState &fs, Params ¶ms, unsigned phaseIdx, bool isGasPhase)
Computes molar volumes where the Peng-Robinson EOS is true.
Definition: PengRobinson.hpp:144
static Evaluation computeFugacity(const Params ¶ms)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:277