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_BLACK_OIL_FLUID_SYSTEM_HPP
00028 #define OPM_BLACK_OIL_FLUID_SYSTEM_HPP
00029
00030 #include "blackoilpvt/OilPvtMultiplexer.hpp"
00031 #include "blackoilpvt/GasPvtMultiplexer.hpp"
00032 #include "blackoilpvt/WaterPvtMultiplexer.hpp"
00033
00034 #include <opm/material/fluidsystems/BaseFluidSystem.hpp>
00035 #include <opm/material/Constants.hpp>
00036
00037 #include <opm/material/common/MathToolbox.hpp>
00038 #include <opm/common/Valgrind.hpp>
00039 #include <opm/material/common/HasMemberGeneratorMacros.hpp>
00040 #include <opm/common/Exceptions.hpp>
00041 #include <opm/common/ErrorMacros.hpp>
00042
00043 #include <memory>
00044 #include <vector>
00045 #include <array>
00046
00047 namespace Opm {
00048 namespace BlackOil {
00049 OPM_GENERATE_HAS_MEMBER(Rs, )
00050 OPM_GENERATE_HAS_MEMBER(Rv, )
00051
00052 template <class FluidSystem, class LhsEval, class FluidState>
00053 LhsEval getRs_(typename std::enable_if<!HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
00054 unsigned regionIdx)
00055 {
00056 const auto& XoG =
00057 Opm::decay<LhsEval>(fluidState.massFraction(FluidSystem::oilPhaseIdx,
00058 FluidSystem::gasCompIdx));
00059 return FluidSystem::convertXoGToRs(XoG, regionIdx);
00060 }
00061
00062 template <class FluidSystem, class LhsEval, class FluidState>
00063 auto getRs_(typename std::enable_if<HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
00064 unsigned regionIdx OPM_UNUSED)
00065 -> decltype(Opm::MathToolbox<typename FluidState::Scalar>
00066 ::template decay<LhsEval>(fluidState.Rs()))
00067 {
00068 return Opm::decay<LhsEval>(fluidState.Rs());
00069 }
00070
00071 template <class FluidSystem, class LhsEval, class FluidState>
00072 LhsEval getRv_(typename std::enable_if<!HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
00073 unsigned regionIdx)
00074 {
00075 const auto& XgO =
00076 Opm::decay<LhsEval>(fluidState.massFraction(FluidSystem::gasPhaseIdx,
00077 FluidSystem::oilCompIdx));
00078 return FluidSystem::convertXgOToRv(XgO, regionIdx);
00079 }
00080
00081 template <class FluidSystem, class LhsEval, class FluidState>
00082 auto getRv_(typename std::enable_if<HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
00083 unsigned regionIdx OPM_UNUSED)
00084 -> decltype(Opm::MathToolbox<typename FluidState::Scalar>
00085 ::template decay<LhsEval>(fluidState.Rv()))
00086 {
00087 return Opm::decay<LhsEval>(fluidState.Rv());
00088 }
00089 }
00090
00091 namespace FluidSystems {
00092
00099 template <class Scalar>
00100 class BlackOil : public BaseFluidSystem<Scalar, BlackOil<Scalar> >
00101 {
00102 typedef BlackOil ThisType;
00103
00104 public:
00105 typedef Opm::GasPvtMultiplexer<Scalar> GasPvt;
00106 typedef Opm::OilPvtMultiplexer<Scalar> OilPvt;
00107 typedef Opm::WaterPvtMultiplexer<Scalar> WaterPvt;
00108
00110 template <class EvaluationT>
00111 struct ParameterCache : public Opm::NullParameterCache<EvaluationT>
00112 {
00113 typedef EvaluationT Evaluation;
00114
00115 public:
00116 ParameterCache(Scalar maxOilSat = 1.0, unsigned regionIdx=0)
00117 {
00118 maxOilSat_ = maxOilSat;
00119 regionIdx_ = regionIdx;
00120 }
00121
00129 template <class OtherCache>
00130 void assignPersistentData(const OtherCache& other)
00131 {
00132 regionIdx_ = other.regionIndex();
00133 maxOilSat_ = other.maxOilSat();
00134 }
00135
00143 unsigned regionIndex() const
00144 { return regionIdx_; }
00145
00153 void setRegionIndex(unsigned val)
00154 { regionIdx_ = val; }
00155
00156 const Scalar maxOilSat() const
00157 { return maxOilSat_; }
00158
00159 void setMaxOilSat(Scalar val)
00160 { maxOilSat_ = val; }
00161
00162 private:
00163 Scalar maxOilSat_;
00164 unsigned regionIdx_;
00165 };
00166
00167
00168
00169
00170 #if HAVE_OPM_PARSER
00171
00174 static void initFromDeck(const Deck& deck, const EclipseState& eclState)
00175 {
00176 auto densityKeyword = deck.getKeyword("DENSITY");
00177 size_t numRegions = densityKeyword.size();
00178 initBegin(numRegions);
00179
00180 numActivePhases_ = 0;
00181 std::fill(&phaseIsActive_[0], &phaseIsActive_[numPhases], false);
00182
00183 if (deck.hasKeyword("OIL")) {
00184 phaseIsActive_[oilPhaseIdx] = true;
00185 ++ numActivePhases_;
00186 }
00187
00188 if (deck.hasKeyword("GAS")) {
00189 phaseIsActive_[gasPhaseIdx] = true;
00190 ++ numActivePhases_;
00191 }
00192
00193 if (deck.hasKeyword("WATER")) {
00194 phaseIsActive_[waterPhaseIdx] = true;
00195 ++ numActivePhases_;
00196 }
00197
00198
00199
00200 setReservoirTemperature(eclState.getTableManager().rtemp());
00201
00202
00203 assert(numActivePhases_ >= 2 && numActivePhases_ <= 3);
00204
00205 setEnableDissolvedGas(deck.hasKeyword("DISGAS"));
00206 setEnableVaporizedOil(deck.hasKeyword("VAPOIL"));
00207
00208
00209 for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
00210 const auto& densityRecord = densityKeyword.getRecord(regionIdx);
00211 setReferenceDensities(densityRecord.getItem("OIL").getSIDouble(0),
00212 densityRecord.getItem("WATER").getSIDouble(0),
00213 densityRecord.getItem("GAS").getSIDouble(0),
00214 regionIdx);
00215 }
00216
00217 if (phaseIsActive(gasPhaseIdx)) {
00218 gasPvt_ = std::make_shared<GasPvt>();
00219 gasPvt_->initFromDeck(deck, eclState);
00220 }
00221
00222 if (phaseIsActive(oilPhaseIdx)) {
00223 oilPvt_ = std::make_shared<OilPvt>();
00224 oilPvt_->initFromDeck(deck, eclState);
00225 }
00226
00227 if (phaseIsActive(waterPhaseIdx)) {
00228 waterPvt_ = std::make_shared<WaterPvt>();
00229 waterPvt_->initFromDeck(deck, eclState);
00230 }
00231
00232 initEnd();
00233 }
00234 #endif // HAVE_OPM_PARSER
00235
00244 static void initBegin(size_t numPvtRegions)
00245 {
00246 enableDissolvedGas_ = true;
00247 enableVaporizedOil_ = false;
00248
00249 numActivePhases_ = numPhases;
00250 std::fill(&phaseIsActive_[0], &phaseIsActive_[numPhases], true);
00251
00252 resizeArrays_(numPvtRegions);
00253
00254 setReservoirTemperature(surfaceTemperature);
00255 }
00256
00263 static void setEnableDissolvedGas(bool yesno)
00264 { enableDissolvedGas_ = yesno; }
00265
00272 static void setEnableVaporizedOil(bool yesno)
00273 { enableVaporizedOil_ = yesno; }
00274
00278 static void setGasPvt(std::shared_ptr<GasPvt> pvtObj)
00279 { gasPvt_ = pvtObj; }
00280
00284 static void setOilPvt(std::shared_ptr<OilPvt> pvtObj)
00285 { oilPvt_ = pvtObj; }
00286
00290 static void setWaterPvt(std::shared_ptr<WaterPvt> pvtObj)
00291 { waterPvt_ = pvtObj; }
00292
00300 static void setReferenceDensities(Scalar rhoOil,
00301 Scalar rhoWater,
00302 Scalar rhoGas,
00303 unsigned regionIdx)
00304 {
00305 referenceDensity_[regionIdx][oilPhaseIdx] = rhoOil;
00306 referenceDensity_[regionIdx][waterPhaseIdx] = rhoWater;
00307 referenceDensity_[regionIdx][gasPhaseIdx] = rhoGas;
00308 }
00309
00313 static void initEnd()
00314 {
00315
00316 size_t numRegions = molarMass_.size();
00317 for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
00318
00319
00320
00321 molarMass_[regionIdx][waterCompIdx] = 18e-3;
00322
00323 if (phaseIsActive(gasPhaseIdx)) {
00324
00325 Scalar p = surfacePressure;
00326 Scalar T = surfaceTemperature;
00327 Scalar rho_g = referenceDensity_[0][gasPhaseIdx];
00328 molarMass_[regionIdx][gasCompIdx] = Opm::Constants<Scalar>::R*T*rho_g / p;
00329 }
00330 else
00331
00332 molarMass_[regionIdx][gasCompIdx] = 2e-3;
00333
00334
00335 molarMass_[regionIdx][oilCompIdx] = 175e-3;
00336 }
00337
00338 isInitialized_ = true;
00339 }
00340
00341 static bool isInitialized()
00342 { return isInitialized_; }
00343
00344
00345
00346
00347
00349 static const unsigned numPhases = 3;
00350
00352 static const unsigned waterPhaseIdx = 0;
00354 static const unsigned oilPhaseIdx = 1;
00356 static const unsigned gasPhaseIdx = 2;
00357
00359 static const Scalar surfacePressure;
00360
00362 static const Scalar surfaceTemperature;
00363
00365 static const char* phaseName(unsigned phaseIdx)
00366 {
00367 static const char* name[] = { "water", "oil", "gas" };
00368
00369 assert(0 <= phaseIdx && phaseIdx < numPhases + 1);
00370 return name[phaseIdx];
00371 }
00372
00374 static bool isLiquid(unsigned phaseIdx)
00375 {
00376 assert(0 <= phaseIdx && phaseIdx < numPhases);
00377 return phaseIdx != gasPhaseIdx;
00378 }
00379
00380
00381
00382
00383
00385 static const unsigned numComponents = 3;
00386
00388 static const unsigned oilCompIdx = 0;
00390 static const unsigned waterCompIdx = 1;
00392 static const unsigned gasCompIdx = 2;
00393
00394 protected:
00395 static const int phaseToSolventCompIdx_[3];
00396 static const int phaseToSoluteCompIdx_[3];
00397
00398 static unsigned char numActivePhases_;
00399 static bool phaseIsActive_[numPhases];
00400
00401 public:
00403 static unsigned numActivePhases()
00404 { return numActivePhases_; }
00405
00407 static unsigned phaseIsActive(unsigned phaseIdx)
00408 {
00409 assert(phaseIdx < numPhases);
00410 return phaseIsActive_[phaseIdx];
00411 }
00412
00414 static constexpr unsigned solventComponentIndex(unsigned phaseIdx)
00415 { return static_cast<unsigned>(phaseToSolventCompIdx_[phaseIdx]); }
00416
00418 static constexpr unsigned soluteComponentIndex(unsigned phaseIdx)
00419 { return static_cast<unsigned>(phaseToSoluteCompIdx_[phaseIdx]); }
00420
00422 static const char* componentName(unsigned compIdx)
00423 {
00424 static const char* name[] = { "Oil", "Water", "Gas" };
00425
00426 assert(0 <= compIdx && compIdx < numComponents);
00427 return name[compIdx];
00428 }
00429
00431 static Scalar molarMass(unsigned compIdx, unsigned regionIdx = 0)
00432 { return molarMass_[regionIdx][compIdx]; }
00433
00435 static bool isIdealMixture(unsigned )
00436 {
00437
00438
00439 return true;
00440 }
00441
00443 static bool isCompressible(unsigned )
00444 { return true; }
00445
00447 static bool isIdealGas(unsigned )
00448 { return false; }
00449
00450
00451
00452
00453
00459 static size_t numRegions()
00460 { return molarMass_.size(); }
00461
00468 static bool enableDissolvedGas()
00469 { return enableDissolvedGas_; }
00470
00477 static bool enableVaporizedOil()
00478 { return enableVaporizedOil_; }
00479
00485 static Scalar referenceDensity(unsigned phaseIdx, unsigned regionIdx)
00486 { return referenceDensity_[regionIdx][phaseIdx]; }
00487
00488
00489
00490
00492 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00493 static LhsEval density(const FluidState& fluidState,
00494 const ParameterCache<ParamCacheEval>& paramCache,
00495 unsigned phaseIdx)
00496 { return density<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
00497
00499 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00500 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00501 const ParameterCache<ParamCacheEval>& paramCache,
00502 unsigned phaseIdx,
00503 unsigned compIdx)
00504 {
00505 return fugacityCoefficient<FluidState, LhsEval>(fluidState,
00506 phaseIdx,
00507 compIdx,
00508 paramCache.regionIndex());
00509 }
00510
00512 template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
00513 static LhsEval viscosity(const FluidState& fluidState,
00514 const ParameterCache<ParamCacheEval>& paramCache,
00515 unsigned phaseIdx)
00516 { return viscosity<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
00517
00518
00519
00520
00521
00522
00524 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00525 static LhsEval density(const FluidState& fluidState,
00526 unsigned phaseIdx,
00527 unsigned regionIdx)
00528 {
00529 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00530 assert(0 <= regionIdx && regionIdx <= numRegions());
00531
00532 switch (phaseIdx) {
00533 case oilPhaseIdx: {
00534 if (!enableDissolvedGas()) {
00535
00536 const auto& bo = inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
00537 return referenceDensity(phaseIdx, regionIdx)*bo;
00538 }
00539
00540
00541 const auto& bo = inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
00542 const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00543
00544 return
00545 bo*referenceDensity(oilPhaseIdx, regionIdx)
00546 + Rs*bo*referenceDensity(gasPhaseIdx, regionIdx);
00547 }
00548
00549 case gasPhaseIdx: {
00550 if (!enableVaporizedOil()) {
00551
00552 const auto& bg = inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
00553 return bg*referenceDensity(phaseIdx, regionIdx);
00554 }
00555
00556
00557 const auto& bg = inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
00558 const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00559
00560 return
00561 bg*referenceDensity(gasPhaseIdx, regionIdx)
00562 + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
00563 }
00564
00565 case waterPhaseIdx:
00566 return
00567 referenceDensity(waterPhaseIdx, regionIdx)
00568 *inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, waterPhaseIdx, regionIdx);
00569 }
00570
00571 OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00572 }
00573
00581 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00582 static LhsEval saturatedDensity(const FluidState& fluidState,
00583 unsigned phaseIdx,
00584 unsigned regionIdx)
00585 {
00586 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00587 assert(0 <= regionIdx && regionIdx <= numRegions());
00588
00589 switch (phaseIdx) {
00590 case oilPhaseIdx: {
00591 if (!enableDissolvedGas()) {
00592
00593 const auto& bo = inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
00594 return referenceDensity(phaseIdx, regionIdx)*bo;
00595 }
00596
00597
00598 const auto& bo = saturatedInverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
00599 const auto& Rs = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
00600
00601 return
00602 bo*referenceDensity(oilPhaseIdx, regionIdx)
00603 + Rs*bo*referenceDensity(gasPhaseIdx, regionIdx);
00604 }
00605
00606 case gasPhaseIdx: {
00607 if (!enableVaporizedOil()) {
00608
00609 const auto& bg = inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
00610 return referenceDensity(phaseIdx, regionIdx)*bg;
00611 }
00612
00613
00614 const auto& bg = saturatedInverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
00615 const auto& Rv = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
00616
00617 return
00618 bg*referenceDensity(gasPhaseIdx, regionIdx)
00619 + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
00620 }
00621
00622 case waterPhaseIdx:
00623 return
00624 referenceDensity(waterPhaseIdx, regionIdx)
00625 *saturatedInverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, waterPhaseIdx, regionIdx);
00626 }
00627
00628 OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00629 }
00630
00639 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00640 static LhsEval inverseFormationVolumeFactor(const FluidState& fluidState,
00641 unsigned phaseIdx,
00642 unsigned regionIdx)
00643 {
00644 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00645 assert(0 <= regionIdx && regionIdx <= numRegions());
00646
00647 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00648 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00649
00650 switch (phaseIdx) {
00651 case oilPhaseIdx: {
00652 if (enableDissolvedGas()) {
00653 if (fluidState.phaseIsPresent(gasPhaseIdx)) {
00654 if (fluidState.saturation(gasPhaseIdx) < 1e-4) {
00655
00656
00657
00658 const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00659 const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
00660 const auto& bSat = oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
00661 const auto& bUndersat = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
00662 return alpha*bSat + (1.0 - alpha)*bUndersat;
00663 }
00664
00665 return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
00666 }
00667
00668 const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00669 return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
00670 }
00671
00672 const LhsEval Rs(0.0);
00673 return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
00674 }
00675 case gasPhaseIdx: {
00676 if (enableVaporizedOil()) {
00677 if (fluidState.phaseIsPresent(oilPhaseIdx)) {
00678 if (fluidState.saturation(oilPhaseIdx) < 1e-4) {
00679
00680
00681
00682 const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00683 const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
00684 const auto& bSat = gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
00685 const auto& bUndersat = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
00686 return alpha*bSat + (1.0 - alpha)*bUndersat;
00687 }
00688
00689 return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
00690 }
00691
00692 const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00693 return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
00694 }
00695
00696 const LhsEval Rv(0.0);
00697 return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
00698 }
00699 case waterPhaseIdx:
00700 return waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p);
00701 default: OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00702 }
00703 }
00704
00712 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00713 static LhsEval saturatedInverseFormationVolumeFactor(const FluidState& fluidState,
00714 unsigned phaseIdx,
00715 unsigned regionIdx)
00716 {
00717 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00718 assert(0 <= regionIdx && regionIdx <= numRegions());
00719
00720 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00721 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00722
00723 switch (phaseIdx) {
00724 case oilPhaseIdx: return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
00725 case gasPhaseIdx: return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
00726 case waterPhaseIdx: return waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p);
00727 default: OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00728 }
00729 }
00730
00732 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00733 static LhsEval fugacityCoefficient(const FluidState& fluidState,
00734 unsigned phaseIdx,
00735 unsigned compIdx,
00736 unsigned regionIdx)
00737 {
00738 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00739 assert(0 <= compIdx && compIdx <= numComponents);
00740 assert(0 <= regionIdx && regionIdx <= numRegions());
00741
00742 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00743 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00744
00745
00746
00747
00748 const LhsEval phi_oO = 20e3/p;
00749
00750
00751 const Scalar phi_gG = 1.0;
00752
00753
00754
00755 const LhsEval phi_wW = 30e3/p;
00756
00757 switch (phaseIdx) {
00758 case gasPhaseIdx:
00759 switch (compIdx) {
00760 case gasCompIdx:
00761 return phi_gG;
00762
00763
00764
00765
00766
00767 case oilCompIdx: {
00768 if (!enableVaporizedOil())
00769
00770
00771 return phi_gG*1e6;
00772
00773 const auto& R_vSat = gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
00774 const auto& X_gOSat = convertRvToXgO(R_vSat, regionIdx);
00775 const auto& x_gOSat = convertXgOToxgO(X_gOSat, regionIdx);
00776
00777 const auto& R_sSat = oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
00778 const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
00779 const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
00780 const auto& x_oOSat = 1.0 - x_oGSat;
00781
00782 const auto& p_o = Opm::decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
00783 const auto& p_g = Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
00784
00785 return phi_oO*p_o*x_oOSat / (p_g*x_gOSat);
00786 }
00787
00788 case waterCompIdx:
00789
00790 return phi_gG*1e6;
00791
00792 default:
00793 OPM_THROW(std::logic_error,
00794 "Invalid component index " << compIdx);
00795 }
00796
00797 case oilPhaseIdx:
00798 switch (compIdx) {
00799 case oilCompIdx:
00800 return phi_oO;
00801
00802
00803
00804 case gasCompIdx: {
00805 if (!enableDissolvedGas())
00806
00807
00808 return phi_oO*1e6;
00809
00810 const auto& R_vSat = gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
00811 const auto& X_gOSat = convertRvToXgO(R_vSat, regionIdx);
00812 const auto& x_gOSat = convertXgOToxgO(X_gOSat, regionIdx);
00813 const auto& x_gGSat = 1.0 - x_gOSat;
00814
00815 const auto& R_sSat = oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
00816 const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
00817 const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
00818
00819 const auto& p_o = Opm::decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
00820 const auto& p_g = Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
00821
00822 return phi_gG*p_g*x_gGSat / (p_o*x_oGSat);
00823 }
00824
00825 case waterCompIdx:
00826 return phi_oO*1e6;
00827
00828 default:
00829 OPM_THROW(std::logic_error,
00830 "Invalid component index " << compIdx);
00831 }
00832
00833 case waterPhaseIdx:
00834
00835
00836
00837
00838
00839
00840 switch (compIdx) {
00841 case waterCompIdx: return phi_wW;
00842 case oilCompIdx: return 1.1e6*phi_wW;
00843 case gasCompIdx: return 1e6*phi_wW;
00844 default:
00845 OPM_THROW(std::logic_error,
00846 "Invalid component index " << compIdx);
00847 }
00848
00849 default:
00850 OPM_THROW(std::logic_error,
00851 "Invalid phase index " << phaseIdx);
00852 }
00853
00854 OPM_THROW(std::logic_error, "Unhandled phase or component index");
00855 }
00856
00858 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00859 static LhsEval viscosity(const FluidState& fluidState,
00860 unsigned phaseIdx,
00861 unsigned regionIdx)
00862 {
00863 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00864 assert(0 <= regionIdx && regionIdx <= numRegions());
00865
00866 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00867 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00868
00869 switch (phaseIdx) {
00870 case oilPhaseIdx: {
00871 if (enableDissolvedGas()) {
00872 if (fluidState.phaseIsPresent(gasPhaseIdx)) {
00873 if (fluidState.saturation(gasPhaseIdx) < 1e-4) {
00874
00875
00876
00877 const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00878 const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
00879 const auto& muSat = oilPvt_->saturatedViscosity(regionIdx, T, p);
00880 const auto& muUndersat = oilPvt_->viscosity(regionIdx, T, p, Rs);
00881 return alpha*muSat + (1.0 - alpha)*muUndersat;
00882 }
00883
00884 return oilPvt_->saturatedViscosity(regionIdx, T, p);
00885 }
00886
00887 const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00888 return oilPvt_->viscosity(regionIdx, T, p, Rs);
00889 }
00890
00891 const LhsEval Rs(0.0);
00892 return oilPvt_->viscosity(regionIdx, T, p, Rs);
00893 }
00894
00895 case gasPhaseIdx: {
00896 if (enableVaporizedOil()) {
00897 if (fluidState.phaseIsPresent(oilPhaseIdx)) {
00898 if (fluidState.saturation(oilPhaseIdx) < 1e-4) {
00899
00900
00901
00902 const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00903 const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
00904 const auto& muSat = gasPvt_->saturatedViscosity(regionIdx, T, p);
00905 const auto& muUndersat = gasPvt_->viscosity(regionIdx, T, p, Rv);
00906 return alpha*muSat + (1.0 - alpha)*muUndersat;
00907 }
00908
00909 return gasPvt_->saturatedViscosity(regionIdx, T, p);
00910 }
00911
00912 const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
00913 return gasPvt_->viscosity(regionIdx, T, p, Rv);
00914 }
00915
00916 const LhsEval Rv(0.0);
00917 return gasPvt_->viscosity(regionIdx, T, p, Rv);
00918 }
00919
00920 case waterPhaseIdx:
00921
00922
00923 return waterPvt_->viscosity(regionIdx, T, p);
00924 }
00925
00926 OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00927 }
00928
00935 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00936 static LhsEval saturatedDissolutionFactor(const FluidState& fluidState,
00937 unsigned phaseIdx,
00938 unsigned regionIdx,
00939 Scalar maxOilSaturation)
00940 {
00941 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00942 assert(0 <= regionIdx && regionIdx <= numRegions());
00943
00944 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00945 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00946 const auto& So = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx));
00947
00948 switch (phaseIdx) {
00949 case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p, So, maxOilSaturation);
00950 case gasPhaseIdx: return gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p, So, maxOilSaturation);
00951 case waterPhaseIdx: return 0.0;
00952 default: OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00953 }
00954 }
00955
00964 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00965 static LhsEval saturatedDissolutionFactor(const FluidState& fluidState,
00966 unsigned phaseIdx,
00967 unsigned regionIdx)
00968 {
00969 assert(0 <= phaseIdx && phaseIdx <= numPhases);
00970 assert(0 <= regionIdx && regionIdx <= numRegions());
00971
00972 const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
00973 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
00974
00975 switch (phaseIdx) {
00976 case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
00977 case gasPhaseIdx: return gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
00978 case waterPhaseIdx: return 0.0;
00979 default: OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
00980 }
00981 }
00982
00986 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00987 static LhsEval bubblePointPressure(const FluidState& fluidState,
00988 unsigned regionIdx)
00989 {
00990 return saturationPressure(fluidState, oilPhaseIdx, regionIdx);
00991 }
00992
00993
00997 template <class FluidState, class LhsEval = typename FluidState::Scalar>
00998 static LhsEval dewPointPressure(const FluidState& fluidState,
00999 unsigned regionIdx)
01000 {
01001 return saturationPressure(fluidState, gasPhaseIdx, regionIdx);
01002 }
01003
01014 template <class FluidState, class LhsEval = typename FluidState::Scalar>
01015 static LhsEval saturationPressure(const FluidState& fluidState,
01016 unsigned phaseIdx,
01017 unsigned regionIdx)
01018 {
01019 assert(0 <= phaseIdx && phaseIdx <= numPhases);
01020 assert(0 <= regionIdx && regionIdx <= numRegions());
01021
01022 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
01023
01024 switch (phaseIdx) {
01025 case oilPhaseIdx: return oilPvt_->saturationPressure(regionIdx, T, Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx));
01026 case gasPhaseIdx: return gasPvt_->saturationPressure(regionIdx, T, Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx));
01027 case waterPhaseIdx: return 0.0;
01028 default: OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
01029 }
01030 }
01031
01032
01033
01034
01039 template <class LhsEval>
01040 static LhsEval convertXoGToRs(const LhsEval& XoG, unsigned regionIdx)
01041 {
01042 Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
01043 Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
01044
01045 return XoG/(1.0 - XoG)*(rho_oRef/rho_gRef);
01046 }
01047
01052 template <class LhsEval>
01053 static LhsEval convertXgOToRv(const LhsEval& XgO, unsigned regionIdx)
01054 {
01055 Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
01056 Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
01057
01058 return XgO/(1.0 - XgO)*(rho_gRef/rho_oRef);
01059 }
01060
01065 template <class LhsEval>
01066 static LhsEval convertRsToXoG(const LhsEval& Rs, unsigned regionIdx)
01067 {
01068 Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
01069 Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
01070
01071 const LhsEval& rho_oG = Rs*rho_gRef;
01072 return rho_oG/(rho_oRef + rho_oG);
01073 }
01074
01079 template <class LhsEval>
01080 static LhsEval convertRvToXgO(const LhsEval& Rv, unsigned regionIdx)
01081 {
01082 Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
01083 Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
01084
01085 const LhsEval& rho_gO = Rv*rho_oRef;
01086 return rho_gO/(rho_gRef + rho_gO);
01087 }
01088
01092 template <class LhsEval>
01093 static LhsEval convertXoGToxoG(const LhsEval& XoG, unsigned regionIdx)
01094 {
01095 Scalar MO = molarMass_[regionIdx][oilCompIdx];
01096 Scalar MG = molarMass_[regionIdx][gasCompIdx];
01097
01098 return XoG*MO / (MG*(1 - XoG) + XoG*MO);
01099 }
01100
01104 template <class LhsEval>
01105 static LhsEval convertxoGToXoG(const LhsEval& xoG, unsigned regionIdx)
01106 {
01107 Scalar MO = molarMass_[regionIdx][oilCompIdx];
01108 Scalar MG = molarMass_[regionIdx][gasCompIdx];
01109
01110 return xoG*MG / (xoG*(MG - MO) + MO);
01111 }
01112
01116 template <class LhsEval>
01117 static LhsEval convertXgOToxgO(const LhsEval& XgO, unsigned regionIdx)
01118 {
01119 Scalar MO = molarMass_[regionIdx][oilCompIdx];
01120 Scalar MG = molarMass_[regionIdx][gasCompIdx];
01121
01122 return XgO*MG / (MO*(1 - XgO) + XgO*MG);
01123 }
01124
01128 template <class LhsEval>
01129 static LhsEval convertxgOToXgO(const LhsEval& xgO, unsigned regionIdx)
01130 {
01131 Scalar MO = molarMass_[regionIdx][oilCompIdx];
01132 Scalar MG = molarMass_[regionIdx][gasCompIdx];
01133
01134 return xgO*MO / (xgO*(MO - MG) + MG);
01135 }
01136
01144 static const GasPvt& gasPvt()
01145 { return *gasPvt_; }
01146
01154 static const OilPvt& oilPvt()
01155 { return *oilPvt_; }
01156
01164 static const WaterPvt& waterPvt()
01165 { return *waterPvt_; }
01166
01172 static Scalar reservoirTemperature(unsigned pvtRegionIdx OPM_UNUSED = 0)
01173 { return reservoirTemperature_; }
01174
01180 static void setReservoirTemperature(Scalar value)
01181 { reservoirTemperature_ = value; }
01182
01183 private:
01184 static void resizeArrays_(size_t numRegions)
01185 {
01186 molarMass_.resize(numRegions);
01187 referenceDensity_.resize(numRegions);
01188 }
01189
01190 static Scalar reservoirTemperature_;
01191
01192 static std::shared_ptr<GasPvt> gasPvt_;
01193 static std::shared_ptr<OilPvt> oilPvt_;
01194 static std::shared_ptr<WaterPvt> waterPvt_;
01195
01196 static bool enableDissolvedGas_;
01197 static bool enableVaporizedOil_;
01198
01199
01200
01201
01202 static std::vector<std::array<Scalar, 3> > referenceDensity_;
01203 static std::vector<std::array<Scalar, 3> > molarMass_;
01204
01205 static bool isInitialized_;
01206 };
01207
01208 template <class Scalar>
01209 const int BlackOil<Scalar>::phaseToSolventCompIdx_[3] =
01210 {
01211 waterCompIdx,
01212 oilCompIdx,
01213 gasCompIdx
01214 };
01215
01216
01217 template <class Scalar>
01218 const int BlackOil<Scalar>::phaseToSoluteCompIdx_[3] =
01219 {
01220 -1,
01221 gasCompIdx,
01222 oilCompIdx
01223 };
01224
01225 template <class Scalar>
01226 unsigned char BlackOil<Scalar>::numActivePhases_;
01227
01228 template <class Scalar>
01229 bool BlackOil<Scalar>::phaseIsActive_[numPhases];
01230
01231 template <class Scalar>
01232 const Scalar
01233 BlackOil<Scalar>::surfaceTemperature = 273.15 + 15.56;
01234
01235 template <class Scalar>
01236 const Scalar
01237 BlackOil<Scalar>::surfacePressure = 101325.0;
01238
01239 template <class Scalar>
01240 Scalar
01241 BlackOil<Scalar>::reservoirTemperature_;
01242
01243 template <class Scalar>
01244 bool BlackOil<Scalar>::enableDissolvedGas_;
01245
01246 template <class Scalar>
01247 bool BlackOil<Scalar>::enableVaporizedOil_;
01248
01249 template <class Scalar>
01250 std::shared_ptr<OilPvtMultiplexer<Scalar> >
01251 BlackOil<Scalar>::oilPvt_;
01252
01253 template <class Scalar>
01254 std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> >
01255 BlackOil<Scalar>::gasPvt_;
01256
01257 template <class Scalar>
01258 std::shared_ptr<WaterPvtMultiplexer<Scalar> >
01259 BlackOil<Scalar>::waterPvt_;
01260
01261 template <class Scalar>
01262 std::vector<std::array<Scalar, 3> >
01263 BlackOil<Scalar>::referenceDensity_;
01264
01265 template <class Scalar>
01266 std::vector<std::array<Scalar, 3> >
01267 BlackOil<Scalar>::molarMass_;
01268
01269 template <class Scalar>
01270 bool BlackOil<Scalar>::isInitialized_ = false;
01271
01272 }}
01273
01274 #endif