28#ifndef OPM_PENG_ROBINSON_HPP
29#define OPM_PENG_ROBINSON_HPP
54template <
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;
141 template <
class Flu
idState,
class Params>
143 typename FluidState::Scalar
152 typedef typename FluidState::Scalar Evaluation;
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);
165 return std::numeric_limits<Scalar>::quiet_NaN();
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};
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;
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 =
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>
284 const Params& params,
288 Evaluation Tcrit, pcrit, Vcrit;
304 template <
class Evaluation>
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));
351 Evaluation fStar = maxVm - minVm;
356 Evaluation fPrime = (fStar - f)/eps;
359 pcrit = (minP + maxP)/2;
360 Vcrit = (maxVm + minVm)/2;
365 Evaluation delta = f/fPrime;
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;
386 if (
findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
402 template <
class Evaluation>
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;
434 Evaluation V = b*1.1;
435 Evaluation delta = 1.0;
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));
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); }
Provides free functions to invert polynomials of degree 1, 2 and 3.
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
Definition: Exceptions.hpp:46
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 void handleCriticalFluid_(Evaluation &Vm, const FluidState &, const Params ¶ms, unsigned phaseIdx, bool isGasPhase)
Definition: PengRobinson.hpp:282
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 void init(Scalar, Scalar, unsigned, Scalar, Scalar, unsigned)
Definition: PengRobinson.hpp:64
static Evaluation ambroseWalton_(const Params &, const Evaluation &T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition: PengRobinson.hpp:497
static bool findExtrema_(Evaluation &Vmin, Evaluation &Vmax, Evaluation &, Evaluation &, const Evaluation &a, const Evaluation &b, const Evaluation &T)
Definition: PengRobinson.hpp:403
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 void findCriticalPoint_(Evaluation &Tcrit, Evaluation &pcrit, Evaluation &Vcrit, const Evaluation &a, const Evaluation &b)
Definition: PengRobinson.hpp:305
static Evaluation computeFugacity(const Params ¶ms)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:277
bool CheckDefined(const T &value)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition: Valgrind.hpp:74
void SetUndefined(const T &value)
Make the memory on which an object resides undefined in valgrind runs.
Definition: Valgrind.hpp:171
Definition: Air_Mesitylene.hpp:34
bool isfinite(const Evaluation &value)
Definition: MathToolbox.hpp:420
Evaluation exp(const Evaluation &value)
Definition: MathToolbox.hpp:403
Evaluation sqrt(const Evaluation &value)
Definition: MathToolbox.hpp:399
ReturnEval_< Evaluation1, Evaluation2 >::type min(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:346
ReturnEval_< Evaluation1, Evaluation2 >::type max(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:341
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:331
auto scalarValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::scalarValue(val))
Definition: MathToolbox.hpp:335
Evaluation abs(const Evaluation &value)
Definition: MathToolbox.hpp:350
ReturnEval_< Evaluation1, Evaluation2 >::type pow(const Evaluation1 &base, const Evaluation2 &exp)
Definition: MathToolbox.hpp:416