PengRobinsonMixture.hpp
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  This file is part of the Open Porous Media project (OPM).
5 
6  OPM 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 2 of the License, or
9  (at your option) any later version.
10 
11  OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_PENG_ROBINSON_MIXTURE_HPP
28 #define OPM_PENG_ROBINSON_MIXTURE_HPP
29 
30 #include "PengRobinson.hpp"
31 
33 
34 #include <iostream>
35 
36 namespace Opm {
41 template <class Scalar, class StaticParameters>
43 {
44  enum { numComponents = StaticParameters::numComponents };
46 
47  // this class cannot be instantiated!
49 
50  // the ideal gas constant
51  static const Scalar R;
52 
53  // the u and w parameters as given by the Peng-Robinson EOS
54  static const Scalar u;
55  static const Scalar w;
56 
57 public:
64  template <class MutableParams, class FluidState>
65  static int computeMolarVolumes(Scalar* Vm,
66  const MutableParams& params,
67  unsigned phaseIdx,
68  const FluidState& fs)
69  {
70  return PengRobinson::computeMolarVolumes(Vm, params, phaseIdx, fs);
71  }
72 
90  template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
91  static LhsEval computeFugacityCoefficient(const FluidState& fs,
92  const Params& params,
93  unsigned phaseIdx,
94  unsigned compIdx)
95  {
96  // note that we normalize the component mole fractions, so
97  // that their sum is 100%. This increases numerical stability
98  // considerably if the fluid state is not physical.
99  LhsEval Vm = params.molarVolume(phaseIdx);
100 
101  // Calculate b_i / b
102  LhsEval bi_b = params.bPure(phaseIdx, compIdx) / params.b(phaseIdx);
103 
104  // Calculate the compressibility factor
105  LhsEval RT = R*fs.temperature(phaseIdx);
106  LhsEval p = fs.pressure(phaseIdx); // molar volume in [bar]
107  LhsEval Z = p*Vm/RT; // compressibility factor
108 
109  // Calculate A^* and B^* (see: Reid, p. 42)
110  LhsEval Astar = params.a(phaseIdx)*p/(RT*RT);
111  LhsEval Bstar = params.b(phaseIdx)*p/(RT);
112 
113  // calculate delta_i (see: Reid, p. 145)
114  LhsEval sumMoleFractions = 0.0;
115  for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx)
116  sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
117  LhsEval deltai = 2*Opm::sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
118  LhsEval tmp = 0;
119  for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
120  tmp +=
121  fs.moleFraction(phaseIdx, compJIdx)
122  / sumMoleFractions
123  * Opm::sqrt(params.aPure(phaseIdx, compJIdx))
124  * (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
125  };
126  deltai *= tmp;
127 
128  LhsEval base =
129  (2*Z + Bstar*(u + std::sqrt(u*u - 4*w))) /
130  (2*Z + Bstar*(u - std::sqrt(u*u - 4*w)));
131  LhsEval expo = Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);
132 
133  LhsEval fugCoeff =
134  Opm::exp(bi_b*(Z - 1))/Opm::max(1e-9, Z - Bstar) *
135  Opm::pow(base, expo);
136 
138  // limit the fugacity coefficient to a reasonable range:
139  //
140  // on one side, we want the mole fraction to be at
141  // least 10^-3 if the fugacity is at the current pressure
142  //
143  fugCoeff = Opm::min(1e10, fugCoeff);
144  //
145  // on the other hand, if the mole fraction of the component is 100%, we want the
146  // fugacity to be at least 10^-3 Pa
147  //
148  fugCoeff = Opm::max(1e-10, fugCoeff);
150 
151  return fugCoeff;
152  }
153 
154 };
155 
156 template <class Scalar, class StaticParameters>
157 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::R = Opm::Constants<Scalar>::R;
158 template<class Scalar, class StaticParameters>
159 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::u = 2.0;
160 template<class Scalar, class StaticParameters>
161 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::w = -1.0;
162 
163 } // namespace Opm
164 
165 #endif
Implements the Peng-Robinson equation of state for liquids and gases.
Definition: PengRobinson.hpp:56
static int computeMolarVolumes(Scalar *Vm, const MutableParams &params, unsigned phaseIdx, const FluidState &fs)
Computes molar volumes where the Peng-Robinson EOS is true.
Definition: PengRobinsonMixture.hpp:65
Definition: Air_Mesitylene.hpp:33
Implements the Peng-Robinson equation of state for liquids and gases.
A central place for various physical constants occuring in some equations.
Definition: Constants.hpp:40
Implements the Peng-Robinson equation of state for a mixture.
Definition: PengRobinsonMixture.hpp:42
A central place for various physical constants occuring in some equations.
static LhsEval computeFugacityCoefficient(const FluidState &fs, const Params &params, unsigned phaseIdx, unsigned compIdx)
Returns the fugacity coefficient of an individual component in the phase.
Definition: PengRobinsonMixture.hpp:91