added Svenns cubic solver
This commit is contained in:
Trine Mykkeltvedt
parent
a003b6c35c
commit
16f7fb8e9d
@ -100,7 +100,7 @@ public:
|
||||
* This basically boils down to creating an uninitialized object of sufficient size.
|
||||
* This is method only non-trivial for dynamically-sized Evaluation objects.
|
||||
*/
|
||||
static Scalar createBlank(Scalar)
|
||||
static Scalar createBlank(Scalar value OPM_UNUSED)
|
||||
{ return Scalar(); }
|
||||
|
||||
/*!
|
||||
@ -136,7 +136,7 @@ public:
|
||||
* function. In general, this returns an evaluation object for which all derivatives
|
||||
* are zero.
|
||||
*/
|
||||
static Scalar createConstant(Scalar, Scalar value)
|
||||
static Scalar createConstant(Scalar x OPM_UNUSED, Scalar value)
|
||||
{ return value; }
|
||||
|
||||
/*!
|
||||
@ -146,7 +146,7 @@ public:
|
||||
* regard to x. For scalars (which do not consider derivatives), this method does
|
||||
* nothing.
|
||||
*/
|
||||
static Scalar createVariable(Scalar, unsigned)
|
||||
static Scalar createVariable(Scalar value OPM_UNUSED, unsigned varIdx OPM_UNUSED)
|
||||
{ throw std::logic_error("Plain floating point objects cannot represent variables"); }
|
||||
|
||||
/*!
|
||||
@ -157,7 +157,7 @@ public:
|
||||
* regard to x. For scalars (which do not consider derivatives), this method does
|
||||
* nothing.
|
||||
*/
|
||||
static Scalar createVariable(Scalar, Scalar, unsigned)
|
||||
static Scalar createVariable(Scalar x OPM_UNUSED, Scalar value OPM_UNUSED, unsigned varIdx OPM_UNUSED)
|
||||
{ throw std::logic_error("Plain floating point objects cannot represent variables"); }
|
||||
|
||||
/*!
|
||||
@ -227,6 +227,14 @@ public:
|
||||
static Scalar asin(Scalar arg)
|
||||
{ return std::asin(arg); }
|
||||
|
||||
//! The sine hyperbolicus of a value
|
||||
static Scalar sinh(Scalar arg)
|
||||
{ return std::sinh(arg); }
|
||||
|
||||
//! The arcus sine hyperbolicus of a value
|
||||
static Scalar asinh(Scalar arg)
|
||||
{ return std::asinh(arg); }
|
||||
|
||||
//! The cosine of a value
|
||||
static Scalar cos(Scalar arg)
|
||||
{ return std::cos(arg); }
|
||||
@ -235,6 +243,14 @@ public:
|
||||
static Scalar acos(Scalar arg)
|
||||
{ return std::acos(arg); }
|
||||
|
||||
//! The cosine hyperbolicus of a value
|
||||
static Scalar cosh(Scalar arg)
|
||||
{ return std::cosh(arg); }
|
||||
|
||||
//! The arcus cosine hyperbolicus of a value
|
||||
static Scalar acosh(Scalar arg)
|
||||
{ return std::acosh(arg); }
|
||||
|
||||
//! The square root of a value
|
||||
static Scalar sqrt(Scalar arg)
|
||||
{ return std::sqrt(arg); }
|
||||
@ -357,6 +373,14 @@ template <class Evaluation>
|
||||
Evaluation asin(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::asin(value); }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation sinh(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::sinh(value); }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation asinh(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::asinh(value); }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation cos(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::cos(value); }
|
||||
@ -365,6 +389,14 @@ template <class Evaluation>
|
||||
Evaluation acos(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::acos(value); }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation cosh(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::cosh(value); }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation acosh(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::acosh(value); }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation sqrt(const Evaluation& value)
|
||||
{ return MathToolbox<Evaluation>::sqrt(value); }
|
||||
|
||||
@ -309,6 +309,106 @@ unsigned invertCubicPolynomial(SolContainer* sol,
|
||||
|
||||
return 3;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Math
|
||||
* \brief Invert a cubic polynomial analytically
|
||||
*
|
||||
* The polynomial is defined as
|
||||
* \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f]
|
||||
*
|
||||
* This method teturns the number of solutions which are in the real
|
||||
* numbers. The "sol" argument contains the real roots of the cubic
|
||||
* polynomial in order with the smallest root first.
|
||||
*
|
||||
* \param sol Container into which the solutions are written
|
||||
* \param a The coefficient for the cubic term
|
||||
* \param b The coefficient for the quadratic term
|
||||
* \param c The coefficient for the linear term
|
||||
* \param d The coefficient for the constant term
|
||||
*/
|
||||
template <class Scalar, class SolContainer>
|
||||
unsigned cubicRoots(SolContainer* sol,
|
||||
Scalar a,
|
||||
Scalar b,
|
||||
Scalar c,
|
||||
Scalar d)
|
||||
{
|
||||
// reduces to a quadratic polynomial
|
||||
if (std::abs(scalarValue(a)) < 1e-30)
|
||||
return invertQuadraticPolynomial(sol, b, c, d);
|
||||
|
||||
// We need to reduce the cubic equation to its "depressed cubic" form (however strange that sounds)
|
||||
// Depressed cubic form: t^3 + p*t + q, where x = t - b/3*a is the transform we use when we have
|
||||
// roots for t. p and q are defined below.
|
||||
// Formula for p and q:
|
||||
Scalar p = (3.0 * a * c - b * b) / (3.0 * a * a);
|
||||
Scalar q = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * d * a * a) / (27.0 * a * a * a);
|
||||
|
||||
// Check if we have three or one real root by looking at the discriminant, and solve accordingly with
|
||||
// correct formula
|
||||
Scalar discr = 4.0 * p * p * p + 27.0 * q * q;
|
||||
if (discr < 0.0) {
|
||||
// Find three real roots of a depressed cubic, using the trigonometric method
|
||||
// Help calculation
|
||||
Scalar theta = (1.0 / 3.0) * acos( ((3.0 * q) / (2.0 * p)) * sqrt(-3.0 / p) );
|
||||
|
||||
// Calculate the three roots
|
||||
sol[0] = 2.0 * sqrt(-p / 3.0) * cos( theta );
|
||||
sol[1] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((2.0 * M_PI) / 3.0) );
|
||||
sol[2] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((4.0 * M_PI) / 3.0) );
|
||||
|
||||
// Sort in ascending order
|
||||
std::sort(sol, sol + 3);
|
||||
|
||||
// Return confirmation of three roots
|
||||
// std::cout << "Z (discr < 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl;
|
||||
return 3;
|
||||
}
|
||||
else if (discr > 0.0) {
|
||||
// Find one real root of a depressed cubic using hyperbolic method. Different solutions depending on
|
||||
// sign of p
|
||||
Scalar t;
|
||||
if (p < 0) {
|
||||
// Help calculation
|
||||
Scalar theta = (1.0 / 3.0) * acosh( ((-3.0 * abs(q)) / (2.0 * p)) * sqrt(-3.0 / p) );
|
||||
|
||||
// Root
|
||||
t = ( (-2.0 * abs(q)) / q ) * sqrt(-p / 3.0) * cosh(theta);
|
||||
}
|
||||
else if (p > 0) {
|
||||
// Help calculation
|
||||
Scalar theta = (1.0 / 3.0) * asinh( ((3.0 * q) / (2.0 * p)) * sqrt(3.0 / p) );
|
||||
// Root
|
||||
t = -2.0 * sqrt(p / 3.0) * sinh(theta);
|
||||
|
||||
}
|
||||
else {
|
||||
std::runtime_error(" p = 0 in cubic root solver!");
|
||||
}
|
||||
|
||||
// Transform t to output solution
|
||||
sol[0] = t - b / (3.0 * a);
|
||||
// std::cout << "Z (discr > 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl;
|
||||
return 1;
|
||||
|
||||
}
|
||||
else {
|
||||
// The discriminant, 4*p^3 + 27*q^2 = 0, thus we have simple (real) roots
|
||||
// If p = 0 then also q = 0, and t = 0 is a triple root
|
||||
if (p == 0) {
|
||||
sol[0] = sol[1] = sol[2] = 0.0 - b / (3.0 * a);
|
||||
}
|
||||
// If p != 0, the we have a simple root and a double root
|
||||
else {
|
||||
sol[0] = (3.0 * q / p) - b / (3.0 * a);
|
||||
sol[1] = sol[2] = (-3.0 * q) / (2.0 * p) - b / (3.0 * a);
|
||||
std::sort(sol, sol + 3);
|
||||
}
|
||||
// std::cout << "Z (disc = 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl;
|
||||
return 3;
|
||||
}
|
||||
}
|
||||
} // end Opm
|
||||
|
||||
#endif
|
||||
|
||||
@ -225,6 +225,40 @@ Evaluation<ValueType, numVars, staticSize> asin(const Evaluation<ValueType, numV
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class ValueType, int numVars, unsigned staticSize>
|
||||
Evaluation<ValueType, numVars, staticSize> sinh(const Evaluation<ValueType, numVars, staticSize>& x)
|
||||
{
|
||||
typedef MathToolbox<ValueType> ValueTypeToolbox;
|
||||
|
||||
Evaluation<ValueType, numVars, staticSize> result(x);
|
||||
|
||||
result.setValue(ValueTypeToolbox::sinh(x.value()));
|
||||
|
||||
// derivatives use the chain rule
|
||||
const ValueType& df_dx = ValueTypeToolbox::cosh(x.value());
|
||||
for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx)
|
||||
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class ValueType, int numVars, unsigned staticSize>
|
||||
Evaluation<ValueType, numVars, staticSize> asinh(const Evaluation<ValueType, numVars, staticSize>& x)
|
||||
{
|
||||
typedef MathToolbox<ValueType> ValueTypeToolbox;
|
||||
|
||||
Evaluation<ValueType, numVars, staticSize> result(x);
|
||||
|
||||
result.setValue(ValueTypeToolbox::asinh(x.value()));
|
||||
|
||||
// derivatives use the chain rule
|
||||
const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(x.value()*x.value() + 1);
|
||||
for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx)
|
||||
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class ValueType, int numVars, unsigned staticSize>
|
||||
Evaluation<ValueType, numVars, staticSize> cos(const Evaluation<ValueType, numVars, staticSize>& x)
|
||||
{
|
||||
@ -259,6 +293,40 @@ Evaluation<ValueType, numVars, staticSize> acos(const Evaluation<ValueType, numV
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class ValueType, int numVars, unsigned staticSize>
|
||||
Evaluation<ValueType, numVars, staticSize> cosh(const Evaluation<ValueType, numVars, staticSize>& x)
|
||||
{
|
||||
typedef MathToolbox<ValueType> ValueTypeToolbox;
|
||||
|
||||
Evaluation<ValueType, numVars, staticSize> result(x);
|
||||
|
||||
result.setValue(ValueTypeToolbox::cosh(x.value()));
|
||||
|
||||
// derivatives use the chain rule
|
||||
const ValueType& df_dx = ValueTypeToolbox::sinh(x.value());
|
||||
for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx)
|
||||
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class ValueType, int numVars, unsigned staticSize>
|
||||
Evaluation<ValueType, numVars, staticSize> acosh(const Evaluation<ValueType, numVars, staticSize>& x)
|
||||
{
|
||||
typedef MathToolbox<ValueType> ValueTypeToolbox;
|
||||
|
||||
Evaluation<ValueType, numVars, staticSize> result(x);
|
||||
|
||||
result.setValue(ValueTypeToolbox::acosh(x.value()));
|
||||
|
||||
// derivatives use the chain rule
|
||||
const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(x.value()*x.value() - 1);
|
||||
for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx)
|
||||
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template <class ValueType, int numVars, unsigned staticSize>
|
||||
Evaluation<ValueType, numVars, staticSize> sqrt(const Evaluation<ValueType, numVars, staticSize>& x)
|
||||
{
|
||||
|
||||
@ -55,6 +55,8 @@ namespace Opm {
|
||||
template <class Scalar>
|
||||
class PengRobinson
|
||||
{
|
||||
//! The ideal gas constant [Pa * m^3/mol/K]
|
||||
static const Scalar R;
|
||||
|
||||
PengRobinson()
|
||||
{ }
|
||||
@ -127,7 +129,7 @@ public:
|
||||
const Evaluation& delta = f/df_dp;
|
||||
pVap = pVap - delta;
|
||||
|
||||
if (std::abs(scalarValue(delta/pVap)) < 1e-10)
|
||||
if (std::abs(Opm::scalarValue(delta/pVap)) < 1e-10)
|
||||
break;
|
||||
}
|
||||
|
||||
@ -160,13 +162,13 @@ public:
|
||||
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)
|
||||
if (!std::isfinite(Opm::scalarValue(a))
|
||||
|| std::abs(Opm::scalarValue(a)) < 1e-30)
|
||||
return std::numeric_limits<Scalar>::quiet_NaN();
|
||||
if (!std::isfinite(scalarValue(b)) || b <= 0)
|
||||
if (!std::isfinite(Opm::scalarValue(b)) || b <= 0)
|
||||
return std::numeric_limits<Scalar>::quiet_NaN();
|
||||
|
||||
const Evaluation& RT= Constants<Scalar>::R*T;
|
||||
const Evaluation& RT= R*T;
|
||||
const Evaluation& Astar = a*p/(RT*RT);
|
||||
const Evaluation& Bstar = b*p/RT;
|
||||
|
||||
@ -184,27 +186,26 @@ public:
|
||||
Valgrind::CheckDefined(a2);
|
||||
Valgrind::CheckDefined(a3);
|
||||
Valgrind::CheckDefined(a4);
|
||||
int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4);
|
||||
// std::cout << "Cubic params : " << a1 << " " << a2 << " " << a3 << " " << a4 << std::endl;
|
||||
// int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4);
|
||||
int numSol = cubicRoots(Z, a1, a2, a3, a4);
|
||||
// std::cout << "Z = " << Z[0] << " " << Z[1] << " " << Z[2] << std::endl;
|
||||
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
|
||||
#warning HACK, should investigate why
|
||||
if (isGasPhase)
|
||||
// Vm = Z[2]*RT/p;
|
||||
Vm = max(1e-7, Z[2]*RT/p);
|
||||
Vm = Opm::max(1e-7, Z[2]*RT/p);
|
||||
else
|
||||
// Vm = Z[2]*RT/p;
|
||||
Vm = max(1e-7, Z[0]*RT/p);
|
||||
Vm = Opm::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 = Z[0]*RT/p;
|
||||
Evaluation VmCubic = Opm::max(1e-7, Z[0]*RT/p);
|
||||
Vm = VmCubic;
|
||||
|
||||
#warning, should investigate why here
|
||||
// find the extrema (if they are present)
|
||||
// Evaluation Vmin, Vmax, pmin, pmax;
|
||||
// if (findExtrema_(Vmin, Vmax,
|
||||
@ -229,7 +230,7 @@ public:
|
||||
}
|
||||
|
||||
Valgrind::CheckDefined(Vm);
|
||||
assert(std::isfinite(scalarValue(Vm)));
|
||||
assert(Opm::isfinite(Vm));
|
||||
assert(Vm > 0);
|
||||
return Vm;
|
||||
}
|
||||
@ -251,7 +252,7 @@ public:
|
||||
const Evaluation& p = params.pressure();
|
||||
const Evaluation& Vm = params.molarVolume();
|
||||
|
||||
const Evaluation& RT = Constants<Scalar>::R*T;
|
||||
const Evaluation& RT = R*T;
|
||||
const Evaluation& Z = p*Vm/RT;
|
||||
const Evaluation& Bstar = p*params.b() / RT;
|
||||
|
||||
@ -260,8 +261,8 @@ public:
|
||||
(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);
|
||||
Opm::exp(Z - 1) / (Z - Bstar) *
|
||||
Opm::pow(tmp, expo);
|
||||
|
||||
return fugCoeff;
|
||||
}
|
||||
@ -299,9 +300,9 @@ protected:
|
||||
//Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx));
|
||||
|
||||
if (isGasPhase)
|
||||
Vm = max(Vm, Vcrit);
|
||||
Vm = Opm::max(Vm, Vcrit);
|
||||
else
|
||||
Vm = min(Vm, Vcrit);
|
||||
Vm = Opm::min(Vm, Vcrit);
|
||||
}
|
||||
|
||||
template <class Evaluation>
|
||||
@ -351,14 +352,14 @@ protected:
|
||||
const Scalar eps = - 1e-11;
|
||||
bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
|
||||
assert(hasExtrema);
|
||||
assert(std::isfinite(scalarValue(maxVm)));
|
||||
assert(std::isfinite(Opm::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) {
|
||||
if (std::abs(Opm::scalarValue(fPrime)) < 1e-40) {
|
||||
Tcrit = T;
|
||||
pcrit = (minP + maxP)/2;
|
||||
Vcrit = (maxVm + minVm)/2;
|
||||
@ -367,7 +368,7 @@ protected:
|
||||
|
||||
// update value for the current iteration
|
||||
Evaluation delta = f/fPrime;
|
||||
assert(std::isfinite(scalarValue(delta)));
|
||||
assert(std::isfinite(Opm::scalarValue(delta)));
|
||||
if (delta > 0)
|
||||
delta = -10;
|
||||
|
||||
@ -415,7 +416,7 @@ protected:
|
||||
Scalar u = 2;
|
||||
Scalar w = -1;
|
||||
|
||||
const Evaluation& RT = Constants<Scalar>::R*T;
|
||||
const Evaluation& RT = R*T;
|
||||
|
||||
// calculate coefficients of the 4th order polynominal in
|
||||
// monomial basis
|
||||
@ -425,11 +426,11 @@ protected:
|
||||
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)));
|
||||
assert(std::isfinite(Opm::scalarValue(a1)));
|
||||
assert(std::isfinite(Opm::scalarValue(a2)));
|
||||
assert(std::isfinite(Opm::scalarValue(a3)));
|
||||
assert(std::isfinite(Opm::scalarValue(a4)));
|
||||
assert(std::isfinite(Opm::scalarValue(a5)));
|
||||
|
||||
// Newton method to find first root
|
||||
|
||||
@ -438,11 +439,11 @@ protected:
|
||||
// above the covolume
|
||||
Evaluation V = b*1.1;
|
||||
Evaluation delta = 1.0;
|
||||
for (unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
|
||||
for (unsigned i = 0; std::abs(Opm::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) {
|
||||
if (std::abs(Opm::scalarValue(fPrime)) < 1e-20) {
|
||||
// give up if the derivative is zero
|
||||
return false;
|
||||
}
|
||||
@ -456,7 +457,7 @@ protected:
|
||||
return false;
|
||||
}
|
||||
}
|
||||
assert(std::isfinite(scalarValue(V)));
|
||||
assert(std::isfinite(Opm::scalarValue(V)));
|
||||
|
||||
// polynomial division
|
||||
Evaluation b1 = a1;
|
||||
@ -467,7 +468,7 @@ protected:
|
||||
// invert resulting cubic polynomial analytically
|
||||
Evaluation allV[4];
|
||||
allV[0] = V;
|
||||
int numSol = 1 + invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
|
||||
int numSol = 1 + Opm::invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
|
||||
|
||||
// sort all roots of the derivative
|
||||
std::sort(allV + 0, allV + numSol);
|
||||
@ -507,9 +508,9 @@ protected:
|
||||
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;
|
||||
const Evaluation& f0 = (tau*(-5.97616 + Opm::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Opm::pow(tau, 5))/Tr;
|
||||
const Evaluation& f1 = (tau*(-5.03365 + Opm::sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*Opm::pow(tau, 5))/Tr;
|
||||
const Evaluation& f2 = (tau*(-0.64771 + Opm::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Opm::pow(tau, 5))/Tr;
|
||||
|
||||
return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
|
||||
}
|
||||
@ -539,6 +540,8 @@ protected:
|
||||
*/
|
||||
};
|
||||
|
||||
template <class Scalar>
|
||||
const Scalar PengRobinson<Scalar>::R = Opm::Constants<Scalar>::R;
|
||||
|
||||
/*
|
||||
template <class Scalar>
|
||||
|
||||
Loading…
Reference in New Issue
Block a user