My Project
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 namespace Opm {
39 template <class Scalar, class StaticParameters>
41 {
42  enum { numComponents = StaticParameters::numComponents };
43  typedef ::Opm::PengRobinson<Scalar> PengRobinson;
44 
45  // this class cannot be instantiated!
46  PengRobinsonMixture() = default;
47 
48  // the ideal gas constant
49  static const Scalar R;
50 
51  // the u and w parameters as given by the Peng-Robinson EOS
52  static const Scalar u;
53  static const Scalar w;
54 
55 public:
62  template <class MutableParams, class FluidState>
63  static int computeMolarVolumes(Scalar* Vm,
64  const MutableParams& params,
65  unsigned phaseIdx,
66  const FluidState& fs)
67  {
68  return PengRobinson::computeMolarVolumes(Vm, params, phaseIdx, fs);
69  }
70 
88  template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
89  static LhsEval computeFugacityCoefficient(const FluidState& fs,
90  const Params& params,
91  unsigned phaseIdx,
92  unsigned compIdx)
93  {
94  // note that we normalize the component mole fractions, so
95  // that their sum is 100%. This increases numerical stability
96  // considerably if the fluid state is not physical.
97  LhsEval Vm = params.molarVolume(phaseIdx);
98 
99  // Calculate b_i / b
100  LhsEval bi_b = params.bPure(phaseIdx, compIdx) / params.b(phaseIdx);
101 
102  // Calculate the compressibility factor
103  LhsEval RT = R*fs.temperature(phaseIdx);
104  LhsEval p = fs.pressure(phaseIdx); // molar volume in [bar]
105  LhsEval Z = p*Vm/RT; // compressibility factor
106 
107  // Calculate A^* and B^* (see: Reid, p. 42)
108  LhsEval Astar = params.a(phaseIdx)*p/(RT*RT);
109  LhsEval Bstar = params.b(phaseIdx)*p/(RT);
110 
111  LhsEval A_s = 0.0;
112  for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
113  A_s += params.aCache(phaseIdx, compIdx, compJIdx) * fs.moleFraction(phaseIdx, compJIdx) * p / (RT * RT);
114  }
115 
116  LhsEval alpha;
117  LhsEval betta;
118  LhsEval gamma;
119  LhsEval ln_phi;
120  LhsEval fugCoeff;
121 
122  Scalar m1;
123  Scalar m2;
124 
125  m1 = 0.5*(u + std::sqrt(u*u - 4*w));
126  m2 = 0.5*(u - std::sqrt(u*u - 4*w));
127 
128  alpha = -log(Z - Bstar) + bi_b * (Z - 1);
129  betta = log((Z + m2 * Bstar) / (Z + m1 * Bstar)) * Astar / ((m1 - m2) * Bstar);
130  gamma = (2 / Astar ) * A_s - bi_b;
131  ln_phi = alpha + (betta * gamma);
132 
133  fugCoeff = exp(ln_phi);
134 
136  // limit the fugacity coefficient to a reasonable range:
137  //
138  // on one side, we want the mole fraction to be at
139  // least 10^-3 if the fugacity is at the current pressure
140  //
141  fugCoeff = min(1e10, fugCoeff);
142  //
143  // on the other hand, if the mole fraction of the component is 100%, we want the
144  // fugacity to be at least 10^-3 Pa
145  //
146  fugCoeff = max(1e-10, fugCoeff);
148 
149  return fugCoeff;
150  }
151 
152 };
153 
154 template <class Scalar, class StaticParameters>
155 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::R = Constants<Scalar>::R;
156 template<class Scalar, class StaticParameters>
157 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::u = 2.0;
158 template<class Scalar, class StaticParameters>
159 const Scalar PengRobinsonMixture<Scalar, StaticParameters>::w = -1.0;
160 
161 } // namespace Opm
162 
163 #endif
A central place for various physical constants occuring in some equations.
Implements the Peng-Robinson equation of state for liquids and gases.
static const Scalar R
The ideal gas constant [J/(mol K)].
Definition: Constants.hpp:45
Implements the Peng-Robinson equation of state for a mixture.
Definition: PengRobinsonMixture.hpp:41
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:63
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:89
Implements the Peng-Robinson equation of state for liquids and gases.
Definition: PengRobinson.hpp:56