Loading...
Searching...
No Matches
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::Scalar
computeMolarVolume(const FluidState& fs,
Params& params,
unsigned phaseIdx,
bool isGasPhase)
{
Valgrind::CheckDefined(fs.temperature(phaseIdx));
Valgrind::CheckDefined(fs.pressure(phaseIdx));
typedef typename FluidState::Scalar 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::Scalar>
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
std::array<Evaluation,4> allV;
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
Exceptions.hpp
Provides the OPM specific exception classes.
IdealGas.hpp
Relations valid for an ideal gas.
PolynomialUtils.hpp
Provides free functions to invert polynomials of degree 1, 2 and 3.
TemperatureOverlayFluidState.hpp
This is a fluid state which allows to set the fluid temperatures and takes all other quantities from ...
UniformTabulated2DFunction.hpp
Implements a scalar function that depends on two variables and which is sampled on an uniform X-Y gri...
Opm::Component
Abstract base class of a pure chemical species.
Definition Component.hpp:44
Opm::Component::acentricFactor
static Scalar acentricFactor()
Returns the acentric factor of the component.
Definition Component.hpp:112
Opm::Component::criticalPressure
static Scalar criticalPressure()
Returns the critical pressure in of the component.
Definition Component.hpp:105
Opm::Component::criticalTemperature
static Scalar criticalTemperature()
Returns the critical temperature in of the component.
Definition Component.hpp:99
Opm::Constants
A central place for various physical constants occuring in some equations.
Definition Constants.hpp:41
Opm::Constants::R
static const Scalar R
The ideal gas constant [J/(mol K)].
Definition Constants.hpp:45
Opm::PengRobinson
Implements the Peng-Robinson equation of state for liquids and gases.
Definition PengRobinson.hpp:60
Opm::PengRobinson::computeVaporPressure
static Evaluation computeVaporPressure(const Params &params, const Evaluation &T)
Predicts the vapor pressure for the temperature given in setTP().
Definition PengRobinson.hpp:105
Opm::PengRobinson::fugacityDifference_
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:530
Opm::PengRobinson::computeMolarVolume
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:148
Opm::PengRobinson::computeFugacity
static Evaluation computeFugacity(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition PengRobinson.hpp:279
Opm::PengRobinson::computeFugacityCoeffient
static Evaluation computeFugacityCoeffient(const Params &params)
Returns the fugacity coefficient for a given pressure and molar volume.
Definition PengRobinson.hpp:247
Opm::PengRobinson::ambroseWalton_
static Evaluation ambroseWalton_(const Params &, const Evaluation &T)
The Ambrose-Walton method to estimate the vapor pressure.
Definition PengRobinson.hpp:504
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
Opm::cubicRoots
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition PolynomialUtils.hpp:331