00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_IAPWS_COMMON_HPP
00028 #define OPM_IAPWS_COMMON_HPP
00029
00030 #include <opm/material/Constants.hpp>
00031 #include <opm/material/common/MathToolbox.hpp>
00032
00033 #include <cmath>
00034
00035 namespace Opm {
00036 namespace IAPWS {
00037
00053 template <class Scalar>
00054 class Common
00055 {
00056 public:
00058 static const Scalar molarMass;
00059
00061 static const Scalar Rs;
00062
00064 static const Scalar criticalTemperature;
00065
00067 static const Scalar criticalPressure;
00068
00070 static const Scalar criticalDensity;
00071
00073 static const Scalar criticalMolarVolume;
00074
00076 static const Scalar acentricFactor;
00077
00079 static const Scalar tripleTemperature;
00080
00082 static const Scalar triplePressure;
00083
00098 template <class Evaluation>
00099 static Evaluation viscosity(const Evaluation& temperature, const Evaluation& rho)
00100 {
00101 Evaluation rhoBar = rho/322.0;
00102 Evaluation TBar = temperature/criticalTemperature;
00103
00104
00105 const Scalar Hij[6][7] = {
00106 { 5.20094e-1, 2.22531e-1,-2.81378e-1, 1.61913e-1,-3.25372e-2, 0, 0 },
00107 { 8.50895e-2, 9.99115e-1,-9.06851e-1, 2.57399e-1, 0, 0, 0 },
00108 { -1.08374, 1.88797 ,-7.72479e-1, 0, 0, 0, 0 },
00109 { -2.89555e-1, 1.26613 ,-4.89837e-1, 0, 6.98452e-2, 0,-4.35673e-3 },
00110 { 0, 0,-2.57040e-1, 0, 0, 8.72102e-3, 0 },
00111 { 0, 1.20573e-1, 0, 0, 0, 0,-5.93264e-4 }
00112 };
00113
00114 Evaluation tmp, tmp2, tmp3 = 1;
00115 Evaluation muBar = 0;
00116 for (int i = 0; i <= 5; ++i) {
00117 tmp = 0;
00118 tmp2 = 1;
00119 for (int j = 0; j <= 6; ++j) {
00120 tmp += Hij[i][j]*tmp2;
00121 tmp2 *= (rhoBar - 1);
00122 };
00123 muBar += tmp3 * tmp;
00124 tmp3 *= 1.0/TBar - 1;
00125 };
00126 muBar *= rhoBar;
00127 muBar = Opm::exp(muBar);
00128
00129
00130 muBar *= 100*Opm::sqrt(TBar);
00131 const Scalar H[4] = {
00132 1.67752, 2.20462, 0.6366564, -0.241605
00133 };
00134
00135 tmp = 0, tmp2 = 1;
00136 for (int i = 0; i < 4; ++i) {
00137 tmp += H[i]/tmp2;
00138 tmp2 *= TBar;
00139 };
00140 muBar /= tmp;
00141
00142 return 1e-6*muBar;
00143 }
00144
00158 template <class Evaluation>
00159 static Evaluation thermalConductivityIAPWS(const Evaluation& T, const Evaluation& rho)
00160 {
00161 static const Scalar thcond_tstar = 647.26 ;
00162 static const Scalar thcond_rhostar = 317.7 ;
00163
00164
00165 static const Scalar thcond_b0 = -0.397070 ;
00166 static const Scalar thcond_b1 = 0.400302 ;
00167 static const Scalar thcond_b2 = 1.060000 ;
00168 static const Scalar thcond_B1 = -0.171587 ;
00169 static const Scalar thcond_B2 = 2.392190 ;
00170
00171 static const Scalar thcond_c1 = 0.642857 ;
00172 static const Scalar thcond_c2 = -4.11717 ;
00173 static const Scalar thcond_c3 = -6.17937 ;
00174 static const Scalar thcond_c4 = 0.00308976 ;
00175 static const Scalar thcond_c5 = 0.0822994 ;
00176 static const Scalar thcond_c6 = 10.0932 ;
00177
00178 static const Scalar thcond_d1 = 0.0701309 ;
00179 static const Scalar thcond_d2 = 0.0118520 ;
00180 static const Scalar thcond_d3 = 0.00169937 ;
00181 static const Scalar thcond_d4 = -1.0200 ;
00182 static const int thcond_a_count = 4;
00183 static const Scalar thcond_a[thcond_a_count] = {
00184 0.0102811
00185 ,0.0299621
00186 ,0.0156146
00187 ,-0.00422464
00188 };
00189
00190 Evaluation Tbar = T / thcond_tstar;
00191 Evaluation rhobar = rho / thcond_rhostar;
00192
00193
00194 Evaluation Troot = Opm::sqrt(Tbar);
00195 Evaluation Tpow = Troot;
00196 Evaluation lam = 0;
00197
00198 for(int k = 0; k < thcond_a_count; ++k) {
00199 lam += thcond_a[k] * Tpow;
00200 Tpow *= Tbar;
00201 }
00202
00203 lam +=
00204 thcond_b0 + thcond_b1
00205 * rhobar + thcond_b2
00206 * Opm::exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
00207
00208 Evaluation DTbar = Opm::abs(Tbar - 1) + thcond_c4;
00209 Evaluation DTbarpow = Opm::pow(DTbar, 3./5);
00210 Evaluation Q = 2. + thcond_c5 / DTbarpow;
00211
00212 Evaluation S;
00213 if(Tbar >= 1)
00214 S = 1. / DTbar;
00215 else
00216 S = thcond_c6 / DTbarpow;
00217
00218 Evaluation rhobar18 = Opm::pow(rhobar, 1.8);
00219 Evaluation rhobarQ = Opm::pow(rhobar, Q);
00220
00221 lam +=
00222 (thcond_d1 / Opm::pow(Tbar,10.0) + thcond_d2) * rhobar18 *
00223 Opm::exp(thcond_c1 * (1 - rhobar * rhobar18))
00224 + thcond_d3 * S * rhobarQ *
00225 Opm::exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
00226 + thcond_d4 *
00227 Opm::exp(thcond_c2 * Opm::pow(Troot,3.0) + thcond_c3 / Opm::pow(rhobar,5.0));
00228 return lam;
00229 }
00230 };
00231
00232 template <class Scalar>
00233 const Scalar Common<Scalar>::molarMass = 18.01518e-3;
00234 template <class Scalar>
00235 const Scalar Common<Scalar>::Rs = Opm::Constants<Scalar>::R/molarMass;
00236 template <class Scalar>
00237 const Scalar Common<Scalar>::criticalTemperature = 647.096;
00238 template <class Scalar>
00239 const Scalar Common<Scalar>::criticalPressure = 22.064e6;
00240 template <class Scalar>
00241 const Scalar Common<Scalar>::criticalDensity = 322.0;
00242 template <class Scalar>
00243 const Scalar Common<Scalar>::criticalMolarVolume = molarMass/criticalDensity;
00244 template <class Scalar>
00245 const Scalar Common<Scalar>::acentricFactor = 0.344;
00246 template <class Scalar>
00247 const Scalar Common<Scalar>::tripleTemperature = 273.16;
00248 template <class Scalar>
00249 const Scalar Common<Scalar>::triplePressure = 611.657;
00250
00251 }
00252 }
00253
00254 #endif