OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make all components local-AD aware

Browse Source
This commit is contained in:
Andreas Lauser 2015-05-21 15:33:07 +02:00
parent 35bbae78ab
commit dadefcd3f3
20 changed files with 1215 additions and 917 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
opm/material/IdealGas.hpp
View File

@ -25,9 +25,7 @@
#include <opm/material/Constants.hpp>
namespace Opm
{
namespace Opm {
/*!
* \brief Relations valid for an ideal gas.
*/
@ -42,25 +40,28 @@ public:
* \brief The density of the gas in \f$\mathrm{[kg/m^3]}\f$, depending on
* pressure, temperature and average molar mass of the gas.
*/
static Scalar density(Scalar avgMolarMass,
Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation density(const Evaluation& avgMolarMass,
const Evaluation& temperature,
const Evaluation& pressure)
{ return pressure*avgMolarMass/(R*temperature); }
/*!
* \brief The pressure of the gas in \f$\mathrm{[N/m^2]}\f$, depending on
* the molar density and temperature.
*/
static Scalar pressure(Scalar temperature,
Scalar rhoMolar)
template <class Evaluation>
static Evaluation pressure(const Evaluation& temperature,
const Evaluation& rhoMolar)
{ return R*temperature*rhoMolar; }
/*!
* \brief The molar density of the gas \f$\mathrm{[mol/m^3]}\f$,
* depending on pressure and temperature.
*/
static Scalar molarDensity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation molarDensity(const Evaluation& temperature,
const Evaluation& pressure)
{ return pressure/(R*temperature); }
};

145
opm/material/components/Air.hpp
View File

@ -26,6 +26,7 @@
#define OPM_AIR_HPP
#include <opm/material/components/Component.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/IdealGas.hpp>
#include <opm/material/common/Exceptions.hpp>
@ -46,12 +47,30 @@ class Air : public Component<Scalar, Air<Scalar> >
typedef Opm::IdealGas<Scalar> IdealGas;
public:
/*!
* \brief Returns true iff the liquid phase is assumed to be compressible
*/
static bool liquidIsCompressible()
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidIsCompressible()"); }
/*!
* \brief A human readable name for the \f$Air\f$.
*/
static const char *name()
{ return "Air"; }
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
*/
static bool gasIsCompressible()
{ return true; }
/*!
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static bool gasIsIdeal()
{ return true; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of \f$AIR\f$.
*
@ -77,24 +96,10 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of phase in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), temperature, pressure);
}
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
*/
static bool gasIsCompressible()
{ return true; }
/*!
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static bool gasIsIdeal()
{ return true; }
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return IdealGas::density(Evaluation(molarMass()), temperature, pressure); }
/*!
* \brief The pressure of gaseous \f$AIR\f$ at a given density and temperature \f$\mathrm{[Pa]}\f$.
@ -102,11 +107,10 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
}
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, Scalar density)
{ return IdealGas::pressure(temperature, density/molarMass()); }
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of \f$AIR\f$ at a given pressure and temperature.
*
@ -128,41 +132,41 @@ public:
* V_c = (R*T_c)/p_c
*
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar Tc = criticalTemperature();
const Scalar Vc = 84.525138; // critical specific volume [cm^3/mol]
const Scalar omega = 0.078; // accentric factor
const Scalar M = molarMass() * 1e3; // molar mas [g/mol]
const Scalar dipole = 0.0; // dipole moment [debye]
Scalar Tc = criticalTemperature();
Scalar Vc = 84.525138; // critical specific volume [cm^3/mol]
Scalar omega = 0.078; // accentric factor
Scalar M = molarMass() * 1e3; // molar mas [g/mol]
Scalar dipole = 0.0; // dipole moment [debye]
Scalar mu_r4 = 131.3 * dipole / std::sqrt(Vc * Tc);
mu_r4 *= mu_r4;
mu_r4 *= mu_r4;
Scalar Fc = 1 - 0.2756*omega + 0.059035*mu_r4;
Scalar Tstar = 1.2593 * temperature/Tc;
Scalar Omega_v =
1.16145*std::pow(Tstar, -0.14874) +
0.52487*std::exp(- 0.77320*Tstar) +
2.16178*std::exp(- 2.43787*Tstar);
Scalar mu = 40.785*Fc*std::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;
Evaluation Tstar = 1.2593 * temperature/Tc;
Evaluation Omega_v =
1.16145*Toolbox::pow(Tstar, -0.14874) +
0.52487*Toolbox::exp(- 0.77320*Tstar) +
2.16178*Toolbox::exp(- 2.43787*Tstar);
return 40.7851e-7*Fc*Toolbox::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
}
// simpler method, from old constrelAir.hh
static Scalar simpleGasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation simpleGasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar r;
if(temperature < 273.15 || temperature > 660.)
OPM_THROW(NumericalIssue, "Air: Temperature out of range at ");
r = 1.496*1.E-6*std::pow(temperature,1.5)/(temperature+120.);
return (r);
typedef MathToolbox<Evaluation> Toolbox;
if(temperature < 273.15 || temperature > 660.) {
OPM_THROW(NumericalIssue,
"Air: Temperature (" << temperature << "K) out of range");
}
return 1.496e-6*Toolbox::pow(temperature, 1.5)/(temperature + 120);
}
/*!
@ -176,9 +180,10 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
return 1005*(temperature-273.15);
return 1005*(temperature - 273.15);
}
/*!
@ -192,14 +197,13 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
gasEnthalpy(temperature, pressure)
-
IdealGas::R * temperature // = pressure * molar volume for an ideal gas
/ molarMass(); // conversion from [J/(mol K)] to [J/(kg K)]
- (IdealGas::R*temperature/molarMass()); // <- specific volume of an ideal gas
}
/*!
@ -213,8 +217,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasThermalConductivity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature,
const Evaluation& pressure)
{
// Isobaric Properties for Nitrogen in: NIST Standard
// see http://webbook.nist.gov/chemistry/fluid/
@ -241,26 +246,30 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{
// scale temperature with referenence temp of 100K
Scalar phi = temperature/100;
typedef MathToolbox<Evaluation> Toolbox;
Scalar c_p =
// scale temperature by reference temp of 100K
Evaluation phi = temperature/100;
Evaluation c_p =
0.661738E+01
-0.105885E+01 * phi
+0.201650E+00 * std::pow(phi,2.)
-0.196930E-01 * std::pow(phi,3.)
+0.106460E-02 * std::pow(phi,4.)
-0.303284E-04 * std::pow(phi,5.)
+0.355861E-06 * std::pow(phi,6.);
+0.201650E+00 * Toolbox::pow(phi,2.)
-0.196930E-01 * Toolbox::pow(phi,3.)
+0.106460E-02 * Toolbox::pow(phi,4.)
-0.303284E-04 * Toolbox::pow(phi,5.)
+0.355861E-06 * Toolbox::pow(phi,6.);
c_p +=
-0.549169E+01 * std::pow(phi,-1.)
+0.585171E+01* std::pow(phi,-2.)
-0.372865E+01* std::pow(phi,-3.)
+0.133981E+01* std::pow(phi,-4.)
-0.233758E+00* std::pow(phi,-5.)
+0.125718E-01* std::pow(phi,-6.);
-0.549169E+01 * Toolbox::pow(phi,-1.)
+0.585171E+01* Toolbox::pow(phi,-2.)
-0.372865E+01* Toolbox::pow(phi,-3.)
+0.133981E+01* Toolbox::pow(phi,-4.)
-0.233758E+00* Toolbox::pow(phi,-5.)
+0.125718E-01* Toolbox::pow(phi,-6.);
c_p *= IdealGas::R / (molarMass() * 1000); // in J/mol/K * mol / kg / 1000 = kJ/kg/K
return c_p;

202
opm/material/components/Brine.hpp
View File

@ -25,8 +25,7 @@
#define OPM_BRINE_HPP
#include <opm/material/components/Component.hpp>
#include <cmath>
#include <opm/material/common/MathToolbox.hpp>
namespace Opm {
@ -51,6 +50,24 @@ public:
static const char *name()
{ return "Brine"; }
/*!
* \copydoc H2O::gasIsIdeal
*/
static bool gasIsIdeal()
{ return H2O::gasIsIdeal(); }
/*!
* \copydoc H2O::gasIsCompressible
*/
static bool gasIsCompressible()
{ return H2O::gasIsCompressible(); }
/*!
* \copydoc H2O::liquidIsCompressible
*/
static bool liquidIsCompressible()
{ return H2O::liquidIsCompressible(); }
/*!
* \copydoc Component::molarMass
*
@ -91,14 +108,16 @@ public:
/*!
* \copydoc H2O::vaporPressure
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{ return H2O::vaporPressure(T); /* [N/m^2] */ }
/*!
* \copydoc Component::gasEnthalpy
*/
static const Scalar gasEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{ return H2O::gasEnthalpy(temperature, pressure); /* [J/kg] */ }
/*!
@ -109,71 +128,71 @@ public:
* - Michaelides 1981
* - Daubert & Danner 1989
*/
static const Scalar liquidEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
/*Numerical coefficents from PALLISER*/
typedef MathToolbox<Evaluation> Toolbox;
// Numerical coefficents from Palliser and McKibbin
static const Scalar f[] = {
2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10
2.63500e-1, 7.48368e-6, 1.44611e-6, -3.80860e-10
};
/*Numerical coefficents from MICHAELIDES for the enthalpy of brine*/
// Numerical coefficents from Michaelides for the enthalpy of brine
static const 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 }
{ -0.90963, -0.36524, +0.249667e-1 },
{ +0.17965e-2, +0.71924e-3, -0.4900e-4 }
};
Scalar theta, h_NaCl;
Scalar m, h_ls, h_ls1, d_h;
Scalar S_lSAT, delta_h;
int i, j;
Scalar hw;
theta = temperature - 273.15;
Evaluation theta = temperature - 273.15;
Scalar S = salinity;
S_lSAT = f[0] + f[1]*theta + f[2]*std::pow(theta,2) + f[3]*std::pow(theta,3);
/*Regularization*/
if (S>S_lSAT) {
Scalar S_lSAT =
f[0]
+ f[1]*Toolbox::value(theta)
+ f[2]*std::pow(Toolbox::value(theta), 2)
+ f[3]*std::pow(Toolbox::value(theta), 3);
// Regularization
if (S > S_lSAT)
S = S_lSAT;
}
hw = H2O::liquidEnthalpy(temperature, pressure)/1E3; /* kJ/kg */
Evaluation hw = H2O::liquidEnthalpy(temperature, pressure)/1e3; // [kJ/kg]
/*DAUBERT and DANNER*/
/*U=*/h_NaCl = (3.6710E4*temperature + 0.5*(6.2770E1)*temperature*temperature - ((6.6670E-2)/3)*temperature*temperature*temperature
+((2.8000E-5)/4)*std::pow(temperature,4))/(58.44E3)- 2.045698e+02; /* kJ/kg */
// From Daubert and Danner
Evaluation h_NaCl =
(3.6710e4*temperature
+ (6.2770e1/2)*temperature*temperature
- (6.6670e-2/3)*temperature*temperature*temperature
+ (2.8000e-5/4)*std::pow(temperature,4))/58.44e3
- 2.045698e+02; // [kJ/kg]
m = (1E3/58.44)*(S/(1-S));
i = 0;
j = 0;
d_h = 0;
Scalar m = S/(1-S)/58.44e-3;
for (i = 0; i<=3; i++) {
for (j=0; j<=2; j++) {
d_h = d_h + a[i][j] * std::pow(theta, i) * std::pow(m, j);
Evaluation d_h = 0;
for (int i = 0; i<=3; ++i) {
for (int j = 0; j <= 2; ++j) {
d_h += a[i][j] * std::pow(theta, i) * std::pow(m, j);
}
}
delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
Evaluation delta_h = 4.184/(1e3 + (58.44 * m))*d_h;
/* Enthalpy of brine */
h_ls1 =(1-S)*hw + S*h_NaCl + S*delta_h; /* kJ/kg */
h_ls = h_ls1*1E3; /*J/kg*/
return (h_ls);
// Enthalpy of brine
Evaluation h_ls = (1-S)*hw + S*h_NaCl + S*delta_h; // [kJ/kg]
return h_ls*1e3; // convert to [J/kg]
}
/*!
* \copydoc H2O::liquidHeatCapacity
*/
static const Scalar liquidHeatCapacity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{
Scalar eps = temperature*1e-8;
return (liquidEnthalpy(temperature + eps, pressure) - liquidEnthalpy(temperature, pressure))/eps;
@ -182,15 +201,17 @@ public:
/*!
* \copydoc H2O::gasHeatCapacity
*/
static const Scalar gasHeatCapacity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{ return H2O::gasHeatCapacity(temperature, pressure); }
/*!
* \copydoc H2O::gasInternalEnergy
*/
static const Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
gasEnthalpy(temperature, pressure) -
@ -200,27 +221,22 @@ public:
/*!
* \copydoc H2O::liquidInternalEnergy
*/
static const Scalar liquidInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
liquidEnthalpy(temperature, pressure) -
pressure/liquidDensity(temperature, pressure);
}
/*!
* \copydoc H2O::gasDensity
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return H2O::gasDensity(temperature, pressure); }
/*!
* \copydoc H2O::gasIsIdeal
*/
static bool gasIsIdeal()
{ return H2O::gasIsIdeal(); }
/*!
* \copydoc Component::liquidDensity
*
@ -228,12 +244,13 @@ public:
* - Batzle & Wang (1992)
* - cited by: Adams & Bachu in Geofluids (2002) 2, 257-271
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar TempC = temperature - 273.15;
Scalar pMPa = pressure/1.0E6;
Evaluation tempC = temperature - 273.15;
Evaluation pMPa = pressure/1.0E6;
Scalar rhow = H2O::liquidDensity(temperature, pressure);
const Evaluation rhow = H2O::liquidDensity(temperature, pressure);
return
rhow +
1000*salinity*(
@ -242,49 +259,44 @@ public:
1.0E-6*(
300*pMPa -
2400*pMPa*salinity +
TempC*(
tempC*(
80.0 -
3*TempC -
3*tempC -
3300*salinity -
13*pMPa +
47*pMPa*salinity)));
}
/*!
* \copydoc H2O::gasIsCompressible
*/
static bool gasIsCompressible()
{ return H2O::gasIsCompressible(); }
/*!
* \copydoc H2O::liquidIsCompressible
*/
static bool liquidIsCompressible()
{ return H2O::liquidIsCompressible(); }
/*!
* \copydoc H2O::gasPressure
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, const Evaluation& density)
{ return H2O::gasPressure(temperature, density); }
/*!
* \copydoc H2O::liquidPressure
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, const Evaluation& density)
{
typedef MathToolbox<Evaluation> Toolbox;
// We use the newton method for this. For the initial value we
// assume the pressure to be 10% higher than the vapor
// pressure
Scalar pressure = 1.1*vaporPressure(temperature);
Scalar eps = pressure*1e-7;
Evaluation pressure = 1.1*vaporPressure(temperature);
Scalar eps = Toolbox::value(pressure)*1e-7;
Scalar deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
Scalar f = liquidDensity(temperature, pressure) - density;
Evaluation deltaP = pressure*2;
for (int i = 0;
i < 5
&& std::abs(Toolbox::value(pressure)*1e-9) < std::abs(Toolbox::value(deltaP));
++i)
{
const Evaluation& f = liquidDensity(temperature, pressure) - density;
Scalar df_dp;
df_dp = liquidDensity(temperature, pressure + eps);
Evaluation df_dp = liquidDensity(temperature, pressure + eps);
df_dp -= liquidDensity(temperature, pressure - eps);
df_dp /= 2*eps;
@ -299,7 +311,8 @@ public:
/*!
* \copydoc H2O::gasViscosity
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{ return H2O::gasViscosity(temperature, pressure); }
/*!
@ -310,16 +323,19 @@ public:
* - cited by: Bachu & Adams (2002)
* "Equations of State for basin geofluids"
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
if(temperature <= 275.) // regularisation
{ temperature = 275; }
Scalar T_C = temperature - 273.15;
typedef MathToolbox<Evaluation> Toolbox;
Scalar A = (0.42*std::pow((std::pow(salinity, 0.8)-0.17), 2) + 0.045)*std::pow(T_C, 0.8);
Scalar mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*std::exp(-A);
Evaluation T_C = temperature - 273.15;
if(temperature <= 275.) // regularization
T_C = Toolbox::createConstant(275.0);
return mu_brine/1000.0; /* unit: Pa s */
Evaluation A = (0.42*std::pow((std::pow(salinity, 0.8)-0.17), 2) + 0.045)*Toolbox::pow(T_C, 0.8);
Evaluation mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*Toolbox::exp(-A);
return mu_brine/1000.0; // convert to [Pa s] (todo: check if correct cP->Pa s is times 10...)
}
};

