PengRobinson.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef OPM_PENG_ROBINSON_HPP
#define OPM_PENG_ROBINSON_HPP

#include <opm/common/Exceptions.hpp>

#include <opm/material/fluidstates/TemperatureOverlayFluidState.hpp>
#include <opm/material/IdealGas.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>

#include <opm/material/common/PolynomialUtils.hpp>

#include <cassert>
#include <csignal>
#include <string>

namespace Opm {

template <class Scalar, bool UseLegacy=true>
class PengRobinson
{
    static const Scalar R;

    PengRobinson()
    { }

public:
    static void  init(Scalar /*aMin*/, Scalar /*aMax*/, unsigned /*na*/,
                      Scalar /*bMin*/, Scalar /*bMax*/, unsigned /*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);

        Scalar VmCrit, pCrit, TCrit;
        for (unsigned i = 0; i < na; ++i) {
            Scalar a = criticalTemperature_.iToX(i);
            assert(std::abs(criticalTemperature_.xToI(criticalTemperature_.iToX(i)) - i) < 1e-10);

            for (unsigned j = 0; j < nb; ++j) {
                Scalar b = criticalTemperature_.jToY(j);
                assert(std::abs(criticalTemperature_.yToJ(criticalTemperature_.jToY(j)) - j) < 1e-10);

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

                criticalTemperature_.setSamplePoint(i, j, TCrit);
                criticalPressure_.setSamplePoint(i, j, pCrit);
                criticalMolarVolume_.setSamplePoint(i, j, VmCrit);
            }
        }
        */
    }

    template <class Evaluation, class Params>
    static Evaluation computeVaporPressure(const Params& params, const Evaluation& T)
    {
        typedef typename Params::Component Component;
        if (T >= Component::criticalTemperature())
            return Component::criticalPressure();

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

        // use the Ambrose-Walton method to get an initial guess of
        // the vapor pressure
        Evaluation pVap = ambroseWalton_(params, T);

        // Newton-Raphson method
        for (unsigned i = 0; i < 5; ++i) {
            // calculate the molar densities
            assert(molarVolumes(Vm, params, T, pVap) == 3);

            const Evaluation& f = fugacityDifference_(params, T, pVap, Vm[0], Vm[2]);
            Evaluation df_dp =
                fugacityDifference_(params, T, pVap  + eps, Vm[0], Vm[2])
                -
                fugacityDifference_(params, T, pVap - eps, Vm[0], Vm[2]);
            df_dp /= 2*eps;

            const Evaluation& delta = f/df_dp;
            pVap = pVap - delta;

            if (std::abs(scalarValue(delta/pVap)) < 1e-10)
                break;
        }

        return pVap;
    }

