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_BRINE_CO2_SYSTEM_HPP
00029 #define OPM_BRINE_CO2_SYSTEM_HPP
00030
00031 #include "BaseFluidSystem.hpp"
00032 #include "NullParameterCache.hpp"
00033
00034 #include <opm/material/IdealGas.hpp>
00035
00036 #include <opm/material/components/Brine.hpp>
00037 #include <opm/material/components/CO2.hpp>
00038 #include <opm/material/components/SimpleCO2.hpp>
00039 #include <opm/material/components/SimpleH2O.hpp>
00040 #include <opm/material/components/TabulatedComponent.hpp>
00041
00042 #include <opm/material/binarycoefficients/H2O_CO2.hpp>
00043 #include <opm/material/binarycoefficients/Brine_CO2.hpp>
00044 #include <opm/material/binarycoefficients/H2O_N2.hpp>
00045
00046 #include <opm/common/Unused.hpp>
00047
00048 #include <iostream>
00049
00050 namespace Opm {
00051 namespace FluidSystems {
00052
00060 template <class Scalar, class CO2Tables>
00061 class BrineCO2
00062 : public BaseFluidSystem<Scalar, BrineCO2<Scalar, CO2Tables> >
00063 {
00064 typedef Opm::H2O<Scalar> H2O_IAPWS;
00065 typedef Opm::Brine<Scalar, H2O_IAPWS> Brine_IAPWS;
00066 typedef Opm::TabulatedComponent<Scalar, H2O_IAPWS> H2O_Tabulated;
00067 typedef Opm::TabulatedComponent<Scalar, Brine_IAPWS> Brine_Tabulated;
00068
00069 typedef H2O_Tabulated H2O;
00070
00071 public:
00072 template <class Evaluation>
00073 struct ParameterCache : public Opm::NullParameterCache<Evaluation>
00074 {};
00075
00077 typedef Opm::BinaryCoeff::Brine_CO2<Scalar, CO2Tables> BinaryCoeffBrineCO2;
00078
00079
00080
00081
00082
00084 static const int numPhases = 2;
00085
00087 static const int liquidPhaseIdx = 0;
00089 static const int gasPhaseIdx = 1;
00090
00094 static const char* phaseName(unsigned phaseIdx)
00095 {
00096 static const char* name[] = {
00097 "liquid",
00098 "gas"
00099 };
00100
00101 assert(0 <= phaseIdx && phaseIdx < numPhases);
00102 return name[phaseIdx];
00103 }
00104
00108 static bool isLiquid(unsigned phaseIdx)
00109 {
00110 assert(0 <= phaseIdx && phaseIdx < numPhases);
00111
00112 return phaseIdx != gasPhaseIdx;
00113 }
00114
00118 static bool isIdealGas(unsigned phaseIdx)
00119 {
00120 assert(0 <= phaseIdx && phaseIdx < numPhases);
00121
00122 if (phaseIdx == gasPhaseIdx)
00123 return CO2::gasIsIdeal();
00124 return false;
00125 }
00126
00130 static bool isIdealMixture(unsigned phaseIdx OPM_OPTIM_UNUSED)
00131 {
00132 assert(0 <= phaseIdx && phaseIdx < numPhases);
00133
00134 return true;
00135 }
00136
00140 static bool isCompressible(unsigned phaseIdx OPM_OPTIM_UNUSED)
00141 {
00142 assert(0 <= phaseIdx && phaseIdx < numPhases);
00143
00144 return true;
00145 }
00146
00147
00148
00149
00151 static const int numComponents = 2;
00152
00154 static const int BrineIdx = 0;
00156 static const int CO2Idx = 1;
00157
00159 typedef Brine_Tabulated Brine;
00161 typedef Opm::CO2<Scalar, CO2Tables> CO2;
00162
00166 static const char* componentName(unsigned compIdx)
00167 {
00168 static const char* name[] = {
00169 Brine::name(),
00170 CO2::name(),
00171 };
00172
00173 assert(0 <= compIdx && compIdx < numComponents);
00174 return name[compIdx];
00175 }
00176
00180 static Scalar molarMass(unsigned compIdx)
00181 {
00182 assert(0 <= compIdx && compIdx < numComponents);
00183 return (compIdx==BrineIdx)
00184 ? Brine::molarMass()
00185 : CO2::molarMass();
00186 }
00187
00188
00189
00190
00191
00195 static void init()
00196 {
00197 init(273.15, 623.15, 50,
00198 1e4, 40e6, 50);
00199 }
00200
00212 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
00213 Scalar pressMin, Scalar pressMax, unsigned nPress)
00214 {
00215 if (H2O::isTabulated) {
00216 H2O_Tabulated::init(tempMin, tempMax, nTemp,
00217 pressMin, pressMax, nPress);
00218 }
00219
00220
00221 Brine_IAPWS::salinity = CO2Tables::brineSalinity;
00222
00223 if (Brine::isTabulated) {
00224 Brine_Tabulated::init(tempMin, tempMax, nTemp,
00225 pressMin, pressMax, nPress);
00226 }
00227 }
00228
00232 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00233 static LhsEval density(const FluidState& fluidState,
00234 const ParameterCache<ParamCacheEval>& ,
00235 unsigned phaseIdx)
00236 {
00237 assert(0 <= phaseIdx && phaseIdx < numPhases);
00238
00239 const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00240 const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00241
00242 if (phaseIdx == liquidPhaseIdx) {
00243
00244
00245
00246 LhsEval xlBrine = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, BrineIdx))));
00247 LhsEval xlCO2 = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, CO2Idx))));
00248 LhsEval sumx = xlBrine + xlCO2;
00249 xlBrine /= sumx;
00250 xlCO2 /= sumx;
00251
00252 LhsEval result = liquidDensity_(temperature,
00253 pressure,
00254 xlBrine,
00255 xlCO2);
00256
00257 Valgrind::CheckDefined(result);
00258 return result;
00259 }
00260
00261 assert(phaseIdx == gasPhaseIdx);
00262
00263
00264
00265
00266 LhsEval xgBrine = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, BrineIdx))));
00267 LhsEval xgCO2 = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, CO2Idx))));
00268 LhsEval sumx = xgBrine + xgCO2;
00269 xgBrine /= sumx;
00270 xgCO2 /= sumx;
00271
00272 LhsEval result = gasDensity_(temperature,
00273 pressure,
00274 xgBrine,
00275 xgCO2);
00276 Valgrind::CheckDefined(result);
00277 return result;
00278 }
00279
00283 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00284 static LhsEval viscosity(const FluidState& fluidState,
00285 const ParameterCache<ParamCacheEval>& ,
00286 unsigned phaseIdx)
00287 {
00288 assert(0 <= phaseIdx && phaseIdx < numPhases);
00289
00290 const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00291 const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00292
00293 if (phaseIdx == liquidPhaseIdx) {
00294
00295
00296 LhsEval result = Brine::liquidViscosity(temperature, pressure);
00297 Valgrind::CheckDefined(result);
00298 return result;
00299 }
00300
00301 assert(phaseIdx == gasPhaseIdx);
00302 LhsEval result = CO2::gasViscosity(temperature, pressure);
00303 Valgrind::CheckDefined(result);
00304 return result;
00305 }
00306
00310 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00311 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00312 const ParameterCache<ParamCacheEval>& ,
00313 unsigned phaseIdx,
00314 unsigned compIdx)
00315 {
00316 assert(0 <= phaseIdx && phaseIdx < numPhases);
00317 assert(0 <= compIdx && compIdx < numComponents);
00318
00319 if (phaseIdx == gasPhaseIdx)
00320
00321
00322
00323 return 1.0;
00324
00325 const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00326 const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00327 assert(temperature > 0);
00328 assert(pressure > 0);
00329
00330
00331
00332
00333 LhsEval xlH2O, xgH2O;
00334 LhsEval xlCO2, xgCO2;
00335 BinaryCoeffBrineCO2::calculateMoleFractions(temperature,
00336 pressure,
00337 Brine_IAPWS::salinity,
00338 -1,
00339 xlCO2,
00340 xgH2O);
00341
00342
00343 xlCO2 = Opm::max(0.0, Opm::min(1.0, xlCO2));
00344 xgH2O = Opm::max(0.0, Opm::min(1.0, xgH2O));
00345
00346 xlH2O = 1.0 - xlCO2;
00347 xgCO2 = 1.0 - xgH2O;
00348
00349 if (compIdx == BrineIdx) {
00350 Scalar phigH2O = 1.0;
00351 return phigH2O * xgH2O / xlH2O;
00352 }
00353 else {
00354 assert(compIdx == CO2Idx);
00355
00356 Scalar phigCO2 = 1.0;
00357 return phigCO2 * xgCO2 / xlCO2;
00358 };
00359 }
00360
00364 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00365 static LhsEval diffusionCoefficient(const FluidState& fluidState,
00366 const ParameterCache<ParamCacheEval>& ,
00367 unsigned phaseIdx,
00368 unsigned )
00369 {
00370 const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00371 const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00372 if (phaseIdx == liquidPhaseIdx)
00373 return BinaryCoeffBrineCO2::liquidDiffCoeff(temperature, pressure);
00374
00375 assert(phaseIdx == gasPhaseIdx);
00376 return BinaryCoeffBrineCO2::gasDiffCoeff(temperature, pressure);
00377 }
00378
00382 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00383 static LhsEval enthalpy(const FluidState& fluidState,
00384 const ParameterCache<ParamCacheEval>& ,
00385 unsigned phaseIdx)
00386 {
00387 assert(0 <= phaseIdx && phaseIdx < numPhases);
00388
00389 const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00390 const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00391
00392 if (phaseIdx == liquidPhaseIdx) {
00393 const LhsEval& XlCO2 = Opm::decay<LhsEval>(fluidState.massFraction(phaseIdx, CO2Idx));
00394 const LhsEval& result = liquidEnthalpyBrineCO2_(temperature,
00395 pressure,
00396 Brine_IAPWS::salinity,
00397 XlCO2);
00398 Valgrind::CheckDefined(result);
00399 return result;
00400 }
00401 else {
00402 const LhsEval& XCO2 = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, CO2Idx));
00403 const LhsEval& XBrine = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, BrineIdx));
00404
00405 LhsEval result = 0;
00406 result += XBrine * Brine::gasEnthalpy(temperature, pressure);
00407 result += XCO2 * CO2::gasEnthalpy(temperature, pressure);
00408 Valgrind::CheckDefined(result);
00409 return result;
00410 }
00411 }
00412
00416 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00417 static LhsEval thermalConductivity(const FluidState& ,
00418 const ParameterCache<ParamCacheEval>& ,
00419 unsigned phaseIdx)
00420 {
00421
00422 if (phaseIdx == liquidPhaseIdx)
00423 return 0.6;
00424
00425
00426 return 0.025;
00427 }
00428
00441 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00442 static LhsEval heatCapacity(const FluidState& fluidState,
00443 const ParameterCache<ParamCacheEval>& ,
00444 unsigned phaseIdx)
00445 {
00446 assert(0 <= phaseIdx && phaseIdx < numPhases);
00447
00448 const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00449 const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00450
00451 if(phaseIdx == liquidPhaseIdx)
00452 return H2O::liquidHeatCapacity(temperature, pressure);
00453 else
00454 return CO2::gasHeatCapacity(temperature, pressure);
00455 }
00456
00457 private:
00458 template <class LhsEval>
00459 static LhsEval gasDensity_(const LhsEval& T,
00460 const LhsEval& pg,
00461 const LhsEval& xgH2O,
00462 const LhsEval& xgCO2)
00463 {
00464 Valgrind::CheckDefined(T);
00465 Valgrind::CheckDefined(pg);
00466 Valgrind::CheckDefined(xgH2O);
00467 Valgrind::CheckDefined(xgCO2);
00468
00469 return CO2::gasDensity(T, pg);
00470 }
00471
00472
00473
00474
00475
00476
00477
00478 template <class LhsEval>
00479 static LhsEval liquidDensity_(const LhsEval& T,
00480 const LhsEval& pl,
00481 const LhsEval& xlH2O,
00482 const LhsEval& xlCO2)
00483 {
00484 Valgrind::CheckDefined(T);
00485 Valgrind::CheckDefined(pl);
00486 Valgrind::CheckDefined(xlH2O);
00487 Valgrind::CheckDefined(xlCO2);
00488
00489 if(T < 273.15) {
00490 OPM_THROW(NumericalProblem,
00491 "Liquid density for Brine and CO2 is only "
00492 "defined above 273.15K (is " << T << "K)");
00493 }
00494 if(pl >= 2.5e8) {
00495 OPM_THROW(NumericalProblem,
00496 "Liquid density for Brine and CO2 is only "
00497 "defined below 250MPa (is " << pl << "Pa)");
00498 }
00499
00500 const LhsEval& rho_brine = Brine::liquidDensity(T, pl);
00501 const LhsEval& rho_pure = H2O::liquidDensity(T, pl);
00502 const LhsEval& rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlH2O, xlCO2);
00503 const LhsEval& contribCO2 = rho_lCO2 - rho_pure;
00504
00505 return rho_brine + contribCO2;
00506 }
00507
00508 template <class LhsEval>
00509 static LhsEval liquidDensityWaterCO2_(const LhsEval& temperature,
00510 const LhsEval& pl,
00511 const LhsEval& ,
00512 const LhsEval& xlCO2)
00513 {
00514 Scalar M_CO2 = CO2::molarMass();
00515 Scalar M_H2O = H2O::molarMass();
00516
00517 const LhsEval& tempC = temperature - 273.15;
00518 const LhsEval& rho_pure = H2O::liquidDensity(temperature, pl);
00519
00520
00521
00522 const LhsEval xlH2O = 1.0 - xlCO2;
00523 const LhsEval& M_T = M_H2O * xlH2O + M_CO2 * xlCO2;
00524 const LhsEval& V_phi =
00525 (37.51 +
00526 tempC*(-9.585e-2 +
00527 tempC*(8.74e-4 -
00528 tempC*5.044e-7))) / 1.0e6;
00529 return 1/ (xlCO2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T));
00530 }
00531
00532 template <class LhsEval>
00533 static LhsEval liquidEnthalpyBrineCO2_(const LhsEval& T,
00534 const LhsEval& p,
00535 Scalar S,
00536 const LhsEval& X_CO2_w)
00537 {
00538
00539
00540
00541
00542
00543 static Scalar f[] = {
00544 2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10
00545 };
00546
00547
00548 static Scalar a[4][3] = {
00549 { 9633.6, -4080.0, +286.49 },
00550 { +166.58, +68.577, -4.6856 },
00551 { -0.90963, -0.36524, +0.249667E-1 },
00552 { +0.17965E-2, +0.71924E-3, -0.4900E-4 }
00553 };
00554
00555 LhsEval theta, h_NaCl;
00556 LhsEval h_ls1, d_h;
00557 LhsEval delta_h;
00558 LhsEval delta_hCO2, hg, hw;
00559
00560 theta = T - 273.15;
00561
00562
00563 Scalar scalarTheta = Opm::scalarValue(theta);
00564 Scalar S_lSAT = f[0] + scalarTheta*(f[1] + scalarTheta*(f[2] + scalarTheta*f[3]));
00565 if (S > S_lSAT)
00566 S = S_lSAT;
00567
00568 hw = H2O::liquidEnthalpy(T, p) /1E3;
00569
00570
00571 h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T
00572 +((2.8000E-5)/4)*(T*T*T*T))/(58.44E3)- 2.045698e+02;
00573
00574 Scalar m = 1E3/58.44 * S/(1-S);
00575 int i = 0;
00576 int j = 0;
00577 d_h = 0;
00578
00579 for (i = 0; i<=3; i++) {
00580 for (j=0; j<=2; j++) {
00581 d_h = d_h + a[i][j] * Opm::pow(theta, static_cast<Scalar>(i)) * std::pow(m, j);
00582 }
00583 }
00584
00585 delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
00586
00587
00588 h_ls1 =(1-S)*hw + S*h_NaCl + S*delta_h;
00589
00590
00591
00592
00593 delta_hCO2 = (-57.4375 + T * 0.1325) * 1000/44;
00594
00595
00596 hg = CO2::gasEnthalpy(T, p)/1E3 + delta_hCO2;
00597
00598
00599 return (h_ls1 - X_CO2_w*hw + hg*X_CO2_w)*1E3;
00600 }
00601 };
00602
00603 }
00604 }
00605
00606 #endif