pengrobinson.hh 19.8 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
6
7
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   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
10
11
12
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
Simon Scholz's avatar
Simon Scholz committed
21
 * \ingroup EOS
22
23
24
25
26
 * \brief Implements the Peng-Robinson equation of state for liquids
 *        and gases.
 *
 * See:
 *
27
 * D.-Y. Peng, D.B. Robinson (1976, pp. 59–64) \cite peng1976 <BR>
28
 * https://doi.org/10.1021/i160057a011
29
 *
30
31
 * R. Reid, et al. "The properties of gases and liquids" (1987, pp. 42-44, 82)
 * \cite reid1987
32
33
34
35
 */
#ifndef DUMUX_PENG_ROBINSON_HH
#define DUMUX_PENG_ROBINSON_HH

36
#include <dune/common/math.hh>
37
#include <dumux/material/fluidstates/temperatureoverlay.hh>
38
#include <dumux/material/idealgas.hh>
Bernd Flemisch's avatar
Bernd Flemisch committed
39
#include <dumux/common/exceptions.hh>
40
#include <dumux/common/math.hh>
41
#include <dumux/common/tabulated2dfunction.hh>
42
43
44

#include <iostream>

45
namespace Dumux {
46
47

/*!
48
 * \ingroup EOS
49
50
51
52
53
 * \brief Implements the Peng-Robinson equation of state for liquids
 *        and gases.
 *
 * See:
 *
54
 * D.-Y. Peng, D.B. Robinson (1976, pp. 59–64) \cite peng1976 <BR>
55
 *
56
 * R. Reid, et al. (1987, pp. 42-44, 82) \cite reid1987
57
58
59
60
61
62
63
64
 */
template <class Scalar>
class PengRobinson
{
    PengRobinson()
    { }

public:
65
66
67
68
69
70
71
72
    static void  init(Scalar aMin, Scalar aMax, int na,
                      Scalar bMin, Scalar bMax, int nb)
    {
        // resize the tabulation for the critical points
        criticalTemperature_.resize(aMin, aMax, na, bMin, bMax, nb);
        criticalPressure_.resize(aMin, aMax, na, bMin, bMax, nb);
        criticalMolarVolume_.resize(aMin, aMax, na, bMin, bMax, nb);

73
74
75
        Scalar VmCrit = 18e-6;
        Scalar pCrit = 1e7;
        Scalar TCrit = 273;
76
77
78
79
80
81
82
83
84
85
86
87
88
        for (int i = 0; i < na; ++i) {
            Scalar a = criticalTemperature_.iToX(i);
            for (int j = 0; j < nb; ++j) {
                Scalar b = criticalTemperature_.jToY(j);

                findCriticalPoint_(TCrit, pCrit, VmCrit,  a, b);

                criticalTemperature_.setSamplePoint(i, j, TCrit);
                criticalPressure_.setSamplePoint(i, j, pCrit);
                criticalMolarVolume_.setSamplePoint(i, j, VmCrit);
            }
        }
    };
Andreas Lauser's avatar
Andreas Lauser committed
89

90
    /*!
91
     * \brief Predicts the vapor pressure \f$\mathrm{[Pa]}\f$ for the temperature given in
92
     *        setTP().
93
94
     * \param T temperature in \f$\mathrm{[K]}\f$
     * \param params Parameters
95
     *
96
97
98
99
100
     * Initially, the vapor pressure is roughly estimated by using the
     * Ambrose-Walton method, then the Newton method is used to make
     * difference between the gas and liquid phase fugacity zero.
     */
    template <class Params>
101
    static Scalar computeVaporPressure(const Params &params, Scalar T)
102
    {
103
        using Component = typename Params::Component;
104
105
106
107
108
109
        if (T >= Component::criticalTemperature())
            return Component::criticalPressure();

        // initial guess of the vapor pressure
        Scalar Vm[3];
        const Scalar eps = Component::criticalPressure()*1e-10;
110

111
112
113
114
115
116
117
        // use the Ambrose-Walton method to get an initial guess of
        // the vapor pressure
        Scalar pVap = ambroseWalton_(params, T);

        // Newton-Raphson method
        for (int i = 0; i < 5; ++i) {
            // calculate the molar densities
Thomas Fetzer's avatar
Thomas Fetzer committed
118
            assert(molarVolumes(Vm, params, T, pVap) == 3);
119
120

            Scalar f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]);
121
            Scalar df_dp =
122
123
124
125
                fugacityDifference_(params, T, pVap  + eps, Vm[0], Vm[2])
                -
                fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]);
            df_dp /= 2*eps;
