| 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 a | ||
| 11 | * mixture. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_PENG_ROBINSON_MIXTURE_HH | ||
| 14 | #define DUMUX_PENG_ROBINSON_MIXTURE_HH | ||
| 15 | |||
| 16 | #include "pengrobinson.hh" | ||
| 17 | |||
| 18 | #include <dumux/material/constants.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup EOS | ||
| 24 | * \brief Implements the Peng-Robinson equation of state for a | ||
| 25 | * mixture. | ||
| 26 | */ | ||
| 27 | template <class Scalar, class StaticParameters> | ||
| 28 | class PengRobinsonMixture | ||
| 29 | { | ||
| 30 | enum { numComponents = StaticParameters::numComponents }; | ||
| 31 | using PengRobinson = Dumux::PengRobinson<Scalar>; | ||
| 32 | |||
| 33 | // this class cannot be instantiated! | ||
| 34 | PengRobinsonMixture() {}; | ||
| 35 | |||
| 36 | // the u and w parameters as given by the Peng-Robinson EOS | ||
| 37 | static const Scalar u; | ||
| 38 | static const Scalar w; | ||
| 39 | |||
| 40 | public: | ||
| 41 | |||
| 42 | /*! | ||
| 43 | * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual | ||
| 44 | * component in the phase. | ||
| 45 | * \param fs Thermodynamic state of the fluids | ||
| 46 | * \param params Parameters | ||
| 47 | * \param phaseIdx The phase index | ||
| 48 | * \param compIdx The index of the component | ||
| 49 | * | ||
| 50 | * The fugacity coefficient \f$\phi_i\f$ of a component \f$i\f$ is | ||
| 51 | * defined as | ||
| 52 | * \f[f_i = \phi_i x_i \;,\f] | ||
| 53 | * where \f$f_i\f$ is the component's fugacity and \f$x_i\f$ is | ||
| 54 | * the component's mole fraction. | ||
| 55 | * | ||
| 56 | * See: | ||
| 57 | * | ||
| 58 | * R. Reid, et al. (1987, pp. 42-44, 143-145) \cite reid1987 | ||
| 59 | */ | ||
| 60 | template <class FluidState, class Params> | ||
| 61 | 265216 | static Scalar computeFugacityCoefficient(const FluidState &fs, | |
| 62 | const Params ¶ms, | ||
| 63 | int phaseIdx, | ||
| 64 | int compIdx) | ||
| 65 | { | ||
| 66 | // note that we normalize the component mole fractions, so | ||
| 67 | // that their sum is 100%. This increases numerical stability | ||
| 68 | // considerably if the fluid state is not physical. | ||
| 69 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 265216 times.
|
265216 | Scalar Vm = params.molarVolume(phaseIdx); |
| 70 | |||
| 71 | // Calculate b_i / b | ||
| 72 | 265216 | Scalar bi_b = params.bPure(phaseIdx, compIdx) / params.b(phaseIdx); | |
| 73 | |||
| 74 | // Calculate the compressibility factor | ||
| 75 | 265216 | Scalar RT = Constants<Scalar>::R*fs.temperature(phaseIdx); | |
| 76 | 265216 | Scalar p = fs.pressure(phaseIdx); // molar volume in [bar] | |
| 77 | 265216 | Scalar Z = p*Vm/RT; // compressibility factor | |
| 78 | |||
| 79 | // Calculate A^* and B^* (see: Reid, p. 42) | ||
| 80 | 265216 | Scalar Astar = params.a(phaseIdx)*p/(RT*RT); | |
| 81 | 265216 | Scalar Bstar = params.b(phaseIdx)*p/(RT); | |
| 82 | |||
| 83 | // calculate delta_i (see: Reid, p. 145) | ||
| 84 | 265216 | Scalar sumMoleFractions = 0.0; | |
| 85 |
2/2✓ Branch 0 taken 1856512 times.
✓ Branch 1 taken 265216 times.
|
2121728 | for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) |
| 86 | 3712926 | sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx); | |
| 87 | |||
| 88 | using std::sqrt; | ||
| 89 | 265216 | Scalar deltai = 2*sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx); | |
| 90 | 265216 | Scalar tmp = 0; | |
| 91 |
2/2✓ Branch 0 taken 1856512 times.
✓ Branch 1 taken 265216 times.
|
2121728 | for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) { |
| 92 | 1856512 | tmp += | |
| 93 | fs.moleFraction(phaseIdx, compJIdx) | ||
| 94 | 1856512 | / sumMoleFractions | |
| 95 | 1856512 | * sqrt(params.aPure(phaseIdx, compJIdx)) | |
| 96 |
2/2✓ Branch 0 taken 795648 times.
✓ Branch 1 taken 1060864 times.
|
3713024 | * (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx)); |
| 97 | } | ||
| 98 | 265216 | deltai *= tmp; | |
| 99 | |||
| 100 | 265216 | Scalar base = | |
| 101 | 265216 | (2*Z + Bstar*(u + sqrt(u*u - 4*w))) / | |
| 102 | 265216 | (2*Z + Bstar*(u - sqrt(u*u - 4*w))); | |
| 103 | 265216 | Scalar expo = Astar/(Bstar*sqrt(u*u - 4*w))*(bi_b - deltai); | |
| 104 | |||
| 105 | using std::exp; | ||
| 106 | using std::max; | ||
| 107 | using std::min; | ||
| 108 | using std::pow; | ||
| 109 | 265216 | Scalar fugCoeff = | |
| 110 |
1/2✓ Branch 0 taken 265216 times.
✗ Branch 1 not taken.
|
265216 | exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) * |
| 111 | 265216 | pow(base, expo); | |
| 112 | |||
| 113 | //////// | ||
| 114 | // limit the fugacity coefficient to a reasonable range: | ||
| 115 | // | ||
| 116 | // on one side, we want the mole fraction to be at | ||
| 117 | // least 10^-3 if the fugacity is at the current pressure | ||
| 118 | // | ||
| 119 | |||
| 120 |
2/2✓ Branch 0 taken 265040 times.
✓ Branch 1 taken 176 times.
|
265216 | fugCoeff = min(1e10, fugCoeff); |
| 121 | // | ||
| 122 | // on the other hand, if the mole fraction of the component is 100%, we want the | ||
| 123 | // fugacity to be at least 10^-3 Pa | ||
| 124 | // | ||
| 125 |
1/2✓ Branch 0 taken 265216 times.
✗ Branch 1 not taken.
|
265216 | fugCoeff = max(1e-10, fugCoeff); |
| 126 | /////////// | ||
| 127 | using std::isfinite; | ||
| 128 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 265216 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 265216 times.
|
530432 | if (!isfinite(fugCoeff)) { |
| 129 | ✗ | std::cout << "Non finite phi: " << fugCoeff << "\n"; | |
| 130 | } | ||
| 131 | |||
| 132 | 265216 | return fugCoeff; | |
| 133 | } | ||
| 134 | |||
| 135 | }; | ||
| 136 | |||
| 137 | template<class Scalar, class StaticParameters> | ||
| 138 | const Scalar PengRobinsonMixture<Scalar, StaticParameters>::u = 2.0; | ||
| 139 | template<class Scalar, class StaticParameters> | ||
| 140 | const Scalar PengRobinsonMixture<Scalar, StaticParameters>::w = -1.0; | ||
| 141 | |||
| 142 | } // end namespace Dumux | ||
| 143 | |||
| 144 | #endif | ||
| 145 |