112
opm/material/components/CO2.hpp
View File

@ -29,6 +29,7 @@
#include <opm/material/Constants.hpp>
#include <opm/material/IdealGas.hpp>
#include <opm/material/components/Component.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <cmath>
#include <iostream>
@ -124,22 +125,24 @@ public:
* Physical and Chemical Reference Data, 25 (6), pp. 1509-1596,
* 1996
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar a[4] =
{ -7.0602087, 1.9391218, -1.6463597, -3.2995634 };
static const Scalar t[4] =
{ 1.0, 1.5, 2.0, 4.0 };
// this is on page 1524 of the reference
Scalar exponent = 0;
Scalar Tred = T/criticalTemperature();
for (int i = 0; i < 5; ++i) {
exponent += a[i]*std::pow(1 - Tred, t[i]);
}
Evaluation exponent = 0;
Evaluation Tred = T/criticalTemperature();
for (int i = 0; i < 5; ++i)
exponent += a[i]*Toolbox::pow(1 - Tred, t[i]);
exponent *= 1.0/Tred;
return std::exp(exponent)*criticalPressure();
return Toolbox::exp(exponent)*criticalPressure();
}
@ -158,8 +161,9 @@ public:
/*!
* \brief Specific enthalpy of gaseous CO2 [J/kg].
*/
static Scalar gasEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
return CO2Tables::tabulatedEnthalpy.eval(temperature, pressure);
}
@ -167,19 +171,21 @@ public:
/*!
* \brief Specific internal energy of CO2 [J/kg].
*/
static Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
Scalar h = gasEnthalpy(temperature, pressure);
Scalar rho = gasDensity(temperature, pressure);
const Evaluation& h = gasEnthalpy(temperature, pressure);
const Evaluation& rho = gasDensity(temperature, pressure);
return h - (pressure / rho);
}
/*!
* \brief The density of CO2 at a given pressure and temperature [kg/m^3].
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
*/
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
return CO2Tables::tabulatedDensity.eval(temperature, pressure);
}
@ -190,57 +196,46 @@ public:
* Equations given in: - Vesovic et al., 1990
* - Fenhour etl al., 1998
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(Evaluation temperature, const Evaluation& pressure)
{
static const double a0 = 0.235156;
static const double a1 = -0.491266;
static const double a2 = 5.211155E-2;
static const double a3 = 5.347906E-2;
static const double a4 = -1.537102E-2;
typedef MathToolbox<Evaluation> Toolbox;
static const double d11 = 0.4071119E-2;
static const double d21 = 0.7198037E-4;
static const double d64 = 0.2411697E-16;
static const double d81 = 0.2971072E-22;
static const double d82 = -0.1627888E-22;
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;
static const double ESP = 251.196;
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;
double mu0, SigmaStar, TStar;
double dmu, rho;
double visco_CO2;
const Scalar ESP = 251.196;
if(temperature < 275.) // regularization
temperature = 275;
temperature = Toolbox::createConstant(275.0);
Evaluation TStar = temperature/ESP;
TStar = temperature/ESP;
// mu0: viscosity in zero-density limit
const Evaluation& logTStar = Toolbox::log(TStar);
Evaluation SigmaStar = Toolbox::exp(a0 + logTStar*(a1 + logTStar*(a2 + logTStar*(a3 + logTStar*a4))));
/* mu0: viscosity in zero-density limit */
SigmaStar = exp(a0 + a1*log(TStar)
+ a2*log(TStar)*log(TStar)
+ a3*log(TStar)*log(TStar)*log(TStar)
+ a4*log(TStar)*log(TStar)*log(TStar)*log(TStar) );
Evaluation mu0 = 1.00697*Toolbox::sqrt(temperature) / SigmaStar;
mu0 = 1.00697*sqrt(temperature) / SigmaStar;
const Evaluation& rho = gasDensity(temperature, pressure); // CO2 mass density [kg/m^3]
/* dmu : excess viscosity at elevated density */
rho = gasDensity(temperature, pressure); /* CO2 mass density [kg/m^3] */
// dmu : excess viscosity at elevated density
Evaluation dmu =
d11*rho
+ d21*rho*rho
+ d64*Toolbox::pow(rho, 6.0)/(TStar*TStar*TStar)
+ d81*Toolbox::pow(rho, 8.0)
+ d82*Toolbox::pow(rho, 8.0)/TStar;
dmu = d11*rho + d21*rho*rho + d64*pow(rho,6)/(TStar*TStar*TStar)
+ d81*pow(rho,8) + d82*pow(rho,8)/TStar;
/* dmucrit : viscosity increase near the critical point */
// False (Lybke 2July2007)
//e1 = 5.5930E-3;
//e2 = 6.1757E-5;
//e4 = 2.6430E-11;
//dmucrit = e1*rho + e2*rho*rho + e4*rho*rho*rho;
//visco_CO2 = (mu0 + dmu + dmucrit)/1.0E6; /* conversion to [Pa s] */
visco_CO2 = (mu0 + dmu)/1.0E6; /* conversion to [Pa s] */
return visco_CO2;
return (mu0 + dmu)/1.0e6; // conversion to [Pa s]
}
/*!
@ -253,15 +248,16 @@ public:
* \param temperature Temperature of component \f$\mathrm{[K]}\f$
* \param pressure Pressure of component \f$\mathrm{[Pa]}\f$
*/
static Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{
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...)
Scalar h1 = gasEnthalpy(temperature - eps, pressure);
Scalar h2 = gasEnthalpy(temperature + eps, pressure);
const Evaluation& h1 = gasEnthalpy(temperature - eps, pressure);
const Evaluation& h2 = gasEnthalpy(temperature + eps, pressure);
return (h2 - h1) / (2*eps) ;
}

45
opm/material/components/Component.hpp
View File