126

127
128
            Scalar delta = f/df_dp;
            pVap = pVap - delta;
129
130
            using std::abs;
            if (abs(delta/pVap) < 1e-10)
131
132
                break;
        }
133

134
135
136
137
        return pVap;
    }

    /*!
138
     * \brief Computes molar volumes \f$\mathrm{[m^3 / mol]}\f$ where the Peng-Robinson EOS is
139
     *        true.
140
141
142
143
     * \param fs Thermodynamic state of the fluids
     * \param params Parameters
     * \param phaseIdx The phase index
     * \param isGasPhase Specifies the phase state
144
145
146
147
148
149
150
151
152
153
     * \param handleUnphysicalPhase Special handling of the case when the EOS has only one
              intersection with the pressure, and the intersection does not correspond to
              the given phase (the phase is thus considered unphysical). If it happens in
              the case of critical fluid, the critical molar volume is returned for the
              unphysical phase. If the fluid is not critical, a proper extremum of the
              EOS is returned for the unphysical phase. If the parameter is false and the
              EOS has only one intersection with the pressure, the molar volume is computed
              from that single intersection, not depending of the given phase (gas or
              fluid). If the EOS has three intersections with the pressure, this parameter
              is ignored.
154
     */
155
156
    template <class FluidState, class Params>
    static Scalar computeMolarVolume(const FluidState &fs,
157
                                     Params &params,
158
                                     int phaseIdx,
159
                                     bool isGasPhase,
160
                                     bool handleUnphysicalPhase = true)
161
    {
162
        Scalar Vm = 0;
163

164
165
        Scalar T = fs.temperature(phaseIdx);
        Scalar p = fs.pressure(phaseIdx);
166
        
167
168
169
        Scalar a = params.a(phaseIdx); // "attractive factor"
        Scalar b = params.b(phaseIdx); // "co-volume"

170
        Scalar RT = Constants<Scalar>::R*T;
171
172
        Scalar Astar = a*p/(RT*RT);
        Scalar Bstar = b*p/RT;
173
174
175
176
177

        Scalar a1 = 1.0;
        Scalar a2 = - (1 - Bstar);
        Scalar a3 = Astar - Bstar*(3*Bstar + 2);
        Scalar a4 = Bstar*(- Astar + Bstar*(1 + Bstar));
178

179
        Scalar Z[3] = {}; // zero-initializer
180
181
        int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4);

182
183
184
185
        // ignore the first two results if the smallest
        // compressibility factor is <= 0.0. (this means that if we
        // would get negative molar volumes for the liquid phase, we
        // consider the liquid phase non-existant.)
Timo Koch's avatar
Timo Koch committed
186
        // Note that invertCubicPolynomial sorts the roots in ascending order
Timo Koch's avatar
Timo Koch committed
187
        if (numSol > 1 && Z[0] <= 0.0) {
188
189
190
191
192
193
194
            Z[0] = Z[numSol - 1];
            numSol = 1;
        }

        if (numSol > 0 && Z[0] <= 0)
            DUNE_THROW(NumericalProblem, "No positive compressibility factors found");
        
195
196
197
198
        if (numSol == 3) {
            // the EOS has three intersections with the pressure,
            // i.e. the molar volume of gas is the largest one and the
            // molar volume of liquid is the smallest one
199
            if (isGasPhase)
200
201
202
203
204
                Vm = Z[2]*RT/p;
            else
                Vm = Z[0]*RT/p;
        }
        else if (numSol == 1) {
205
            // the EOS has only one intersection with the pressure
206
            
207
            Vm = Z[0]*RT/p;
208

209
            // Handle single root case specially unless told otherwise
210
            if (handleUnphysicalPhase)
211
                Vm = handleSingleRoot_(Vm, fs, params, phaseIdx, isGasPhase);
212
        }
