Merge pull request #521 from akva2/component_fixes

Fixes in some component classes
This commit is contained in:
Markus Blatt
2022-08-02 14:21:11 +02:00
committed by GitHub
3 changed files with 42 additions and 49 deletions

View File

@@ -80,7 +80,7 @@ public:
static Scalar molarMass()
{
const Scalar M1 = H2O::molarMass();
const Scalar M2 = 58e-3; // molar mass of NaCl [kg/mol]
constexpr Scalar M2 = 58e-3; // molar mass of NaCl [kg/mol]
const Scalar X2 = salinity; // mass fraction of salt in brine
return M1*M2/(M2 + X2*(M1 - M2));
}
@@ -149,22 +149,22 @@ public:
const Evaluation& pressure)
{
// Numerical coefficents from Palliser and McKibbin
static const Scalar f[] = {
static constexpr Scalar f[] = {
2.63500e-1, 7.48368e-6, 1.44611e-6, -3.80860e-10
};
// Numerical coefficents from Michaelides for the enthalpy of brine
static const Scalar a[4][3] = {
static constexpr Scalar a[4][3] = {
{ -9633.6, -4080.0, +286.49 },
{ +166.58, +68.577, -4.6856 },
{ -0.90963, -0.36524, +0.249667e-1 },
{ +0.17965e-2, +0.71924e-3, -0.4900e-4 }
};
const Evaluation& theta = temperature - 273.15;
const Evaluation theta = temperature - 273.15;
Evaluation S = salinity;
const Evaluation& S_lSAT =
const Evaluation S_lSAT =
f[0]
+ f[1]*theta
+ f[2]*pow(theta, 2)
@@ -174,17 +174,17 @@ public:
if (S > S_lSAT)
S = S_lSAT;
const Evaluation& hw = H2O::liquidEnthalpy(temperature, pressure)/1e3; // [kJ/kg]
const Evaluation hw = H2O::liquidEnthalpy(temperature, pressure)/1e3; // [kJ/kg]
// From Daubert and Danner
const Evaluation& h_NaCl =
const Evaluation h_NaCl =
(3.6710e4*temperature
+ (6.2770e1/2)*temperature*temperature
- (6.6670e-2/3)*temperature*temperature*temperature
+ (2.8000e-5/4)*pow(temperature, 4.0))/58.44e3
- 2.045698e+02; // [kJ/kg]
const Evaluation& m = S/(1-S)/58.44e-3;
const Evaluation m = S/(1-S)/58.44e-3;
Evaluation d_h = 0;
for (int i = 0; i<=3; ++i) {
@@ -193,10 +193,10 @@ public:
}
}
const Evaluation& delta_h = 4.184/(1e3 + (58.44 * m))*d_h;
const Evaluation delta_h = 4.184/(1e3 + (58.44 * m))*d_h;
// Enthalpy of brine
const Evaluation& h_ls = (1-S)*hw + S*h_NaCl + S*delta_h; // [kJ/kg]
const Evaluation h_ls = (1-S)*hw + S*h_NaCl + S*delta_h; // [kJ/kg]
return h_ls*1e3; // convert to [J/kg]
}
@@ -306,7 +306,7 @@ public:
&& std::abs(scalarValue(pressure)*1e-9) < std::abs(scalarValue(deltaP));
++i)
{
const Evaluation& f = liquidDensity(temperature, pressure) - density;
const Evaluation f = liquidDensity(temperature, pressure) - density;
Evaluation df_dp = liquidDensity(temperature, pressure + eps);
df_dp -= liquidDensity(temperature, pressure - eps);

View File

@@ -51,7 +51,7 @@ namespace Opm {
template <class Scalar, class CO2Tables>
class CO2 : public Component<Scalar, CO2<Scalar, CO2Tables> >
{
static const Scalar R;
static constexpr Scalar R = Constants<Scalar>::R;
public:
/*!
@@ -129,9 +129,9 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{
static const Scalar a[4] =
static constexpr Scalar a[4] =
{ -7.0602087, 1.9391218, -1.6463597, -3.2995634 };
static const Scalar t[4] =
static constexpr Scalar t[4] =
{ 1.0, 1.5, 2.0, 4.0 };
// this is on page 1524 of the reference
@@ -176,8 +176,8 @@ public:
const Evaluation& pressure,
bool extrapolate = false)
{
const Evaluation& h = gasEnthalpy(temperature, pressure, extrapolate);
const Evaluation& rho = gasDensity(temperature, pressure, extrapolate);
const Evaluation h = gasEnthalpy(temperature, pressure, extrapolate);
const Evaluation rho = gasDensity(temperature, pressure, extrapolate);
return h - (pressure / rho);
}
@@ -204,31 +204,31 @@ public:
const Evaluation& pressure,
bool extrapolate = false)
{
const Scalar a0 = 0.235156;
const Scalar a1 = -0.491266;
const Scalar a2 = 5.211155e-2;
const Scalar a3 = 5.347906e-2;
const Scalar a4 = -1.537102e-2;
constexpr Scalar a0 = 0.235156;
constexpr Scalar a1 = -0.491266;
constexpr Scalar a2 = 5.211155e-2;
constexpr Scalar a3 = 5.347906e-2;
constexpr Scalar a4 = -1.537102e-2;
const Scalar d11 = 0.4071119e-2;
const Scalar d21 = 0.7198037e-4;
const Scalar d64 = 0.2411697e-16;
const Scalar d81 = 0.2971072e-22;
const Scalar d82 = -0.1627888e-22;
constexpr Scalar d11 = 0.4071119e-2;
constexpr Scalar d21 = 0.7198037e-4;
constexpr Scalar d64 = 0.2411697e-16;
constexpr Scalar d81 = 0.2971072e-22;
constexpr Scalar d82 = -0.1627888e-22;
const Scalar ESP = 251.196;
constexpr Scalar ESP = 251.196;
if(temperature < 275.) // regularization
temperature = 275.0;
Evaluation TStar = temperature/ESP;
// mu0: viscosity in zero-density limit
const Evaluation& logTStar = log(TStar);
const Evaluation logTStar = log(TStar);
Evaluation SigmaStar = exp(a0 + logTStar*(a1 + logTStar*(a2 + logTStar*(a3 + logTStar*a4))));
Evaluation mu0 = 1.00697*sqrt(temperature) / SigmaStar;
const Evaluation& rho = gasDensity(temperature, pressure, extrapolate); // CO2 mass density [kg/m^3]
const Evaluation rho = gasDensity(temperature, pressure, extrapolate); // CO2 mass density [kg/m^3]
// dmu : excess viscosity at elevated density
Evaluation dmu =
@@ -254,22 +254,18 @@ public:
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar eps = 1e-6;
constexpr Scalar eps = 1e-6;
// use central differences here because one-sided methods do
// not come with a performance improvement. (central ones are
// more accurate, though...)
const Evaluation& h1 = gasEnthalpy(temperature - eps, pressure);
const Evaluation& h2 = gasEnthalpy(temperature + eps, pressure);
const Evaluation h1 = gasEnthalpy(temperature - eps, pressure);
const Evaluation h2 = gasEnthalpy(temperature + eps, pressure);
return (h2 - h1) / (2*eps) ;
}
};
template <class Scalar, class CO2Tables>
const Scalar CO2<Scalar, CO2Tables>::R = Constants<Scalar>::R;
} // namespace Opm
#endif

View File

@@ -64,12 +64,12 @@ namespace Opm {
* \tparam Scalar The type used for representing scalar values
*/
template <class Scalar>
class SimpleHuDuanH2O : public Component<Scalar, SimpleHuDuanH2O<Scalar> >
class SimpleHuDuanH2O : public Component<Scalar, SimpleHuDuanH2O<Scalar>>
{
typedef ::Opm::IdealGas<Scalar> IdealGas;
typedef IAPWS::Common<Scalar> Common;
using IdealGas = ::Opm::IdealGas<Scalar>;
using Common = IAPWS::Common<Scalar>;
static const Scalar R; // specific gas constant of water
static constexpr Scalar R = Constants<Scalar>::R / 18e-3; // specific gas constant of water
public:
/*!
@@ -146,7 +146,7 @@ public:
if (T < tripleTemperature())
return 0; // water is solid: We don't take sublimation into account
static const Scalar n[10] = {
static constexpr Scalar n[10] = {
0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
-0.48232657361591e4, 0.40511340542057e6, -0.23855557567849,
@@ -364,7 +364,7 @@ public:
throw NumericalIssue(oss.str());
}
const Evaluation& rho = liquidDensity(temperature, pressure, extrapolate);
const Evaluation rho = liquidDensity(temperature, pressure, extrapolate);
return Common::viscosity(temperature, rho);
}
@@ -402,10 +402,10 @@ private:
Evaluation p = pressure / 1e6; // to MPa
Scalar Mw = molarMass() * 1e3; //kg/kmol
static const Scalar k0[5] = { 3.27225e-07, -4.20950e-04, 2.32594e-01, -4.16920e+01, 5.71292e+03 };
static const Scalar k1[5] = { -2.32306e-10, 2.91138e-07, -1.49662e-04, 3.59860e-02, -3.55071 };
static const Scalar k2[3] = { 2.57241e-14, -1.24336e-11, 5.42707e-07 };
static const Scalar k3[3] = { -4.42028e-18, 2.10007e-15, -8.11491e-11 };
static constexpr Scalar k0[5] = { 3.27225e-07, -4.20950e-04, 2.32594e-01, -4.16920e+01, 5.71292e+03 };
static constexpr Scalar k1[5] = { -2.32306e-10, 2.91138e-07, -1.49662e-04, 3.59860e-02, -3.55071 };
static constexpr Scalar k2[3] = { 2.57241e-14, -1.24336e-11, 5.42707e-07 };
static constexpr Scalar k3[3] = { -4.42028e-18, 2.10007e-15, -8.11491e-11 };
Evaluation k0_eval = 1e-3 * (((k0[0]*T + k0[1])*T + k0[2])*T + k0[3] + k0[4]/T);
Evaluation k1_eval = 1e-2 * (((k1[0]*T + k1[1])*T + k1[2])*T + k1[3] + k1[4]/T);
Evaluation k2_eval = 1e-1 * ((k2[0]*T + k2[1])*T*T + k2[2]);
@@ -421,9 +421,6 @@ private:
};
template <class Scalar>
const Scalar SimpleHuDuanH2O<Scalar>::R = Constants<Scalar>::R / 18e-3;
} // namespace Opm
#endif