    template <class FluidState, class Params>
    static
    typename FluidState::ValueType
    computeMolarVolume(const FluidState& fs,
                       Params& params,
                       unsigned phaseIdx,
                       bool isGasPhase)
    {
        Valgrind::CheckDefined(fs.temperature(phaseIdx));
        Valgrind::CheckDefined(fs.pressure(phaseIdx));

        typedef typename FluidState::ValueType Evaluation;

        Evaluation Vm = 0;
        Valgrind::SetUndefined(Vm);

        const Evaluation& T = fs.temperature(phaseIdx);
        const Evaluation& p = fs.pressure(phaseIdx);

        const Evaluation& a = params.a(phaseIdx); // "attractive factor"
        const Evaluation& b = params.b(phaseIdx); // "co-volume"

        if (!std::isfinite(scalarValue(a))
            || std::abs(scalarValue(a)) < 1e-30)
            return std::numeric_limits<Scalar>::quiet_NaN();
        if (!std::isfinite(scalarValue(b)) || b <= 0)
            return std::numeric_limits<Scalar>::quiet_NaN();

        const Evaluation& RT= Constants<Scalar>::R*T;
        const Evaluation& Astar = a*p/(RT*RT);
        const Evaluation& Bstar = b*p/RT;

        const Evaluation& a1 = 1.0;
        const Evaluation& a2 = - (1 - Bstar);
        const Evaluation& a3 = Astar - Bstar*(3*Bstar + 2);
        const Evaluation& a4 = Bstar*(- Astar + Bstar*(1 + Bstar));

        // 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.)
        Evaluation Z[3] = {0.0,0.0,0.0};
        Valgrind::CheckDefined(a1);
        Valgrind::CheckDefined(a2);
        Valgrind::CheckDefined(a3);
        Valgrind::CheckDefined(a4);

        int numSol = cubicRoots(Z, a1, a2, a3, a4);
        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
            if (isGasPhase)
                Vm = max(1e-7, Z[2]*RT/p);
            else
                Vm = max(1e-7, Z[0]*RT/p);
        }
        else if (numSol == 1) {
            // the EOS only has one intersection with the pressure,
            // for the other phase, we take the extremum of the EOS
            // with the largest distance from the intersection.
            Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
            Vm = VmCubic;

            if (UseLegacy) {
                // find the extrema (if they are present)
                Evaluation Vmin, Vmax, pmin, pmax;
                if (findExtrema_(Vmin, Vmax, pmin, pmax, a, b, T)) {
                    if (isGasPhase)
                        Vm = std::max(Vmax, VmCubic);
                    else {
                        if (Vmin > 0)
                            Vm = std::min(Vmin, VmCubic);
                        else
                            Vm = VmCubic;
                    }
                } else {
                    // the EOS does not exhibit any physically meaningful
                    // extrema, and the fluid is critical...
                    Vm = VmCubic;
                    handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
                }
            }
        }

        Valgrind::CheckDefined(Vm);
        assert(std::isfinite(scalarValue(Vm)));
        assert(Vm > 0);
        return Vm;
    }

    template <class Evaluation, class Params>
    static Evaluation computeFugacityCoeffient(const Params& params)
    {
        const Evaluation& T = params.temperature();
        const Evaluation& p = params.pressure();
        const Evaluation& Vm = params.molarVolume();

        const Evaluation& RT = Constants<Scalar>::R*T;
        const Evaluation& Z = p*Vm/RT;
        const Evaluation& Bstar = p*params.b() / RT;

        const Evaluation& tmp =
            (Vm + params.b()*(1 + std::sqrt(2))) /
            (Vm + params.b()*(1 - std::sqrt(2)));
        const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
        const Evaluation& fugCoeff =
            exp(Z - 1) / (Z - Bstar) *
            pow(tmp, expo);

        return fugCoeff;
    }

    template <class Evaluation, class Params>
    static Evaluation computeFugacity(const Params& params)
    { return params.pressure()*computeFugacityCoeff(params); }

protected:
    template <class FluidState, class Params, class Evaluation = typename FluidState::ValueType>
    static void handleCriticalFluid_(Evaluation& Vm,
                                     const FluidState& /*fs*/,
                                     const Params& params,
                                     unsigned phaseIdx,
                                     bool isGasPhase)
    {
        Evaluation Tcrit, pcrit, Vcrit;
        findCriticalPoint_(Tcrit,
                           pcrit,
                           Vcrit,
                           Evaluation(params.a(phaseIdx)),
                           Evaluation(params.b(phaseIdx)) );


        //Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx));

        if (isGasPhase)
            Vm = max(Vm, Vcrit);
        else
            Vm = min(Vm, Vcrit);
    }

