3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
pengrobinsonmixture.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_PENG_ROBINSON_MIXTURE_HH
26#define DUMUX_PENG_ROBINSON_MIXTURE_HH
27
28#include "pengrobinson.hh"
29
31
32namespace Dumux {
33
39template <class Scalar, class StaticParameters>
41{
42 enum { numComponents = StaticParameters::numComponents };
44
45 // this class cannot be instantiated!
47
48 // the u and w parameters as given by the Peng-Robinson EOS
49 static const Scalar u;
50 static const Scalar w;
51
52public:
62 template <class MutableParams, class FluidState>
63 static int computeMolarVolumes(Scalar *Vm,
64 const MutableParams &params,
65 int phaseIdx,
66 const FluidState &fs)
67 {
68 return PengRobinson::computeMolarVolumes(Vm, params, phaseIdx, fs);
69 }
70
89 template <class FluidState, class Params>
90 static Scalar computeFugacityCoefficient(const FluidState &fs,
91 const Params &params,
92 int phaseIdx,
93 int compIdx)
94 {
95 // note that we normalize the component mole fractions, so
96 // that their sum is 100%. This increases numerical stability
97 // considerably if the fluid state is not physical.
98 Scalar Vm = params.molarVolume(phaseIdx);
99
100 // Calculate b_i / b
101 Scalar bi_b = params.bPure(phaseIdx, compIdx) / params.b(phaseIdx);
102
103 // Calculate the compressibility factor
104 Scalar RT = Constants<Scalar>::R*fs.temperature(phaseIdx);
105 Scalar p = fs.pressure(phaseIdx); // molar volume in [bar]
106 Scalar Z = p*Vm/RT; // compressibility factor
107
108 // Calculate A^* and B^* (see: Reid, p. 42)
109 Scalar Astar = params.a(phaseIdx)*p/(RT*RT);
110 Scalar Bstar = params.b(phaseIdx)*p/(RT);
111
112 // calculate delta_i (see: Reid, p. 145)
113 Scalar sumMoleFractions = 0.0;
114 for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
115 sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
116
117 using std::sqrt;
118 Scalar deltai = 2*sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
119 Scalar tmp = 0;
120 for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
121 tmp +=
122 fs.moleFraction(phaseIdx, compJIdx)
123 / sumMoleFractions
124 * sqrt(params.aPure(phaseIdx, compJIdx))
125 * (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
126 }
127 deltai *= tmp;
128
129 Scalar base =
130 (2*Z + Bstar*(u + sqrt(u*u - 4*w))) /
131 (2*Z + Bstar*(u - sqrt(u*u - 4*w)));
132 Scalar expo = Astar/(Bstar*sqrt(u*u - 4*w))*(bi_b - deltai);
133
134 using std::exp;
135 using std::max;
136 using std::min;
137 using std::pow;
138 Scalar fugCoeff =
139 exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) *
140 pow(base, expo);
141
143 // limit the fugacity coefficient to a reasonable range:
144 //
145 // on one side, we want the mole fraction to be at
146 // least 10^-3 if the fugacity is at the current pressure
147 //
148
149 fugCoeff = min(1e10, fugCoeff);
150 //
151 // on the other hand, if the mole fraction of the component is 100%, we want the
152 // fugacity to be at least 10^-3 Pa
153 //
154 fugCoeff = max(1e-10, fugCoeff);
156 using std::isfinite;
157 if (!isfinite(fugCoeff)) {
158 std::cout << "Non finite phi: " << fugCoeff << "\n";
159 }
160
161 return fugCoeff;
162 }
163
164};
165
166template<class Scalar, class StaticParameters>
167const Scalar PengRobinsonMixture<Scalar, StaticParameters>::u = 2.0;
168template<class Scalar, class StaticParameters>
169const Scalar PengRobinsonMixture<Scalar, StaticParameters>::w = -1.0;
170
171} // end namespace Dumux
172
173#endif
Implements the Peng-Robinson equation of state for liquids and gases.
A central place for various physical constants occuring in some equations.
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
A central place for various physical constants occuring in some equations.
Definition: constants.hh:39
Implements the Peng-Robinson equation of state for liquids and gases.
Definition: pengrobinson.hh:57
Implements the Peng-Robinson equation of state for a mixture.
Definition: pengrobinsonmixture.hh:41
static int computeMolarVolumes(Scalar *Vm, const MutableParams &params, int phaseIdx, const FluidState &fs)
Computes molar volumes where the Peng-Robinson EOS is true.
Definition: pengrobinsonmixture.hh:63
static Scalar computeFugacityCoefficient(const FluidState &fs, const Params &params, int phaseIdx, int compIdx)
Returns the fugacity coefficient of an individual component in the phase.
Definition: pengrobinsonmixture.hh:90