My Project
PengRobinsonParamsMixture.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_PARAMS_MIXTURE_HPP
28 #define OPM_PENG_ROBINSON_PARAMS_MIXTURE_HPP
29 
30 #include "PengRobinsonParams.hpp"
31 
34 
35 #include <algorithm>
36 
37 namespace Opm
38 {
39 
57 template <class Scalar, class FluidSystem, unsigned phaseIdx, bool useSpe5Relations=false>
59  : public PengRobinsonParams<Scalar>
60 {
61  enum { numComponents = FluidSystem::numComponents };
62 
63  // Peng-Robinson parameters for pure substances
65 
67 
68 public:
69 
73  Scalar getaCache(unsigned compIIdx, unsigned compJIdx ) const
74  {
75  return aCache_[compIIdx][compJIdx];
76  }
77 
81  template <class FluidState>
82  void updatePure(const FluidState& fluidState)
83  {
84  updatePure(fluidState.temperature(phaseIdx),
85  fluidState.pressure(phaseIdx));
86  }
87 
93  void updatePure(Scalar temperature, Scalar pressure)
94  {
95  Valgrind::CheckDefined(temperature);
96  Valgrind::CheckDefined(pressure);
97 
98  // Calculate the Peng-Robinson parameters of the pure
99  // components
100  //
101  // See: R. Reid, et al.: The Properties of Gases and Liquids,
102  // 4th edition, McGraw-Hill, 1987, p. 43
103  for (unsigned i = 0; i < numComponents; ++i) {
104  Scalar pc = FluidSystem::criticalPressure(i);
105  Scalar omega = FluidSystem::acentricFactor(i);
106  Scalar Tr = temperature/FluidSystem::criticalTemperature(i);
107  Scalar RTc = Constants<Scalar>::R*FluidSystem::criticalTemperature(i);
108 
109  Scalar f_omega;
110 
111  if (omega < 0.49)
112  f_omega = 0.37464 + omega*(1.54226 + omega*(-0.26992));
113  else
114  f_omega = 0.379642 + omega*(1.48503 + omega*(-0.164423 + omega*0.016666));
115 
116  Valgrind::CheckDefined(f_omega);
117 
118  Scalar tmp = 1 + f_omega*(1 - sqrt(Tr));
119  tmp = tmp*tmp;
120 
121  Scalar newA = 0.4572355*RTc*RTc/pc * tmp;
122  Scalar newB = 0.0777961 * RTc / pc;
123  assert(std::isfinite(scalarValue(newA)));
124  assert(std::isfinite(scalarValue(newB)));
125 
126  this->pureParams_[i].setA(newA);
127  this->pureParams_[i].setB(newB);
128  Valgrind::CheckDefined(this->pureParams_[i].a());
129  Valgrind::CheckDefined(this->pureParams_[i].b());
130  }
131 
132  updateACache_();
133  }
134 
142  template <class FluidState>
143  void updateMix(const FluidState& fs)
144  {
145  using FlashEval = typename FluidState::Scalar;
146  Scalar sumx = 0.0;
147  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
148  sumx += fs.moleFraction(phaseIdx, compIdx);
149  sumx = std::max(Scalar(1e-10), sumx);
150 
151  // Calculate the Peng-Robinson parameters of the mixture
152  //
153  // See: R. Reid, et al.: The Properties of Gases and Liquids,
154  // 4th edition, McGraw-Hill, 1987, p. 82
155  FlashEval newA = 0;
156  FlashEval newB = 0;
157  for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
158  const FlashEval moleFracI = fs.moleFraction(phaseIdx, compIIdx);
159  FlashEval xi = max(0.0, min(1.0, moleFracI));
160  Valgrind::CheckDefined(xi);
161 
162  for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
163  const FlashEval moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
164  FlashEval xj = max(0.0, min(1.0, moleFracJ));
165  Valgrind::CheckDefined(xj);
166 
167  // mixing rule from Reid, page 82
168  newA += xi * xj * aCache_[compIIdx][compJIdx];
169 
170  assert(std::isfinite(scalarValue(newA)));
171  }
172 
173  // mixing rule from Reid, page 82
174  newB += max(0.0, xi) * this->pureParams_[compIIdx].b();
175  assert(std::isfinite(scalarValue(newB)));
176  }
177 
178  // assert(newB > 0);
179  this->setA(newA);
180  this->setB(newB);
181 
182  Valgrind::CheckDefined(this->a());
183  Valgrind::CheckDefined(this->b());
184 
185  }
186 
195  template <class FluidState>
196  void updateSingleMoleFraction(const FluidState& fs,
197  unsigned /*compIdx*/)
198  {
199  updateMix(fs);
200  }
201 
205  const PureParams& pureParams(unsigned compIdx) const
206  { return pureParams_[compIdx]; }
207 
211  const PureParams& operator[](unsigned compIdx) const
212  {
213  assert(0 <= compIdx && compIdx < numComponents);
214  return pureParams_[compIdx];
215  }
216 
221  void checkDefined() const
222  {
223 #ifndef NDEBUG
224  for (unsigned i = 0; i < numComponents; ++i)
225  pureParams_[i].checkDefined();
226 
227  Valgrind::CheckDefined(this->a());
228  Valgrind::CheckDefined(this->b());
229 #endif
230  }
231 
232 protected:
233  PureParams pureParams_[numComponents];
234 
235 private:
236  void updateACache_()
237  {
238  for (unsigned compIIdx = 0; compIIdx < numComponents; ++ compIIdx) {
239  for (unsigned compJIdx = 0; compJIdx < numComponents; ++ compJIdx) {
240  // interaction coefficient as given in SPE5
241  Scalar Psi = FluidSystem::interactionCoefficient(compIIdx, compJIdx);
242 
243  aCache_[compIIdx][compJIdx] =
244  sqrt(this->pureParams_[compIIdx].a()
245  * this->pureParams_[compJIdx].a())
246  * (1 - Psi);
247  }
248  }
249  }
250 
251  Scalar aCache_[numComponents][numComponents];
252 };
253 
254 
255 } // namespace Opm
256 
257 #endif
A central place for various physical constants occuring in some equations.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Stores and provides access to the Peng-Robinson parameters.
A central place for various physical constants occuring in some equations.
Definition: Constants.hpp:41
The mixing rule for the oil and the gas phases of the SPE5 problem.
Definition: PengRobinsonParamsMixture.hpp:60
void updateSingleMoleFraction(const FluidState &fs, unsigned)
Calculates the "a" and "b" Peng-Robinson parameters for the mixture provided that only a single mole ...
Definition: PengRobinsonParamsMixture.hpp:196
void updatePure(Scalar temperature, Scalar pressure)
Peng-Robinson parameters for the pure components.
Definition: PengRobinsonParamsMixture.hpp:93
void updatePure(const FluidState &fluidState)
Update Peng-Robinson parameters for the pure components.
Definition: PengRobinsonParamsMixture.hpp:82
const PureParams & operator[](unsigned compIdx) const
Returns the Peng-Robinson parameters for a pure component.
Definition: PengRobinsonParamsMixture.hpp:211
Scalar getaCache(unsigned compIIdx, unsigned compJIdx) const
TODO.
Definition: PengRobinsonParamsMixture.hpp:73
void checkDefined() const
If run under valgrind, this method produces an warning if the parameters where not determined correct...
Definition: PengRobinsonParamsMixture.hpp:221
void updateMix(const FluidState &fs)
Calculates the "a" and "b" Peng-Robinson parameters for the mixture.
Definition: PengRobinsonParamsMixture.hpp:143
const PureParams & pureParams(unsigned compIdx) const
Return the Peng-Robinson parameters of a pure substance,.
Definition: PengRobinsonParamsMixture.hpp:205
Stores and provides access to the Peng-Robinson parameters.
Definition: PengRobinsonParams.hpp:44
Scalar a() const
Returns the attractive parameter 'a' of the Peng-Robinson fluid.
Definition: PengRobinsonParams.hpp:50
void setB(Scalar value)
Set the repulsive parameter 'b' of the Peng-Robinson fluid.
Definition: PengRobinsonParams.hpp:83
Scalar b() const
Returns the repulsive parameter 'b' of the Peng-Robinson fluid.
Definition: PengRobinsonParams.hpp:57
void setA(Scalar value)
Set the attractive parameter 'a' of the Peng-Robinson fluid.
Definition: PengRobinsonParams.hpp:76
Definition: MathToolbox.hpp:50