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