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_WET_GAS_PVT_HPP
00028 #define OPM_WET_GAS_PVT_HPP
00029
00030 #include <opm/material/Constants.hpp>
00031 #include <opm/material/common/OpmFinal.hpp>
00032 #include <opm/material/common/UniformXTabulated2DFunction.hpp>
00033 #include <opm/material/common/Tabulated1DFunction.hpp>
00034
00035 #if HAVE_OPM_PARSER
00036 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00037 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
00038 #endif
00039
00040 namespace Opm {
00041
00046 template <class Scalar>
00047 class WetGasPvt
00048 {
00049 typedef Opm::UniformXTabulated2DFunction<Scalar> TabulatedTwoDFunction;
00050 typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
00051 typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
00052
00053 public:
00054 WetGasPvt()
00055 {
00056 vapPar1_ = 0.0;
00057 }
00058
00059 #if HAVE_OPM_PARSER
00060
00065 void initFromDeck(const Deck& deck, const EclipseState& eclState)
00066 {
00067 const auto& pvtgTables = eclState.getTableManager().getPvtgTables();
00068 const auto& densityKeyword = deck.getKeyword("DENSITY");
00069
00070 assert(pvtgTables.size() == densityKeyword.size());
00071
00072 size_t numRegions = pvtgTables.size();
00073 setNumRegions(numRegions);
00074
00075 for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
00076 Scalar rhoRefO = densityKeyword.getRecord(regionIdx).getItem("OIL").getSIDouble(0);
00077 Scalar rhoRefG = densityKeyword.getRecord(regionIdx).getItem("GAS").getSIDouble(0);
00078 Scalar rhoRefW = densityKeyword.getRecord(regionIdx).getItem("WATER").getSIDouble(0);
00079
00080 setReferenceDensities(regionIdx, rhoRefO, rhoRefG, rhoRefW);
00081 }
00082
00083 for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
00084 const auto& pvtgTable = pvtgTables[regionIdx];
00085
00086 const auto& saturatedTable = pvtgTable.getSaturatedTable();
00087 assert(saturatedTable.numRows() > 1);
00088
00089 auto& gasMu = gasMu_[regionIdx];
00090 auto& invGasB = inverseGasB_[regionIdx];
00091 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
00092 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
00093 auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
00094
00095 oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
00096 saturatedTable.getColumn("PG"),
00097 saturatedTable.getColumn("RV"));
00098
00099 std::vector<Scalar> invSatGasBArray;
00100 std::vector<Scalar> invSatGasBMuArray;
00101
00102
00103 for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
00104 Scalar pg = saturatedTable.get("PG" , outerIdx);
00105 Scalar B = saturatedTable.get("BG" , outerIdx);
00106 Scalar mu = saturatedTable.get("MUG" , outerIdx);
00107
00108 invGasB.appendXPos(pg);
00109 gasMu.appendXPos(pg);
00110
00111 invSatGasBArray.push_back(1.0/B);
00112 invSatGasBMuArray.push_back(1.0/(mu*B));
00113
00114 assert(invGasB.numX() == outerIdx + 1);
00115 assert(gasMu.numX() == outerIdx + 1);
00116
00117 const auto& underSaturatedTable = pvtgTable.getUnderSaturatedTable(outerIdx);
00118 size_t numRows = underSaturatedTable.numRows();
00119 for (size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
00120 Scalar Rv = underSaturatedTable.get("RV" , innerIdx);
00121 Scalar Bg = underSaturatedTable.get("BG" , innerIdx);
00122 Scalar mug = underSaturatedTable.get("MUG" , innerIdx);
00123
00124 invGasB.appendSamplePoint(outerIdx, Rv, 1.0/Bg);
00125 gasMu.appendSamplePoint(outerIdx, Rv, mug);
00126 }
00127 }
00128
00129 {
00130 std::vector<double> tmpPressure = saturatedTable.getColumn("PG").vectorCopy( );
00131
00132 invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
00133 invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
00134 }
00135
00136
00137 for (unsigned xIdx = 0; xIdx < invGasB.numX(); ++xIdx) {
00138
00139 assert(invGasB.numY(xIdx) > 0);
00140
00141
00142
00143 if (invGasB.numY(xIdx) > 1)
00144 continue;
00145
00146
00147
00148
00149 size_t masterTableIdx = xIdx + 1;
00150 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
00151 {
00152 if (pvtgTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
00153 break;
00154 }
00155
00156 if (masterTableIdx >= saturatedTable.numRows())
00157 OPM_THROW(std::runtime_error,
00158 "PVTG tables are invalid: The last table must exhibit at least one "
00159 "entry for undersaturated gas!");
00160
00161
00162
00163 extendPvtgTable_(regionIdx,
00164 xIdx,
00165 pvtgTable.getUnderSaturatedTable(xIdx),
00166 pvtgTable.getUnderSaturatedTable(masterTableIdx));
00167 }
00168 }
00169
00170 vapPar1_ = 0.0;
00171 if (deck.hasKeyword("VAPPARS")) {
00172 const auto& vapParsKeyword = deck.getKeyword("VAPPARS");
00173 vapPar1_ = vapParsKeyword.getRecord(0).getItem("OIL_VAP_PROPENSITY").template get<double>(0);
00174 }
00175
00176 initEnd();
00177 }
00178
00179 private:
00180 void extendPvtgTable_(unsigned regionIdx,
00181 unsigned xIdx,
00182 const SimpleTable& curTable,
00183 const SimpleTable& masterTable)
00184 {
00185 std::vector<double> RvArray = curTable.getColumn("RV").vectorCopy();
00186 std::vector<double> gasBArray = curTable.getColumn("BG").vectorCopy();
00187 std::vector<double> gasMuArray = curTable.getColumn("MUG").vectorCopy();
00188
00189 auto& invGasB = inverseGasB_[regionIdx];
00190 auto& gasMu = gasMu_[regionIdx];
00191
00192 for (size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
00193 const auto& RVColumn = masterTable.getColumn("RV");
00194 const auto& BGColumn = masterTable.getColumn("BG");
00195 const auto& viscosityColumn = masterTable.getColumn("MUG");
00196
00197
00198 Scalar diffRv = RVColumn[newRowIdx] - RVColumn[newRowIdx - 1];
00199 Scalar newRv = RvArray.back() + diffRv;
00200
00201
00202 Scalar B1 = BGColumn[newRowIdx];
00203 Scalar B2 = BGColumn[newRowIdx - 1];
00204 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
00205
00206
00207
00208 Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
00209
00210
00211 Scalar mu1 = viscosityColumn[newRowIdx];
00212 Scalar mu2 = viscosityColumn[newRowIdx - 1];
00213 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
00214
00215
00216
00217 Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
00218
00219
00220
00221 RvArray.push_back(newRv);
00222 gasBArray.push_back(newBg);
00223 gasMuArray.push_back(newMug);
00224
00225
00226 invGasB.appendSamplePoint(xIdx, newRv, 1.0/newBg);
00227 gasMu.appendSamplePoint(xIdx, newRv, newMug);
00228 }
00229 }
00230
00231 public:
00232 #endif // HAVE_OPM_PARSER
00233
00234 void setNumRegions(size_t numRegions)
00235 {
00236 oilReferenceDensity_.resize(numRegions);
00237 gasReferenceDensity_.resize(numRegions);
00238 inverseGasB_.resize(numRegions);
00239 inverseGasBMu_.resize(numRegions);
00240 inverseSaturatedGasB_.resize(numRegions);
00241 inverseSaturatedGasBMu_.resize(numRegions);
00242 gasMu_.resize(numRegions);
00243 saturatedOilVaporizationFactorTable_.resize(numRegions);
00244 saturationPressure_.resize(numRegions);
00245 }
00246
00250 void setReferenceDensities(unsigned regionIdx,
00251 Scalar rhoRefOil,
00252 Scalar rhoRefGas,
00253 Scalar )
00254 {
00255 oilReferenceDensity_[regionIdx] = rhoRefOil;
00256 gasReferenceDensity_[regionIdx] = rhoRefGas;
00257 }
00258
00264 void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
00265 { saturatedOilVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
00266
00276 void setSaturatedGasFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
00277 {
00278 auto& invGasB = inverseGasB_[regionIdx];
00279
00280 const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
00281
00282 Scalar T = 273.15 + 15.56;
00283
00284 Scalar RvMin = 0.0;
00285 Scalar RvMax = RvTable.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), true);
00286
00287 Scalar poMin = samplePoints.front().first;
00288 Scalar poMax = samplePoints.back().first;
00289
00290 size_t nRv = 20;
00291 size_t nP = samplePoints.size()*2;
00292
00293 Scalar rhogRef = gasReferenceDensity_[regionIdx];
00294 Scalar rhooRef = oilReferenceDensity_[regionIdx];
00295
00296 TabulatedOneDFunction gasFormationVolumeFactor;
00297 gasFormationVolumeFactor.setContainerOfTuples(samplePoints);
00298
00299 updateSaturationPressure_(regionIdx);
00300
00301
00302
00303
00304
00305 for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
00306 Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
00307
00308 invGasB.appendXPos(Rv);
00309
00310 for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
00311 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
00312
00313 Scalar poSat = saturationPressure(regionIdx, T, Rv);
00314 Scalar BgSat = gasFormationVolumeFactor.eval(poSat, true);
00315 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
00316 Scalar rhoo = rhooRef/BgSat*(1 + drhoo_dp*(pg - poSat));
00317
00318 Scalar Bg = rhooRef/rhoo;
00319
00320 invGasB.appendSamplePoint(RvIdx, pg, 1.0/Bg);
00321 }
00322 }
00323 }
00324
00337 void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction& invBg)
00338 { inverseGasB_[regionIdx] = invBg; }
00339
00345 void setGasViscosity(unsigned regionIdx, const TabulatedTwoDFunction& mug)
00346 { gasMu_[regionIdx] = mug; }
00347
00355 void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints& samplePoints )
00356 {
00357 auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
00358
00359 Scalar RvMin = 0.0;
00360 Scalar RvMax = oilVaporizationFac.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(), true);
00361
00362 Scalar poMin = samplePoints.front().first;
00363 Scalar poMax = samplePoints.back().first;
00364
00365 size_t nRv = 20;
00366 size_t nP = samplePoints.size()*2;
00367
00368 TabulatedOneDFunction mugTable;
00369 mugTable.setContainerOfTuples(samplePoints);
00370
00371
00372
00373 for (size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
00374 Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
00375
00376 gasMu_[regionIdx].appendXPos(Rv);
00377
00378 for (size_t pIdx = 0; pIdx < nP; ++pIdx) {
00379 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
00380 Scalar mug = mugTable.eval(pg, true);
00381
00382 gasMu_[regionIdx].appendSamplePoint(RvIdx, pg, mug);
00383 }
00384 }
00385 }
00386
00390 void initEnd()
00391 {
00392
00393 size_t numRegions = gasMu_.size();
00394 for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
00395
00396
00397 const auto& gasMu = gasMu_[regionIdx];
00398 const auto& invGasB = inverseGasB_[regionIdx];
00399 assert(gasMu.numX() == invGasB.numX());
00400
00401 auto& invGasBMu = inverseGasBMu_[regionIdx];
00402 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
00403 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
00404
00405 std::vector<Scalar> satPressuresArray;
00406 std::vector<Scalar> invSatGasBArray;
00407 std::vector<Scalar> invSatGasBMuArray;
00408 for (size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
00409 invGasBMu.appendXPos(gasMu.xAt(pIdx));
00410
00411 assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
00412
00413 size_t numRv = gasMu.numY(pIdx);
00414 for (size_t rvIdx = 0; rvIdx < numRv; ++rvIdx)
00415 invGasBMu.appendSamplePoint(pIdx,
00416 gasMu.yAt(pIdx, rvIdx),
00417 invGasB.valueAt(pIdx, rvIdx)
00418 / gasMu.valueAt(pIdx, rvIdx));
00419
00420
00421
00422
00423 satPressuresArray.push_back(gasMu.xAt(pIdx));
00424 invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRv - 1));
00425 invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRv - 1));
00426 }
00427
00428 invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
00429 invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
00430
00431 updateSaturationPressure_(regionIdx);
00432 }
00433 }
00434
00438 unsigned numRegions() const
00439 { return gasReferenceDensity_.size(); }
00440
00444 template <class Evaluation>
00445 Evaluation viscosity(unsigned regionIdx,
00446 const Evaluation& ,
00447 const Evaluation& pressure,
00448 const Evaluation& Rv) const
00449 {
00450 const Evaluation& invBg = inverseGasB_[regionIdx].eval(pressure, Rv, true);
00451 const Evaluation& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rv, true);
00452
00453 return invBg/invMugBg;
00454 }
00455
00459 template <class Evaluation>
00460 Evaluation saturatedViscosity(unsigned regionIdx,
00461 const Evaluation& ,
00462 const Evaluation& pressure) const
00463 {
00464 const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure, true);
00465 const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure, true);
00466
00467 return invBg/invMugBg;
00468 }
00469
00473 template <class Evaluation>
00474 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
00475 const Evaluation& ,
00476 const Evaluation& pressure,
00477 const Evaluation& Rv) const
00478 { return inverseGasB_[regionIdx].eval(pressure, Rv, true); }
00479
00483 template <class Evaluation>
00484 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
00485 const Evaluation& ,
00486 const Evaluation& pressure) const
00487 { return inverseSaturatedGasB_[regionIdx].eval(pressure, true); }
00488
00492 template <class Evaluation>
00493 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
00494 const Evaluation& ,
00495 const Evaluation& pressure) const
00496 {
00497 return saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, true);
00498 }
00499
00507 template <class Evaluation>
00508 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
00509 const Evaluation& ,
00510 const Evaluation& pressure,
00511 const Evaluation& oilSaturation,
00512 Scalar maxOilSaturation) const
00513 {
00514 Evaluation tmp =
00515 saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, true);
00516
00517
00518
00519 maxOilSaturation = std::min(maxOilSaturation, Scalar(1.0));
00520 if (vapPar1_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
00521 static const Scalar eps = 0.001;
00522 const Evaluation& So = Opm::max(oilSaturation, eps);
00523 tmp *= Opm::max(1e-3, Opm::pow(So/maxOilSaturation, vapPar1_));
00524 }
00525
00526 return tmp;
00527 }
00528
00539 template <class Evaluation>
00540 Evaluation saturationPressure(unsigned regionIdx,
00541 const Evaluation& temperature OPM_UNUSED,
00542 const Evaluation& Rv) const
00543 {
00544 typedef Opm::MathToolbox<Evaluation> Toolbox;
00545
00546 const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
00547 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
00548
00549
00550 Evaluation pSat = saturationPressure_[regionIdx].eval(Rv, true);
00551
00552
00553
00554
00555 bool onProbation = false;
00556 for (unsigned i = 0; i < 20; ++i) {
00557 const Evaluation& f = RvTable.eval(pSat, true) - Rv;
00558 const Evaluation& fPrime = RvTable.evalDerivative(pSat, true);
00559
00560
00561
00562 if (std::abs(Opm::scalarValue(fPrime)) < 1.0e-30) {
00563 return pSat;
00564 }
00565
00566 const Evaluation& delta = f/fPrime;
00567
00568 pSat -= delta;
00569
00570 if (pSat < 0.0) {
00571
00572
00573 if (onProbation)
00574 return 0.0;
00575
00576 onProbation = true;
00577 pSat = 0.0;
00578 }
00579
00580 if (std::abs(Opm::scalarValue(delta)) < std::abs(Opm::scalarValue(pSat))*eps)
00581 return pSat;
00582 }
00583
00584 std::stringstream errlog;
00585 errlog << "Finding saturation pressure did not converge:"
00586 << " pSat = " << pSat
00587 << ", Rv = " << Rv;
00588 OpmLog::debug("Wet gas saturation pressure", errlog.str());
00589 OPM_THROW_NOLOG(NumericalProblem, errlog.str());
00590 }
00591
00592 private:
00593 void updateSaturationPressure_(unsigned regionIdx)
00594 {
00595 typedef std::pair<Scalar, Scalar> Pair;
00596 const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
00597
00598
00599
00600 size_t n = oilVaporizationFac.numSamples();
00601 Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
00602
00603 SamplingPoints pSatSamplePoints;
00604 Scalar Rv = 0;
00605 for (size_t i = 0; i <= n; ++ i) {
00606 Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
00607 Rv = saturatedOilVaporizationFactor(regionIdx, Scalar(1e30), pSat);
00608
00609 Pair val(Rv, pSat);
00610 pSatSamplePoints.push_back(val);
00611 }
00612
00613
00614 auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
00615 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
00616 pSatSamplePoints.erase(last, pSatSamplePoints.end());
00617
00618 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
00619 }
00620
00621 std::vector<Scalar> gasReferenceDensity_;
00622 std::vector<Scalar> oilReferenceDensity_;
00623 std::vector<TabulatedTwoDFunction> inverseGasB_;
00624 std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
00625 std::vector<TabulatedTwoDFunction> gasMu_;
00626 std::vector<TabulatedTwoDFunction> inverseGasBMu_;
00627 std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
00628 std::vector<TabulatedOneDFunction> saturatedOilVaporizationFactorTable_;
00629 std::vector<TabulatedOneDFunction> saturationPressure_;
00630
00631 Scalar vapPar1_;
00632 };
00633
00634 }
00635
00636 #endif