@ -36,10 +36,12 @@ namespace Opm {
* \tparam Scalar The type used for scalar values
* \tparam Implementation Necessary for static polymorphism
*/
template <class Scalar, class Implementation>
template <class ScalarT, class Implementation>
class Component
{
public:
typedef ScalarT Scalar;
static const bool isTabulated = false;
/*!
@ -116,9 +118,10 @@ public:
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
* temperature in \f$\mathrm{[K]}\f$.
*
* \param T temperature of the component in \f$\mathrm{[K]}\f$
* \param temperature temperature of the component in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::vaporPressure()"); }
/*!
@ -127,7 +130,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::gasDensity()"); }
/*!
@ -136,7 +140,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidDensity()"); }
/*!
@ -145,7 +150,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::gasEnthalpy()"); }
/*!
@ -154,7 +160,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidEnthalpy()"); }
/*!
@ -163,7 +170,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::gasInternalEnergy()"); }
/*!
@ -172,7 +180,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidInternalEnergy()"); }
/*!
@ -182,7 +191,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::gasViscosity()"); }
/*!
@ -191,31 +201,36 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidViscosity()"); }
/*!
* \brief Thermal conductivity of the component [W/(m^2 K/m)] as a gas.
*/
static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::gasThermalConductivity()"); }
/*!
* \brief Thermal conductivity of the component [W/(m^2 K/m)] as a liquid.
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidThermalConductivity()"); }
/*!
* \brief Specific isobaric heat capacity of the component [J/kg] as a gas.
*/
static Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::gasHeatCapacity()"); }
/*!
* \brief Specific isobaric heat capacity of the component [J/kg] as a liquid.
*/
static Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{ OPM_THROW(std::runtime_error, "Not implemented: Component::liquidHeatCapacity()"); }
};

90
opm/material/components/Dnapl.hpp
View File

@ -24,9 +24,11 @@
#ifndef OPM_DNAPL_HPP
#define OPM_DNAPL_HPP
#include <opm/material/IdealGas.hpp>
#include "Component.hpp"
#include <opm/material/IdealGas.hpp>
#include <opm/material/common/MathToolbox.hpp>
namespace Opm {
/*!
* \ingroup Components
@ -49,19 +51,10 @@ public:
{ return "DNAPL"; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of TCE.
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static Scalar molarMass()
{ return 131.39e-3; /* [kg/mol] */ }
/*!
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure TCE
* at a given temperature.
*
* \param T temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar T)
{ return 3900; /* [Pa] (at 20C) */ }
static bool gasIsIdeal()
{ return true; }
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
@ -75,33 +68,55 @@ public:
static bool liquidIsCompressible()
{ return false; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of TCE.
*/
static Scalar molarMass()
{
return 131.39e-3; // [kg/mol]
}
/*!
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure TCE
* at a given temperature.
*
* \param T temperature of component in \f$\mathrm{[K]}\f$
*/
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(3900); // [Pa] (at 20C)
}
/*!
* \brief The density of steam at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
*/
template <class Evaluation>
static Evaluation gasDensity(Evaluation temperature, Evaluation pressure)
{
return IdealGas<Scalar>::density(molarMass(),
return IdealGas<Scalar>::density(Evaluation(molarMass()),
temperature,
pressure);
}
/*!
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static bool gasIsIdeal()
{ return true; }
/*!
* \brief The density of pure TCE at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{ return 1460.0; /* [kg/m^3] */ }
template <class Evaluation>
static Evaluation liquidDensity(Evaluation temperature, Evaluation pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(1460.0); // [kg/m^3]
}
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure TCE.
@ -109,8 +124,13 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{ return 5.7e-4; /* [Pa s] */ }
template <class Evaluation>
static Evaluation liquidViscosity(Evaluation temperature, Evaluation pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(5.7e-4); // [Pa s]
}
/*!
* \brief The enthalpy of pure TCE at a given pressure and temperature \f$\mathrm{[J/kg]}\f$.
@ -118,19 +138,27 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
{ return 120.0/molarMass() * temperature; /* [J/kg] */ }
template <class Evaluation>
static Evaluation liquidEnthalpy(Evaluation temperature, Evaluation pressure)
{
return 120.0/molarMass() * temperature; // [J/kg]
}
/*!
* \brief Specific heat conductivity of liquid TCE \f$\mathrm{[W/(m K)]}\f$.
*
* \todo The value returned here is a guess which does not necessarily correspond to reality in any way!
* \todo The value returned here is a guess which does not necessarily correspond to reality in any way!
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
{ return 0.3; }
template <class Evaluation>
static Evaluation liquidThermalConductivity(Evaluation temperature, Evaluation pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(0.3);
}
};
} // namespace Opm

291
opm/material/components/H2O.hpp
View File

@ -125,14 +125,15 @@ public:
*
* \param T Absolute temperature of the system in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(Evaluation temperature)
{
if (T > criticalTemperature())
T = criticalTemperature();
if (T < tripleTemperature())
T = tripleTemperature();
if (temperature > criticalTemperature())
temperature = criticalTemperature();
if (temperature < tripleTemperature())
temperature = tripleTemperature();
return Region4::saturationPressure(T);
return Region4::saturationPressure(temperature);
}
/*!
* \brief The vapor temperature in \f$\mathrm{[Ka]}\f$ of pure water
@ -146,7 +147,8 @@ public:
*
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar vaporTemperature(Scalar pressure)
template <class Evaluation>
static Evaluation vaporTemperature(const Evaluation& pressure)
{
if (pressure > criticalPressure())
pressure = criticalPressure();
@ -168,14 +170,17 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Enthalpy of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Enthalpy of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
@ -186,13 +191,13 @@ public:
// dependence on pressure, so we can just return the
// specific enthalpy at the point of regularization, i.e.
// the triple pressure - 100Pa
return enthalpyRegion2_(temperature, triplePressure() - 100);
return enthalpyRegion2_(temperature, Toolbox::createConstant(triplePressure() - 100));
}
Scalar pv = vaporPressure(temperature);
Evaluation pv = vaporPressure(temperature);
if (pressure > pv) {
// the pressure is too high, in this case we use the slope
// of the enthalpy at the vapor pressure to regularize
Scalar dh_dp =
Evaluation dh_dp =
Rs*temperature*
Region2::tau(temperature)*
Region2::dpi_dp(pv)*
@ -218,25 +223,28 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Enthalpy of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Enthalpy of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
Scalar pv = vaporPressure(temperature);
const Evaluation& pv = vaporPressure(temperature);
if (pressure < pv) {
// the pressure is too low, in this case we use the slope
// of the enthalpy at the vapor pressure to regularize
Scalar dh_dp =
const Evaluation& dh_dp =
Rs * temperature*
Region1::tau(temperature)*
Region1::dpi_dp(pv)*
Region1::dpi_dp(Toolbox::value(pv))*
Region1::ddgamma_dtaudpi(temperature, pv);
return
@ -259,27 +267,26 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasHeatCapacity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Heat capacity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Heat capacity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
if (pressure < triplePressure() - 100) {
return heatCap_p_Region2_(temperature, triplePressure() - 100);
}
Scalar pv = vaporPressure(temperature);
if (pressure > pv) {
if (pressure < triplePressure() - 100)
return heatCap_p_Region2_(temperature, Evaluation(triplePressure() - 100));
const Evaluation& pv = vaporPressure(temperature);
if (pressure > pv)
// the pressure is too high, in this case we use the heat
// cap at the vapor pressure to regularize
return
heatCap_p_Region2_(temperature, pv);
};
return heatCap_p_Region2_(temperature, pv);
return heatCap_p_Region2_(temperature, pressure);
}
@ -295,22 +302,23 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidHeatCapacity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"heat Capacity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"heat Capacity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
Scalar pv = vaporPressure(temperature);
const Evaluation& pv = vaporPressure(temperature);
if (pressure < pv) {
// the pressure is too low, in this case we use the heat cap at the vapor pressure to regularize
return
heatCap_p_Region1_(temperature, pv);
// the pressure is too low, in this case we use the heat capacity at the
// vapor pressure to regularize
return heatCap_p_Region1_(temperature, pv);
};
return heatCap_p_Region1_(temperature, pressure);
@ -328,14 +336,15 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Internal Energy of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Internal Energy of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@ -385,13 +394,14 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Internal Energy of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Internal Energy of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
@ -458,14 +468,15 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidHeatCapacityConstVolume(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacityConstVolume(const Evaluation& temperature,
const Evaluation& pressure)
{
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Heat capacity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Heat capacity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@ -492,13 +503,14 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasHeatCapacityConstVolume(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacityConstVolume(const Evaluation& temperature, const Evaluation& pressure)
{
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Heat capacity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Heat capacity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
@ -538,31 +550,36 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Density of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Density of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
if (pressure < triplePressure() - 100) {
// We assume an ideal gas for low pressures to avoid the
// 0/0 for the internal energy and enthalpy.
Scalar rho0IAPWS = 1.0/volumeRegion2_(temperature,
triplePressure() - 100);
Scalar rho0Id = IdealGas<Scalar>::density(molarMass(),
temperature,
triplePressure() - 100);
return
rho0IAPWS/rho0Id *
IdealGas<Scalar>::density(molarMass(),
const Evaluation& rho0IAPWS =
1.0/volumeRegion2_(temperature,
Evaluation(triplePressure() - 100));
const Evaluation& rho0Id =
IdealGas<Scalar>::density(Evaluation(molarMass()),
temperature,
pressure);
Evaluation(triplePressure() - 100));
return
rho0IAPWS/rho0Id
*IdealGas<Scalar>::density(Evaluation(molarMass()),
temperature,
pressure);
}
Scalar pv = vaporPressure(temperature);
Evaluation pv = vaporPressure(temperature);
if (pressure > pv) {
// the pressure is too high, in this case we use the slope
// of the density energy at the vapor pressure to
@ -570,10 +587,10 @@ public:
// calculate the partial derivative of the specific volume
// to the pressure at the vapor pressure.
Scalar eps = pv*1e-8;
Scalar v0 = volumeRegion2_(temperature, pv);
Scalar v1 = volumeRegion2_(temperature, pv + eps);
Scalar dv_dp = (v1 - v0)/eps;
Scalar eps = Toolbox::value(pv)*1e-8;
Evaluation v0 = volumeRegion2_(temperature, pv);
Evaluation v1 = volumeRegion2_(temperature, pv + eps);
Evaluation dv_dp = (v1 - v0)/eps;
/*
Scalar pi = Region2::pi(pv);
Scalar dp_dpi = Region2::dp_dpi(pv);
@ -589,7 +606,7 @@ public:
// calculate the partial derivative of the density to the
// pressure at vapor pressure
Scalar drho_dp = - 1/(v0*v0)*dv_dp;
Evaluation drho_dp = - 1/(v0*v0)*dv_dp;
// use a straight line for extrapolation
return 1.0/v0 + (pressure - pv)*drho_dp;
@ -616,23 +633,26 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param density Density in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, Scalar density)
{
typedef MathToolbox<Evaluation> Toolbox;
Valgrind::CheckDefined(temperature);
Valgrind::CheckDefined(density);
// We use the newton method for this. For the initial value we
// assume steam to be an ideal gas
Scalar pressure = IdealGas<Scalar>::pressure(temperature, density/molarMass());
Evaluation pressure = IdealGas<Scalar>::pressure(temperature, density/molarMass());
Scalar eps = pressure*1e-7;
Scalar deltaP = pressure*2;
Evaluation deltaP = pressure*2;
Valgrind::CheckDefined(pressure);
Valgrind::CheckDefined(deltaP);
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
Scalar f = gasDensity(temperature, pressure) - density;
for (int i = 0; i < 5 && std::abs(Toolbox::value(pressure)*1e-9) < std::abs(Toolbox::value(deltaP)); ++i) {
Evaluation f = gasDensity(temperature, pressure) - density;
Scalar df_dp;
Evaluation df_dp;
df_dp = gasDensity(temperature, pressure + eps);
df_dp -= gasDensity(temperature, pressure - eps);
df_dp /= 2*eps;
@ -659,27 +679,30 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Density of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Density of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
// regularization
Scalar pv = vaporPressure(temperature);
Evaluation pv = vaporPressure(temperature);
if (pressure < pv) {
// the pressure is too low, in this case we use the slope
// of the density at the vapor pressure to regularize
// calculate the partial derivative of the specific volume
// to the pressure at the vapor pressure.
Scalar eps = pv*1e-8;
Scalar v0 = volumeRegion1_(temperature, pv);
Scalar v1 = volumeRegion1_(temperature, pv + eps);
Scalar dv_dp = (v1 - v0)/eps;
Scalar eps = Toolbox::value(pv)*1e-8;
Evaluation v0 = volumeRegion1_(temperature, pv);
Evaluation v1 = volumeRegion1_(temperature, pv + eps);
Evaluation dv_dp = (v1 - v0)/eps;
/*
Scalar v0 = volumeRegion1_(temperature, pv);
@ -697,7 +720,7 @@ public:
// calculate the partial derivative of the density to the
// pressure at vapor pressure
Scalar drho_dp = - 1/(v0*v0)*dv_dp;
Evaluation drho_dp = - 1/(v0*v0)*dv_dp;
// use a straight line for extrapolation
return 1.0/v0 + (pressure - pv)*drho_dp;
@ -719,19 +742,22 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param density Density of the fluid in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, Scalar density)
{
// We use the newton method for this. For the initial value we
typedef MathToolbox<Evaluation> Toolbox;
// We use the Newton method for this. For the initial value we
// assume the pressure to be 10% higher than the vapor
// pressure
Scalar pressure = 1.1*vaporPressure(temperature);
Scalar eps = pressure*1e-7;
Evaluation pressure = 1.1*vaporPressure(temperature);
Scalar eps = Toolbox::value(pressure)*1e-7;
Scalar deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
Scalar f = liquidDensity(temperature, pressure) - density;
Evaluation deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(Toolbox::value(pressure)*1e-9) < std::abs(Toolbox::value(deltaP)); ++i) {
Evaluation f = liquidDensity(temperature, pressure) - density;
Scalar df_dp;
Evaluation df_dp;
df_dp = liquidDensity(temperature, pressure + eps);
df_dp -= liquidDensity(temperature, pressure - eps);
df_dp /= 2*eps;
@ -758,16 +784,17 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Viscosity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Viscosity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
Scalar rho = gasDensity(temperature, pressure);
Evaluation rho = gasDensity(temperature, pressure);
return Common::viscosity(temperature, rho);
}
@ -782,16 +809,17 @@ public:
* \param temperature Absolute temperature of the fluid in \f$\mathrm{[K]}\f$
* \param pressure Phase pressure in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalIssue,
"Viscosity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
"Viscosity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
};
Scalar rho = liquidDensity(temperature, pressure);
const Evaluation& rho = liquidDensity(temperature, pressure);
return Common::viscosity(temperature, rho);
}
@ -808,9 +836,10 @@ public:
* \param temperature Absolute temperature in K
* \param pressure Phase pressure of the phase in Pa
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar rho = liquidDensity(temperature, pressure);
const Evaluation& rho = liquidDensity(temperature, pressure);
return Common::thermalConductivityIAPWS(temperature, rho);
}
@ -827,15 +856,17 @@ public:
* \param temperature Absolute temperature in K
* \param pressure Phase pressure of the phase in Pa
*/
static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar rho = gasDensity(temperature, pressure);
const Evaluation& rho = gasDensity(temperature, pressure);
return Common::thermalConductivityIAPWS(temperature, rho);
}
private:
// the unregularized specific enthalpy for liquid water
static Scalar enthalpyRegion1_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation enthalpyRegion1_(const Evaluation& temperature, const Evaluation& pressure)
{
return
Region1::tau(temperature) *
@ -844,16 +875,19 @@ private:
}
// the unregularized specific isobaric heat capacity
static Scalar heatCap_p_Region1_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation heatCap_p_Region1_(const Evaluation& temperature, const Evaluation& pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
return
- std::pow(Region1::tau(temperature), 2 ) *
- Toolbox::pow(Region1::tau(temperature), 2.0) *
Region1::ddgamma_ddtau(temperature, pressure) *
Rs;
}
// the unregularized specific isochoric heat capacity
static Scalar heatCap_v_Region1_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation heatCap_v_Region1_(const Evaluation& temperature, const Evaluation& pressure)
{
double tau = Region1::tau(temperature);
double num = Region1::dgamma_dpi(temperature, pressure) - tau * Region1::ddgamma_dtaudpi(temperature, pressure);
@ -866,7 +900,8 @@ private:
}
// the unregularized specific internal energy for liquid water
static Scalar internalEnergyRegion1_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation internalEnergyRegion1_(const Evaluation& temperature, const Evaluation& pressure)
{
return
Rs * temperature *
@ -875,7 +910,8 @@ private:
}
// the unregularized specific volume for liquid water
static Scalar volumeRegion1_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation volumeRegion1_(const Evaluation& temperature, const Evaluation& pressure)
{
return
Region1::pi(pressure)*
@ -884,7 +920,8 @@ private:
}
// the unregularized specific enthalpy for steam
static Scalar enthalpyRegion2_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation enthalpyRegion2_(const Evaluation& temperature, const Evaluation& pressure)
{
return
Region2::tau(temperature) *
@ -893,7 +930,8 @@ private:
}
// the unregularized specific internal energy for steam
static Scalar internalEnergyRegion2_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation internalEnergyRegion2_(const Evaluation& temperature, const Evaluation& pressure)
{
return
Rs * temperature *
@ -902,21 +940,25 @@ private:
}
// the unregularized specific isobaric heat capacity
static Scalar heatCap_p_Region2_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation heatCap_p_Region2_(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return
- std::pow(Region2::tau(temperature), 2 ) *
- Toolbox::pow(Region2::tau(temperature), 2 ) *
Region2::ddgamma_ddtau(temperature, pressure) *
Rs;
}
// the unregularized specific isochoric heat capacity
static Scalar heatCap_v_Region2_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation heatCap_v_Region2_(const Evaluation& temperature, const Evaluation& pressure)
{
double tau = Region2::tau(temperature);
double pi = Region2::pi(pressure);
double num = 1 + pi * Region2::dgamma_dpi(temperature, pressure) + tau * pi * Region2::ddgamma_dtaudpi(temperature, pressure);
double diff = num * num / (1 - pi * pi * Region2::ddgamma_ddpi(temperature, pressure));
const Evaluation& tau = Region2::tau(temperature);
const Evaluation& pi = Region2::pi(pressure);
const Evaluation& num = 1 + pi * Region2::dgamma_dpi(temperature, pressure) + tau * pi * Region2::ddgamma_dtaudpi(temperature, pressure);
const Evaluation& diff = num * num / (1 - pi * pi * Region2::ddgamma_ddpi(temperature, pressure));
return
- std::pow(tau, 2 ) *
Region2::ddgamma_ddtau(temperature, pressure) * Rs
@ -924,7 +966,8 @@ private:
}
// the unregularized specific volume for steam
static Scalar volumeRegion2_(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation volumeRegion2_(const Evaluation& temperature, const Evaluation& pressure)
{
return
Region2::pi(pressure)*

6
opm/material/components/Lnapl.hpp
View File

@ -57,7 +57,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return 890; }
/*!
@ -66,7 +67,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{ return 8e-3; }
};

114
opm/material/components/Mesitylene.hpp
View File

@ -28,8 +28,9 @@
#include <opm/material/components/Component.hpp>
#include <opm/material/Constants.hpp>
namespace Opm
{
#include <opm/material/common/MathToolbox.hpp>
namespace Opm {
/*!
* \ingroup Components
* \brief Component for Mesitylene
@ -91,15 +92,18 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar temperature)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar A = 7.07638;
const Scalar B = 1571.005;
const Scalar C = 209.728;
const Scalar T = temperature - 273.15;
const Evaluation& T = temperature - 273.15;
return 100 * 1.334 * std::pow(Scalar(10.0), Scalar(A - (B / (T + C))));
return 100 * 1.334 * Toolbox::pow(10.0, A - (B / (T + C)));
}
@ -109,7 +113,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
// Gauss quadrature rule:
// Interval: [0K; temperature (K)]
@ -133,11 +138,13 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar heatVap(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation heatVap(Evaluation temperature, const Evaluation& pressure)
{
temperature = std::min(temperature, criticalTemperature()); // regularization
temperature = std::max(temperature, 0.0); // regularization
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, criticalTemperature()); // regularization
temperature = Toolbox::max(temperature, 0.0); // regularization
const Scalar T_crit = criticalTemperature();
const Scalar Tr1 = boilingTemperature()/criticalTemperature();
@ -150,9 +157,9 @@ public:
/ (1.07 - Tr1); /* [J/mol] */
/* Variation with temp according to Watson relation eq 7-12.1*/
const Scalar Tr2 = temperature/criticalTemperature();
const Evaluation& Tr2 = temperature/criticalTemperature();
const Scalar n = 0.375;
const Scalar DH_vap = DH_v_boil * std::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
const Evaluation& DH_vap = DH_v_boil * Toolbox::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
return (DH_vap/molarMass()); // we need [J/kg]
}
@ -167,7 +174,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
return liquidEnthalpy(temperature,pressure) + heatVap(temperature, pressure);
}
@ -178,12 +186,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
return IdealGas<Scalar>::density(molarMass(),
temperature,
pressure);
}
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return IdealGas<Scalar>::density(Evaluation(molarMass()), temperature, pressure); }
/*!
* \brief The density of pure mesitylene at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
@ -191,7 +196,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return molarLiquidDensity_(temperature)*molarMass(); }
/*!
@ -219,20 +225,23 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
* \param regularize defines, if the functions is regularized or not, set to true by default
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure, bool regularize=true)
template <class Evaluation>
static Evaluation gasViscosity(Evaluation temperature, const Evaluation& pressure, bool regularize=true)
{
temperature = std::min(temperature, 500.0); // regularization
temperature = std::max(temperature, 250.0);
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0);
// reduced temperature
Scalar Tr = temperature/criticalTemperature();
const Evaluation& Tr = temperature/criticalTemperature();
Scalar Fp0 = 1.0;
Scalar xi = 0.00474;
Scalar eta_xi =
Fp0*(0.807*std::pow(Tr,0.618)
- 0.357*std::exp(-0.449*Tr)
+ 0.34*std::exp(-4.058*Tr)
const Evaluation& eta_xi =
Fp0*(0.807*Toolbox::pow(Tr,0.618)
- 0.357*Toolbox::exp(-0.449*Tr)
+ 0.34*Toolbox::exp(-4.058*Tr)
+ 0.018);
return eta_xi/xi/1e7; // [Pa s]
@ -244,15 +253,18 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(Evaluation temperature, const Evaluation& pressure)
{
temperature = std::min(temperature, 500.0); // regularization
temperature = std::max(temperature, 250.0);
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0);
const Scalar A = -6.749;
const Scalar B = 2010.0;
return std::exp(A + B/temperature)*1e-3; // [Pa s]
return Toolbox::exp(A + B/temperature)*1e-3; // [Pa s]
}
/*!
@ -264,28 +276,29 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
*/
static Scalar liquidHeatCapacity(const Scalar temperature,
const Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{
/* according Reid et al. : Missenard group contrib. method (s. example 5-8) */
/* Mesitylen: C9H12 : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was to much!) */
/* linear interpolation between table values [J/(mol K)]*/
Scalar H, CH3, C6H5;
Evaluation H, CH3, C6H5;
if(temperature<298.) {
// extrapolation for temperature < 273K
H = 13.4+1.2*(temperature-273.0)/25.; // 13.4 + 1.2 = 14.6 = H(T=298K) i.e. interpolation of table values 273<T<298
CH3 = 40.0+1.6*(temperature-273.0)/25.; // 40 + 1.6 = 41.6 = CH3(T=298K)
C6H5 = 113.0+4.2*(temperature-273.0)/25.; // 113 + 4.2 =117.2 = C6H5(T=298K)
H = 13.4 + 1.2*(temperature-273.0)/25.; // 13.4 + 1.2 = 14.6 = H(T=298K) i.e. interpolation of table values 273<T<298
CH3 = 40.0 + 1.6*(temperature-273.0)/25.; // 40 + 1.6 = 41.6 = CH3(T=298K)
C6H5 = 113.0 + 4.2*(temperature-273.0)/25.; // 113 + 4.2 =117.2 = C6H5(T=298K)
}
else if((temperature>=298.0)&&(temperature<323.)){ // i.e. interpolation of table values 298<T<323
H = 14.6+0.9*(temperature-298.0)/25.;
CH3 = 41.6+1.9*(temperature-298.0)/25.;
C6H5 = 117.2+6.2*(temperature-298.0)/25.;
H = 14.6 + 0.9*(temperature-298.0)/25.;
CH3 = 41.6 + 1.9*(temperature-298.0)/25.;
C6H5 = 117.2 + 6.2*(temperature-298.0)/25.;
}
else if((temperature>=323.0)&&(temperature<348.)){// i.e. interpolation of table values 323<T<348
H = 15.5+1.2*(temperature-323.0)/25.;
CH3 = 43.5+2.3*(temperature-323.0)/25.;
C6H5 = 123.4+6.3*(temperature-323.0)/25.;
H = 15.5 + 1.2*(temperature-323.0)/25.;
CH3 = 43.5 + 2.3*(temperature-323.0)/25.;
C6H5 = 123.4 + 6.3*(temperature-323.0)/25.;
}
else {
assert(temperature>=348.0);
@ -308,14 +321,17 @@ protected:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar molarLiquidDensity_(Scalar temperature)
template <class Evaluation>
static Evaluation molarLiquidDensity_(Evaluation temperature)
{
temperature = std::min(temperature, 500.0); // regularization
temperature = std::max(temperature, 250.0);
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0);
const Scalar Z_RA = 0.2556; // from equation
const Scalar expo = 1.0 + std::pow(1.0 - temperature/criticalTemperature(), 2.0/7.0);
Scalar V = Consts::R*criticalTemperature()/criticalPressure()*std::pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
const Evaluation& expo = 1.0 + Toolbox::pow(1.0 - temperature/criticalTemperature(), 2.0/7.0);
const Evaluation& V = Consts::R*criticalTemperature()/criticalPressure()*Toolbox::pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
return 1.0/V; // molar density [mol/m^3]
}

89
opm/material/components/N2.hpp
View File

@ -27,6 +27,7 @@
#define OPM_N2_HPP
#include <opm/material/IdealGas.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include "Component.hpp"
@ -88,7 +89,7 @@ public:
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure molecular nitrogen
* at a given temperature.
*
*\param T temperature of component in \f$\mathrm{[K]}\f$
*\param temperature temperature of component in \f$\mathrm{[K]}\f$
*
* Taken from:
*
@ -98,28 +99,31 @@ public:
* Physical and Chemical Refefence Data, Vol. 29, No. 6,
* pp. 1361-1433
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
if (T > criticalTemperature())
typedef MathToolbox<Evaluation> Toolbox;
if (temperature > criticalTemperature())
return criticalPressure();
if (T < tripleTemperature())
if (temperature < tripleTemperature())
return 0; // N2 is solid: We don't take sublimation into
// account
// note: this is the ancillary equation given on page 1368
Scalar sigma = Scalar(1.0) - T/criticalTemperature();
Scalar sqrtSigma = std::sqrt(sigma);
const Evaluation& sigma = 1.0 - temperature/criticalTemperature();
const Evaluation& sqrtSigma = Toolbox::sqrt(sigma);
const Scalar N1 = -6.12445284;
const Scalar N2 = 1.26327220;
const Scalar N3 = -0.765910082;
const Scalar N4 = -1.77570564;
return
criticalPressure() *
std::exp(criticalTemperature()/T*
(sigma*(N1 +
sqrtSigma*N2 +
sigma*(sqrtSigma*N3 +
sigma*sigma*sigma*N4))));
Toolbox::exp(criticalTemperature()/temperature*
(sigma*(N1 +
sqrtSigma*N2 +
sigma*(sqrtSigma*N3 +
sigma*sigma*sigma*N4))));
}
/*!
@ -128,13 +132,14 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), temperature, pressure);
return IdealGas::density(Evaluation(molarMass()), temperature, pressure);
}
/*!
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
*/
static bool gasIsCompressible()
@ -152,7 +157,8 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, const Evaluation& density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
@ -161,14 +167,15 @@ public:
/*!
* \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure nitrogen gas.
*
* \param T temperature of component in \f$\mathrm{[K]}\f$
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* See: R. Reid, et al.: The Properties of Gases and Liquids, 4th
* edition, McGraw-Hill, 1987, pp 154, 657, 665
*/
static const Scalar gasEnthalpy(Scalar T,
Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
// method of Joback
const Scalar cpVapA = 31.15;
@ -180,10 +187,10 @@ public:
return
1/molarMass()* // conversion from [J/(mol K)] to [J/(kg K)]
T*(cpVapA + T*
(cpVapB/2 + T*
(cpVapC/3 + T*
(cpVapD/4))));
temperature*(cpVapA + temperature*
(cpVapB/2 + temperature*
(cpVapC/3 + temperature*
(cpVapD/4))));
//#warning NIST DATA STUPID INTERPOLATION
// Scalar T2 = 300.;
@ -207,8 +214,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
gasEnthalpy(temperature, pressure) -
@ -223,8 +231,9 @@ public:
* This is equivalent to the partial derivative of the specific
* enthalpy to the temperature.
*/
static const Scalar gasHeatCapacity(Scalar T,
Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature,
const Evaluation& pressure)
{
// method of Joback
const Scalar cpVapA = 31.15;
@ -235,9 +244,9 @@ public:
return
1/molarMass()* // conversion from [J/(mol K)] to [J/(kg K)]
cpVapA + T*
(cpVapB + T*
(cpVapC + T*
cpVapA + temperature*
(cpVapB + temperature*
(cpVapC + temperature*
(cpVapD)));
}
@ -254,8 +263,11 @@ public:
* 5th edition, McGraw-Hill, 2001 pp 9.7-9.8 (omega and V_c taken from p. A.19)
*
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar Tc = criticalTemperature();
const Scalar Vc = 90.1; // critical specific volume [cm^3/mol]
const Scalar omega = 0.037; // accentric factor
@ -267,12 +279,12 @@ public:
mu_r4 *= mu_r4;
Scalar Fc = 1 - 0.2756*omega + 0.059035*mu_r4;
Scalar Tstar = 1.2593 * temperature/Tc;
Scalar Omega_v =
1.16145*std::pow(Tstar, -0.14874) +
0.52487*std::exp(- 0.77320*Tstar) +
2.16178*std::exp(- 2.43787*Tstar);
Scalar mu = 40.785*Fc*std::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
const Evaluation& Tstar = 1.2593 * temperature/Tc;
const Evaluation& Omega_v =
1.16145*Toolbox::pow(Tstar, -0.14874) +
0.52487*Toolbox::exp(- 0.77320*Tstar) +
2.16178*Toolbox::exp(- 2.43787*Tstar);
const Evaluation& mu = 40.785*Fc*Toolbox::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;
@ -289,8 +301,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasThermalConductivity(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature,
const Evaluation& pressure)
{ return 0.024572; }
};

2
opm/material/components/NullComponent.hpp
View File

@ -35,7 +35,7 @@ namespace Opm
* Its main purpose is to make things compile and give runtime errors
* if it is actually used.
*/
template <class Scalar, int isLiquidV=-1>
template <class Scalar>
class NullComponent : public Component<Scalar, NullComponent<Scalar> >
{
};

23
opm/material/components/SimpleCO2.hpp
View File

@ -97,22 +97,25 @@ public:
/*!
* \copydoc Component::gasEnthalpy
*/
static const Scalar gasEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{ return 571.3e3 + (temperature - 298.15)*0.85e3; }
/*!
* \copydoc Component::liquidEnthalpy
*/
static const Scalar liquidEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{ return (temperature - 298.15)*5e3; }
/*!
* \copydoc Component::gasInternalEnergy
*/
static const Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
gasEnthalpy(temperature, pressure) -
@ -123,10 +126,11 @@ public:
/*!
* \copydoc Component::gasDensity
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), temperature, pressure);
return IdealGas::density(Evaluation(molarMass()), temperature, pressure);
}
/*!
@ -137,7 +141,8 @@ public:
* See: R. Reid, et al.: The Properties of Gases and Liquids, 4th
* edition, McGraw-Hill, 1987, pp 396-397, 667
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
const Scalar Tc = criticalTemperature();
const Scalar Vc = 93.9; // critical specific volume [cm^3/mol]

131
opm/material/components/SimpleH2O.hpp
View File

@ -24,6 +24,7 @@
#define OPM_SIMPLE_H2O_HPP
#include <opm/material/IdealGas.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include "Component.hpp"
@ -59,6 +60,24 @@ public:
static const char *name()
{ return "H2O"; }
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
*/
static bool gasIsCompressible()
{ return true; }
/*!
* \brief Returns true iff the liquid phase is assumed to be compressible
*/
static bool liquidIsCompressible()
{ return false; }
/*!
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static bool gasIsIdeal()
{ return true; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of water.
*/
@ -101,8 +120,11 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (T > criticalTemperature())
return criticalPressure();
if (T < tripleTemperature())
@ -115,17 +137,17 @@ public:
0.65017534844798e3
};
Scalar sigma = T + n[8]/(T - n[9]);
Evaluation sigma = T + n[8]/(T - n[9]);
Scalar A = (sigma + n[0])*sigma + n[1];
Scalar B = (n[2]*sigma + n[3])*sigma + n[4];
Scalar C = (n[5]*sigma + n[6])*sigma + n[7];
Evaluation A = (sigma + n[0])*sigma + n[1];
Evaluation B = (n[2]*sigma + n[3])*sigma + n[4];
Evaluation C = (n[5]*sigma + n[6])*sigma + n[7];
Scalar tmp = Scalar(2.0)*C/(std::sqrt(B*B - Scalar(4.0)*A*C) - B);
Evaluation tmp = 2.0*C/(Toolbox::sqrt(B*B - 4.0*A*C) - B);
tmp *= tmp;
tmp *= tmp;
return Scalar(1e6)*tmp;
return 1e6*tmp;
}
/*!
@ -134,8 +156,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{ return 1976*(temperature - 293.15) + 2.45e6; }
/*!
@ -144,8 +167,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidEnthalpy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{ return 4180*(temperature - 293.15); }
/*!
@ -161,8 +185,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
gasEnthalpy(temperature, pressure) -
@ -176,12 +201,14 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidInternalEnergy(Scalar temperature,
Scalar pressure)
{ return
template <class Evaluation>
static Evaluation liquidInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
return
liquidEnthalpy(temperature, pressure) -
pressure/liquidDensity(temperature, pressure); }
pressure/liquidDensity(temperature, pressure);
}
/*!
* \brief Specific heat conductivity of liquid water \f$\mathrm{[W/(m K)]}\f$.
@ -189,9 +216,13 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidThermalConductivity(Scalar temperature,
Scalar pressure)
{ return 0.578078; } // conductivity of liquid water [W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C
template <class Evaluation>
static Evaluation liquidThermalConductivity(const Evaluation& temperature,
const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(0.578078); // conductivity of liquid water [W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C
}
/*!
* \brief Specific heat conductivity of steam \f$\mathrm{[W/(m K)]}\f$.
@ -199,21 +230,13 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasThermalConductivity(Scalar temperature,
Scalar pressure)
{ return 0.028224; } // conductivity of steam [W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
*/
static bool gasIsCompressible()
{ return true; }
/*!
* \brief Returns true iff the liquid phase is assumed to be compressible
*/
static bool liquidIsCompressible()
{ return false; }
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature,
const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(0.028224); // conductivity of steam [W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C
}
/*!
* \brief The density \f$\mathrm{[kg/m^3]}\f$ of steam at a given pressure and temperature.
@ -221,25 +244,21 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
// Assume an ideal gas
return molarMass()*IdealGas::molarDensity(temperature, pressure);
}
/*!
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static bool gasIsIdeal()
{ return true; }
/*!
* \brief The pressure of steam in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, const Evaluation& density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
@ -251,9 +270,11 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{
return 1000;
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(1000);
}
/*!
@ -262,10 +283,11 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, const Evaluation& density)
{
OPM_THROW(std::logic_error,
"The liquid pressure is undefined for incompressible fluids");
"The liquid pressure is undefined for incompressible fluids");
}
/*!
@ -275,9 +297,12 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
* \param regularize defines, if the functions is regularized or not, set to true by default
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure, bool regularize=true)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature,
const Evaluation& pressure)
{
return 1e-05;
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(1e-05);
}
/*!
@ -286,9 +311,11 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
return 1e-03;
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::createConstant(1e-03);
}
};

298
opm/material/components/TabulatedComponent.hpp
View File

@ -32,9 +32,9 @@
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/ErrorMacros.hpp>
namespace Opm
{
#include <opm/material/common/MathToolbox.hpp>
namespace Opm {
/*!
* \ingroup Components
*
@ -50,10 +50,12 @@ namespace Opm
* vapor pressure curve, if false use the
* pressure range [p_min, p_max]
*/
template <class Scalar, class RawComponent, bool useVaporPressure=true>
template <class ScalarT, class RawComponent, bool useVaporPressure=true>
class TabulatedComponent
{
public:
typedef ScalarT Scalar;
static const bool isTabulated = true;
/*!
@ -248,12 +250,14 @@ public:
*
* \param T temperature of component
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
Scalar result = interpolateT_(vaporPressure_, T);
if (std::isnan(result)) {
return RawComponent::vaporPressure(T);
}
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateT_(vaporPressure_, temperature);
if (std::isnan(Toolbox::value(result)))
return RawComponent::vaporPressure(temperature);
return result;
}
@ -263,14 +267,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateGasTP_(gasEnthalpy_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasEnthalpy_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::gasEnthalpy(temperature, pressure);
}
return result;
}
@ -280,14 +286,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateLiquidTP_(liquidEnthalpy_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidEnthalpy_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::liquidEnthalpy(temperature, pressure);
}
return result;
}
@ -297,14 +305,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateGasTP_(gasHeatCapacity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasHeatCapacity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::gasHeatCapacity(temperature, pressure);
}
return result;
}
@ -314,14 +324,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateLiquidTP_(liquidHeatCapacity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidHeatCapacity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::liquidHeatCapacity(temperature, pressure);
}
return result;
}
@ -331,12 +343,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
{
Scalar result =
gasEnthalpy(temperature, pressure) - pressure/gasDensity(temperature, pressure);
return result;
}
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{ return gasEnthalpy(temperature, pressure) - pressure/gasDensity(temperature, pressure); }
/*!
* \brief Specific internal energy of the liquid \f$\mathrm{[J/kg]}\f$.
@ -344,12 +353,9 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
{
Scalar result =
liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure);
return result;
}
template <class Evaluation>
static Evaluation liquidInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{ return liquidEnthalpy(temperature, pressure) - pressure/liquidDensity(temperature, pressure); }
/*!
* \brief The pressure of gas in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
@ -357,15 +363,17 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar gasPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, Scalar density)
{
Scalar result = interpolateGasTRho_(gasPressure_,
temperature,
density);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTRho_(gasPressure_,
temperature,
density);
if (std::isnan(Toolbox::value(result)))
return RawComponent::gasPressure(temperature,
density);
}
return result;
}
@ -375,15 +383,17 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
static Scalar liquidPressure(Scalar temperature, Scalar density)
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, Scalar density)
{
Scalar result = interpolateLiquidTRho_(liquidPressure_,
temperature,
density);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTRho_(liquidPressure_,
temperature,
density);
if (std::isnan(Toolbox::value(result)))
return RawComponent::liquidPressure(temperature,
density);
}
return result;
}
@ -413,14 +423,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateGasTP_(gasDensity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasDensity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::gasDensity(temperature, pressure);
}
return result;
}
@ -431,14 +443,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateLiquidTP_(liquidDensity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidDensity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::liquidDensity(temperature, pressure);
}
return result;
}
@ -448,14 +462,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateGasTP_(gasViscosity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasViscosity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::gasViscosity(temperature, pressure);
}
return result;
}
@ -465,14 +481,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateLiquidTP_(liquidViscosity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidViscosity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::liquidViscosity(temperature, pressure);
}
return result;
}
@ -482,14 +500,16 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateGasTP_(gasThermalConductivity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasThermalConductivity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::gasThermalConductivity(temperature, pressure);
}
return result;
}
@ -499,26 +519,31 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar result = interpolateLiquidTP_(liquidThermalConductivity_,
temperature,
pressure);
if (std::isnan(result)) {
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidThermalConductivity_,
temperature,
pressure);
if (std::isnan(Toolbox::value(result)))
return RawComponent::liquidThermalConductivity(temperature, pressure);
}
return result;
}
private:
// returns an interpolated value depending on temperature
static Scalar interpolateT_(const Scalar *values, Scalar T)
template <class Evaluation>
static Evaluation interpolateT_(const Scalar *values, const Evaluation& T)
{
Scalar alphaT = tempIdx_(T);
typedef Opm::MathToolbox<Evaluation> Toolbox;
Evaluation alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1)
return std::numeric_limits<Scalar>::quiet_NaN();
unsigned iT = (unsigned) alphaT;
unsigned iT = (unsigned) Toolbox::value(alphaT);
alphaT -= iT;
return
@ -528,34 +553,37 @@ private:
// returns an interpolated value for liquid depending on
// temperature and pressure
static Scalar interpolateLiquidTP_(const Scalar *values, Scalar T, Scalar p)
template <class Evaluation>
static Evaluation interpolateLiquidTP_(const Scalar *values, const Evaluation& T, const Evaluation& p)
{
Scalar alphaT = tempIdx_(T);
typedef MathToolbox<Evaluation> Toolbox;
Evaluation alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1) {
return std::numeric_limits<Scalar>::quiet_NaN();
return Toolbox::createConstant(std::numeric_limits<Scalar>::quiet_NaN());
}
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, Toolbox::value(alphaT)));
alphaT -= iT;
Scalar alphaP1 = pressLiquidIdx_(p, iT);
Scalar alphaP2 = pressLiquidIdx_(p, iT + 1);
Evaluation alphaP1 = pressLiquidIdx_(p, iT);
Evaluation alphaP2 = pressLiquidIdx_(p, iT + 1);
unsigned iP1 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP1));
unsigned iP2 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP2));
unsigned iP1 = std::max<int>(0, std::min<int>(nPress_ - 2, Toolbox::value(alphaP1)));
unsigned iP2 = std::max<int>(0, std::min<int>(nPress_ - 2, Toolbox::value(alphaP2)));
alphaP1 -= iP1;
alphaP2 -= iP2;
#if 0 && !defined NDEBUG
if(!(0 <= alphaT && alphaT <= 1.0))
OPM_THROW(NumericalIssue, "Temperature out of range: "
<< "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
<< "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
if(!(0 <= alphaP1 && alphaP1 <= 1.0))
OPM_THROW(NumericalIssue, "First liquid pressure out of range: "
<< "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T)) << ", " << maxLiquidPressure_(tempIdx_(T)) << "]");
<< "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T)) << ", " << maxLiquidPressure_(tempIdx_(T)) << "]");
if(!(0 <= alphaP2 && alphaP2 <= 1.0))
OPM_THROW(NumericalIssue, "Second liquid pressure out of range: "
<< "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T) + 1) << ", " << maxLiquidPressure_(tempIdx_(T) + 1) << "]");
<< "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T) + 1) << ", " << maxLiquidPressure_(tempIdx_(T) + 1) << "]");
#endif
return
@ -567,33 +595,36 @@ private:
// returns an interpolated value for gas depending on
// temperature and pressure
static Scalar interpolateGasTP_(const Scalar *values, Scalar T, Scalar p)
template <class Evaluation>
static Evaluation interpolateGasTP_(const Scalar *values, const Evaluation& T, const Evaluation& p)
{
Scalar alphaT = tempIdx_(T);
typedef MathToolbox<Evaluation> Toolbox;
Evaluation alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1) {
return std::numeric_limits<Scalar>::quiet_NaN();
return Toolbox::createConstant(std::numeric_limits<Scalar>::quiet_NaN());
}
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, Toolbox::value(alphaT)));
alphaT -= iT;
Scalar alphaP1 = pressGasIdx_(p, iT);
Scalar alphaP2 = pressGasIdx_(p, iT + 1);
unsigned iP1 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP1));
unsigned iP2 = std::max<int>(0, std::min<int>(nPress_ - 2, (int) alphaP2));
Evaluation alphaP1 = pressGasIdx_(p, iT);
Evaluation alphaP2 = pressGasIdx_(p, iT + 1);
unsigned iP1 = std::max<int>(0, std::min<int>(nPress_ - 2, Toolbox::value(alphaP1)));
unsigned iP2 = std::max<int>(0, std::min<int>(nPress_ - 2, Toolbox::value(alphaP2)));
alphaP1 -= iP1;
alphaP2 -= iP2;
#if 0 && !defined NDEBUG
if(!(0 <= alphaT && alphaT <= 1.0))
OPM_THROW(NumericalIssue, "Temperature out of range: "
<< "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
<< "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
if(!(0 <= alphaP1 && alphaP1 <= 1.0))
OPM_THROW(NumericalIssue, "First gas pressure out of range: "
<< "p=" << p << " range: [" << minGasPressure_(tempIdx_(T)) << ", " << maxGasPressure_(tempIdx_(T)) << "]");
<< "p=" << p << " range: [" << minGasPressure_(tempIdx_(T)) << ", " << maxGasPressure_(tempIdx_(T)) << "]");
if(!(0 <= alphaP2 && alphaP2 <= 1.0))
OPM_THROW(NumericalIssue, "Second gas pressure out of range: "
<< "p=" << p << " range: [" << minGasPressure_(tempIdx_(T) + 1) << ", " << maxGasPressure_(tempIdx_(T) + 1) << "]");
<< "p=" << p << " range: [" << minGasPressure_(tempIdx_(T) + 1) << ", " << maxGasPressure_(tempIdx_(T) + 1) << "]");
#endif
return
@ -605,14 +636,15 @@ private:
// returns an interpolated value for gas depending on
// temperature and density
static Scalar interpolateGasTRho_(const Scalar *values, Scalar T, Scalar rho)
template <class Evaluation>
static Evaluation interpolateGasTRho_(const Scalar *values, const Evaluation& T, const Evaluation& rho)
{
Scalar alphaT = tempIdx_(T);
Evaluation alphaT = tempIdx_(T);
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
alphaT -= iT;
Scalar alphaP1 = densityGasIdx_(rho, iT);
Scalar alphaP2 = densityGasIdx_(rho, iT + 1);
Evaluation alphaP1 = densityGasIdx_(rho, iT);
Evaluation alphaP2 = densityGasIdx_(rho, iT + 1);
unsigned iP1 = std::max<int>(0, std::min<int>(nDensity_ - 2, (int) alphaP1));
unsigned iP2 = std::max<int>(0, std::min<int>(nDensity_ - 2, (int) alphaP2));
alphaP1 -= iP1;
@ -627,14 +659,15 @@ private:
// returns an interpolated value for liquid depending on
// temperature and density
static Scalar interpolateLiquidTRho_(const Scalar *values, Scalar T, Scalar rho)
template <class Evaluation>
static Evaluation interpolateLiquidTRho_(const Scalar *values, const Evaluation& T, const Evaluation& rho)
{
Scalar alphaT = tempIdx_(T);
Evaluation alphaT = tempIdx_(T);
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
alphaT -= iT;
Scalar alphaP1 = densityLiquidIdx_(rho, iT);
Scalar alphaP2 = densityLiquidIdx_(rho, iT + 1);
Evaluation alphaP1 = densityLiquidIdx_(rho, iT);
Evaluation alphaP2 = densityLiquidIdx_(rho, iT + 1);
unsigned iP1 = std::max<int>(0, std::min<int>(nDensity_ - 2, (int) alphaP1));
unsigned iP2 = std::max<int>(0, std::min<int>(nDensity_ - 2, (int) alphaP2));
alphaP1 -= iP1;
@ -649,13 +682,15 @@ private:
// returns the index of an entry in a temperature field
static Scalar tempIdx_(Scalar temperature)
template <class Evaluation>
static Evaluation tempIdx_(const Evaluation& temperature)
{
return (nTemp_ - 1)*(temperature - tempMin_)/(tempMax_ - tempMin_);
}
// returns the index of an entry in a pressure field
static Scalar pressLiquidIdx_(Scalar pressure, unsigned tempIdx)
template <class Evaluation>
static Evaluation pressLiquidIdx_(const Evaluation& pressure, unsigned tempIdx)
{
Scalar plMin = minLiquidPressure_(tempIdx);
Scalar plMax = maxLiquidPressure_(tempIdx);
@ -664,7 +699,8 @@ private:
}
// returns the index of an entry in a temperature field
static Scalar pressGasIdx_(Scalar pressure, unsigned tempIdx)
template <class Evaluation>
static Evaluation pressGasIdx_(const Evaluation& pressure, unsigned tempIdx)
{
Scalar pgMin = minGasPressure_(tempIdx);
Scalar pgMax = maxGasPressure_(tempIdx);
@ -673,7 +709,8 @@ private:
}
// returns the index of an entry in a density field
static Scalar densityLiquidIdx_(Scalar density, unsigned tempIdx)
template <class Evaluation>
static Evaluation densityLiquidIdx_(const Evaluation& density, unsigned tempIdx)
{
Scalar densityMin = minLiquidDensity_(tempIdx);
Scalar densityMax = maxLiquidDensity_(tempIdx);
@ -681,7 +718,8 @@ private:
}
// returns the index of an entry in a density field
static Scalar densityGasIdx_(Scalar density, unsigned tempIdx)
template <class Evaluation>
static Evaluation densityGasIdx_(const Evaluation& density, unsigned tempIdx)
{
Scalar densityMin = minGasDensity_(tempIdx);
Scalar densityMax = maxGasDensity_(tempIdx);

42
opm/material/components/Unit.hpp
View File

@ -79,10 +79,12 @@ public:
{ return 1.0; }
/*!
* \copydoc Component::liquidIsCompressible
* \copydoc Component::vaporPressure
*/
static Scalar vaporPressure(Scalar T)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{ return 1.0; }
/*!
* \copydoc Component::liquidIsCompressible
*/
@ -104,74 +106,86 @@ public:
/*!
* \copydoc Component::liquidDensity
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::liquidViscosity
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::gasDensity
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::gasViscosity
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::gasEnthalpy
*/
static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::liquidEnthalpy
*/
static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::gasInternalEnergy
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::liquidInternalEnergy
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidInternalEnergy(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::gasThermalConductivity
*/
static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::liquidThermalConductivity
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::gasHeatCapacity
*/
static Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
/*!
* \copydoc Component::liquidHeatCapacity
*/
static Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{ return 1.0; }
};

164
opm/material/components/Xylene.hpp
View File

@ -25,14 +25,15 @@
#ifndef OPM_XYLENE_HPP
#define OPM_XYLENE_HPP
#include <cmath>
#include <opm/material/IdealGas.hpp>
#include <opm/material/components/Component.hpp>
#include <opm/material/Constants.hpp>
#include <opm/material/common/MathToolbox.hpp>
namespace Opm
{
#include <cmath>
namespace Opm {
/*!
* \ingroup Components
* \brief Component for Xylene
@ -43,6 +44,7 @@ template <class Scalar>
class Xylene : public Component<Scalar, Xylene<Scalar> >
{
typedef Constants<Scalar> Consts;
typedef Opm::IdealGas<Scalar> IdealGas;
public:
/*!
@ -72,13 +74,13 @@ public:
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at xylene's triple point.
*/
static const Scalar tripleTemperature()
static Scalar tripleTemperature()
{ OPM_THROW(std::runtime_error, "Not implemented: tripleTemperature for xylene"); }
/*!
* \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at xylene's triple point.
*/
static const Scalar triplePressure()
static Scalar triplePressure()
{ OPM_THROW(std::runtime_error, "Not implemented: triplePressure for xylene"); }
/*!
@ -87,17 +89,16 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar temperature)
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar A = 7.00909;;
const Scalar B = 1462.266;;
const Scalar C = 215.110;;
Scalar T = temperature - 273.15;
Scalar psat = 1.334*std::pow(10.0, (A - (B/(T + C)))); // in [mbar]
psat *= 100.0; // in [Pa] (0.001*1.E5)
return psat;
return 100*1.334*Toolbox::pow(10.0, (A - (B/(temperature - 273.15 + C))));
}
/*!
@ -105,35 +106,36 @@ public:
*
* source : Reid et al. (fourth edition): Missenard group contrib. method (chap 5-7, Table 5-11, s. example 5-8)
*
* \param temp temperature of component in \f$\mathrm{[K]}\f$
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar spHeatCapLiquidPhase(Scalar temp, Scalar pressure)
template <class Evaluation>
static Evaluation spHeatCapLiquidPhase(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar CH3,C6H5,H;
Evaluation CH3,C6H5,H;
// after Reid et al. : Missenard group contrib. method (s. example 5-8)
// Xylene: C9H12 : 3* CH3 ; 1* C6H5 (phenyl-ring) ; -2* H (this was too much!)
// linear interpolation between table values [J/(mol K)]
if(temp < 298.0){ // take care: extrapolation for Temp<273
H = 13.4 + 1.2*(temp - 273.0)/25.0; // 13.4 + 1.2 = 14.6 = H(T=298K) i.e. interpolation of table values 273<T<298
CH3 = 40.0 + 1.6*(temp - 273.0)/25.0; // 40 + 1.6 = 41.6 = CH3(T=298K)
C6H5 = 113.0 + 4.2*(temp - 273.0)/25.0; // 113 + 4.2 = 117.2 = C6H5(T=298K)
if(temperature < 298.0){ // take care: extrapolation for Temperature<273
H = 13.4 + 1.2*(temperature - 273.0)/25.0; // 13.4 + 1.2 = 14.6 = H(T=298K) i.e. interpolation of table values 273<T<298
CH3 = 40.0 + 1.6*(temperature - 273.0)/25.0; // 40 + 1.6 = 41.6 = CH3(T=298K)
C6H5 = 113.0 + 4.2*(temperature - 273.0)/25.0; // 113 + 4.2 = 117.2 = C6H5(T=298K)
}
else if(temp < 323.0){
H = 14.6 + 0.9*(temp - 298.0)/25.0; // i.e. interpolation of table values 298<T<323
CH3 = 41.6 + 1.9*(temp - 298.0)/25.0;
C6H5 = 117.2 + 6.2*(temp - 298.0)/25.0;
else if(temperature < 323.0){
H = 14.6 + 0.9*(temperature - 298.0)/25.0; // i.e. interpolation of table values 298<T<323
CH3 = 41.6 + 1.9*(temperature - 298.0)/25.0;
C6H5 = 117.2 + 6.2*(temperature - 298.0)/25.0;
}
else if(temp < 348.0){
H = 15.5 + 1.2*(temp - 323.0)/25.0; // i.e. interpolation of table values 323<T<348
CH3 = 43.5 + 2.3*(temp - 323.0)/25.0;
C6H5 = 123.4 + 6.3*(temp - 323.0)/25.0;
else if(temperature < 348.0){
H = 15.5 + 1.2*(temperature - 323.0)/25.0; // i.e. interpolation of table values 323<T<348
CH3 = 43.5 + 2.3*(temperature - 323.0)/25.0;
C6H5 = 123.4 + 6.3*(temperature - 323.0)/25.0;
}
else {
H = 16.7 + 2.1*(temp - 348.0)/25.0; // i.e. interpolation of table values 348<T<373
CH3 = 45.8 + 2.5*(temp - 348.0)/25.0; // take care: extrapolation for Temp>373
C6H5 = 129.7 + 6.3*(temp - 348.0)/25.0; // most likely leads to underestimation
H = 16.7 + 2.1*(temperature - 348.0)/25.0; // i.e. interpolation of table values 348<T<373
CH3 = 45.8 + 2.5*(temperature - 348.0)/25.0; // take care: extrapolation for Temperature>373
C6H5 = 129.7 + 6.3*(temperature - 348.0)/25.0; // most likely leads to underestimation
}
return (C6H5 + 2*CH3 - H)/molarMass();// J/(mol K) -> J/(kg K)
@ -143,7 +145,8 @@ public:
/*!
* \copydoc Component::liquidEnthalpy
*/
static Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
// Gauss quadrature rule:
// Interval: [0K; temperature (K)]
@ -157,10 +160,10 @@ public:
// enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input
return 0.5*temperature*(spHeatCapLiquidPhase(0.2113*temperature,pressure)
+ spHeatCapLiquidPhase(0.7887*temperature,pressure));
+ spHeatCapLiquidPhase(0.7887*temperature,pressure));
}
/*!
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at xylene's boiling point (1 atm).
*/
static Scalar boilingTemperature()
@ -174,11 +177,14 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar heatVap(Scalar temperature,
Scalar pressure)
template <class Evaluation>
static Evaluation heatVap(Evaluation temperature,
const Evaluation& pressure)
{
temperature = std::min(temperature, criticalTemperature()); // regularization
temperature = std::max(temperature, 0.0); // regularization
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, criticalTemperature()); // regularization
temperature = Toolbox::max(temperature, 0.0); // regularization
const Scalar T_crit = criticalTemperature();
const Scalar Tr1 = boilingTemperature()/criticalTemperature();
@ -189,10 +195,10 @@ public:
* (3.978 * Tr1 - 3.958 + 1.555*std::log(p_crit * 1e-5 /*Pa->bar*/ ) )
/ (1.07 - Tr1); /* [J/mol] */
/* Variation with temp according to Watson relation eq 7-12.1*/
const Scalar Tr2 = temperature/criticalTemperature();
/* Variation with temperature according to Watson relation eq 7-12.1*/
const Evaluation& Tr2 = temperature/criticalTemperature();
const Scalar n = 0.375;
const Scalar DH_vap = DH_v_boil * std::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
const Evaluation& DH_vap = DH_v_boil * Toolbox::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
return (DH_vap/molarMass()); // we need [J/kg]
}
@ -203,7 +209,8 @@ public:
* The relation used here is true on the vapor pressure curve, i.e. as long
* as there is a liquid phase present.
*/
static Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
return liquidEnthalpy(temperature, pressure) + heatVap(temperature, pressure);
}
@ -211,11 +218,10 @@ public:
/*!
* \copydoc Component::gasDensity
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
return IdealGas<Scalar>::density(molarMass(),
temperature,
pressure);
return IdealGas::density(Evaluation(molarMass()), temperature, pressure);
}
/*!
@ -223,8 +229,9 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar molarGasDensity(Scalar temperature, Scalar pressure)
*/
template <class Evaluation>
static Evaluation molarGasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
return gasDensity(temperature, pressure) / molarMass();
}
@ -235,32 +242,32 @@ public:
*
* source : Reid et al. (fourth edition): Modified Racket technique (chap. 3-11, eq. 3-11.9)
*
* \param temp temperature of component in \f$\mathrm{[K]}\f$
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar molarLiquidDensity(Scalar temp, Scalar pressure)
template <class Evaluation>
static Evaluation molarLiquidDensity(Evaluation temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
// saturated molar volume according to Lide, CRC Handbook of
// Thermophysical and Thermochemical Data, CRC Press, 1994
// valid for 245 < Temp < 600
temp = std::min(temp, 500.0); // regularization
temp = std::max(temp, 250.0); // regularization
// valid for 245 < Temperature < 600
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0); // regularization
const Scalar A1 = 0.25919; // from table
const Scalar A2 = 0.0014569; // from table
const Scalar expo = 1.0 + std::pow((1.0 - temp/criticalTemperature()), (2.0/7.0));
const Scalar V = A2*std::pow(A1, expo); // liquid molar volume [m^3/mol]
return 1.0/V; // molar density [mol/m^3]
const Evaluation& expo = 1.0 + Toolbox::pow((1.0 - temperature/criticalTemperature()), (2.0/7.0));
return 1.0/(A2*Toolbox::pow(A1, expo));
}
/*!
* \copydoc Component::liquidDensity
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{
return molarLiquidDensity(temperature, pressure)*molarMass(); // [kg/m^3]
}
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{ return molarLiquidDensity(temperature, pressure)*molarMass(); }
/*!
* \copydoc Component::gasIsCompressible
@ -283,41 +290,44 @@ public:
/*!
* \copydoc Component::gasViscosity
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gasViscosity(Evaluation temperature, const Evaluation& pressure)
{
temperature = std::min(temperature, 500.0); // regularization
temperature = std::max(temperature, 250.0); // regularization
typedef MathToolbox<Evaluation> Toolbox;
const Scalar Tr = std::max(temperature/criticalTemperature(), 1e-10);
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0); // regularization
const Evaluation& Tr = Toolbox::max(temperature/criticalTemperature(), 1e-10);
const Scalar Fp0 = 1.0;
const Scalar xi = 0.004623;
const Scalar eta_xi = Fp0*(0.807*std::pow(Tr, 0.618)
- 0.357*std::exp(-0.449*Tr)
+ 0.34*std::exp(-4.058*Tr)
const Evaluation& eta_xi = Fp0*(0.807*Toolbox::pow(Tr, 0.618)
- 0.357*Toolbox::exp(-0.449*Tr)
+ 0.34*Toolbox::exp(-4.058*Tr)
+ 0.018);
Scalar r = eta_xi/xi; // [1e-6 P]
r /= 1.0e7; // [Pa s]
return r;
return eta_xi/xi / 1e7; // [Pa s]
}
/*!
* \copydoc Component::liquidViscosity
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation liquidViscosity(Evaluation temperature, const Evaluation& pressure)
{
temperature = std::min(temperature, 500.0); // regularization
temperature = std::max(temperature, 250.0); // regularization
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0); // regularization
const Scalar A = -3.82;
const Scalar B = 1027.0;
const Scalar C = -6.38e-4;
const Scalar D = 4.52e-7;
Scalar r = std::exp(A + B/temperature + C*temperature + D*temperature*temperature); // in [cP]
r *= 1.0e-3; // in [Pa s]
return r; // [Pa s]
return 1e-3*Toolbox::exp(A
+ B/temperature
+ C*temperature
+ D*temperature*temperature); // in [cP]
}
};

60
opm/material/components/iapws/Common.hpp
View File

@ -26,6 +26,7 @@
#define OPM_IAPWS_COMMON_HPP
#include <opm/material/Constants.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <cmath>
@ -92,10 +93,13 @@ public:
* IAPWS: "Release on the IAPWS Formulation 2008 for the Viscosity
* of Ordinary Water Substance", http://www.iapws.org/relguide/visc.pdf
*/
static Scalar viscosity(Scalar temperature, Scalar rho)
template <class Evaluation>
static Evaluation viscosity(const Evaluation& temperature, const Evaluation& rho)
{
Scalar rhoBar = rho/322.0;
Scalar TBar = temperature/criticalTemperature;
typedef MathToolbox<Evaluation> Toolbox;
Evaluation rhoBar = rho/322.0;
Evaluation TBar = temperature/criticalTemperature;
// muBar = muBar_1
const Scalar Hij[6][7] = {
@ -107,8 +111,8 @@ public:
{ 0, 1.20573e-1, 0, 0, 0, 0,-5.93264e-4 }
};
Scalar tmp, tmp2, tmp3 = 1;
Scalar muBar = 0;
Evaluation tmp, tmp2, tmp3 = 1;
Evaluation muBar = 0;
for (int i = 0; i <= 5; ++i) {
tmp = 0;
tmp2 = 1;
@ -120,10 +124,10 @@ public:
tmp3 *= 1.0/TBar - 1;
};
muBar *= rhoBar;
muBar = std::exp(muBar);
muBar = Toolbox::exp(muBar);
// muBar *= muBar_0
muBar *= 100*std::sqrt(TBar);
muBar *= 100*Toolbox::sqrt(TBar);
const Scalar H[4] = {
1.67752, 2.20462, 0.6366564, -0.241605
};
@ -151,8 +155,11 @@ public:
* \param T absolute temperature in K
* \param rho density of water in kg/m^3
*/
static Scalar thermalConductivityIAPWS(Scalar T, Scalar rho)
template <class Evaluation>
static Evaluation thermalConductivityIAPWS(const Evaluation& T, const Evaluation& rho)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar thcond_tstar = 647.26 ;
static const Scalar thcond_rhostar = 317.7 ;
/*static const Scalar thcond_kstar = 1.0 ;*/
@ -182,13 +189,13 @@ public:
,-0.00422464
};
Scalar Tbar = T / thcond_tstar;
Scalar rhobar = rho / thcond_rhostar;
Evaluation Tbar = T / thcond_tstar;
Evaluation rhobar = rho / thcond_rhostar;
/* fast implementation... minimised calls to 'pow' routine... */
Scalar Troot = std::sqrt(Tbar);
Scalar Tpow = Troot;
Scalar lam = 0;
Evaluation Troot = Toolbox::sqrt(Tbar);
Evaluation Tpow = Troot;
Evaluation lam = 0;
for(int k = 0; k < thcond_a_count; ++k) {
lam += thcond_a[k] * Tpow;
@ -198,29 +205,28 @@ public:
lam +=
thcond_b0 + thcond_b1
* rhobar + thcond_b2
* std::exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
* Toolbox::exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
Scalar DTbar = std::abs(Tbar - 1) + thcond_c4;
Scalar DTbarpow = std::pow(DTbar, 3./5);
Scalar Q = 2. + thcond_c5 / DTbarpow;
Evaluation DTbar = Toolbox::abs(Tbar - 1) + thcond_c4;
Evaluation DTbarpow = Toolbox::pow(DTbar, 3./5);
Evaluation Q = 2. + thcond_c5 / DTbarpow;
Scalar S;
if(Tbar >= 1){
Evaluation S;
if(Tbar >= 1)
S = 1. / DTbar;
}else{
else
S = thcond_c6 / DTbarpow;
}
Scalar rhobar18 = std::pow(rhobar, 1.8);
Scalar rhobarQ = std::pow(rhobar, Q);
Evaluation rhobar18 = Toolbox::pow(rhobar, 1.8);
Evaluation rhobarQ = Toolbox::pow(rhobar, Q);
lam +=
(thcond_d1 / std::pow(Tbar,10) + thcond_d2) * rhobar18 *
std::exp(thcond_c1 * (1 - rhobar * rhobar18))
(thcond_d1 / Toolbox::pow(Tbar,10.0) + thcond_d2) * rhobar18 *
Toolbox::exp(thcond_c1 * (1 - rhobar * rhobar18))
+ thcond_d3 * S * rhobarQ *
std::exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
Toolbox::exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
+ thcond_d4 *
std::exp(thcond_c2 * std::pow(Troot,3) + thcond_c3 / std::pow(rhobar,5));
Toolbox::exp(thcond_c2 * Toolbox::pow(Troot,3.0) + thcond_c3 / Toolbox::pow(rhobar,5.0));
return /*thcond_kstar * */ lam;
}
};

