00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00028 #ifndef OPM_BINARY_COEFF_BRINE_CO2_HPP
00029 #define OPM_BINARY_COEFF_BRINE_CO2_HPP
00030
00031 #include <opm/material/components/Brine.hpp>
00032 #include <opm/material/components/H2O.hpp>
00033 #include <opm/material/components/CO2.hpp>
00034 #include <opm/material/IdealGas.hpp>
00035
00036 namespace Opm {
00037 namespace BinaryCoeff {
00038
00043 template<class Scalar, class CO2Tables, bool verbose = true>
00044 class Brine_CO2 {
00045 typedef Opm::H2O<Scalar> H2O;
00046 typedef Opm::CO2<Scalar, CO2Tables> CO2;
00047 typedef Opm::IdealGas<Scalar> IdealGas;
00048 static const int liquidPhaseIdx = 0;
00049 static const int gasPhaseIdx = 1;
00050
00051 public:
00059 template <class Evaluation>
00060 static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
00061 {
00062
00063 Scalar k = 1.3806504e-23;
00064 Scalar c = 4;
00065 Scalar R_h = 1.72e-10;
00066 const Evaluation& mu = CO2::gasViscosity(temperature, pressure);
00067 return k / (c * M_PI * R_h) * (temperature / mu);
00068 }
00069
00076 template <class Evaluation>
00077 static Evaluation liquidDiffCoeff(const Evaluation& , const Evaluation& )
00078 {
00079
00080 return 2e-9;
00081 }
00082
00100 template <class Evaluation>
00101 static void calculateMoleFractions(const Evaluation& temperature,
00102 const Evaluation& pg,
00103 Scalar salinity,
00104 const int knownPhaseIdx,
00105 Evaluation& xlCO2,
00106 Evaluation& ygH2O)
00107 {
00108 Evaluation A = computeA_(temperature, pg);
00109
00110
00111 Scalar x_NaCl = salinityToMolFrac_(salinity);
00112
00113
00114
00115 if (knownPhaseIdx < 0) {
00116 Scalar molalityNaCl = moleFracToMolality_(x_NaCl);
00117 Evaluation m0_CO2 = molalityCO2inPureWater_(temperature, pg);
00118 Evaluation gammaStar = activityCoefficient_(temperature, pg, molalityNaCl);
00119 Evaluation m_CO2 = m0_CO2 / gammaStar;
00120 xlCO2 = m_CO2 / (molalityNaCl + 55.508 + m_CO2);
00121 ygH2O = A * (1 - xlCO2 - x_NaCl);
00122 }
00123
00124
00125
00126
00127 if (knownPhaseIdx == liquidPhaseIdx)
00128 ygH2O = A * (1 - xlCO2 - x_NaCl);
00129
00130
00131
00132
00133 if (knownPhaseIdx == gasPhaseIdx)
00134
00135 xlCO2 = 1 - x_NaCl - ygH2O / A;
00136 }
00137
00141 template <class Evaluation>
00142 static Evaluation henry(const Evaluation& temperature)
00143 { return fugacityCoefficientCO2(temperature, 1e5)*1e5; }
00144
00153 template <class Evaluation>
00154 static Evaluation fugacityCoefficientCO2(const Evaluation& temperature, const Evaluation& pg)
00155 {
00156 Valgrind::CheckDefined(temperature);
00157 Valgrind::CheckDefined(pg);
00158
00159 Evaluation V = 1 / (CO2::gasDensity(temperature, pg) / CO2::molarMass()) * 1.e6;
00160 Evaluation pg_bar = pg / 1.e5;
00161 Evaluation a_CO2 = (7.54e7 - 4.13e4 * temperature);
00162 Scalar b_CO2 = 27.8;
00163 Scalar R = IdealGas::R * 10.;
00164 Evaluation lnPhiCO2;
00165
00166 lnPhiCO2 = Opm::log(V / (V - b_CO2));
00167 lnPhiCO2 += b_CO2 / (V - b_CO2);
00168 lnPhiCO2 -= 2 * a_CO2 / (R * Opm::pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
00169 lnPhiCO2 +=
00170 a_CO2 * b_CO2
00171 / (R
00172 * Opm::pow(temperature, 1.5)
00173 * b_CO2
00174 * b_CO2)
00175 * (Opm::log((V + b_CO2) / V)
00176 - b_CO2 / (V + b_CO2));
00177 lnPhiCO2 -= Opm::log(pg_bar * V / (R * temperature));
00178
00179 return Opm::exp(lnPhiCO2);
00180 }
00181
00190 template <class Evaluation>
00191 static Evaluation fugacityCoefficientH2O(const Evaluation& temperature, const Evaluation& pg)
00192 {
00193 const Evaluation& V = 1 / (CO2::gasDensity(temperature, pg) / CO2::molarMass()) * 1.e6;
00194 const Evaluation& pg_bar = pg / 1.e5;
00195 const Evaluation& a_CO2 = (7.54e7 - 4.13e4 * temperature);
00196 Scalar a_CO2_H2O = 7.89e7;
00197 Scalar b_CO2 = 27.8;
00198 Scalar b_H2O = 18.18;
00199 Scalar R = IdealGas::R * 10.;
00200 Evaluation lnPhiH2O;
00201
00202 lnPhiH2O =
00203 Opm::log(V/(V - b_CO2))
00204 + b_H2O/(V - b_CO2) - 2*a_CO2_H2O
00205 / (R*Opm::pow(temperature, 1.5)*b_CO2)*Opm::log((V + b_CO2)/V)
00206 + a_CO2*b_H2O/(R*Opm::pow(temperature, 1.5)*b_CO2*b_CO2)
00207 *(Opm::log((V + b_CO2)/V) - b_CO2/(V + b_CO2))
00208 - Opm::log(pg_bar*V/(R*temperature));
00209 return Opm::exp(lnPhiH2O);
00210 }
00211
00212 private:
00218 static Scalar salinityToMolFrac_(Scalar salinity) {
00219
00220 const Scalar Mw = H2O::molarMass();
00221 const Scalar Ms = 58.8e-3;
00222
00223 const Scalar X_NaCl = salinity;
00224
00225 const Scalar x_NaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
00226 return x_NaCl;
00227 }
00228
00234 static Scalar moleFracToMolality_(Scalar x_NaCl)
00235 {
00236
00237 return 55.508 * x_NaCl / (1 - x_NaCl);
00238 }
00239
00247 template <class Evaluation>
00248 static Evaluation molalityCO2inPureWater_(const Evaluation& temperature, const Evaluation& pg)
00249 {
00250 const Evaluation& A = computeA_(temperature, pg);
00251 const Evaluation& B = computeB_(temperature, pg);
00252 const Evaluation& yH2OinGas = (1 - B) / (1. / A - B);
00253 const Evaluation& xCO2inWater = B * (1 - yH2OinGas);
00254 return (xCO2inWater * 55.508) / (1 - xCO2inWater);
00255 }
00256
00266 template <class Evaluation>
00267 static Evaluation activityCoefficient_(const Evaluation& temperature,
00268 const Evaluation& pg,
00269 Scalar molalityNaCl)
00270 {
00271 const Evaluation& lambda = computeLambda_(temperature, pg);
00272 const Evaluation& xi = computeXi_(temperature, pg);
00273 const Evaluation& lnGammaStar =
00274 2*molalityNaCl*lambda + xi*molalityNaCl*molalityNaCl;
00275 return Opm::exp(lnGammaStar);
00276 }
00277
00286 template <class Evaluation>
00287 static Evaluation computeA_(const Evaluation& temperature, const Evaluation& pg)
00288 {
00289 const Evaluation& deltaP = pg / 1e5 - 1;
00290 Scalar v_av_H2O = 18.1;
00291 Scalar R = IdealGas::R * 10;
00292 const Evaluation& k0_H2O = equilibriumConstantH2O_(temperature);
00293 const Evaluation& phi_H2O = fugacityCoefficientH2O(temperature, pg);
00294 const Evaluation& pg_bar = pg / 1.e5;
00295 return k0_H2O/(phi_H2O*pg_bar)*Opm::exp(deltaP*v_av_H2O/(R*temperature));
00296 }
00297
00306 template <class Evaluation>
00307 static Evaluation computeB_(const Evaluation& temperature, const Evaluation& pg)
00308 {
00309 const Evaluation& deltaP = pg / 1e5 - 1;
00310 const Scalar v_av_CO2 = 32.6;
00311 const Scalar R = IdealGas::R * 10;
00312 const Evaluation& k0_CO2 = equilibriumConstantCO2_(temperature);
00313 const Evaluation& phi_CO2 = fugacityCoefficientCO2(temperature, pg);
00314 const Evaluation& pg_bar = pg / 1.e5;
00315 return phi_CO2*pg_bar/(55.508*k0_CO2)*Opm::exp(-(deltaP*v_av_CO2)/(R*temperature));
00316 }
00317
00325 template <class Evaluation>
00326 static Evaluation computeLambda_(const Evaluation& temperature, const Evaluation& pg)
00327 {
00328 static const Scalar c[6] =
00329 { -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
00330
00331 Evaluation pg_bar = pg / 1.0E5;
00332 return
00333 c[0]
00334 + c[1]*temperature
00335 + c[2]/temperature
00336 + c[3]*pg_bar/temperature
00337 + c[4]*pg_bar/(630.0 - temperature)
00338 + c[5]*temperature*Opm::log(pg_bar);
00339 }
00340
00348 template <class Evaluation>
00349 static Evaluation computeXi_(const Evaluation& temperature, const Evaluation& pg)
00350 {
00351 static const Scalar c[4] =
00352 { 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
00353
00354 Evaluation pg_bar = pg / 1.0E5;
00355 return c[0] + c[1]*temperature + c[2]*pg_bar/temperature + c[3]*pg_bar/(630.0 - temperature);
00356 }
00357
00364 template <class Evaluation>
00365 static Evaluation equilibriumConstantCO2_(const Evaluation& temperature)
00366 {
00367 Evaluation temperatureCelcius = temperature - 273.15;
00368 static const Scalar c[3] = { 1.189, 1.304e-2, -5.446e-5 };
00369 Evaluation logk0_CO2 = c[0] + temperatureCelcius*(c[1] + temperatureCelcius*c[2]);
00370 Evaluation k0_CO2 = Opm::pow(10.0, logk0_CO2);
00371 return k0_CO2;
00372 }
00373
00380 template <class Evaluation>
00381 static Evaluation equilibriumConstantH2O_(const Evaluation& temperature)
00382 {
00383 Evaluation temperatureCelcius = temperature - 273.15;
00384 static const Scalar c[4] = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7 };
00385 Evaluation logk0_H2O =
00386 c[0] + temperatureCelcius*(c[1] + temperatureCelcius*(c[2] + temperatureCelcius*c[3]));
00387 return Opm::pow(10.0, logk0_H2O);
00388 }
00389
00390 };
00391
00392 }
00393 }
00394
00395 #endif