213
214
        else
            DUNE_THROW(NumericalProblem, "Unexpected number of roots in the equation for compressibility factor: " << numSol);
Andreas Lauser's avatar
Andreas Lauser committed
215

216
        using std::isfinite;
217
218
219
220
        
        if (!isfinite(Vm) || Vm <= 0)
            DUNE_THROW(NumericalProblem, "Bad molar volume");

221
222
223
        return Vm;
    }

224
225
226
227
228
229
230
231
232
233
234
235
236
237
    /*!
     * \brief Returns the critical temperature for a given mix
     *
     * \param params EOS (equation of state) parameters of a single-component fluid \
     *        (usually PengRobinsonParms) or a mixture (usually PengRobinsonMixtureParams)
     */
    template <class Params>
    static Scalar criticalTemperature(const Params &params)
    {
        // Here, a() is the attractive force parameter and b()
        // is the repulsive force parameter of the EOS
        return criticalTemperature_(params.a(), params.b());
    }

238
    /*!
239
     * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ for a given pressure
240
241
242
243
244
245
     *        and molar volume.
     *
     * This is the same value as computeFugacity() because the mole
     * fraction of a component in a pure fluid is obviously always
     * 100%.
     *
246
     * \param params Parameters
247
248
     */
    template <class Params>
249
    static Scalar computeFugacityCoeffient(const Params &params)
250
251
252
253
254
    {
        Scalar T = params.temperature();
        Scalar p = params.pressure();
        Scalar Vm = params.molarVolume();

255
        Scalar RT = Constants<Scalar>::R*T;
256
257
        Scalar Z = p*Vm/RT;
        Scalar Bstar = p*params.b() / RT;
258
        using std::sqrt;
259
        Scalar tmp =
260
261
262
263
264
            (Vm + params.b()*(1 + sqrt(2))) /
            (Vm + params.b()*(1 - sqrt(2)));
        Scalar expo = - params.a()/(RT * 2 * params.b() * sqrt(2));
        using std::exp;
        using std::pow;
265
        Scalar fugCoeff =
266
267
            exp(Z - 1) / (Z - Bstar) *
            pow(tmp, expo);
268
269

        return fugCoeff;
270
    }
271
272

    /*!
273
     * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ for a given pressure
274
275
276
277
278
279
     *        and molar volume.
     *
     * This is the fugacity coefficient times the pressure. The mole
     * fraction of a component in a pure fluid is obviously always
     * 100%, so it is not required.
     *
280
     * \param params Parameters
281
282
283
284
285
286
     */
    template <class Params>
    static Scalar computeFugacity(const Params &params)
    { return params.pressure()*computeFugacityCoeff(params); }

protected:
287

288
289
290
291
292
293
    /*
     * Handle single-root case in the equation of state. The problem here is
     * that only one phase exists in this case, while we should also figure out
     * something to return for the other phase. For reference, see Andreas
     * Lauser's thesis: http://dx.doi.org/10.18419/opus-516 (page 32).
     */
294
    template <class FluidState, class Params>
295
    static Scalar handleSingleRoot_(Scalar Vm,
296
297
298
299
300
301
302
303
304
                                  const FluidState &fs,
                                  const Params &params,
                                  int phaseIdx,
                                  bool isGasPhase)
    {
        Scalar T = fs.temperature(phaseIdx);
        
        // for the other phase, we take the extremum of the EOS
        // with the largest distance from the intersection.
305
        if (T > criticalTemperature_(params.a(phaseIdx), params.b(phaseIdx))) {
306
307
            // if the EOS does not exhibit any extrema, the fluid
            // is critical...
308
            return handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
309
310
311
312
313
314
        }
        else {
            // find the extrema (if they are present)
            Scalar Vmin, Vmax, pmin, pmax;
            using std::min;
            using std::max;
315
316
317
            if (findExtrema_(Vmin, Vmax, pmin, pmax,
                             params.a(phaseIdx), params.b(phaseIdx), T))
                return isGasPhase ? max(Vmax, Vm) : min(Vmin, Vm);
318
        }
319
        return Vm;
320
321
    }
    
322
    template <class FluidState, class Params>
323
    static Scalar handleCriticalFluid_(Scalar Vm,
324
                                     const FluidState &fs,
325
                                     const Params &params,
326
                                     int phaseIdx,
327
                                     bool isGasPhase)
328
    {
329
330
        Scalar Vcrit = criticalMolarVolume_(params.a(phaseIdx), params.b(phaseIdx));

331
332
        using std::max;
        using std::min;
333
        return isGasPhase ? max(Vm, Vcrit) : min(Vm, Vcrit);
334
335
    }

336
337
338
339
340
    static void findCriticalPoint_(Scalar &Tcrit,
                                   Scalar &pcrit,
                                   Scalar &Vcrit,
                                   Scalar a,
                                   Scalar b)
Andreas Lauser's avatar
Andreas Lauser committed
341
    {
342
343
        Scalar minVm;
        Scalar maxVm;
344

345
346
        Scalar minP(0);
        Scalar maxP(1e100);
347
348
349
350
351
352
353

        // first, we need to find an isotherm where the EOS exhibits
        // a maximum and a minimum
        Scalar Tmin = 250; // [K]
        for (int i = 0; i < 30; ++i) {
            bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, Tmin);
            if (hasExtrema)
354
                break;
355
            Tmin /= 2;
356
        }
357

358
        Scalar T = Tmin;
359
360
361
362
363
364
365
366

        // Newton's method: Start at minimum temperature and minimize
        // the "gap" between the extrema of the EOS
        for (int i = 0; i < 25; ++i) {
            // calculate function we would like to minimize
            Scalar f = maxVm - minVm;

            // check if we're converged
367
368
369
            if (f < 1e-10) {
                Tcrit = T;
                pcrit = (minP + maxP)/2;
370
371
372
373
374
375
376
377
378
                Vcrit = (maxVm + minVm)/2;
                return;
            }

            // backward differences. Forward differences are not
            // robust, since the isotherm could be critical if an
            // epsilon was added to the temperature. (this is case
            // rarely happens, though)
            const Scalar eps = - 1e-8;
379
            [[maybe_unused]] bool hasExtrema = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
380
            assert(hasExtrema);
381
            Scalar fStar = maxVm - minVm;
Andreas Lauser's avatar
Andreas Lauser committed
382

383
384
385
386
            // derivative of the difference between the maximum's
            // molar volume and the minimum's molar volume regarding
            // temperature
            Scalar fPrime = (fStar - f)/eps;
387

388
389
390
391
392
393
394
395
            // update value for the current iteration
            Scalar delta = f/fPrime;
            if (delta > 0)
                delta = -10;

            // line search (we have to make sure that both extrema
            // still exist after the update)
            for (int j = 0; ; ++j) {
396
                if (j >= 20) {
397
                    DUNE_THROW(NumericalProblem,
398
                               "Could not determine the critical point for a=" << a << ", b=" << b);
399
                }
400

401
402
403
404
                if (findExtrema_(minVm, maxVm, minP, maxP, a, b, T - delta)) {
                    // if the isotherm for T - delta exhibits two
                    // extrema the update is finished
                    T -= delta;
405
406
                    break;
                }
407
                else
408
                    delta /= 2;
409
            }
410
        }
411
    }
412
413
414

    // find the two molar volumes where the EOS exhibits extrema and
    // which are larger than the covolume of the phase
415
    static bool findExtrema_(Scalar &Vmin,
416
                             Scalar &Vmax,
417
418
419
420
421
                             Scalar &pMin,
                             Scalar &pMax,
                             Scalar a,
                             Scalar b,
                             Scalar T)