    template <class Evaluation>
    static void findCriticalPoint_(Evaluation& Tcrit,
                                   Evaluation& pcrit,
                                   Evaluation& Vcrit,
                                   const Evaluation& a,
                                   const Evaluation& b)
    {
        Evaluation minVm(0);
        Evaluation maxVm(1e30);

        Evaluation minP(0);
        Evaluation maxP(1e30);

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

        Evaluation T = Tmin;

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

            // check if we're converged
            if (f < 1e-10 || (i == iMax - 1 && f < 1e-8)) {
                Tcrit = T;
                pcrit = (minP + maxP)/2;
                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-11;
            [[maybe_unused]] bool found = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
            assert(found);
            assert(std::isfinite(scalarValue(maxVm)));
            Evaluation fStar = maxVm - minVm;

            // derivative of the difference between the maximum's
            // molar volume and the minimum's molar volume regarding
            // temperature
            Evaluation fPrime = (fStar - f)/eps;
            if (std::abs(scalarValue(fPrime)) < 1e-40) {
                Tcrit = T;
                pcrit = (minP + maxP)/2;
                Vcrit = (maxVm + minVm)/2;
                return;
            }

            // update value for the current iteration
            Evaluation delta = f/fPrime;
            assert(std::isfinite(scalarValue(delta)));
            if (delta > 0)
                delta = -10;

            // line search (we have to make sure that both extrema
            // still exist after the update)
            for (unsigned j = 0; ; ++j) {
                if (j >= 20) {
                    if (f < 1e-8) {
                        Tcrit = T;
                        pcrit = (minP + maxP)/2;
                        Vcrit = (maxVm + minVm)/2;
                        return;
                    }

                    const std::string msg =
                        "Could not determine the critical point for a="
                            + std::to_string(getValue(a))
                        + ", b=" + std::to_string(getValue(b));
                    throw NumericalProblem(msg);
                }

                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;
                    break;
                }
                else
                    delta /= 2;
            };
        }

        assert(false);
    }

    // find the two molar volumes where the EOS exhibits extrema and
    // which are larger than the covolume of the phase
    template <class Evaluation>
    static bool findExtrema_(Evaluation& Vmin,
                             Evaluation& Vmax,
                             Evaluation& /*pMin*/,
                             Evaluation& /*pMax*/,
                             const Evaluation& a,
                             const Evaluation& b,
                             const Evaluation& T)
    {
        Scalar u = 2;
        Scalar w = -1;

        const Evaluation& RT = Constants<Scalar>::R*T;
        // calculate coefficients of the 4th order polynominal in
        // monomial basis
        const Evaluation& a1 = RT;
        const Evaluation& a2 = 2*RT*u*b - 2*a;
        const Evaluation& a3 = 2*RT*w*b*b + RT*u*u*b*b  + 4*a*b - u*a*b;
        const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
        const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;

        assert(std::isfinite(scalarValue(a1)));
        assert(std::isfinite(scalarValue(a2)));
        assert(std::isfinite(scalarValue(a3)));
        assert(std::isfinite(scalarValue(a4)));
        assert(std::isfinite(scalarValue(a5)));

        // 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
        Evaluation V = b*1.1;
        Evaluation delta = 1.0;
        for (unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
            const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
            const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));

            if (std::abs(scalarValue(fPrime)) < 1e-20) {
                // give up if the derivative is zero
                return false;
            }


            delta = f/fPrime;
            V -= delta;

            if (i > 200) {
                // give up after 200 iterations...
                return false;
            }
        }
        assert(std::isfinite(scalarValue(V)));

        // polynomial division
        Evaluation b1 = a1;
        Evaluation b2 = a2 + V*b1;
        Evaluation b3 = a3 + V*b2;
        Evaluation b4 = a4 + V*b3;

        // invert resulting cubic polynomial analytically
        // a vector is used instead of an array due to faulty
        // compiler diagnostics with gcc 14.2
        std::vector<Evaluation> allV(4);
        allV[0] = V;
        const decltype(allV.size()) numSol = 1 + invertCubicPolynomial<Evaluation>(allV.data() + 1, b1, b2, b3, b4);

        assert (numSol <= allV.size());

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

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


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

    template <class Evaluation, class Params>
    static Evaluation ambroseWalton_(const Params& /*params*/, const Evaluation& T)
    {
        typedef typename Params::Component Component;

        const Evaluation& Tr = T / Component::criticalTemperature();
        const Evaluation& tau = 1 - Tr;
        const Evaluation& omega = Component::acentricFactor();

        const Evaluation& f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
        const Evaluation& f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
        const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;

        return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
    }

    template <class Evaluation, class Params>
    static Evaluation fugacityDifference_(const Params& params,
                                          const Evaluation& T,
                                          const Evaluation& p,
                                          const Evaluation& VmLiquid,
                                          const Evaluation& VmGas)
    { return fugacity(params, T, p, VmLiquid) - fugacity(params, T, p, VmGas); }


};

} // namespace Opm

#endif