Commit 1bf769d3 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[material] allow adl for math functions in components

parent d48e2b8a
......@@ -219,7 +219,9 @@ public:
{
const Scalar epsk = 103.3; // [K]
using namespace std;
using std::log;
using std::exp;
using std::sqrt;
const Scalar logTstar = log(temperature/epsk);
const Scalar Omega = exp(0.431
- 0.4623*logTstar
......@@ -230,6 +232,7 @@ public:
const Scalar sigma = 0.36; // [nm]
const Scalar eta0 = 0.0266958*sqrt(1000.0*molarMass()*temperature)/(sigma*sigma*Omega);
using std::pow;
const Scalar tau = criticalTemperature()/temperature;
const Scalar rhoc = 10.4477; // [mol/m^3]
const Scalar delta = 0.001*pressure/(temperature*8.3144598)/rhoc;
......
......@@ -150,7 +150,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 */
......@@ -160,6 +162,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++) {
......@@ -280,25 +283,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;
}
......@@ -329,7 +333,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;
......@@ -369,9 +375,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*temperature); // [Pa s]
//using std::exp;
//return 1027919.422*exp(-0.04862*temperature); // [Pa s]
//according to http://www.ecltechnology.com/subsur/reports/pvt_tgb.pdf [Page 10]
Scalar temperatureFahrenheit = (9/5)*(temperature-273.15)+32;
Scalar API = 9;
return ((std::pow(10,0.10231*std::pow(API,2)-3.9464*API+46.5037))*(std::pow(temperatureFahrenheit,-0.04542*std::pow(API,2)+1.70405*API-19.18)))*0.001;
using std::pow;
return ((pow(10,0.10231*pow(API,2)-3.9464*API+46.5037))*(pow(temperatureFahrenheit,-0.04542*pow(API,2)+1.70405*API-19.18)))*0.001;
}
/*!
......@@ -461,12 +473,15 @@ protected:
*/
static Scalar molarLiquidDensity_(Scalar temperature)
{
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);
using std::pow;
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 Scalar expo = 1.0 + pow(1.0 - temperature/criticalTemperature(), 2.0/7.0);
Scalar V = Consts::R*criticalTemperature()/criticalPressure()*pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
return 1.0/V; // molar density [mol/m^3]
}
......
......@@ -130,11 +130,13 @@ public:
muBar += tmp3 * tmp;
tmp3 *= 1.0/TBar - 1;
}
using std::exp;
muBar *= rhoBar;
muBar = std::exp(muBar);
muBar = exp(muBar);
// muBar *= muBar_0
muBar *= 100*std::sqrt(TBar);
using std::sqrt;
muBar *= 100*sqrt(TBar);
constexpr Scalar H[4] = {
1.67752, 2.20462, 0.6366564, -0.241605
};
......@@ -197,6 +199,7 @@ public: