00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00028 #ifndef OPM_PENG_ROBINSON_HPP
00029 #define OPM_PENG_ROBINSON_HPP
00030
00031 #include <opm/material/fluidstates/TemperatureOverlayFluidState.hpp>
00032 #include <opm/material/IdealGas.hpp>
00033 #include <opm/material/common/UniformTabulated2DFunction.hpp>
00034
00035 #include <opm/common/Unused.hpp>
00036 #include <opm/material/common/PolynomialUtils.hpp>
00037
00038 #include <csignal>
00039
00040 namespace Opm {
00041
00055 template <class Scalar>
00056 class PengRobinson
00057 {
00059 static const Scalar R;
00060
00061 PengRobinson()
00062 { }
00063
00064 public:
00065 static void init(Scalar , Scalar , unsigned ,
00066 Scalar , Scalar , unsigned )
00067 {
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 }
00092
00101 template <class Evaluation, class Params>
00102 static Evaluation computeVaporPressure(const Params& params, const Evaluation& T)
00103 {
00104 typedef typename Params::Component Component;
00105 if (T >= Component::criticalTemperature())
00106 return Component::criticalPressure();
00107
00108
00109 Evaluation Vm[3];
00110 const Scalar eps = Component::criticalPressure()*1e-10;
00111
00112
00113
00114 Evaluation pVap = ambroseWalton_(params, T);
00115
00116
00117 for (unsigned i = 0; i < 5; ++i) {
00118
00119 int numSol OPM_OPTIM_UNUSED = molarVolumes(Vm, params, T, pVap);
00120 assert(numSol == 3);
00121
00122 const Evaluation& f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]);
00123 Evaluation df_dp =
00124 fugacityDifference_(params, T, pVap + eps, Vm[0], Vm[2])
00125 -
00126 fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]);
00127 df_dp /= 2*eps;
00128
00129 const Evaluation& delta = f/df_dp;
00130 pVap = pVap - delta;
00131
00132 if (std::abs(Opm::scalarValue(delta/pVap)) < 1e-10)
00133 break;
00134 }
00135
00136 return pVap;
00137 }
00138
00143 template <class FluidState, class Params>
00144 static
00145 typename FluidState::Scalar
00146 computeMolarVolume(const FluidState& fs,
00147 Params& params,
00148 unsigned phaseIdx,
00149 bool isGasPhase)
00150 {
00151 Valgrind::CheckDefined(fs.temperature(phaseIdx));
00152 Valgrind::CheckDefined(fs.pressure(phaseIdx));
00153
00154 typedef typename FluidState::Scalar Evaluation;
00155
00156 Evaluation Vm = 0;
00157 Valgrind::SetUndefined(Vm);
00158
00159 const Evaluation& T = fs.temperature(phaseIdx);
00160 const Evaluation& p = fs.pressure(phaseIdx);
00161
00162 const Evaluation& a = params.a(phaseIdx);
00163 const Evaluation& b = params.b(phaseIdx);
00164
00165 if (!std::isfinite(Opm::scalarValue(a))
00166 || std::abs(Opm::scalarValue(a)) < 1e-30)
00167 return std::numeric_limits<Scalar>::quiet_NaN();
00168 if (!std::isfinite(Opm::scalarValue(b)) || b <= 0)
00169 return std::numeric_limits<Scalar>::quiet_NaN();
00170
00171 const Evaluation& RT= R*T;
00172 const Evaluation& Astar = a*p/(RT*RT);
00173 const Evaluation& Bstar = b*p/RT;
00174
00175 const Evaluation& a1 = 1.0;
00176 const Evaluation& a2 = - (1 - Bstar);
00177 const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
00178 const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
00179
00180
00181
00182
00183
00184 Evaluation Z[3] = {0.0,0.0,0.0};
00185 Valgrind::CheckDefined(a1);
00186 Valgrind::CheckDefined(a2);
00187 Valgrind::CheckDefined(a3);
00188 Valgrind::CheckDefined(a4);
00189 int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4);
00190 if (numSol == 3) {
00191
00192
00193
00194 if (isGasPhase)
00195 Vm = Z[2]*RT/p;
00196 else
00197 Vm = Z[0]*RT/p;
00198 }
00199 else if (numSol == 1) {
00200
00201
00202
00203 Evaluation VmCubic = Z[0]*RT/p;
00204 Vm = VmCubic;
00205
00206
00207 Evaluation Vmin, Vmax, pmin, pmax;
00208 if (findExtrema_(Vmin, Vmax,
00209 pmin, pmax,
00210 a, b, T))
00211 {
00212 if (isGasPhase)
00213 Vm = std::max(Vmax, VmCubic);
00214 else {
00215 if (Vmin > 0)
00216 Vm = std::min(Vmin, VmCubic);
00217 else
00218 Vm = VmCubic;
00219 }
00220 }
00221 else {
00222
00223
00224 Vm = VmCubic;
00225 handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
00226 }
00227 }
00228
00229 Valgrind::CheckDefined(Vm);
00230 assert(std::isfinite(Opm::scalarValue(Vm)));
00231 assert(Vm > 0);
00232 return Vm;
00233 }
00234
00245 template <class Evaluation, class Params>
00246 static Evaluation computeFugacityCoeffient(const Params& params)
00247 {
00248 const Evaluation& T = params.temperature();
00249 const Evaluation& p = params.pressure();
00250 const Evaluation& Vm = params.molarVolume();
00251
00252 const Evaluation& RT = R*T;
00253 const Evaluation& Z = p*Vm/RT;
00254 const Evaluation& Bstar = p*params.b() / RT;
00255
00256 const Evaluation& tmp =
00257 (Vm + params.b()*(1 + std::sqrt(2))) /
00258 (Vm + params.b()*(1 - std::sqrt(2)));
00259 const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
00260 const Evaluation& fugCoeff =
00261 Opm::exp(Z - 1) / (Z - Bstar) *
00262 Opm::pow(tmp, expo);
00263
00264 return fugCoeff;
00265 }
00266
00277 template <class Evaluation, class Params>
00278 static Evaluation computeFugacity(const Params& params)
00279 { return params.pressure()*computeFugacityCoeff(params); }
00280
00281 protected:
00282 template <class FluidState, class Params, class Evaluation = typename FluidState::Scalar>
00283 static void handleCriticalFluid_(Evaluation& Vm,
00284 const FluidState& ,
00285 const Params& params,
00286 unsigned phaseIdx,
00287 bool isGasPhase)
00288 {
00289 Evaluation Tcrit, pcrit, Vcrit;
00290 findCriticalPoint_(Tcrit,
00291 pcrit,
00292 Vcrit,
00293 params.a(phaseIdx),
00294 params.b(phaseIdx));
00295
00296
00297
00298
00299 if (isGasPhase)
00300 Vm = Opm::max(Vm, Vcrit);
00301 else
00302 Vm = Opm::min(Vm, Vcrit);
00303 }
00304
00305 template <class Evaluation>
00306 static void findCriticalPoint_(Evaluation& Tcrit,
00307 Evaluation& pcrit,
00308 Evaluation& Vcrit,
00309 const Evaluation& a,
00310 const Evaluation& b)
00311 {
00312 Evaluation minVm(0);
00313 Evaluation maxVm(1e30);
00314
00315 Evaluation minP(0);
00316 Evaluation maxP(1e30);
00317
00318
00319
00320 Evaluation Tmin = 250;
00321 for (unsigned i = 0; i < 30; ++i) {
00322 bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
00323 if (hasExtrema)
00324 break;
00325 Tmin /= 2;
00326 };
00327
00328 Evaluation T = Tmin;
00329
00330
00331
00332 unsigned iMax = 100;
00333 for (unsigned i = 0; i < iMax; ++i) {
00334
00335 Evaluation f = maxVm - minVm;
00336
00337
00338 if (f < 1e-10 || (i == iMax - 1 && f < 1e-8)) {
00339 Tcrit = T;
00340 pcrit = (minP + maxP)/2;
00341 Vcrit = (maxVm + minVm)/2;
00342 return;
00343 }
00344
00345
00346
00347
00348
00349 const Scalar eps = - 1e-11;
00350 bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
00351 assert(hasExtrema);
00352 assert(std::isfinite(Opm::scalarValue(maxVm)));
00353 Evaluation fStar = maxVm - minVm;
00354
00355
00356
00357
00358 Evaluation fPrime = (fStar - f)/eps;
00359 if (std::abs(Opm::scalarValue(fPrime)) < 1e-40) {
00360 Tcrit = T;
00361 pcrit = (minP + maxP)/2;
00362 Vcrit = (maxVm + minVm)/2;
00363 return;
00364 }
00365
00366
00367 Evaluation delta = f/fPrime;
00368 assert(std::isfinite(Opm::scalarValue(delta)));
00369 if (delta > 0)
00370 delta = -10;
00371
00372
00373
00374 for (unsigned j = 0; ; ++j) {
00375 if (j >= 20) {
00376 if (f < 1e-8) {
00377 Tcrit = T;
00378 pcrit = (minP + maxP)/2;
00379 Vcrit = (maxVm + minVm)/2;
00380 return;
00381 }
00382 OPM_THROW(NumericalProblem,
00383 "Could not determine the critical point for a=" << a << ", b=" << b);
00384 }
00385
00386 if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
00387
00388
00389 T -= delta;
00390 break;
00391 }
00392 else
00393 delta /= 2;
00394 };
00395 }
00396
00397 assert(false);
00398 }
00399
00400
00401
00402 template <class Evaluation>
00403 static bool findExtrema_(Evaluation& Vmin,
00404 Evaluation& Vmax,
00405 Evaluation& ,
00406 Evaluation& ,
00407 const Evaluation& a,
00408 const Evaluation& b,
00409 const Evaluation& T)
00410 {
00411 Scalar u = 2;
00412 Scalar w = -1;
00413
00414 const Evaluation& RT = R*T;
00415
00416
00417
00418 const Evaluation& a1 = RT;
00419 const Evaluation& a2 = 2*RT*u*b - 2*a;
00420 const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
00421 const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
00422 const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
00423
00424 assert(std::isfinite(Opm::scalarValue(a1)));
00425 assert(std::isfinite(Opm::scalarValue(a2)));
00426 assert(std::isfinite(Opm::scalarValue(a3)));
00427 assert(std::isfinite(Opm::scalarValue(a4)));
00428 assert(std::isfinite(Opm::scalarValue(a5)));
00429
00430
00431
00432
00433
00434
00435 Evaluation V = b*1.1;
00436 Evaluation delta = 1.0;
00437 for (unsigned i = 0; std::abs(Opm::scalarValue(delta)) > 1e-12; ++i) {
00438 const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
00439 const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
00440
00441 if (std::abs(Opm::scalarValue(fPrime)) < 1e-20) {
00442
00443 return false;
00444 }
00445
00446
00447 delta = f/fPrime;
00448 V -= delta;
00449
00450 if (i > 200) {
00451
00452 return false;
00453 }
00454 }
00455 assert(std::isfinite(Opm::scalarValue(V)));
00456
00457
00458 Evaluation b1 = a1;
00459 Evaluation b2 = a2 + V*b1;
00460 Evaluation b3 = a3 + V*b2;
00461 Evaluation b4 = a4 + V*b3;
00462
00463
00464 Evaluation allV[4];
00465 allV[0] = V;
00466 int numSol = 1 + Opm::invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
00467
00468
00469 std::sort(allV + 0, allV + numSol);
00470
00471
00472 if (numSol != 4 || allV[numSol - 2] < b) {
00473
00474
00475 return false;
00476 }
00477
00478
00479
00480 Vmin = allV[numSol - 2];
00481 Vmax = allV[numSol - 1];
00482 return true;
00483 }
00484
00497 template <class Evaluation, class Params>
00498 static Evaluation ambroseWalton_(const Params& , const Evaluation& T)
00499 {
00500 typedef typename Params::Component Component;
00501
00502 const Evaluation& Tr = T / Component::criticalTemperature();
00503 const Evaluation& tau = 1 - Tr;
00504 const Evaluation& omega = Component::acentricFactor();
00505
00506 const Evaluation& f0 = (tau*(-5.97616 + Opm::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Opm::pow(tau, 5))/Tr;
00507 const Evaluation& f1 = (tau*(-5.03365 + Opm::sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*Opm::pow(tau, 5))/Tr;
00508 const Evaluation& f2 = (tau*(-0.64771 + Opm::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Opm::pow(tau, 5))/Tr;
00509
00510 return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
00511 }
00512
00523 template <class Evaluation, class Params>
00524 static Evaluation fugacityDifference_(const Params& params,
00525 const Evaluation& T,
00526 const Evaluation& p,
00527 const Evaluation& VmLiquid,
00528 const Evaluation& VmGas)
00529 { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
00530
00531
00532
00533
00534
00535
00536 };
00537
00538 template <class Scalar>
00539 const Scalar PengRobinson<Scalar>::R = Opm::Constants<Scalar>::R;
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552 }
00553
00554 #endif