| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup EOS | ||
| 10 | * \brief Implements the Peng-Robinson equation of state for liquids | ||
| 11 | * and gases. | ||
| 12 | * | ||
| 13 | * See: | ||
| 14 | * | ||
| 15 | * D.-Y. Peng, D.B. Robinson (1976, pp. 59–64) \cite peng1976 <BR> | ||
| 16 | * https://doi.org/10.1021/i160057a011 | ||
| 17 | * | ||
| 18 | * R. Reid, et al. "The properties of gases and liquids" (1987, pp. 42-44, 82) | ||
| 19 | * \cite reid1987 | ||
| 20 | */ | ||
| 21 | #ifndef DUMUX_PENG_ROBINSON_HH | ||
| 22 | #define DUMUX_PENG_ROBINSON_HH | ||
| 23 | |||
| 24 | #include <dune/common/math.hh> | ||
| 25 | #include <dumux/material/fluidstates/temperatureoverlay.hh> | ||
| 26 | #include <dumux/material/idealgas.hh> | ||
| 27 | #include <dumux/common/exceptions.hh> | ||
| 28 | #include <dumux/common/math.hh> | ||
| 29 | #include <dumux/common/tabulated2dfunction.hh> | ||
| 30 | |||
| 31 | #include <iostream> | ||
| 32 | |||
| 33 | namespace Dumux { | ||
| 34 | |||
| 35 | /*! | ||
| 36 | * \ingroup EOS | ||
| 37 | * \brief Implements the Peng-Robinson equation of state for liquids | ||
| 38 | * and gases. | ||
| 39 | * | ||
| 40 | * See: | ||
| 41 | * | ||
| 42 | * D.-Y. Peng, D.B. Robinson (1976, pp. 59–64) \cite peng1976 <BR> | ||
| 43 | * | ||
| 44 | * R. Reid, et al. (1987, pp. 42-44, 82) \cite reid1987 | ||
| 45 | */ | ||
| 46 | template <class Scalar> | ||
| 47 | class PengRobinson | ||
| 48 | { | ||
| 49 | PengRobinson() | ||
| 50 | { } | ||
| 51 | |||
| 52 | public: | ||
| 53 | 2 | static void init(Scalar aMin, Scalar aMax, int na, | |
| 54 | Scalar bMin, Scalar bMax, int nb) | ||
| 55 | { | ||
| 56 | // resize the tabulation for the critical points | ||
| 57 | 2 | criticalTemperature_.resize(aMin, aMax, na, bMin, bMax, nb); | |
| 58 | 2 | criticalPressure_.resize(aMin, aMax, na, bMin, bMax, nb); | |
| 59 | 2 | criticalMolarVolume_.resize(aMin, aMax, na, bMin, bMax, nb); | |
| 60 | |||
| 61 | 2 | Scalar VmCrit = 18e-6; | |
| 62 | 2 | Scalar pCrit = 1e7; | |
| 63 | 2 | Scalar TCrit = 273; | |
| 64 |
2/2✓ Branch 0 taken 200 times.
✓ Branch 1 taken 2 times.
|
202 | for (int i = 0; i < na; ++i) { |
| 65 | 200 | Scalar a = criticalTemperature_.iToX(i); | |
| 66 |
2/2✓ Branch 1 taken 40000 times.
✓ Branch 2 taken 200 times.
|
40200 | for (int j = 0; j < nb; ++j) { |
| 67 | 40000 | Scalar b = criticalTemperature_.jToY(j); | |
| 68 | |||
| 69 | 40000 | findCriticalPoint_(TCrit, pCrit, VmCrit, a, b); | |
| 70 | |||
| 71 | 40000 | criticalTemperature_.setSamplePoint(i, j, TCrit); | |
| 72 | 40000 | criticalPressure_.setSamplePoint(i, j, pCrit); | |
| 73 | 40000 | criticalMolarVolume_.setSamplePoint(i, j, VmCrit); | |
| 74 | } | ||
| 75 | } | ||
| 76 | 2 | }; | |
| 77 | |||
| 78 | /*! | ||
| 79 | * \brief Predicts the vapor pressure \f$\mathrm{[Pa]}\f$ for the temperature given in | ||
| 80 | * setTP(). | ||
| 81 | * \param T temperature in \f$\mathrm{[K]}\f$ | ||
| 82 | * \param params Parameters | ||
| 83 | * | ||
| 84 | * Initially, the vapor pressure is roughly estimated by using the | ||
| 85 | * Ambrose-Walton method, then the Newton method is used to make | ||
| 86 | * difference between the gas and liquid phase fugacity zero. | ||
| 87 | */ | ||
| 88 | template <class Params> | ||
| 89 | static Scalar computeVaporPressure(const Params ¶ms, Scalar T) | ||
| 90 | { | ||
| 91 | using Component = typename Params::Component; | ||
| 92 | if (T >= Component::criticalTemperature()) | ||
| 93 | return Component::criticalPressure(); | ||
| 94 | |||
| 95 | // initial guess of the vapor pressure | ||
| 96 | Scalar Vm[3]; | ||
| 97 | const Scalar eps = Component::criticalPressure()*1e-10; | ||
| 98 | |||
| 99 | // use the Ambrose-Walton method to get an initial guess of | ||
| 100 | // the vapor pressure | ||
| 101 | Scalar pVap = ambroseWalton_(params, T); | ||
| 102 | |||
| 103 | // Newton-Raphson method | ||
| 104 | for (int i = 0; i < 5; ++i) { | ||
| 105 | // calculate the molar densities | ||
| 106 | assert(molarVolumes(Vm, params, T, pVap) == 3); | ||
| 107 | |||
| 108 | Scalar f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]); | ||
| 109 | Scalar df_dp = | ||
| 110 | fugacityDifference_(params, T, pVap + eps, Vm[0], Vm[2]) | ||
| 111 | - | ||
| 112 | fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]); | ||
| 113 | df_dp /= 2*eps; | ||
| 114 | |||
| 115 | Scalar delta = f/df_dp; | ||
| 116 | pVap = pVap - delta; | ||
| 117 | using std::abs; | ||
| 118 | if (abs(delta/pVap) < 1e-10) | ||
| 119 | break; | ||
| 120 | } | ||
| 121 | |||
| 122 | return pVap; | ||
| 123 | } | ||
| 124 | |||
| 125 | /*! | ||
| 126 | * \brief Computes molar volumes \f$\mathrm{[m^3 / mol]}\f$ where the Peng-Robinson EOS is | ||
| 127 | * true. | ||
| 128 | * \param fs Thermodynamic state of the fluids | ||
| 129 | * \param params Parameters | ||
| 130 | * \param phaseIdx The phase index | ||
| 131 | * \param isGasPhase Specifies the phase state | ||
| 132 | * \param handleUnphysicalPhase Special handling of the case when the EOS has only one | ||
| 133 | intersection with the pressure, and the intersection does not correspond to | ||
| 134 | the given phase (the phase is thus considered unphysical). If it happens in | ||
| 135 | the case of critical fluid, the critical molar volume is returned for the | ||
| 136 | unphysical phase. If the fluid is not critical, a proper extremum of the | ||
| 137 | EOS is returned for the unphysical phase. If the parameter is false and the | ||
| 138 | EOS has only one intersection with the pressure, the molar volume is computed | ||
| 139 | from that single intersection, not depending of the given phase (gas or | ||
| 140 | fluid). If the EOS has three intersections with the pressure, this parameter | ||
| 141 | is ignored. | ||
| 142 | */ | ||
| 143 | template <class FluidState, class Params> | ||
| 144 | 38009 | static Scalar computeMolarVolume(const FluidState &fs, | |
| 145 | Params ¶ms, | ||
| 146 | int phaseIdx, | ||
| 147 | bool isGasPhase, | ||
| 148 | bool handleUnphysicalPhase = true) | ||
| 149 | { | ||
| 150 | 38009 | Scalar Vm = 0; | |
| 151 | |||
| 152 | 38009 | Scalar T = fs.temperature(phaseIdx); | |
| 153 | 38009 | Scalar p = fs.pressure(phaseIdx); | |
| 154 | |||
| 155 | 38009 | Scalar a = params.a(phaseIdx); // "attractive factor" | |
| 156 | 38009 | Scalar b = params.b(phaseIdx); // "co-volume" | |
| 157 | |||
| 158 | 38009 | Scalar RT = Constants<Scalar>::R*T; | |
| 159 | 38009 | Scalar Astar = a*p/(RT*RT); | |
| 160 | 38009 | Scalar Bstar = b*p/RT; | |
| 161 | |||
| 162 | 38009 | Scalar a1 = 1.0; | |
| 163 | 38009 | Scalar a2 = - (1 - Bstar); | |
| 164 | 38009 | Scalar a3 = Astar - Bstar*(3*Bstar + 2); | |
| 165 | 38009 | Scalar a4 = Bstar*(- Astar + Bstar*(1 + Bstar)); | |
| 166 | |||
| 167 | 38009 | Scalar Z[3] = {}; // zero-initializer | |
| 168 | 38009 | int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4); | |
| 169 | |||
| 170 | // ignore the first two results if the smallest | ||
| 171 | // compressibility factor is <= 0.0. (this means that if we | ||
| 172 | // would get negative molar volumes for the liquid phase, we | ||
| 173 | // consider the liquid phase non-existent.) | ||
| 174 | // Note that invertCubicPolynomial sorts the roots in ascending order | ||
| 175 |
3/4✓ Branch 0 taken 12270 times.
✓ Branch 1 taken 25739 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12270 times.
|
38009 | if (numSol > 1 && Z[0] <= 0.0) { |
| 176 | ✗ | Z[0] = Z[numSol - 1]; | |
| 177 | ✗ | numSol = 1; | |
| 178 | } | ||
| 179 | |||
| 180 |
2/4✓ Branch 0 taken 38009 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 38009 times.
|
38009 | if (numSol > 0 && Z[0] <= 0) |
| 181 | ✗ | DUNE_THROW(NumericalProblem, "No positive compressibility factors found"); | |
| 182 | |||
| 183 |
2/2✓ Branch 0 taken 12270 times.
✓ Branch 1 taken 25739 times.
|
38009 | if (numSol == 3) { |
| 184 | // the EOS has three intersections with the pressure, | ||
| 185 | // i.e. the molar volume of gas is the largest one and the | ||
| 186 | // molar volume of liquid is the smallest one | ||
| 187 |
2/2✓ Branch 0 taken 77 times.
✓ Branch 1 taken 12193 times.
|
12270 | if (isGasPhase) |
| 188 | 77 | Vm = Z[2]*RT/p; | |
| 189 | else | ||
| 190 | 12193 | Vm = Z[0]*RT/p; | |
| 191 | } | ||
| 192 |
1/2✓ Branch 0 taken 25739 times.
✗ Branch 1 not taken.
|
25739 | else if (numSol == 1) { |
| 193 | // the EOS has only one intersection with the pressure | ||
| 194 | |||
| 195 | 25739 | Vm = Z[0]*RT/p; | |
| 196 | |||
| 197 | // Handle single root case specially unless told otherwise | ||
| 198 |
1/2✓ Branch 0 taken 25739 times.
✗ Branch 1 not taken.
|
25739 | if (handleUnphysicalPhase) |
| 199 | 25739 | Vm = handleSingleRoot_(Vm, fs, params, phaseIdx, isGasPhase); | |
| 200 | } | ||
| 201 | else | ||
| 202 | ✗ | DUNE_THROW(NumericalProblem, "Unexpected number of roots in the equation for compressibility factor: " << numSol); | |
| 203 | |||
| 204 | using std::isfinite; | ||
| 205 | |||
| 206 |
3/6✓ Branch 0 taken 38009 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 38009 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 38009 times.
✗ Branch 5 not taken.
|
76018 | if (!isfinite(Vm) || Vm <= 0) |
| 207 | ✗ | DUNE_THROW(NumericalProblem, "Bad molar volume"); | |
| 208 | |||
| 209 | 38009 | return Vm; | |
| 210 | } | ||
| 211 | |||
| 212 | /*! | ||
| 213 | * \brief Returns the critical temperature for a given mix | ||
| 214 | * | ||
| 215 | * \param params EOS (equation of state) parameters of a single-component fluid \ | ||
| 216 | * (usually PengRobinsonParms) or a mixture (usually PengRobinsonMixtureParams) | ||
| 217 | */ | ||
| 218 | template <class Params> | ||
| 219 | static Scalar criticalTemperature(const Params ¶ms) | ||
| 220 | { | ||
| 221 | // Here, a() is the attractive force parameter and b() | ||
| 222 | // is the repulsive force parameter of the EOS | ||
| 223 | return criticalTemperature_(params.a(), params.b()); | ||
| 224 | } | ||
| 225 | |||
| 226 | /*! | ||
| 227 | * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ for a given pressure | ||
| 228 | * and molar volume. | ||
| 229 | * | ||
| 230 | * This is the same value as computeFugacity() because the mole | ||
| 231 | * fraction of a component in a pure fluid is obviously always | ||
| 232 | * 100%. | ||
| 233 | * | ||
| 234 | * \param params Parameters | ||
| 235 | */ | ||
| 236 | template <class Params> | ||
| 237 | static Scalar computeFugacityCoeffient(const Params ¶ms) | ||
| 238 | { | ||
| 239 | Scalar T = params.temperature(); | ||
| 240 | Scalar p = params.pressure(); | ||
| 241 | Scalar Vm = params.molarVolume(); | ||
| 242 | |||
| 243 | Scalar RT = Constants<Scalar>::R*T; | ||
| 244 | Scalar Z = p*Vm/RT; | ||
| 245 | Scalar Bstar = p*params.b() / RT; | ||
| 246 | using std::sqrt; | ||
| 247 | Scalar tmp = | ||
| 248 | (Vm + params.b()*(1 + sqrt(2))) / | ||
| 249 | (Vm + params.b()*(1 - sqrt(2))); | ||
| 250 | Scalar expo = - params.a()/(RT * 2 * params.b() * sqrt(2)); | ||
| 251 | using std::exp; | ||
| 252 | using std::pow; | ||
| 253 | Scalar fugCoeff = | ||
| 254 | exp(Z - 1) / (Z - Bstar) * | ||
| 255 | pow(tmp, expo); | ||
| 256 | |||
| 257 | return fugCoeff; | ||
| 258 | } | ||
| 259 | |||
| 260 | /*! | ||
| 261 | * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ for a given pressure | ||
| 262 | * and molar volume. | ||
| 263 | * | ||
| 264 | * This is the fugacity coefficient times the pressure. The mole | ||
| 265 | * fraction of a component in a pure fluid is obviously always | ||
| 266 | * 100%, so it is not required. | ||
| 267 | * | ||
| 268 | * \param params Parameters | ||
| 269 | */ | ||
| 270 | template <class Params> | ||
| 271 | static Scalar computeFugacity(const Params ¶ms) | ||
| 272 | { return params.pressure()*computeFugacityCoeff(params); } | ||
| 273 | |||
| 274 | protected: | ||
| 275 | |||
| 276 | /* | ||
| 277 | * Handle single-root case in the equation of state. The problem here is | ||
| 278 | * that only one phase exists in this case, while we should also figure out | ||
| 279 | * something to return for the other phase. For reference, see Andreas | ||
| 280 | * Lauser's thesis: http://dx.doi.org/10.18419/opus-516 (page 32). | ||
| 281 | */ | ||
| 282 | template <class FluidState, class Params> | ||
| 283 | 25739 | static Scalar handleSingleRoot_(Scalar Vm, | |
| 284 | const FluidState &fs, | ||
| 285 | const Params ¶ms, | ||
| 286 | int phaseIdx, | ||
| 287 | bool isGasPhase) | ||
| 288 | { | ||
| 289 | 25739 | Scalar T = fs.temperature(phaseIdx); | |
| 290 | |||
| 291 | // for the other phase, we take the extremum of the EOS | ||
| 292 | // with the largest distance from the intersection. | ||
| 293 |
2/2✓ Branch 3 taken 18553 times.
✓ Branch 4 taken 7186 times.
|
25739 | if (T > criticalTemperature_(params.a(phaseIdx), params.b(phaseIdx))) { |
| 294 | // if the EOS does not exhibit any extrema, the fluid | ||
| 295 | // is critical... | ||
| 296 | 18553 | return handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase); | |
| 297 | } | ||
| 298 | else { | ||
| 299 | // find the extrema (if they are present) | ||
| 300 | Scalar Vmin, Vmax, pmin, pmax; | ||
| 301 | using std::min; | ||
| 302 | using std::max; | ||
| 303 |
1/2✓ Branch 3 taken 7186 times.
✗ Branch 4 not taken.
|
7186 | if (findExtrema_(Vmin, Vmax, pmin, pmax, |
| 304 | params.a(phaseIdx), params.b(phaseIdx), T)) | ||
| 305 |
5/6✓ Branch 0 taken 396 times.
✓ Branch 1 taken 6790 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 385 times.
✓ Branch 4 taken 6790 times.
✗ Branch 5 not taken.
|
13987 | return isGasPhase ? max(Vmax, Vm) : min(Vmin, Vm); |
| 306 | } | ||
| 307 | ✗ | return Vm; | |
| 308 | } | ||
| 309 | |||
| 310 | template <class FluidState, class Params> | ||
| 311 | 18553 | static Scalar handleCriticalFluid_(Scalar Vm, | |
| 312 | const FluidState &fs, | ||
| 313 | const Params ¶ms, | ||
| 314 | int phaseIdx, | ||
| 315 | bool isGasPhase) | ||
| 316 | { | ||
| 317 | 18553 | Scalar Vcrit = criticalMolarVolume_(params.a(phaseIdx), params.b(phaseIdx)); | |
| 318 | |||
| 319 | using std::max; | ||
| 320 | using std::min; | ||
| 321 |
5/6✓ Branch 0 taken 18531 times.
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 334 times.
✓ Branch 3 taken 18197 times.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
|
18909 | return isGasPhase ? max(Vm, Vcrit) : min(Vm, Vcrit); |
| 322 | } | ||
| 323 | |||
| 324 | 40000 | static void findCriticalPoint_(Scalar &Tcrit, | |
| 325 | Scalar &pcrit, | ||
| 326 | Scalar &Vcrit, | ||
| 327 | Scalar a, | ||
| 328 | Scalar b) | ||
| 329 | { | ||
| 330 | Scalar minVm; | ||
| 331 | Scalar maxVm; | ||
| 332 | |||
| 333 | 40000 | Scalar minP(0); | |
| 334 | 40000 | Scalar maxP(1e100); | |
| 335 | |||
| 336 | // first, we need to find an isotherm where the EOS exhibits | ||
| 337 | // a maximum and a minimum | ||
| 338 | 40000 | Scalar Tmin = 250; // [K] | |
| 339 |
1/2✓ Branch 0 taken 45194 times.
✗ Branch 1 not taken.
|
45194 | for (int i = 0; i < 30; ++i) { |
| 340 | 45194 | bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin); | |
| 341 |
2/2✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 5194 times.
|
45194 | if (hasExtrema) |
| 342 | break; | ||
| 343 | 5194 | Tmin /= 2; | |
| 344 | } | ||
| 345 | |||
| 346 | Scalar T = Tmin; | ||
| 347 | |||
| 348 | // Newton's method: Start at minimum temperature and minimize | ||
| 349 | // the "gap" between the extrema of the EOS | ||
| 350 |
2/2✓ Branch 0 taken 811658 times.
✓ Branch 1 taken 7920 times.
|
819578 | for (int i = 0; i < 25; ++i) { |
| 351 | // calculate function we would like to minimize | ||
| 352 | 811658 | Scalar f = maxVm - minVm; | |
| 353 | |||
| 354 | // check if we're converged | ||
| 355 |
2/2✓ Branch 0 taken 32080 times.
✓ Branch 1 taken 779578 times.
|
811658 | if (f < 1e-10) { |
| 356 | 32080 | Tcrit = T; | |
| 357 | 32080 | pcrit = (minP + maxP)/2; | |
| 358 | 32080 | Vcrit = (maxVm + minVm)/2; | |
| 359 | 32080 | return; | |
| 360 | } | ||
| 361 | |||
| 362 | // backward differences. Forward differences are not | ||
| 363 | // robust, since the isotherm could be critical if an | ||
| 364 | // epsilon was added to the temperature. (this is case | ||
| 365 | // rarely happens, though) | ||
| 366 | 779578 | const Scalar eps = - 1e-8; | |
| 367 | 779578 | [[maybe_unused]] bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps); | |
| 368 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 779578 times.
|
779578 | assert(hasExtrema); |
| 369 | 779578 | Scalar fStar = maxVm - minVm; | |
| 370 | |||
| 371 | // derivative of the difference between the maximum's | ||
| 372 | // molar volume and the minimum's molar volume regarding | ||
| 373 | // temperature | ||
| 374 | 779578 | Scalar fPrime = (fStar - f)/eps; | |
| 375 | |||
| 376 | // update value for the current iteration | ||
| 377 | 779578 | Scalar delta = f/fPrime; | |
| 378 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 779578 times.
|
779578 | if (delta > 0) |
| 379 | ✗ | delta = -10; | |
| 380 | |||
| 381 | // line search (we have to make sure that both extrema | ||
| 382 | // still exist after the update) | ||
| 383 | 1496318 | for (int j = 0; ; ++j) { | |
| 384 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2275896 times.
|
2275896 | if (j >= 20) { |
| 385 | ✗ | DUNE_THROW(NumericalProblem, | |
| 386 | "Could not determine the critical point for a=" << a << ", b=" << b); | ||
| 387 | } | ||
| 388 | |||
| 389 |
2/2✓ Branch 1 taken 779578 times.
✓ Branch 2 taken 1496318 times.
|
2275896 | if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) { |
| 390 | // if the isotherm for T - delta exhibits two | ||
| 391 | // extrema the update is finished | ||
| 392 | 779578 | T -= delta; | |
| 393 | break; | ||
| 394 | } | ||
| 395 | else | ||
| 396 | 1496318 | delta /= 2; | |
| 397 | } | ||
| 398 | } | ||
| 399 | } | ||
| 400 | |||
| 401 | // find the two molar volumes where the EOS exhibits extrema and | ||
| 402 | // which are larger than the covolume of the phase | ||
| 403 | 3107854 | static bool findExtrema_(Scalar &Vmin, | |
| 404 | Scalar &Vmax, | ||
| 405 | Scalar &pMin, | ||
| 406 | Scalar &pMax, | ||
| 407 | Scalar a, | ||
| 408 | Scalar b, | ||
| 409 | Scalar T) | ||
| 410 | { | ||
| 411 | 3107854 | Scalar u = 2; | |
| 412 | 3107854 | Scalar w = -1; | |
| 413 | |||
| 414 | 3107854 | Scalar RT = Constants<Scalar>::R*T; | |
| 415 | |||
| 416 | // calculate coefficients of the 4th order polynominal in | ||
| 417 | // monomial basis | ||
| 418 | 3107854 | Scalar a1 = RT; | |
| 419 | 3107854 | Scalar a2 = 2*RT*u*b - 2*a; | |
| 420 | 3107854 | Scalar a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b; | |
| 421 | 3107854 | Scalar a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b; | |
| 422 | 3107854 | Scalar a5 = RT*w*w*b*b*b*b - u*a*b*b*b; | |
| 423 | |||
| 424 | // Newton method to find first root | ||
| 425 | |||
| 426 | // if the values which we got on Vmin and Vmax are useful, we | ||
| 427 | // will reuse them as initial value, else we will start 10% | ||
| 428 | // above the covolume | ||
| 429 | 3107854 | Scalar V = b*1.1; | |
| 430 | 3107854 | Scalar delta = 1.0; | |
| 431 | using std::abs; | ||
| 432 |
4/4✓ Branch 0 taken 12360375 times.
✓ Branch 1 taken 3107796 times.
✓ Branch 2 taken 12360375 times.
✓ Branch 3 taken 3107796 times.
|
15468171 | for (int i = 0; abs(delta) > 1e-9; ++i) { |
| 433 | 12360375 | Scalar f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1))); | |
| 434 | 12360375 | Scalar fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1)); | |
| 435 | |||
| 436 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 12360375 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12360375 times.
|
24720750 | if (abs(fPrime) < 1e-20) { |
| 437 | // give up if the derivative is zero | ||
| 438 | ✗ | Vmin = 0; | |
| 439 | ✗ | Vmax = 0; | |
| 440 | ✗ | return false; | |
| 441 | } | ||
| 442 | |||
| 443 | |||
| 444 | 12360375 | delta = f/fPrime; | |
| 445 | 12360375 | V -= delta; | |
| 446 | |||
| 447 |
2/2✓ Branch 0 taken 58 times.
✓ Branch 1 taken 12360317 times.
|
12360375 | if (i > 20) { |
| 448 | // give up after 20 iterations... | ||
| 449 | 58 | Vmin = 0; | |
| 450 | 58 | Vmax = 0; | |
| 451 | 58 | return false; | |
| 452 | } | ||
| 453 | } | ||
| 454 | |||
| 455 | // polynomial division | ||
| 456 | 3107796 | Scalar b1 = a1; | |
| 457 | 3107796 | Scalar b2 = a2 + V*b1; | |
| 458 | 3107796 | Scalar b3 = a3 + V*b2; | |
| 459 | 3107796 | Scalar b4 = a4 + V*b3; | |
| 460 | |||
| 461 | // invert resulting cubic polynomial analytically | ||
| 462 | Scalar allV[4]; | ||
| 463 | 3107796 | allV[0] = V; | |
| 464 | 3107796 | int numSol = 1 + invertCubicPolynomial(allV + 1, b1, b2, b3, b4); | |
| 465 | |||
| 466 | // sort all roots of the derivative | ||
| 467 | 3107796 | std::sort(allV + 0, allV + numSol); | |
| 468 | |||
| 469 | // check whether the result is physical | ||
| 470 |
2/2✓ Branch 0 taken 1501454 times.
✓ Branch 1 taken 1606342 times.
|
3107796 | if (allV[numSol - 2] < b) { |
| 471 | // the second largest extremum is smaller than the phase's | ||
| 472 | // covolume which is physically impossible | ||
| 473 | 1501454 | Vmin = Vmax = 0; | |
| 474 | 1501454 | return false; | |
| 475 | } | ||
| 476 | |||
| 477 | // it seems that everything is okay... | ||
| 478 | 1606342 | Vmin = allV[numSol - 2]; | |
| 479 | 1606342 | Vmax = allV[numSol - 1]; | |
| 480 | 1606342 | return true; | |
| 481 | } | ||
| 482 | |||
| 483 | /*! | ||
| 484 | * \brief The Ambrose-Walton method to estimate the vapor | ||
| 485 | * pressure. | ||
| 486 | * | ||
| 487 | * \return Vapor pressure estimate in bar | ||
| 488 | * | ||
| 489 | * See: | ||
| 490 | * | ||
| 491 | * D. Ambrose, J. Walton (1989, pp. 1395-1403) \cite ambrose1989 | ||
| 492 | */ | ||
| 493 | template <class Params> | ||
| 494 | static Scalar ambroseWalton_(const Params ¶ms, Scalar T) | ||
| 495 | { | ||
| 496 | using Component = typename Params::Component; | ||
| 497 | |||
| 498 | Scalar Tr = T / Component::criticalTemperature(); | ||
| 499 | Scalar tau = 1 - Tr; | ||
| 500 | Scalar omega = Component::acentricFactor(); | ||
| 501 | using std::sqrt; | ||
| 502 | using Dune::power; | ||
| 503 | Scalar f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*power(tau, 5))/Tr; | ||
| 504 | Scalar f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*power(tau, 5))/Tr; | ||
| 505 | Scalar f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*power(tau, 5))/Tr; | ||
| 506 | using std::exp; | ||
| 507 | return Component::criticalPressure()*exp(f0 + omega * (f1 + omega*f2)); | ||
| 508 | } | ||
| 509 | |||
| 510 | /*! | ||
| 511 | * \brief Returns the difference between the liquid and the gas phase | ||
| 512 | * fugacities in [bar] | ||
| 513 | * | ||
| 514 | * \param params Parameters | ||
| 515 | * \param T Temperature [K] | ||
| 516 | * \param p Pressure [bar] | ||
| 517 | * \param VmLiquid Molar volume of the liquid phase [cm^3/mol] | ||
| 518 | * \param VmGas Molar volume of the gas phase [cm^3/mol] | ||
| 519 | */ | ||
| 520 | template <class Params> | ||
| 521 | static Scalar fugacityDifference_(const Params ¶ms, | ||
| 522 | Scalar T, | ||
| 523 | Scalar p, | ||
| 524 | Scalar VmLiquid, | ||
| 525 | Scalar VmGas) | ||
| 526 | { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); } | ||
| 527 | |||
| 528 | static Tabulated2DFunction<Scalar> criticalTemperature_; | ||
| 529 | static Tabulated2DFunction<Scalar> criticalPressure_; | ||
| 530 | static Tabulated2DFunction<Scalar> criticalMolarVolume_; | ||
| 531 | }; | ||
| 532 | |||
| 533 | template <class Scalar> | ||
| 534 | Tabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalTemperature_; | ||
| 535 | |||
| 536 | template <class Scalar> | ||
| 537 | Tabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalPressure_; | ||
| 538 | |||
| 539 | template <class Scalar> | ||
| 540 | Tabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalMolarVolume_; | ||
| 541 | |||
| 542 | } // end namespace | ||
| 543 | |||
| 544 | #endif | ||
| 545 |