My Project
PengRobinson.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 */
28 #ifndef OPM_PENG_ROBINSON_HPP
29 #define OPM_PENG_ROBINSON_HPP
30 
34 
36 
37 #include <csignal>
38 
39 namespace Opm {
40 
54 template <class Scalar>
56 {
58  static const Scalar R;
59 
60  PengRobinson()
61  { }
62 
63 public:
64  static void init(Scalar /*aMin*/, Scalar /*aMax*/, unsigned /*na*/,
65  Scalar /*bMin*/, Scalar /*bMax*/, unsigned /*nb*/)
66  {
67 /*
68  // resize the tabulation for the critical points
69  criticalTemperature_.resize(aMin, aMax, na, bMin, bMax, nb);
70  criticalPressure_.resize(aMin, aMax, na, bMin, bMax, nb);
71  criticalMolarVolume_.resize(aMin, aMax, na, bMin, bMax, nb);
72 
73  Scalar VmCrit, pCrit, TCrit;
74  for (unsigned i = 0; i < na; ++i) {
75  Scalar a = criticalTemperature_.iToX(i);
76  assert(std::abs(criticalTemperature_.xToI(criticalTemperature_.iToX(i)) - i) < 1e-10);
77 
78  for (unsigned j = 0; j < nb; ++j) {
79  Scalar b = criticalTemperature_.jToY(j);
80  assert(std::abs(criticalTemperature_.yToJ(criticalTemperature_.jToY(j)) - j) < 1e-10);
81 
82  findCriticalPoint_(TCrit, pCrit, VmCrit, a, b);
83 
84  criticalTemperature_.setSamplePoint(i, j, TCrit);
85  criticalPressure_.setSamplePoint(i, j, pCrit);
86  criticalMolarVolume_.setSamplePoint(i, j, VmCrit);
87  }
88  }
89  */
90  }
91 
100  template <class Evaluation, class Params>
101  static Evaluation computeVaporPressure(const Params& params, const Evaluation& T)
102  {
103  typedef typename Params::Component Component;
106 
107  // initial guess of the vapor pressure
108  Evaluation Vm[3];
109  const Scalar eps = Component::criticalPressure()*1e-10;
110 
111  // use the Ambrose-Walton method to get an initial guess of
112  // the vapor pressure
113  Evaluation pVap = ambroseWalton_(params, T);
114 
115  // Newton-Raphson method
116  for (unsigned i = 0; i < 5; ++i) {
117  // calculate the molar densities
118  assert(molarVolumes(Vm, params, T, pVap) == 3);
119 
120  const Evaluation& f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]);
121  Evaluation df_dp =
122  fugacityDifference_(params, T, pVap + eps, Vm[0], Vm[2])
123  -
124  fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]);
125  df_dp /= 2*eps;
126 
127  const Evaluation& delta = f/df_dp;
128  pVap = pVap - delta;
129 
130  if (std::abs(scalarValue(delta/pVap)) < 1e-10)
131  break;
132  }
133 
134  return pVap;
135  }
136 
141  template <class FluidState, class Params>
142  static
143  typename FluidState::Scalar
144  computeMolarVolume(const FluidState& fs,
145  Params& params,
146  unsigned phaseIdx,
147  bool isGasPhase)
148  {
149  Valgrind::CheckDefined(fs.temperature(phaseIdx));
150  Valgrind::CheckDefined(fs.pressure(phaseIdx));
151 
152  typedef typename FluidState::Scalar Evaluation;
153 
154  Evaluation Vm = 0;
155  Valgrind::SetUndefined(Vm);
156 
157  const Evaluation& T = fs.temperature(phaseIdx);
158  const Evaluation& p = fs.pressure(phaseIdx);
159 
160  const Evaluation& a = params.a(phaseIdx); // "attractive factor"
161  const Evaluation& b = params.b(phaseIdx); // "co-volume"
162 
163  if (!std::isfinite(scalarValue(a))
164  || std::abs(scalarValue(a)) < 1e-30)
165  return std::numeric_limits<Scalar>::quiet_NaN();
166  if (!std::isfinite(scalarValue(b)) || b <= 0)
167  return std::numeric_limits<Scalar>::quiet_NaN();
168 
169  const Evaluation& RT= Constants<Scalar>::R*T;
170  const Evaluation& Astar = a*p/(RT*RT);
171  const Evaluation& Bstar = b*p/RT;
172 
173  const Evaluation& a1 = 1.0;
174  const Evaluation& a2 = - (1 - Bstar);
175  const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
176  const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
177 
178  // ignore the first two results if the smallest
179  // compressibility factor is <= 0.0. (this means that if we
180  // would get negative molar volumes for the liquid phase, we
181  // consider the liquid phase non-existant.)
182  Evaluation Z[3] = {0.0,0.0,0.0};
183  Valgrind::CheckDefined(a1);
184  Valgrind::CheckDefined(a2);
185  Valgrind::CheckDefined(a3);
186  Valgrind::CheckDefined(a4);
187 
188  int numSol = cubicRoots(Z, a1, a2, a3, a4);
189  if (numSol == 3) {
190  // the EOS has three intersections with the pressure,
191  // i.e. the molar volume of gas is the largest one and the
192  // molar volume of liquid is the smallest one
193  if (isGasPhase)
194  Vm = max(1e-7, Z[2]*RT/p);
195  else
196  Vm = max(1e-7, Z[0]*RT/p);
197  }
198  else if (numSol == 1) {
199  // the EOS only has one intersection with the pressure,
200  // for the other phase, we take the extremum of the EOS
201  // with the largest distance from the intersection.
202  Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
203  Vm = VmCubic;
204 
205  // find the extrema (if they are present)
206  Evaluation Vmin, Vmax, pmin, pmax;
207  if (findExtrema_(Vmin, Vmax,
208  pmin, pmax,
209  a, b, T))
210  {
211  if (isGasPhase)
212  Vm = std::max(Vmax, VmCubic);
213  else {
214  if (Vmin > 0)
215  Vm = std::min(Vmin, VmCubic);
216  else
217  Vm = VmCubic;
218  }
219  }
220  else {
221  // the EOS does not exhibit any physically meaningful
222  // extrema, and the fluid is critical...
223  Vm = VmCubic;
224  handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
225  }
226  }
227 
228  Valgrind::CheckDefined(Vm);
229  assert(std::isfinite(scalarValue(Vm)));
230  assert(Vm > 0);
231  return Vm;
232  }
233 
244  template <class Evaluation, class Params>
245  static Evaluation computeFugacityCoeffient(const Params& params)
246  {
247  const Evaluation& T = params.temperature();
248  const Evaluation& p = params.pressure();
249  const Evaluation& Vm = params.molarVolume();
250 
251  const Evaluation& RT = Constants<Scalar>::R*T;
252  const Evaluation& Z = p*Vm/RT;
253  const Evaluation& Bstar = p*params.b() / RT;
254 
255  const Evaluation& tmp =
256  (Vm + params.b()*(1 + std::sqrt(2))) /
257  (Vm + params.b()*(1 - std::sqrt(2)));
258  const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
259  const Evaluation& fugCoeff =
260  exp(Z - 1) / (Z - Bstar) *
261  pow(tmp, expo);
262 
263  return fugCoeff;
264  }
265 
276  template <class Evaluation, class Params>
277  static Evaluation computeFugacity(const Params& params)
278  { return params.pressure()*computeFugacityCoeff(params); }
279 
280 protected:
281  template <class FluidState, class Params, class Evaluation = typename FluidState::Scalar>
282  static void handleCriticalFluid_(Evaluation& Vm,
283  const FluidState& /*fs*/,
284  const Params& params,
285  unsigned phaseIdx,
286  bool isGasPhase)
287  {
288  Evaluation Tcrit, pcrit, Vcrit;
289  findCriticalPoint_(Tcrit,
290  pcrit,
291  Vcrit,
292  params.a(phaseIdx),
293  params.b(phaseIdx));
294 
295 
296  //Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx));
297 
298  if (isGasPhase)
299  Vm = max(Vm, Vcrit);
300  else
301  Vm = min(Vm, Vcrit);
302  }
303 
304  template <class Evaluation>
305  static void findCriticalPoint_(Evaluation& Tcrit,
306  Evaluation& pcrit,
307  Evaluation& Vcrit,
308  const Evaluation& a,
309  const Evaluation& b)
310  {
311  Evaluation minVm(0);
312  Evaluation maxVm(1e30);
313 
314  Evaluation minP(0);
315  Evaluation maxP(1e30);
316 
317  // first, we need to find an isotherm where the EOS exhibits
318  // a maximum and a minimum
319  Evaluation Tmin = 250; // [K]
320  for (unsigned i = 0; i < 30; ++i) {
321  bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
322  if (hasExtrema)
323  break;
324  Tmin /= 2;
325  };
326 
327  Evaluation T = Tmin;
328 
329  // Newton's method: Start at minimum temperature and minimize
330  // the "gap" between the extrema of the EOS
331  unsigned iMax = 100;
332  for (unsigned i = 0; i < iMax; ++i) {
333  // calculate function we would like to minimize
334  Evaluation f = maxVm - minVm;
335 
336  // check if we're converged
337  if (f < 1e-10 || (i == iMax - 1 && f < 1e-8)) {
338  Tcrit = T;
339  pcrit = (minP + maxP)/2;
340  Vcrit = (maxVm + minVm)/2;
341  return;
342  }
343 
344  // backward differences. Forward differences are not
345  // robust, since the isotherm could be critical if an
346  // epsilon was added to the temperature. (this is case
347  // rarely happens, though)
348  const Scalar eps = - 1e-11;
349  assert(findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps));
350  assert(std::isfinite(scalarValue(maxVm)));
351  Evaluation fStar = maxVm - minVm;
352 
353  // derivative of the difference between the maximum's
354  // molar volume and the minimum's molar volume regarding
355  // temperature
356  Evaluation fPrime = (fStar - f)/eps;
357  if (std::abs(scalarValue(fPrime)) < 1e-40) {
358  Tcrit = T;
359  pcrit = (minP + maxP)/2;
360  Vcrit = (maxVm + minVm)/2;
361  return;
362  }
363 
364  // update value for the current iteration
365  Evaluation delta = f/fPrime;
366  assert(std::isfinite(scalarValue(delta)));
367  if (delta > 0)
368  delta = -10;
369 
370  // line search (we have to make sure that both extrema
371  // still exist after the update)
372  for (unsigned j = 0; ; ++j) {
373  if (j >= 20) {
374  if (f < 1e-8) {
375  Tcrit = T;
376  pcrit = (minP + maxP)/2;
377  Vcrit = (maxVm + minVm)/2;
378  return;
379  }
380 
381  std::ostringstream oss;
382  oss << "Could not determine the critical point for a=" << a << ", b=" << b;
383  throw NumericalIssue(oss.str());
384  }
385 
386  if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
387  // if the isotherm for T - delta exhibits two
388  // extrema the update is finished
389  T -= delta;
390  break;
391  }
392  else
393  delta /= 2;
394  };
395  }
396 
397  assert(false);
398  }
399 
400  // find the two molar volumes where the EOS exhibits extrema and
401  // which are larger than the covolume of the phase
402  template <class Evaluation>
403  static bool findExtrema_(Evaluation& Vmin,
404  Evaluation& Vmax,
405  Evaluation& /*pMin*/,
406  Evaluation& /*pMax*/,
407  const Evaluation& a,
408  const Evaluation& b,
409  const Evaluation& T)
410  {
411  Scalar u = 2;
412  Scalar w = -1;
413 
414  const Evaluation& RT = Constants<Scalar>::R*T;
415  // calculate coefficients of the 4th order polynominal in
416  // monomial basis
417  const Evaluation& a1 = RT;
418  const Evaluation& a2 = 2*RT*u*b - 2*a;
419  const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
420  const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
421  const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
422 
423  assert(std::isfinite(scalarValue(a1)));
424  assert(std::isfinite(scalarValue(a2)));
425  assert(std::isfinite(scalarValue(a3)));
426  assert(std::isfinite(scalarValue(a4)));
427  assert(std::isfinite(scalarValue(a5)));
428 
429  // Newton method to find first root
430 
431  // if the values which we got on Vmin and Vmax are usefull, we
432  // will reuse them as initial value, else we will start 10%
433  // above the covolume
434  Evaluation V = b*1.1;
435  Evaluation delta = 1.0;
436  for (unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
437  const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
438  const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
439 
440  if (std::abs(scalarValue(fPrime)) < 1e-20) {
441  // give up if the derivative is zero
442  return false;
443  }
444 
445 
446  delta = f/fPrime;
447  V -= delta;
448 
449  if (i > 200) {
450  // give up after 200 iterations...
451  return false;
452  }
453  }
454  assert(std::isfinite(scalarValue(V)));
455 
456  // polynomial division
457  Evaluation b1 = a1;
458  Evaluation b2 = a2 + V*b1;
459  Evaluation b3 = a3 + V*b2;
460  Evaluation b4 = a4 + V*b3;
461 
462  // invert resulting cubic polynomial analytically
463  Evaluation allV[4];
464  allV[0] = V;
465  int numSol = 1 + invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
466 
467  // sort all roots of the derivative
468  std::sort(allV + 0, allV + numSol);
469 
470  // check whether the result is physical
471  if (numSol != 4 || allV[numSol - 2] < b) {
472  // the second largest extremum is smaller than the phase's
473  // covolume which is physically impossible
474  return false;
475  }
476 
477 
478  // it seems that everything is okay...
479  Vmin = allV[numSol - 2];
480  Vmax = allV[numSol - 1];
481  return true;
482  }
483 
496  template <class Evaluation, class Params>
497  static Evaluation ambroseWalton_(const Params& /*params*/, const Evaluation& T)
498  {
499  typedef typename Params::Component Component;
500 
501  const Evaluation& Tr = T / Component::criticalTemperature();
502  const Evaluation& tau = 1 - Tr;
503  const Evaluation& omega = Component::acentricFactor();
504 
505  const Evaluation& f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
506  const Evaluation& f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
507  const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
508 
509  return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
510  }
511 
522  template <class Evaluation, class Params>
523  static Evaluation fugacityDifference_(const Params& params,
524  const Evaluation& T,
525  const Evaluation& p,
526  const Evaluation& VmLiquid,
527  const Evaluation& VmGas)
528  { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
529 
530 
531 };
532 
533 } // namespace Opm
534 
535 #endif
Relations valid for an ideal gas.
Provides free functions to invert polynomials of degree 1, 2 and 3.
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:331
This is a fluid state which allows to set the fluid temperatures and takes all other quantities from ...
Implements a scalar function that depends on two variables and which is sampled on an uniform X-Y gri...
Abstract base class of a pure chemical species.
Definition: Component.hpp:42
static Scalar acentricFactor()
Returns the acentric factor of the component.
Definition: Component.hpp:110
static Scalar criticalPressure()
Returns the critical pressure in of the component.
Definition: Component.hpp:103
static Scalar criticalTemperature()
Returns the critical temperature in of the component.
Definition: Component.hpp:97
A central place for various physical constants occuring in some equations.
Definition: Constants.hpp:41
static const Scalar R
The ideal gas constant [J/(mol K)].
Definition: Constants.hpp:45
Implements the Peng-Robinson equation of state for liquids and gases.
Definition: PengRobinson.hpp:56
static Evaluation computeFugacityCoeffient(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:245
static Evaluation fugacityDifference_(const Params &params, const Evaluation &T, const Evaluation &p, const Evaluation &VmLiquid, const Evaluation &VmGas)
Returns the difference between the liquid and the gas phase fugacities in [bar].
Definition: PengRobinson.hpp:523
static Evaluation computeVaporPressure(const Params &params, const Evaluation &T)
Predicts the vapor pressure for the temperature given in setTP().
Definition: PengRobinson.hpp:101
static Evaluation ambroseWalton_(const Params &, const Evaluation &T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition: PengRobinson.hpp:497
static FluidState::Scalar computeMolarVolume(const FluidState &fs, Params &params, unsigned phaseIdx, bool isGasPhase)
Computes molar volumes where the Peng-Robinson EOS is true.
Definition: PengRobinson.hpp:144
static Evaluation computeFugacity(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: PengRobinson.hpp:277