61 static const Scalar R;
67 static void init(Scalar , Scalar ,
unsigned ,
68 Scalar , Scalar ,
unsigned )
103 template <
class Evaluation,
class Params>
106 typedef typename Params::Component
Component;
119 for (
unsigned i = 0; i < 5; ++i) {
121 assert(molarVolumes(Vm, params, T, pVap) == 3);
130 const Evaluation& delta = f/df_dp;
133 if (std::abs(scalarValue(delta/pVap)) < 1e-10)
144 template <
class Flu
idState,
class Params>
146 typename FluidState::Scalar
152 Valgrind::CheckDefined(fs.temperature(phaseIdx));
153 Valgrind::CheckDefined(fs.pressure(phaseIdx));
155 typedef typename FluidState::Scalar Evaluation;
158 Valgrind::SetUndefined(Vm);
160 const Evaluation& T = fs.temperature(phaseIdx);
161 const Evaluation& p = fs.pressure(phaseIdx);
163 const Evaluation& a = params.a(phaseIdx);
164 const Evaluation& b = params.b(phaseIdx);
166 if (!std::isfinite(scalarValue(a))
167 || std::abs(scalarValue(a)) < 1e-30)
168 return std::numeric_limits<Scalar>::quiet_NaN();
169 if (!std::isfinite(scalarValue(b)) || b <= 0)
170 return std::numeric_limits<Scalar>::quiet_NaN();
173 const Evaluation& Astar = a*p/(RT*RT);
174 const Evaluation& Bstar = b*p/RT;
176 const Evaluation& a1 = 1.0;
177 const Evaluation& a2 = - (1 - Bstar);
178 const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
179 const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
185 Evaluation Z[3] = {0.0,0.0,0.0};
186 Valgrind::CheckDefined(a1);
187 Valgrind::CheckDefined(a2);
188 Valgrind::CheckDefined(a3);
189 Valgrind::CheckDefined(a4);
197 Vm = max(1e-7, Z[2]*RT/p);
199 Vm = max(1e-7, Z[0]*RT/p);
201 else if (numSol == 1) {
205 Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
210 Evaluation Vmin, Vmax, pmin, pmax;
211 if (findExtrema_(Vmin, Vmax, pmin, pmax, a, b, T)) {
213 Vm = std::max(Vmax, VmCubic);
216 Vm = std::min(Vmin, VmCubic);
224 handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
229 Valgrind::CheckDefined(Vm);
230 assert(std::isfinite(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();
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 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,
293 Evaluation(params.a(phaseIdx)),
294 Evaluation(params.b(phaseIdx)) );
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 assert(findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps));
351 assert(std::isfinite(scalarValue(maxVm)));
352 Evaluation fStar = maxVm - minVm;
357 Evaluation fPrime = (fStar - f)/eps;
358 if (std::abs(scalarValue(fPrime)) < 1e-40) {
360 pcrit = (minP + maxP)/2;
361 Vcrit = (maxVm + minVm)/2;
366 Evaluation delta = f/fPrime;
367 assert(std::isfinite(scalarValue(delta)));
373 for (
unsigned j = 0; ; ++j) {
377 pcrit = (minP + maxP)/2;
378 Vcrit = (maxVm + minVm)/2;
382 const std::string msg =
383 "Could not determine the critical point for a="
384 + std::to_string(getValue(a))
385 +
", b=" + std::to_string(getValue(b));
386 throw NumericalProblem(msg);
389 if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
405 template <
class Evaluation>
406 static bool findExtrema_(Evaluation& Vmin,
420 const Evaluation& a1 = RT;
421 const Evaluation& a2 = 2*RT*u*b - 2*a;
422 const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
423 const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
424 const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
426 assert(std::isfinite(scalarValue(a1)));
427 assert(std::isfinite(scalarValue(a2)));
428 assert(std::isfinite(scalarValue(a3)));
429 assert(std::isfinite(scalarValue(a4)));
430 assert(std::isfinite(scalarValue(a5)));
437 Evaluation V = b*1.1;
438 Evaluation delta = 1.0;
439 for (
unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
440 const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
441 const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
443 if (std::abs(scalarValue(fPrime)) < 1e-20) {
457 assert(std::isfinite(scalarValue(V)));
461 Evaluation b2 = a2 + V*b1;
462 Evaluation b3 = a3 + V*b2;
463 Evaluation b4 = a4 + V*b3;
466 std::array<Evaluation,4> allV;
468 int numSol = 1 + invertCubicPolynomial<Evaluation>(allV.data() + 1, b1, b2, b3, b4);
471 std::sort(allV.begin(), allV.begin() + numSol);
474 if (numSol != 4 || allV[numSol - 2] < b) {
482 Vmin = allV[numSol - 2];
483 Vmax = allV[numSol - 1];
499 template <
class Evaluation,
class Params>
502 typedef typename Params::Component
Component;
505 const Evaluation& tau = 1 - Tr;
508 const Evaluation& f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
509 const Evaluation& f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
510 const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
525 template <
class Evaluation,
class Params>
529 const Evaluation& VmLiquid,
530 const Evaluation& VmGas)
531 {
return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
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:526
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition PolynomialUtils.hpp:331