3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
pengrobinson.hh
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program 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 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
31#ifndef DUMUX_PENG_ROBINSON_HH
32#define DUMUX_PENG_ROBINSON_HH
33
37#include <dumux/common/math.hh>
39
40#include <iostream>
41
42namespace Dumux {
43
55template <class Scalar>
57{
59 { }
60
61public:
62 static void init(Scalar aMin, Scalar aMax, int na,
63 Scalar bMin, Scalar bMax, int nb)
64 {
65 // resize the tabulation for the critical points
66 criticalTemperature_.resize(aMin, aMax, na, bMin, bMax, nb);
67 criticalPressure_.resize(aMin, aMax, na, bMin, bMax, nb);
68 criticalMolarVolume_.resize(aMin, aMax, na, bMin, bMax, nb);
69
70 Scalar VmCrit = 18e-6;
71 Scalar pCrit = 1e7;
72 Scalar TCrit = 273;
73 for (int i = 0; i < na; ++i) {
74 Scalar a = criticalTemperature_.iToX(i);
75 for (int j = 0; j < nb; ++j) {
76 Scalar b = criticalTemperature_.jToY(j);
77
78 findCriticalPoint_(TCrit, pCrit, VmCrit, a, b);
79
80 criticalTemperature_.setSamplePoint(i, j, TCrit);
81 criticalPressure_.setSamplePoint(i, j, pCrit);
82 criticalMolarVolume_.setSamplePoint(i, j, VmCrit);
83 }
84 }
85 };
86
97 template <class Params>
98 static Scalar computeVaporPressure(const Params &params, Scalar T)
99 {
100 using Component = typename Params::Component;
101 if (T >= Component::criticalTemperature())
102 return Component::criticalPressure();
103
104 // initial guess of the vapor pressure
105 Scalar Vm[3];
106 const Scalar eps = Component::criticalPressure()*1e-10;
107
108 // use the Ambrose-Walton method to get an initial guess of
109 // the vapor pressure
110 Scalar pVap = ambroseWalton_(params, T);
111
112 // Newton-Raphson method
113 for (int i = 0; i < 5; ++i) {
114 // calculate the molar densities
115 assert(molarVolumes(Vm, params, T, pVap) == 3);
116
117 Scalar f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]);
118 Scalar df_dp =
119 fugacityDifference_(params, T, pVap + eps, Vm[0], Vm[2])
120 -
121 fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]);
122 df_dp /= 2*eps;
123
124 Scalar delta = f/df_dp;
125 pVap = pVap - delta;
126 using std::abs;
127 if (abs(delta/pVap) < 1e-10)
128 break;
129 }
130
131 return pVap;
132 }
133
142 template <class FluidState, class Params>
143 static Scalar computeMolarVolume(const FluidState &fs,
144 Params &params,
145 int phaseIdx,
146 bool isGasPhase)
147 {
148 Scalar Vm = 0;
149
150 Scalar T = fs.temperature(phaseIdx);
151 Scalar p = fs.pressure(phaseIdx);
152
153 Scalar a = params.a(phaseIdx); // "attractive factor"
154 Scalar b = params.b(phaseIdx); // "co-volume"
155
156 Scalar RT = Constants<Scalar>::R*T;
157 Scalar Astar = a*p/(RT*RT);
158 Scalar Bstar = b*p/RT;
159
160 Scalar a1 = 1.0;
161 Scalar a2 = - (1 - Bstar);
162 Scalar a3 = Astar - Bstar*(3*Bstar + 2);
163 Scalar a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
164
165 // ignore the first two results if the smallest
166 // compressibility factor is <= 0.0. (this means that if we
167 // would get negative molar volumes for the liquid phase, we
168 // consider the liquid phase non-existant.)
169 Scalar Z[3] = {0.0,0.0,0.0};
170 int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4);
171 if (numSol == 3) {
172 // the EOS has three intersections with the pressure,
173 // i.e. the molar volume of gas is the largest one and the
174 // molar volume of liquid is the smallest one
175 if (isGasPhase)
176 Vm = Z[2]*RT/p;
177 else
178 Vm = Z[0]*RT/p;
179 }
180 else if (numSol == 1) {
181 // the EOS only has one intersection with the pressure,
182 // for the other phase, we take the extremum of the EOS
183 // with the largest distance from the intersection.
184 Scalar VmCubic = Z[0]*RT/p;
185
186 if (T > criticalTemperature_(a, b)) {
187 // if the EOS does not exhibit any extrema, the fluid
188 // is critical...
189 Vm = VmCubic;
190 handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
191 }
192 else {
193 // find the extrema (if they are present)
194 Scalar Vmin, Vmax, pmin, pmax;
195 using std::min;
196 using std::max;
197 if (findExtrema_(Vmin, Vmax,
198 pmin, pmax,
199 a, b, T))
200 {
201 if (isGasPhase)
202 Vm = max(Vmax, VmCubic);
203 else
204 Vm = min(Vmin, VmCubic);
205 }
206 else
207 Vm = VmCubic;
208 }
209 }
210
211 using std::isfinite;
212 assert(isfinite(Vm) && Vm > 0);
213 return Vm;
214 }
215
226 template <class Params>
227 static Scalar computeFugacityCoeffient(const Params &params)
228 {
229 Scalar T = params.temperature();
230 Scalar p = params.pressure();
231 Scalar Vm = params.molarVolume();
232
233 Scalar RT = Constants<Scalar>::R*T;
234 Scalar Z = p*Vm/RT;
235 Scalar Bstar = p*params.b() / RT;
236 using std::sqrt;
237 Scalar tmp =
238 (Vm + params.b()*(1 + sqrt(2))) /
239 (Vm + params.b()*(1 - sqrt(2)));
240 Scalar expo = - params.a()/(RT * 2 * params.b() * sqrt(2));
241 using std::exp;
242 using std::pow;
243 Scalar fugCoeff =
244 exp(Z - 1) / (Z - Bstar) *
245 pow(tmp, expo);
246
247 return fugCoeff;
248 }
249
260 template <class Params>
261 static Scalar computeFugacity(const Params &params)
262 { return params.pressure()*computeFugacityCoeff(params); }
263
264protected:
265 template <class FluidState, class Params>
266 static void handleCriticalFluid_(Scalar &Vm,
267 const FluidState &fs,
268 const Params &params,
269 int phaseIdx,
270 bool isGasPhase)
271 {
272 Scalar Vcrit = criticalMolarVolume_(params.a(phaseIdx), params.b(phaseIdx));
273
274 using std::max;
275 using std::min;
276 if (isGasPhase)
277 Vm = max(Vm, Vcrit);
278 else
279 Vm = min(Vm, Vcrit);
280 }
281
282 static void findCriticalPoint_(Scalar &Tcrit,
283 Scalar &pcrit,
284 Scalar &Vcrit,
285 Scalar a,
286 Scalar b)
287 {
288 Scalar minVm;
289 Scalar maxVm;
290
291 Scalar minP(0);
292 Scalar maxP(1e100);
293
294 // first, we need to find an isotherm where the EOS exhibits
295 // a maximum and a minimum
296 Scalar Tmin = 250; // [K]
297 for (int i = 0; i < 30; ++i) {
298 bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
299 if (hasExtrema)
300 break;
301 Tmin /= 2;
302 }
303
304 Scalar T = Tmin;
305
306 // Newton's method: Start at minimum temperature and minimize
307 // the "gap" between the extrema of the EOS
308 for (int i = 0; i < 25; ++i) {
309 // calculate function we would like to minimize
310 Scalar f = maxVm - minVm;
311
312 // check if we're converged
313 if (f < 1e-10) {
314 Tcrit = T;
315 pcrit = (minP + maxP)/2;
316 Vcrit = (maxVm + minVm)/2;
317 return;
318 }
319
320 // backward differences. Forward differences are not
321 // robust, since the isotherm could be critical if an
322 // epsilon was added to the temperature. (this is case
323 // rarely happens, though)
324 const Scalar eps = - 1e-8;
325 bool DUNE_UNUSED hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
326 assert(hasExtrema);
327 Scalar fStar = maxVm - minVm;
328
329 // derivative of the difference between the maximum's
330 // molar volume and the minimum's molar volume regarding
331 // temperature
332 Scalar fPrime = (fStar - f)/eps;
333
334 // update value for the current iteration
335 Scalar delta = f/fPrime;
336 if (delta > 0)
337 delta = -10;
338
339 // line search (we have to make sure that both extrema
340 // still exist after the update)
341 for (int j = 0; ; ++j) {
342 if (j >= 20) {
343 DUNE_THROW(NumericalProblem,
344 "Could not determine the critical point for a=" << a << ", b=" << b);
345 }
346
347 if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
348 // if the isotherm for T - delta exhibits two
349 // extrema the update is finished
350 T -= delta;
351 break;
352 }
353 else
354 delta /= 2;
355 }
356 }
357 }
358
359 // find the two molar volumes where the EOS exhibits extrema and
360 // which are larger than the covolume of the phase
361 static bool findExtrema_(Scalar &Vmin,
362 Scalar &Vmax,
363 Scalar &pMin,
364 Scalar &pMax,
365 Scalar a,
366 Scalar b,
367 Scalar T)
368 {
369 Scalar u = 2;
370 Scalar w = -1;
371
372 Scalar RT = Constants<Scalar>::R*T;
373
374 // calculate coefficients of the 4th order polynominal in
375 // monomial basis
376 Scalar a1 = RT;
377 Scalar a2 = 2*RT*u*b - 2*a;
378 Scalar a3 = 2*RT*w*b*b + RT*u*u*b*b + 4*a*b - u*a*b;
379 Scalar a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
380 Scalar a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
381
382 // Newton method to find first root
383
384 // if the values which we got on Vmin and Vmax are usefull, we
385 // will reuse them as initial value, else we will start 10%
386 // above the covolume
387 Scalar V = b*1.1;
388 Scalar delta = 1.0;
389 using std::abs;
390 for (int i = 0; abs(delta) > 1e-9; ++i) {
391 Scalar f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
392 Scalar fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
393
394 if (abs(fPrime) < 1e-20) {
395 // give up if the derivative is zero
396 Vmin = 0;
397 Vmax = 0;
398 return false;
399 }
400
401
402 delta = f/fPrime;
403 V -= delta;
404
405 if (i > 20) {
406 // give up after 20 iterations...
407 Vmin = 0;
408 Vmax = 0;
409 return false;
410 }
411 }
412
413 // polynomial division
414 Scalar b1 = a1;
415 Scalar b2 = a2 + V*b1;
416 Scalar b3 = a3 + V*b2;
417 Scalar b4 = a4 + V*b3;
418
419 // invert resulting cubic polynomial analytically
420 Scalar allV[4];
421 allV[0] = V;
422 int numSol = 1 + invertCubicPolynomial(allV + 1, b1, b2, b3, b4);
423
424 // sort all roots of the derivative
425 std::sort(allV + 0, allV + numSol);
426
427 // check whether the result is physical
428 if (allV[numSol - 2] < b) {
429 // the second largest extremum is smaller than the phase's
430 // covolume which is physically impossible
431 Vmin = Vmax = 0;
432 return false;
433 }
434
435 // it seems that everything is okay...
436 Vmin = allV[numSol - 2];
437 Vmax = allV[numSol - 1];
438 return true;
439 }
440
451 template <class Params>
452 static Scalar ambroseWalton_(const Params &params, Scalar T)
453 {
454 using Component = typename Params::Component;
455
456 Scalar Tr = T / Component::criticalTemperature();
457 Scalar tau = 1 - Tr;
458 Scalar omega = Component::acentricFactor();
459 using std::sqrt;
460 using std::pow;
461 Scalar f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
462 Scalar f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
463 Scalar f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
464 using std::exp;
465 return Component::criticalPressure()*exp(f0 + omega * (f1 + omega*f2));
466 }
467
478 template <class Params>
479 static Scalar fugacityDifference_(const Params &params,
480 Scalar T,
481 Scalar p,
482 Scalar VmLiquid,
483 Scalar VmGas)
484 { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
485
489};
490
491template <class Scalar>
493
494template <class Scalar>
496
497template <class Scalar>
499
500} // end namespace
501
502#endif
Some exceptions thrown in DuMux
Define some often used mathematical functions.
Implements tabulation for a two-dimensional function.
This is a fluid state which allows to set the fluid temperatures and takes all other quantities from ...
Relations valid for an ideal gas.
int invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: math.hh:242
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
constexpr double eps
epsilon for checking direction of scvf normals
Definition: test_tpfafvgeometry_nonconforming.cc:44
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
Implements tabulation for a two-dimensional function.
Definition: tabulated2dfunction.hh:47
A central place for various physical constants occuring in some equations.
Definition: constants.hh:39
Implements the Peng-Robinson equation of state for liquids and gases.
Definition: pengrobinson.hh:57
static void handleCriticalFluid_(Scalar &Vm, const FluidState &fs, const Params &params, int phaseIdx, bool isGasPhase)
Definition: pengrobinson.hh:266
static Scalar fugacityDifference_(const Params &params, Scalar T, Scalar p, Scalar VmLiquid, Scalar VmGas)
Returns the difference between the liquid and the gas phase fugacities in [bar].
Definition: pengrobinson.hh:479
static Tabulated2DFunction< Scalar > criticalPressure_
Definition: pengrobinson.hh:487
static void init(Scalar aMin, Scalar aMax, int na, Scalar bMin, Scalar bMax, int nb)
Definition: pengrobinson.hh:62
static Tabulated2DFunction< Scalar > criticalMolarVolume_
Definition: pengrobinson.hh:488
static bool findExtrema_(Scalar &Vmin, Scalar &Vmax, Scalar &pMin, Scalar &pMax, Scalar a, Scalar b, Scalar T)
Definition: pengrobinson.hh:361
static void findCriticalPoint_(Scalar &Tcrit, Scalar &pcrit, Scalar &Vcrit, Scalar a, Scalar b)
Definition: pengrobinson.hh:282
static Scalar computeFugacity(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: pengrobinson.hh:261
static Scalar ambroseWalton_(const Params &params, Scalar T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition: pengrobinson.hh:452
static Tabulated2DFunction< Scalar > criticalTemperature_
Definition: pengrobinson.hh:486
static Scalar computeMolarVolume(const FluidState &fs, Params &params, int phaseIdx, bool isGasPhase)
Computes molar volumes where the Peng-Robinson EOS is true.
Definition: pengrobinson.hh:143
static Scalar computeVaporPressure(const Params &params, Scalar T)
Predicts the vapor pressure for the temperature given in setTP().
Definition: pengrobinson.hh:98
static Scalar computeFugacityCoeffient(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition: pengrobinson.hh:227