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_SPE5_FLUID_SYSTEM_HPP
00028 #define OPM_SPE5_FLUID_SYSTEM_HPP
00029
00030 #include "BaseFluidSystem.hpp"
00031 #include "Spe5ParameterCache.hpp"
00032
00033 #include <opm/material/Constants.hpp>
00034 #include <opm/material/eos/PengRobinsonMixture.hpp>
00035
00036 #include <opm/material/common/Spline.hpp>
00037
00038 namespace Opm {
00039 namespace FluidSystems {
00054 template <class Scalar>
00055 class Spe5
00056 : public BaseFluidSystem<Scalar, Spe5<Scalar> >
00057 {
00058 typedef Opm::FluidSystems::Spe5<Scalar> ThisType;
00059
00060 typedef typename Opm::PengRobinsonMixture<Scalar, ThisType> PengRobinsonMixture;
00061 typedef typename Opm::PengRobinson<Scalar> PengRobinson;
00062
00063 static const Scalar R;
00064
00065 public:
00067 template <class Evaluation>
00068 struct ParameterCache : public Opm::Spe5ParameterCache<Evaluation, ThisType>
00069 {};
00070
00071
00072
00073
00074
00076 static const int numPhases = 3;
00077
00079 static const int gasPhaseIdx = 0;
00081 static const int waterPhaseIdx = 1;
00083 static const int oilPhaseIdx = 2;
00084
00086 typedef Opm::H2O<Scalar> H2O;
00087
00089 static const char* phaseName(unsigned phaseIdx)
00090 {
00091 static const char* name[] = {
00092 "gas",
00093 "water",
00094 "oil",
00095 };
00096
00097 assert(0 <= phaseIdx && phaseIdx < numPhases);
00098 return name[phaseIdx];
00099 }
00100
00102 static bool isLiquid(unsigned phaseIdx)
00103 {
00104
00105 return phaseIdx != gasPhaseIdx;
00106 }
00107
00113 static bool isCompressible(unsigned )
00114 {
00115
00116 return true;
00117 }
00118
00120 static bool isIdealGas(unsigned )
00121 {
00122
00123 return false;
00124 }
00125
00127 static bool isIdealMixture(unsigned phaseIdx)
00128 {
00129
00130
00131
00132 return phaseIdx == waterPhaseIdx;
00133 }
00134
00135
00136
00137
00138
00140 static const int numComponents = 7;
00141
00142 static const int H2OIdx = 0;
00143 static const int C1Idx = 1;
00144 static const int C3Idx = 2;
00145 static const int C6Idx = 3;
00146 static const int C10Idx = 4;
00147 static const int C15Idx = 5;
00148 static const int C20Idx = 6;
00149
00151 static const char* componentName(unsigned compIdx)
00152 {
00153 static const char* name[] = {
00154 H2O::name(),
00155 "C1",
00156 "C3",
00157 "C6",
00158 "C10",
00159 "C15",
00160 "C20"
00161 };
00162
00163 assert(0 <= compIdx && compIdx < numComponents);
00164 return name[compIdx];
00165 }
00166
00168 static Scalar molarMass(unsigned compIdx)
00169 {
00170 return
00171 (compIdx == H2OIdx)
00172 ? H2O::molarMass()
00173 : (compIdx == C1Idx)
00174 ? 16.04e-3
00175 : (compIdx == C3Idx)
00176 ? 44.10e-3
00177 : (compIdx == C6Idx)
00178 ? 86.18e-3
00179 : (compIdx == C10Idx)
00180 ? 142.29e-3
00181 : (compIdx == C15Idx)
00182 ? 206.00e-3
00183 : (compIdx == C20Idx)
00184 ? 282.00e-3
00185 : 1e30;
00186 }
00187
00191 static Scalar criticalTemperature(unsigned compIdx)
00192 {
00193 return
00194 (compIdx == H2OIdx)
00195 ? H2O::criticalTemperature()
00196 : (compIdx == C1Idx)
00197 ? 343.0*5/9
00198 : (compIdx == C3Idx)
00199 ? 665.7*5/9
00200 : (compIdx == C6Idx)
00201 ? 913.4*5/9
00202 : (compIdx == C10Idx)
00203 ? 1111.8*5/9
00204 : (compIdx == C15Idx)
00205 ? 1270.0*5/9
00206 : (compIdx == C20Idx)
00207 ? 1380.0*5/9
00208 : 1e30;
00209 }
00210
00214 static Scalar criticalPressure(unsigned compIdx)
00215 {
00216 return
00217 (compIdx == H2OIdx)
00218 ? H2O::criticalPressure()
00219 : (compIdx == C1Idx)
00220 ? 667.8 * 6894.7573
00221 : (compIdx == C3Idx)
00222 ? 616.3 * 6894.7573
00223 : (compIdx == C6Idx)
00224 ? 436.9 * 6894.7573
00225 : (compIdx == C10Idx)
00226 ? 304.0 * 6894.7573
00227 : (compIdx == C15Idx)
00228 ? 200.0 * 6894.7573
00229 : (compIdx == C20Idx)
00230 ? 162.0 * 6894.7573
00231 : 1e30;
00232 }
00233
00237 static Scalar criticalMolarVolume(unsigned compIdx)
00238 {
00239 return
00240 (compIdx == H2OIdx)
00241 ? H2O::criticalMolarVolume()
00242 : (compIdx == C1Idx)
00243 ? 0.290*R*criticalTemperature(C1Idx)/criticalPressure(C1Idx)
00244 : (compIdx == C3Idx)
00245 ? 0.277*R*criticalTemperature(C3Idx)/criticalPressure(C3Idx)
00246 : (compIdx == C6Idx)
00247 ? 0.264*R*criticalTemperature(C6Idx)/criticalPressure(C6Idx)
00248 : (compIdx == C10Idx)
00249 ? 0.257*R*criticalTemperature(C10Idx)/criticalPressure(C10Idx)
00250 : (compIdx == C15Idx)
00251 ? 0.245*R*criticalTemperature(C15Idx)/criticalPressure(C15Idx)
00252 : (compIdx == C20Idx)
00253 ? 0.235*R*criticalTemperature(C20Idx)/criticalPressure(C20Idx)
00254 : 1e30;
00255 }
00256
00260 static Scalar acentricFactor(unsigned compIdx)
00261 {
00262 return
00263 (compIdx == H2OIdx)
00264 ? H2O::acentricFactor()
00265 : (compIdx == C1Idx)
00266 ? 0.0130
00267 : (compIdx == C3Idx)
00268 ? 0.1524
00269 : (compIdx == C6Idx)
00270 ? 0.3007
00271 : (compIdx == C10Idx)
00272 ? 0.4885
00273 : (compIdx == C15Idx)
00274 ? 0.6500
00275 : (compIdx == C20Idx)
00276 ? 0.8500
00277 : 1e30;
00278 }
00279
00285 static Scalar interactionCoefficient(unsigned comp1Idx, unsigned comp2Idx)
00286 {
00287 unsigned i = std::min(comp1Idx, comp2Idx);
00288 unsigned j = std::max(comp1Idx, comp2Idx);
00289 if (i == C1Idx && (j == C15Idx || j == C20Idx))
00290 return 0.05;
00291 else if (i == C3Idx && (j == C15Idx || j == C20Idx))
00292 return 0.005;
00293 return 0;
00294 }
00295
00296
00297
00298
00299
00308 static void init(Scalar minT = 273.15,
00309 Scalar maxT = 373.15,
00310 Scalar minP = 1e4,
00311 Scalar maxP = 100e6)
00312 {
00313 Opm::PengRobinsonParamsMixture<Scalar, ThisType, gasPhaseIdx, true> prParams;
00314
00315
00316
00317
00318
00319
00320
00321 Scalar minA = 1e30, maxA = -1e30;
00322 Scalar minB = 1e30, maxB = -1e30;
00323
00324 prParams.updatePure(minT, minP);
00325 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00326 minA = std::min(prParams.pureParams(compIdx).a(), minA);
00327 maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
00328 minB = std::min(prParams.pureParams(compIdx).b(), minB);
00329 maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
00330 };
00331
00332 prParams.updatePure(maxT, minP);
00333 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00334 minA = std::min(prParams.pureParams(compIdx).a(), minA);
00335 maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
00336 minB = std::min(prParams.pureParams(compIdx).b(), minB);
00337 maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
00338 };
00339
00340 prParams.updatePure(minT, maxP);
00341 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00342 minA = std::min(prParams.pureParams(compIdx).a(), minA);
00343 maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
00344 minB = std::min(prParams.pureParams(compIdx).b(), minB);
00345 maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
00346 };
00347
00348 prParams.updatePure(maxT, maxP);
00349 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00350 minA = std::min(prParams.pureParams(compIdx).a(), minA);
00351 maxA = std::max(prParams.pureParams(compIdx).a(), maxA);
00352 minB = std::min(prParams.pureParams(compIdx).b(), minB);
00353 maxB = std::max(prParams.pureParams(compIdx).b(), maxB);
00354 };
00355
00356 PengRobinson::init(minA, maxA, 100,
00357 minB, maxB, 200);
00358 }
00359
00361 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00362 static LhsEval density(const FluidState& fluidState,
00363 const ParameterCache<ParamCacheEval>& paramCache,
00364 unsigned phaseIdx)
00365 {
00366 assert(0 <= phaseIdx && phaseIdx < numPhases);
00367
00368 return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx);
00369 }
00370
00372 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00373 static LhsEval viscosity(const FluidState& ,
00374 const ParameterCache<ParamCacheEval>& ,
00375 unsigned phaseIdx)
00376 {
00377 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00378
00379 if (phaseIdx == gasPhaseIdx) {
00380
00381
00382 return 0.0170e-2 * 0.1;
00383 }
00384 else if (phaseIdx == waterPhaseIdx)
00385
00386 return 0.7e-2 * 0.1;
00387 else {
00388 assert(phaseIdx == oilPhaseIdx);
00389
00390
00391 return 0.208e-2 * 0.1;
00392 }
00393 }
00394
00396 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00397 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00398 const ParameterCache<ParamCacheEval>& paramCache,
00399 unsigned phaseIdx,
00400 unsigned compIdx)
00401 {
00402 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00403 assert(0 <= compIdx && compIdx <= numComponents);
00404
00405 if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx)
00406 return PengRobinsonMixture::computeFugacityCoefficient(fluidState,
00407 paramCache,
00408 phaseIdx,
00409 compIdx);
00410 else {
00411 assert(phaseIdx == waterPhaseIdx);
00412 return
00413 henryCoeffWater_(compIdx, fluidState.temperature(waterPhaseIdx))
00414 / fluidState.pressure(waterPhaseIdx);
00415 }
00416 }
00417
00418 protected:
00419 template <class LhsEval>
00420 static LhsEval henryCoeffWater_(unsigned compIdx, const LhsEval& temperature)
00421 {
00422
00423
00424 switch (compIdx) {
00425 case H2OIdx: return H2O::vaporPressure(temperature);
00426
00427
00428
00429 case C1Idx: return 5.57601e+09;
00430 case C3Idx: return 1e10;
00431 case C6Idx: return 1e10;
00432 case C10Idx: return 1e10;
00433 case C15Idx: return 1e10;
00434 case C20Idx: return 1e10;
00435 default: OPM_THROW(std::logic_error, "Unknown component index " << compIdx);
00436 }
00437 }
00438 };
00439
00440 template <class Scalar>
00441 const Scalar Spe5<Scalar>::R = Opm::Constants<Scalar>::R;
00442
00443 }
00444 }
00445
00446 #endif