127
opm/material/components/iapws/Region1.hpp
View File

@ -24,6 +24,8 @@
#ifndef OPM_IAPWS_REGION1_HPP
#define OPM_IAPWS_REGION1_HPP
#include <opm/material/common/MathToolbox.hpp>
#include <cmath>
namespace Opm {
@ -52,21 +54,24 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static bool isValid(Scalar temperature, Scalar pressure)
template <class Evaluation>
static bool isValid(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return
temperature <= 623.15 &&
pressure <= 100e6;
Toolbox::value(temperature) <= 623.15 &&
Toolbox::value(pressure) <= 100e6;
// actually this is:
/*
return
(
273.15 <= temperature &&
temperature <= 623.15 &&
pressure >= vaporPressure(temperature) &&
pressure <= 100e6
);
return
(
273.15 <= temperature &&
temperature <= 623.15 &&
pressure >= vaporPressure(temperature) &&
pressure <= 100e6
);
*/
}
@ -75,7 +80,8 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar tau(Scalar temperature)
template <class Evaluation>
static Evaluation tau(const Evaluation& temperature)
{ return 1386.0 / temperature; }
/*!
@ -84,7 +90,8 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar dtau_dT(Scalar temperature)
template <class Evaluation>
static Evaluation dtau_dT(const Evaluation& temperature)
{ return - 1386.0 / (temperature*temperature); }
/*!
@ -92,7 +99,8 @@ public:
*
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar pi(Scalar pressure)
template <class Evaluation>
static Evaluation pi(const Evaluation& pressure)
{ return pressure / 16.53e6; }
/*!
@ -113,7 +121,6 @@ public:
static Scalar dp_dpi(Scalar pressure)
{ return 16.53e6; }
/*!
* \brief The Gibbs free energy (dimensionless) for IAPWS region 1 (i.e. liquid)
*
@ -124,14 +131,17 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar gamma(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gamma(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result = 0;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = 0;
for (int i = 0; i < 34; ++i) {
result += n(i)*std::pow(7.1 - pi_, I(i))*std::pow(tau_ - 1.222, J(i));
result += n(i)*Toolbox::pow(7.1 - pi_, I(i))*Toolbox::pow(tau_ - 1.222, J(i));
}
return result;
@ -149,17 +159,20 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar dgamma_dtau(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation dgamma_dtau(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result = 0.0;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 34; i++) {
result +=
n(i) *
std::pow(7.1 - pi_, I(i)) *
std::pow(tau_ - 1.222, J(i)-1) *
Toolbox::pow(7.1 - pi_, static_cast<Scalar>(I(i))) *
Toolbox::pow(tau_ - 1.222, static_cast<Scalar>(J(i)-1)) *
J(i);
}
@ -177,18 +190,21 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar dgamma_dpi(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation dgamma_dpi(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result = 0.0;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 34; i++) {
result +=
-n(i) *
I(i) *
std::pow(7.1 - pi_, I(i) - 1) *
std::pow(tau_ - 1.222, J(i));
Toolbox::pow(7.1 - pi_, static_cast<Scalar>(I(i) - 1)) *
Toolbox::pow(tau_ - 1.222, static_cast<Scalar>(J(i)));
}
return result;
@ -206,19 +222,22 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar ddgamma_dtaudpi(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation ddgamma_dtaudpi(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result = 0.0;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 34; i++) {
result +=
-n(i) *
I(i) *
J(i) *
std::pow(7.1 - pi_, I(i) - 1) *
std::pow(tau_ - 1.222, J(i) - 1);
Toolbox::pow(7.1 - pi_, static_cast<Scalar>(I(i) - 1)) *
Toolbox::pow(tau_ - 1.222, static_cast<Scalar>(J(i) - 1));
}
return result;
@ -236,19 +255,22 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar ddgamma_ddpi(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation ddgamma_ddpi(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result = 0.0;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 34; i++) {
result +=
n(i) *
I(i) *
(I(i) - 1) *
std::pow(7.1 - pi_, I(i) - 2) *
std::pow(tau_ - 1.222, J(i));
Toolbox::pow(7.1 - pi_, I(i) - 2) *
Toolbox::pow(tau_ - 1.222, J(i));
}
return result;
@ -265,19 +287,22 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar ddgamma_ddtau(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation ddgamma_ddtau(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result = 0.0;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 34; i++) {
result +=
n(i) *
std::pow(7.1 - pi_, I(i)) *
Toolbox::pow(7.1 - pi_, I(i)) *
J(i) *
(J(i) - 1) *
std::pow(tau_ - 1.222, J(i) - 2);
Toolbox::pow(tau_ - 1.222, J(i) - 2);
}
return result;
@ -303,7 +328,7 @@ private:
return n[i];
}
static short int I(int i)
static Scalar I(int i)
{
static const short int I[34] = {
0, 0, 0,
@ -322,7 +347,7 @@ private:
return I[i];
}
static short int J(int i)
static Scalar J(int i)
{
static const short int J[34] = {
-2, -1, 0,

130
opm/material/components/iapws/Region2.hpp
View File

@ -24,6 +24,8 @@
#ifndef OPM_IAPWS_REGION2_HPP
#define OPM_IAPWS_REGION2_HPP
#include <opm/material/common/MathToolbox.hpp>
#include <cmath>
namespace Opm {
@ -53,17 +55,18 @@ public:
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static bool isValid(Scalar temperature, Scalar pressure)
template <class Evaluation>
static bool isValid(const Evaluation& temperature, const Evaluation& pressure)
{
return
temperature <= 623.15 && pressure <= 100e6;
// actually this is:
/*
return
(273.15 <= temperature && temperature <= 623.15 && pressure <= vaporPressure(temperature)) ||
(623.15 < temperature && temperature <= 863.15 && pressure <= auxPressure(temperature)) ||
(863.15 < temperature && temperature <= 1073.15 && pressure < 100e6);
return
(273.15 <= temperature && temperature <= 623.15 && pressure <= vaporPressure(temperature)) ||
(623.15 < temperature && temperature <= 863.15 && pressure <= auxPressure(temperature)) ||
(863.15 < temperature && temperature <= 1073.15 && pressure < 100e6);
*/
}
@ -72,7 +75,8 @@ public:
*
* \param temperature temperature of component
*/
static Scalar tau(Scalar temperature)
template <class Evaluation>
static Evaluation tau(const Evaluation& temperature)
{ return 540.0 / temperature; }
/*!
@ -81,7 +85,8 @@ public:
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*/
static Scalar dtau_dT(Scalar temperature)
template <class Evaluation>
static Evaluation dtau_dT(const Evaluation& temperature)
{ return - 540.0 / (temperature*temperature); }
/*!
@ -89,7 +94,8 @@ public:
*
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar pi(Scalar pressure)
template <class Evaluation>
static Evaluation pi(const Evaluation& pressure)
{ return pressure / 1e6; }
/*!
@ -98,7 +104,8 @@ public:
*
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar dpi_dp(Scalar pressure)
template <class Evaluation>
static Scalar dpi_dp(const Evaluation& pressure)
{ return 1.0 / 1e6; }
/*!
@ -107,7 +114,8 @@ public:
*
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar dp_dpi(Scalar pressure)
template <class Evaluation>
static Evaluation dp_dpi(const Evaluation& pressure)
{ return 1e6; }
/*!
@ -121,24 +129,27 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar gamma(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation gamma(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau = tau(temperature); /* reduced temperature */
Scalar pi = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
Scalar result;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
Evaluation result;
// ideal gas part
result = ln(pi);
result = Toolbox::ln(pi_);
for (int i = 0; i < 9; ++i)
result += n_g(i)*std::pow(tau, J_g(i));
result += n_g(i)*Toolbox::pow(tau_, J_g(i));
// residual part
for (int i = 0; i < 43; ++i)
result +=
n_r(i)*
std::pow(pi, I_r(i))*
std::pow(tau - 0.5, J_r(i));
Toolbox::pow(pi_, I_r(i))*
Toolbox::pow(tau_ - 0.5, J_r(i));
return result;
}
@ -154,27 +165,30 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar dgamma_dtau(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation dgamma_dtau(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Scalar result = 0;
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 9; i++) {
result +=
n_g(i) *
J_g(i) *
std::pow(tau_, J_g(i) - 1);
Toolbox::pow(tau_, static_cast<Scalar>(J_g(i) - 1));
}
// residual part
for (int i = 0; i < 43; i++) {
result +=
n_r(i) *
std::pow(pi_, I_r(i)) *
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i))) *
J_r(i) *
std::pow(tau_ - 0.5, J_r(i) - 1);
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 1));
}
return result;
@ -192,21 +206,24 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar dgamma_dpi(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation dgamma_dpi(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Scalar result = 1/pi_;
Evaluation result = 1/pi_;
// residual part
for (int i = 0; i < 43; i++) {
result +=
n_r(i) *
I_r(i) *
std::pow(pi_, I_r(i) - 1) *
std::pow(tau_ - 0.5, J_r(i));
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i) - 1)) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i)));
}
return result;
@ -224,13 +241,16 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar ddgamma_dtaudpi(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation ddgamma_dtaudpi(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Scalar result = 0;
Evaluation result = Toolbox::createConstant(0.0);
// residual part
for (int i = 0; i < 43; i++) {
@ -238,8 +258,8 @@ public:
n_r(i) *
I_r(i) *
J_r(i) *
std::pow(pi_, I_r(i) - 1) *
std::pow(tau_ - 0.5, J_r(i) - 1);
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i) - 1)) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 1));
}
return result;
@ -257,13 +277,16 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar ddgamma_ddpi(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation ddgamma_ddpi(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Scalar result = -1/(pi_*pi_);
Evaluation result = -1/(pi_*pi_);
// residual part
for (int i = 0; i < 43; i++) {
@ -271,8 +294,8 @@ public:
n_r(i) *
I_r(i) *
(I_r(i) - 1) *
std::pow(pi_, I_r(i) - 2) *
std::pow(tau_ - 0.5, J_r(i));
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i) - 2)) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i)));
}
return result;
@ -290,29 +313,32 @@ public:
* 1997 for the Thermodynamic Properties of Water and Steam",
* http://www.iapws.org/relguide/IF97-Rev.pdf
*/
static Scalar ddgamma_ddtau(Scalar temperature, Scalar pressure)
template <class Evaluation>
static Evaluation ddgamma_ddtau(const Evaluation& temperature, const Evaluation& pressure)
{
Scalar tau_ = tau(temperature); /* reduced temperature */
Scalar pi_ = pi(pressure); /* reduced pressure */
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Scalar result = 0;
Evaluation result = Toolbox::createConstant(0.0);
for (int i = 0; i < 9; i++) {
result +=
n_g(i) *
J_g(i) *
(J_g(i) - 1) *
std::pow(tau_, J_g(i) - 2);
Toolbox::pow(tau_, static_cast<Scalar>(J_g(i) - 2));
}
// residual part
for (int i = 0; i < 43; i++) {
result +=
n_r(i) *
std::pow(pi_, I_r(i)) *
Toolbox::pow(pi_, I_r(i)) *
J_r(i) *
(J_r(i) - 1.) *
std::pow(tau_ - 0.5, J_r(i) - 2.);
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 2));
}
return result;
@ -339,10 +365,10 @@ private:
-0.26674547914087e-4, 0.20481737692309e-7, 0.43870667284435e-6,
-0.32277677238570e-4, -0.15033924542148e-2, -0.40668253562649e-1,
-0.78847309559367e-9, 0.12790717852285e-7, 0.48225372718507e-6,
0.22922076337661e-5, -0.16714766451061e-10, -0.21171472321355e-2,
0.22922076337661e-5, -0.16714766451061e-10, -0.21171472321355e-2,
-0.23895741934104e2, -0.59059564324270e-17, -0.12621808899101e-5,
-0.38946842435739e-1, 0.11256211360459e-10, -0.82311340897998e1,
0.19809712802088e-7, 0.10406965210174e-18, -0.10234747095929e-12,
0.19809712802088e-7, 0.10406965210174e-18, -0.10234747095929e-12,
-0.10018179379511e-8, -0.80882908646985e-10, 0.10693031879409,
-0.33662250574171, 0.89185845355421e-24, 0.30629316876232e-12,
-0.42002467698208e-5, -0.59056029685639e-25, 0.37826947613457e-5,