422
423
424
425
    {
        Scalar u = 2;
        Scalar w = -1;

426
        Scalar RT = Constants<Scalar>::R*T;
427

428
429
430
431
432
433
434
        // calculate coefficients of the 4th order polynominal in
        // monomial basis
        Scalar a1 = RT;
        Scalar a2 = 2*RT*u*b - 2*a;
        Scalar a3 = 2*RT*w*b*b + RT*u*u*b*b  + 4*a*b - u*a*b;
        Scalar a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
        Scalar a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
435

436
437
438
439
440
        // Newton method to find first root

        // if the values which we got on Vmin and Vmax are usefull, we
        // will reuse them as initial value, else we will start 10%
        // above the covolume
441
        Scalar V = b*1.1;
442
        Scalar delta = 1.0;
443
444
        using std::abs;
        for (int i = 0; abs(delta) > 1e-9; ++i) {
445
446
            Scalar f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
            Scalar fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
447

448
            if (abs(fPrime) < 1e-20) {
449
450
451
452
453
                // give up if the derivative is zero
                Vmin = 0;
                Vmax = 0;
                return false;
            }
454

455
456
457
458
459
460
461
462
463
464
465

            delta = f/fPrime;
            V -= delta;

            if (i > 20) {
                // give up after 20 iterations...
                Vmin = 0;
                Vmax = 0;
                return false;
            }
        }
466

467
468
469
470
471
        // polynomial division
        Scalar b1 = a1;
        Scalar b2 = a2 + V*b1;
        Scalar b3 = a3 + V*b2;
        Scalar b4 = a4 + V*b3;
472
473

        // invert resulting cubic polynomial analytically
474
475
        Scalar allV[4];
        allV[0] = V;
476
        int numSol = 1 + invertCubicPolynomial(allV + 1, b1, b2, b3, b4);
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501

        // sort all roots of the derivative
        std::sort(allV + 0, allV + numSol);

        // check whether the result is physical
        if (allV[numSol - 2] < b) {
            // the second largest extremum is smaller than the phase's
            // covolume which is physically impossible
            Vmin = Vmax = 0;
            return false;
        }

        // it seems that everything is okay...
        Vmin = allV[numSol - 2];
        Vmax = allV[numSol - 1];
        return true;
    }

    /*!
     * \brief The Ambrose-Walton method to estimate the vapor
     *        pressure.
     *
     * \return Vapor pressure estimate in bar
     *
     * See:
502
     *
503
     * D. Ambrose, J. Walton (1989, pp. 1395-1403) \cite ambrose1989
504
505
506
507
     */
    template <class Params>
    static Scalar ambroseWalton_(const Params &params, Scalar T)
    {
508
        using Component = typename Params::Component;
509

510
511
512
        Scalar Tr = T / Component::criticalTemperature();
        Scalar tau = 1 - Tr;
        Scalar omega = Component::acentricFactor();
513
        using std::sqrt;
514
515
516
517
        using Dune::power;
        Scalar f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*power(tau, 5))/Tr;
        Scalar f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*power(tau, 5))/Tr;
        Scalar f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*power(tau, 5))/Tr;
518
519
        using std::exp;
        return Component::criticalPressure()*exp(f0 + omega * (f1 + omega*f2));
520
    }
521
522
523
524
525

    /*!
     * \brief Returns the difference between the liquid and the gas phase
     *        fugacities in [bar]
     *
526
     * \param params Parameters
527
     * \param T Temperature [K]
528
529
530
531
532
533
534
535
536
537
     * \param p Pressure [bar]
     * \param VmLiquid Molar volume of the liquid phase [cm^3/mol]
     * \param VmGas Molar volume of the gas phase [cm^3/mol]
     */
    template <class Params>
    static Scalar fugacityDifference_(const Params &params,
                                      Scalar T,
                                      Scalar p,
                                      Scalar VmLiquid,
                                      Scalar VmGas)
538
    { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }
539
540
541
542

    static Tabulated2DFunction<Scalar> criticalTemperature_;
    static Tabulated2DFunction<Scalar> criticalPressure_;
    static Tabulated2DFunction<Scalar> criticalMolarVolume_;
543
544
};

545
546
547
548
549
550
551
552
553
template <class Scalar>
Tabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalTemperature_;

template <class Scalar>
Tabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalPressure_;

template <class Scalar>
Tabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalMolarVolume_;

554
} // end namespace
555
556

#endif