28 #ifndef OPM_PENG_ROBINSON_HPP 29 #define OPM_PENG_ROBINSON_HPP 35 #include <opm/common/Unused.hpp> 55 template <
class Scalar>
59 static const Scalar R;
65 static void init(Scalar , Scalar ,
unsigned ,
66 Scalar , Scalar ,
unsigned )
101 template <
class Evaluation,
class Params>
104 typedef typename Params::Component
Component;
117 for (
unsigned i = 0; i < 5; ++i) {
119 int numSol OPM_OPTIM_UNUSED = molarVolumes(Vm, params, T, pVap);
129 const Evaluation& delta = f/df_dp;
132 if (std::abs(Opm::scalarValue(delta/pVap)) < 1e-10)
143 template <
class Flu
idState,
class Params>
145 typename FluidState::Scalar
151 Valgrind::CheckDefined(fs.temperature(phaseIdx));
152 Valgrind::CheckDefined(fs.pressure(phaseIdx));
154 typedef typename FluidState::Scalar Evaluation;
157 Valgrind::SetUndefined(Vm);
159 const Evaluation& T = fs.temperature(phaseIdx);
160 const Evaluation& p = fs.pressure(phaseIdx);
162 const Evaluation& a = params.a(phaseIdx);
163 const Evaluation& b = params.b(phaseIdx);
165 if (!std::isfinite(Opm::scalarValue(a))
166 || std::abs(Opm::scalarValue(a)) < 1e-30)
167 return std::numeric_limits<Scalar>::quiet_NaN();
168 if (!std::isfinite(Opm::scalarValue(b)) || b <= 0)
169 return std::numeric_limits<Scalar>::quiet_NaN();
171 const Evaluation& RT= R*T;
172 const Evaluation& Astar = a*p/(RT*RT);
173 const Evaluation& Bstar = b*p/RT;
175 const Evaluation& a1 = 1.0;
176 const Evaluation& a2 = - (1 - Bstar);
177 const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
178 const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
184 Evaluation Z[3] = {0.0,0.0,0.0};
185 Valgrind::CheckDefined(a1);
186 Valgrind::CheckDefined(a2);
187 Valgrind::CheckDefined(a3);
188 Valgrind::CheckDefined(a4);
199 else if (numSol == 1) {
203 Evaluation VmCubic = Z[0]*RT/p;
207 Evaluation Vmin, Vmax, pmin, pmax;
208 if (findExtrema_(Vmin, Vmax,
213 Vm = std::max(Vmax, VmCubic);
216 Vm = std::min(Vmin, VmCubic);
225 handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
229 Valgrind::CheckDefined(Vm);
230 assert(std::isfinite(Opm::scalarValue(Vm)));
245 template <
class Evaluation,
class Params>
248 const Evaluation& T = params.temperature();
249 const Evaluation& p = params.pressure();
250 const Evaluation& Vm = params.molarVolume();
252 const Evaluation& RT = R*T;
253 const Evaluation& Z = p*Vm/RT;
254 const Evaluation& Bstar = p*params.b() / RT;
256 const Evaluation& tmp =
257 (Vm + params.b()*(1 + std::sqrt(2))) /
258 (Vm + params.b()*(1 - std::sqrt(2)));
259 const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
260 const Evaluation& fugCoeff =
261 Opm::exp(Z - 1) / (Z - Bstar) *
277 template <
class Evaluation,
class Params>
279 {
return params.pressure()*computeFugacityCoeff(params); }
282 template <
class Flu
idState,
class Params,
class Evaluation =
typename Flu
idState::Scalar>
283 static void handleCriticalFluid_(Evaluation& Vm,
285 const Params& params,
289 Evaluation Tcrit, pcrit, Vcrit;
290 findCriticalPoint_(Tcrit,
300 Vm = Opm::max(Vm, Vcrit);
302 Vm = Opm::min(Vm, Vcrit);
305 template <
class Evaluation>
306 static void findCriticalPoint_(Evaluation& Tcrit,
313 Evaluation maxVm(1e30);
316 Evaluation maxP(1e30);
320 Evaluation Tmin = 250;
321 for (
unsigned i = 0; i < 30; ++i) {
322 bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
333 for (
unsigned i = 0; i < iMax; ++i) {
335 Evaluation f = maxVm - minVm;
338 if (f < 1e-10 || (i == iMax - 1 && f < 1e-8)) {
340 pcrit = (minP + maxP)/2;
341 Vcrit = (maxVm + minVm)/2;
349 const Scalar eps = - 1e-11;
350 bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
352 assert(std::isfinite(Opm::scalarValue(maxVm)));
353 Evaluation fStar = maxVm - minVm;
358 Evaluation fPrime = (fStar - f)/eps;
359 if (std::abs(Opm::scalarValue(fPrime)) < 1e-40) {
361 pcrit = (minP + maxP)/2;
362 Vcrit = (maxVm + minVm)/2;
367 Evaluation delta = f/fPrime;
368 assert(std::isfinite(Opm::scalarValue(delta)));
374 for (
unsigned j = 0; ; ++j) {
378 pcrit = (minP + maxP)/2;
379 Vcrit = (maxVm + minVm)/2;
382 OPM_THROW(NumericalProblem,
383 "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>
403 static bool findExtrema_(Evaluation& Vmin,
414 const Evaluation& RT = R*T;
418 const Evaluation& a1 = RT;
419 const Evaluation& a2 = 2*RT*u*b - 2*a;
420 const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
421 const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
422 const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
424 assert(std::isfinite(Opm::scalarValue(a1)));
425 assert(std::isfinite(Opm::scalarValue(a2)));
426 assert(std::isfinite(Opm::scalarValue(a3)));
427 assert(std::isfinite(Opm::scalarValue(a4)));
428 assert(std::isfinite(Opm::scalarValue(a5)));
435 Evaluation V = b*1.1;
436 Evaluation delta = 1.0;
437 for (
unsigned i = 0; std::abs(Opm::scalarValue(delta)) > 1e-12; ++i) {
438 const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
439 const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
441 if (std::abs(Opm::scalarValue(fPrime)) < 1e-20) {
455 assert(std::isfinite(Opm::scalarValue(V)));
459 Evaluation b2 = a2 + V*b1;
460 Evaluation b3 = a3 + V*b2;
461 Evaluation b4 = a4 + V*b3;
466 int numSol = 1 + Opm::invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
469 std::sort(allV + 0, allV + numSol);
472 if (numSol != 4 || allV[numSol - 2] < b) {
480 Vmin = allV[numSol - 2];
481 Vmax = allV[numSol - 1];
497 template <
class Evaluation,
class Params>
500 typedef typename Params::Component
Component;
503 const Evaluation& tau = 1 - Tr;
504 const Evaluation& omega = Component::acentricFactor();
506 const Evaluation& f0 = (tau*(-5.97616 + Opm::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Opm::pow(tau, 5))/Tr;
507 const Evaluation& f1 = (tau*(-5.03365 + Opm::sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*Opm::pow(tau, 5))/Tr;
508 const Evaluation& f2 = (tau*(-0.64771 + Opm::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Opm::pow(tau, 5))/Tr;
523 template <
class Evaluation,
class Params>
527 const Evaluation& VmLiquid,
528 const Evaluation& VmGas)
529 {
return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
538 template <
class Scalar>
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:524
Implements the Peng-Robinson equation of state for liquids and gases.
Definition: PengRobinson.hpp:56
This is a fluid state which allows to set the fluid temperatures and takes all other quantities from ...
Relations valid for an ideal gas.
Abstract base class of a pure chemical species.
Definition: Component.hpp:43
static Scalar criticalPressure()
Returns the critical pressure in of the component.
Definition: Component.hpp:105
Provides free functions to invert polynomials of degree 1, 2 and 3.
Definition: Air_Mesitylene.hpp:33
static Evaluation computeVaporPressure(const Params ¶ms, const Evaluation &T)
Predicts the vapor pressure for the temperature given in setTP().
Definition: PengRobinson.hpp:102
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
A central place for various physical constants occuring in some equations.
Definition: Constants.hpp:40
static Evaluation computeFugacity(const Params ¶ms)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:278
static Evaluation ambroseWalton_(const Params &, const Evaluation &T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition: PengRobinson.hpp:498
static Scalar criticalTemperature()
Returns the critical temperature in of the component.
Definition: Component.hpp:99
static Evaluation computeFugacityCoeffient(const Params ¶ms)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:246
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:146