40
opm/material/components/iapws/Region4.hpp
View File

@ -25,6 +25,8 @@
#ifndef OPM_IAPWS_REGION4_HPP
#define OPM_IAPWS_REGION4_HPP
#include <opm/material/common/MathToolbox.hpp>
#include <cmath>
namespace Opm {
@ -51,12 +53,15 @@ public:
* \brief Returns the saturation pressure in \f$\mathrm{[Pa]}\f$ of pure water at a given
* temperature.
*
*\param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
*
* The saturation pressure is often also called vapor pressure.
*/
static Scalar saturationPressure(Scalar temperature)
template <class Evaluation>
static Evaluation saturationPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar n[10] = {
0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
@ -64,13 +69,13 @@ public:
0.65017534844798e3
};
Scalar sigma = temperature + n[8]/(temperature - n[9]);
const Evaluation& sigma = temperature + n[8]/(temperature - n[9]);
Scalar A = (sigma + n[0])*sigma + n[1];
Scalar B = (n[2]*sigma + n[3])*sigma + n[4];
Scalar C = (n[5]*sigma + n[6])*sigma + n[7];
const Evaluation& A = (sigma + n[0])*sigma + n[1];
const Evaluation& B = (n[2]*sigma + n[3])*sigma + n[4];
const Evaluation& C = (n[5]*sigma + n[6])*sigma + n[7];
Scalar tmp = 2*C/(std::sqrt(B*B - 4*A*C) - B);
Evaluation tmp = 2*C/(Toolbox::sqrt(B*B - 4*A*C) - B);
tmp *= tmp;
tmp *= tmp;
@ -81,27 +86,30 @@ public:
* \brief Returns the saturation temperature in \f$\mathrm{[K]}\f$ of pure water at a given
* pressure.
*
*\param pressure pressure of component in \f$\mathrm{[Pa]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* The saturation pressure is often also called vapor pressure.
*/
static Scalar vaporTemperature(Scalar pressure)
template <class Evaluation>
static Evaluation vaporTemperature(const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar n[10] = {
0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
-0.48232657361591e4, 0.40511340542057e6, -0.23855557567849,
0.65017534844798e3
};
Scalar beta = pow((pressure/1e6 /*from Pa to MPa*/), (1./4.));
Scalar beta2 = pow(beta, 2.);
Scalar E = beta2 + n[2] * beta + n[5];
Scalar F = n[0]*beta2 + n[3]*beta + n[6];
Scalar G = n[1]*beta2 + n[4]*beta + n[7];
const Evaluation& beta = pow((pressure/1e6 /*from Pa to MPa*/), (1./4.));
const Evaluation& beta2 = pow(beta, 2.);
const Evaluation& E = beta2 + n[2] * beta + n[5];
const Evaluation& F = n[0]*beta2 + n[3]*beta + n[6];
const Evaluation& G = n[1]*beta2 + n[4]*beta + n[7];
Scalar D = ( 2.*G)/(-F -std::sqrt(pow(F,2.) - 4.*E*G));
const Evaluation& D = ( 2.*G)/(-F -Toolbox::sqrt(pow(F,2.) - 4.*E*G));
Scalar temperature = (n[9] + D - std::sqrt(pow(n[9]+D , 2.) - 4.* (n[8] + n[9]*D)) ) * 0.5;
const Evaluation& temperature = (n[9] + D - Toolbox::sqrt(pow(n[9]+D , 2.) - 4.* (n[8] + n[9]*D)) ) * 0.5;
return temperature;
}