Commit 48d30a37 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[material] allow adl for math functions in components

parent e77fc6f9
......@@ -137,14 +137,17 @@ public:
const Scalar omega = 0.078; // accentric factor
const Scalar M = molarMass() * 1e3; // molar mas [g/mol]
using std::pow;
using std::exp;
using std::cbrt;
Scalar Fc = 1 - 0.2756*omega;
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::cbrt(Vc * Vc) * Omega_v);
1.16145*pow(Tstar, -0.14874) +
0.52487*exp(- 0.77320*Tstar) +
2.16178*exp(- 2.43787*Tstar);
Scalar mu = 40.785 * Fc * sqrt(M * temperature)
/ (cbrt(Vc * Vc) * Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;
......@@ -158,8 +161,9 @@ public:
DUNE_THROW(NumericalProblem,
"simpleGasViscosity: Temperature out of range! (T = " << temperature << " K)");
}
return 1.496e-6 * std::sqrt(temperature * temperature * temperature) / (temperature + 120.);
using std::sqrt;
return 1.496e-6 * sqrt(temperature * temperature * temperature) / (temperature + 120.);
}
/*!
......@@ -216,22 +220,21 @@ public:
// scale temperature with reference temp of 100K
Scalar phi = temperature/100;
using std::pow;
Scalar 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);
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.201650E+00 * pow(phi,2)
-0.196930E-01 * pow(phi,3)
+0.106460E-02 * pow(phi,4)
-0.303284E-04 * pow(phi,5)
+0.355861E-06 * pow(phi,6);
c_p += -0.549169E+01 * pow(phi,-1)
+0.585171E+01 * pow(phi,-2)
-0.372865E+01 * pow(phi,-3)
+0.133981E+01 * pow(phi,-4)
-0.233758E+00 * pow(phi,-5)
+0.125718E-01 * pow(phi,-6);
c_p *= IdealGas::R / (molarMass() * 1000); // in J/mol/K * mol / kg / 1000 = kJ/kg/K
return c_p;
}
......
......@@ -151,7 +151,9 @@ public:
const Scalar salSat = f[0] + f[1]*theta + f[2]*theta*theta + f[3]*theta*theta*theta;
/*Regularization*/
salinity = std::min(std::max(salinity,0.0), salSat);
using std::min;
using std::max;
salinity = min(max(salinity,0.0), salSat);
const Scalar hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
......@@ -161,6 +163,7 @@ public:
const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
using std::pow;
Scalar d_h = 0;
for (int i = 0; i<=3; i++) {
for (int j=0; j<=2; j++) {
......@@ -281,25 +284,26 @@ public:
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure, Scalar salinity = constantSalinity)
{
using std::max;
const Scalar TempC = temperature - 273.15;
const Scalar pMPa = pressure/1.0E6;
salinity = std::max(0.0, salinity);
salinity = max(0.0, salinity);
const Scalar rhow = H2O::liquidDensity(temperature, pressure);
const Scalar density = rhow +
1000*salinity*(
0.668 +
0.44*salinity +
1.0E-6*(
300*pMPa -
2400*pMPa*salinity +
TempC*(
80.0 +
3*TempC -
3300*salinity -
13*pMPa +
47*pMPa*salinity)));
const Scalar density = rhow +
1000*salinity*(
0.668 +
0.44*salinity +
1.0E-6*(
300*pMPa -
2400*pMPa*salinity +
TempC*(
80.0 +
3*TempC -
3300*salinity -
13*pMPa +
47*pMPa*salinity)));
assert(density > 0.0);
return density;
}
......@@ -330,7 +334,9 @@ public:
const Scalar eps = pressure*1e-7;
Scalar deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
using std::abs;
for (int i = 0; i < 5 && abs(pressure*1e-9) < abs(deltaP); ++i) {
Scalar f = liquidDensity(temperature, pressure) - density;
Scalar df_dp;
......@@ -370,9 +376,12 @@ public:
static Scalar liquidViscosity(Scalar temperature, Scalar pressure, Scalar salinity = constantSalinity)
{
// regularisation
temperature = std::max(temperature, 275.0);
salinity = std::max(0.0, salinity);
using std::max;
temperature = max(temperature, 275.0);
salinity = max(0.0, salinity);
using std::pow;
using std::exp;
const Scalar T_C = temperature - 273.15;
const Scalar A = (0.42*pow((pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
const Scalar mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*exp(-A);
......
......@@ -152,7 +152,6 @@ public:
// calculate: \int_0^T c_p dT
return
1/molarMass()* // conversion from [J/mol] to [J/kg]
T*(cpVapA + T*
(cpVapB/2 + T*
(cpVapC/3 + T*
......@@ -204,17 +203,20 @@ public:
const Scalar M = molarMass() * 1e3; // molar mas [g/mol]
const Scalar dipole = 0.0; // dipole moment [debye]
Scalar mu_r4 = 131.3 * dipole / std::sqrt(Vc * Tc);
using std::sqrt;
Scalar mu_r4 = 131.3 * dipole / sqrt(Vc * Tc);
mu_r4 *= mu_r4;
mu_r4 *= mu_r4;
using std::exp;
using std::pow;
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);
1.16145*pow(Tstar, -0.14874) +
0.52487*exp(- 0.77320*Tstar) +
2.16178*exp(- 2.43787*Tstar);
Scalar mu = 40.785*Fc*sqrt(M*temperature)/(pow(Vc, 2./3)*Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;
......
......@@ -133,12 +133,14 @@ public:
// 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]);
}
using std::pow;
for (int i = 0; i < 5; ++i)
exponent += a[i]*pow(1 - Tred, t[i]);
exponent *= 1.0/Tred;
return std::exp(exponent)*criticalPressure();
using std::exp;
return exp(exponent)*criticalPressure();
}
/*!
......@@ -310,16 +312,19 @@ public:
TStar = temperature/ESP;
/* mu0: viscosity in zero-density limit */
using std::exp;
using std::log;
using std::sqrt;
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) );
mu0 = 1.00697*sqrt(temperature) / SigmaStar;
/* dmu : excess viscosity at elevated density */
rho = gasDensity(temperature, pressure); /* CO2 mass density [kg/m^3] */
using std::pow;
dmu = d11*rho + d21*rho*rho + d64*pow(rho,6)/(TStar*TStar*TStar)
+ d81*pow(rho,8) + d82*pow(rho,8)/TStar;
......
......@@ -122,7 +122,10 @@ protected:
if (Dune::FloatCmp::eq<Scalar>(temperature, maxTemp()))
return numTempSteps - 2;
const int result = static_cast<int>((temperature - minTemp())/(maxTemp() - minTemp())*(numTempSteps - 1));
return std::max(0, std::min(result, numTempSteps - 2));
using std::min;
using std::max;
return max(0, min(result, numTempSteps - 2));
}
int findPressIdx_(Scalar pressure) const
......@@ -130,7 +133,10 @@ protected:
if (Dune::FloatCmp::eq<Scalar>(pressure, maxPress()))
return numPressSteps - 2;
const int result = static_cast<int>((pressure - minPress())/(maxPress() - minPress())*(numPressSteps - 1));
return std::max(0, std::min(result, numPressSteps - 2));
using std::min;
using std::max;
return max(0, min(result, numPressSteps - 2));
}
Scalar temperatureAt_(int i) const
......
......@@ -103,7 +103,8 @@ public:
const Scalar B = 1.45838;
const Scalar C = -2.77580;
return 1e5 * std::exp(A - B/(temperature + C));
using std::exp;
return 1e5 * exp(A - B/(temperature + C));
}
/*!
......@@ -165,7 +166,6 @@ public:
// calculate: \int_0^T c_p dT
return
1/molarMass()* // conversion from [J/mol] to [J/kg]
T*(cpVapA + T*
(cpVapB/2 + T*
(cpVapC/3 + T*
......@@ -192,17 +192,20 @@ public:
const Scalar M = molarMass() * 1e3; // molar mas [g/mol]
const Scalar dipole = 0.0; // dipole moment [debye]
Scalar mu_r4 = 131.3 * dipole / std::sqrt(Vc * Tc);
using std::sqrt;
Scalar mu_r4 = 131.3 * dipole / sqrt(Vc * Tc);
mu_r4 *= mu_r4;
mu_r4 *= mu_r4;
using std::pow;
using std::exp;
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);
1.16145*pow(Tstar, -0.14874) +
0.52487*exp(- 0.77320*Tstar) +
2.16178*exp(- 2.43787*Tstar);
Scalar mu = 40.785*Fc*sqrt(M*temperature)/(pow(Vc, 2./3)*Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;
......
......@@ -130,8 +130,10 @@ public:
*/
static Scalar vaporPressure(Scalar T)
{
T = std::min(T, criticalTemperature());
T = std::max(T,tripleTemperature());
using std::min;
using std::max;
T = min(T, criticalTemperature());
T = max(T,tripleTemperature());
return Region4::saturationPressure(T);
}
......@@ -150,8 +152,10 @@ public:
*/
static Scalar vaporTemperature(Scalar pressure)
{
pressure = std::min(pressure, criticalPressure());
pressure = std::max(pressure, triplePressure());
using std::min;
using std::max;
pressure = min(pressure, criticalPressure());
pressure = max(pressure, triplePressure());
return Region4::vaporTemperature(pressure);
}
......@@ -580,7 +584,9 @@ public:
Scalar deltaP = pressure*2;
Valgrind::CheckDefined(pressure);
Valgrind::CheckDefined(deltaP);
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
using std::abs;
for (int i = 0; i < 5 && abs(pressure*1e-9) < abs(deltaP); ++i) {
Scalar f = gasDensity(temperature, pressure) - density;
Scalar df_dp;
......@@ -674,7 +680,9 @@ public:
Scalar eps = pressure*1e-7;
Scalar deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(pressure*1e-9) < std::abs(deltaP); ++i) {
using std::abs;
for (int i = 0; i < 5 && abs(pressure*1e-9) < abs(deltaP); ++i) {
Scalar f = liquidDensity(temperature, pressure) - density;
Scalar df_dp;
......
......@@ -102,8 +102,8 @@ public:
const Scalar mW = refComponentMolecularWeight() *1000. ; // in [g/mol];
return A+(B/mW)-(C/std::pow((mW+D),E));
using std::pow;
return A+(B/mW)-(C/pow((mW+D),E));
}
static Scalar perbutationFactorBoilingTemperature()
......@@ -114,12 +114,13 @@ public:
const Scalar D = 9.0180e-2;
const Scalar E = -1.0482;
Scalar deltaSpecificGravity = std::log(refComponentSpecificGravity()/specificGravity());
Scalar deltaMolecularWeight = std::log(refComponentMolecularWeight()/molecularWeight());
using std::log;
Scalar deltaSpecificGravity = log(refComponentSpecificGravity()/specificGravity());
Scalar deltaMolecularWeight = log(refComponentMolecularWeight()/molecularWeight());
return A*std::pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*std::pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
using std::pow;
return A*pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+ E*deltaSpecificGravity*deltaMolecularWeight;
}
static Scalar perbutationFactorCriticalTemperature()
......@@ -130,12 +131,13 @@ public:
const Scalar D = -5.7090e-2;
const Scalar E = -8.4583e-2;
Scalar deltaSpecificGravity = std::log(refComponentSpecificGravity()/specificGravity());
Scalar deltaMolecularWeight = std::log(refComponentMolecularWeight()/molecularWeight());
using std::log;
Scalar deltaSpecificGravity = log(refComponentSpecificGravity()/specificGravity());
Scalar deltaMolecularWeight = log(refComponentMolecularWeight()/molecularWeight());
return A*std::pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*std::pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
using std::pow;
return A*pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+ E*deltaSpecificGravity*deltaMolecularWeight;
}
static Scalar perbutationFactorCriticalPressure()
......@@ -146,75 +148,73 @@ public:
const Scalar D = -2.2389e-1;
const Scalar E = 2.6984;
Scalar deltaSpecificGravity = std::log(refComponentSpecificGravity()/specificGravity());
Scalar deltaMolecularWeight = std::log(refComponentMolecularWeight()/molecularWeight());
using std::log;
Scalar deltaSpecificGravity = log(refComponentSpecificGravity()/specificGravity());
Scalar deltaMolecularWeight = log(refComponentMolecularWeight()/molecularWeight());
return A*std::pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*std::pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
using std::pow;
return A*pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+ E*deltaSpecificGravity*deltaMolecularWeight;
}
static Scalar refComponentBoilingTemperature()
static Scalar refComponentBoilingTemperature()
{
const Scalar A = 477.63; //All factors for 1 atm / 101325 pascals [760 mmHg]
const Scalar B = 88.51;
const Scalar C = 1007;
const Scalar D = 1214.40;
return A*std::log((1000.*refComponentMolecularWeight() + B)/(1000.*refComponentMolecularWeight()+C)) + D;
using std::log;
return A*log((1000.*refComponentMolecularWeight() + B)/(1000.*refComponentMolecularWeight()+C)) + D;
}
static Scalar refComponentCriticalTemperature()
static Scalar refComponentCriticalTemperature()
{
const Scalar A = 226.50;
const Scalar B = 6.78;
const Scalar C = 1.282e6;
const Scalar D = 2668;
return A*std::log((1000.*refComponentMolecularWeight() + B)/(1000.*refComponentMolecularWeight()+C)) + D ;
using std::log;
return A*log((1000.*refComponentMolecularWeight() + B)/(1000.*refComponentMolecularWeight()+C)) + D ;
}
static Scalar refComponentCriticalPressure()
static Scalar refComponentCriticalPressure()
{
const Scalar A = 141.20;
const Scalar B = 45.66e-2;
const Scalar C = 16.59e-3;
const Scalar D = 2.19;
return (A*1000.*molecularWeight())/(std::pow(B + (C*1000.*molecularWeight()),D)) ;
using std::pow;
return (A*1000.*molecularWeight())/(pow(B + (C*1000.*molecularWeight()),D)) ;
}
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at heavyoil's boiling point (1 atm)
*/
static Scalar boilingTemperature()
static Scalar boilingTemperature()
{
return refComponentBoilingTemperature() * std::pow((1 + 2*perbutationFactorBoilingTemperature())/(1 - 2*perbutationFactorBoilingTemperature()),2);
using std::pow;
return refComponentBoilingTemperature() * pow((1 + 2*perbutationFactorBoilingTemperature())/(1 - 2*perbutationFactorBoilingTemperature()),2);
}
/*!
* \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of heavyoil
*/
static Scalar criticalTemperature()
static Scalar criticalTemperature()
{
return refComponentCriticalTemperature() * std::pow((1 + 2*perbutationFactorCriticalTemperature())/(1 - 2*perbutationFactorCriticalTemperature()),2);
using std::pow;
return refComponentCriticalTemperature() * pow((1 + 2*perbutationFactorCriticalTemperature())/(1 - 2*perbutationFactorCriticalTemperature()),2);
}
/*!
* \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of heavyoil
*/
static Scalar criticalPressure()
static Scalar criticalPressure()
{
return refComponentCriticalPressure() * std::pow((1 + 2*perbutationFactorCriticalPressure())/(1 - 2*perbutationFactorCriticalPressure()),2);
using std::pow;
return refComponentCriticalPressure() * pow((1 + 2*perbutationFactorCriticalPressure())/(1 - 2*perbutationFactorCriticalPressure()),2);
}
/*!
......@@ -231,8 +231,9 @@ public:
const Scalar C = 42.95101;
Scalar T = temperature - 273.15;
return 100*1.334*std::pow(10.0, (A - (B/(T + C)))); // in [Pa]
using std::pow;
return 100*1.334*pow(10.0, (A - (B/(T + C)))); // in [Pa]
}
static Scalar vaporTemperature(Scalar pressure)
......@@ -243,8 +244,8 @@ public:
const Scalar P = pressure;
return Scalar ((B/(A-std::log10(P/100*1.334)))-C);
using std::log10;
return Scalar ((B/(A-log10(P/100*1.334)))-C);
}
/*!
......@@ -268,8 +269,8 @@ public:
// Enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input and is
// fit over a certain temperature range. This suggests choosing an interval of integration being in the actual fit range.
// I.e. choosing T=273.15K as reference point for liquid enthalpy.
const Scalar sqrt1over3 = std::sqrt(1./3.);
using std::sqrt;
const Scalar sqrt1over3 = sqrt(1./3.);
const Scalar TEval1 = 0.5*(temperature-273.15)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration
const Scalar TEval2 = 0.5*(temperature-273.15)* (-1)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration
......@@ -289,22 +290,26 @@ public:
static Scalar heatVap(Scalar temperature,
const Scalar pressure)
{
temperature = std::min(temperature, criticalTemperature()); // regularization
temperature = std::max(temperature, 0.0); // regularization
using std::min;
using std::max;
temperature = min(temperature, criticalTemperature()); // regularization
temperature = max(temperature, 0.0); // regularization
Scalar T_crit = criticalTemperature();
Scalar Tr1 = boilingTemperature()/criticalTemperature();
Scalar p_crit = criticalPressure();
// Chen method, eq. 7-11.4 (at boiling)
using std::log;
const Scalar DH_v_boil = Consts::R * T_crit * Tr1
* (3.978 * Tr1 - 3.958 + 1.555*std::log(p_crit * 1e-5 /*Pa->bar*/ ) )
* (3.978 * Tr1 - 3.958 + 1.555*log(p_crit * 1e-5 /*Pa->bar*/ ) )
/ (1.07 - Tr1); /* [J/mol] */
/* Variation with temp according to Watson relation eq 7-12.1*/
using std::pow;
const Scalar 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 Scalar DH_vap = DH_v_boil * pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
return (DH_vap/molarMass()); // we need [J/kg]
}
......@@ -383,18 +388,22 @@ public:
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure, bool regularize=true)
{
temperature = std::min(temperature, 500.0); // regularization
temperature = std::max(temperature, 250.0);
using std::min;
using std::max;
temperature = min(temperature, 500.0); // regularization
temperature = max(temperature, 250.0);
// reduced temperature
Scalar Tr = temperature/criticalTemperature();
using std::pow;
using std::exp;
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)
Fp0*(0.807*pow(Tr,0.618)
- 0.357*exp(-0.449*Tr)
+ 0.34*exp(-4.058*Tr)
+ 0.018);
return eta_xi/xi/1e7; // [Pa s]
......@@ -414,12 +423,15 @@ public:
/* according to Lashanizadegan et al (2008) in Chemical Engineering Communications: */
/* Simultaneous Heat and Fluid Flow in Porous Media: Case Study: Steam Injection for Tertiary Oil Recovery */
//return 1027919.422*std::exp(-0.04862*tempera