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_LIVE_OIL_PVT_HPP
00028 #define OPM_LIVE_OIL_PVT_HPP
00029
00030 #include <opm/material/Constants.hpp>
00031
00032 #include <opm/material/common/OpmFinal.hpp>
00033 #include <opm/material/common/UniformXTabulated2DFunction.hpp>
00034 #include <opm/material/common/Tabulated1DFunction.hpp>
00035
00036 #if HAVE_OPM_PARSER
00037 #include <opm/parser/eclipse/Deck/Deck.hpp>
00038 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00039 #include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
00040 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
00041 #endif
00042
00043 namespace Opm {
00048 template <class Scalar>
00049 class LiveOilPvt
00050 {
00051 typedef Opm::UniformXTabulated2DFunction<Scalar> TabulatedTwoDFunction;
00052 typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
00053 typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
00054
00055 public:
00056 LiveOilPvt()
00057 {
00058 vapPar2_ = 0.0;
00059 }
00060
00061 #if HAVE_OPM_PARSER
00062
00065 void initFromDeck(const Deck& deck, const EclipseState& eclState)
00066 {
00067 const auto& pvtoTables = eclState.getTableManager().getPvtoTables();
00068 const auto& densityKeyword = deck.getKeyword("DENSITY");
00069
00070 assert(pvtoTables.size() == densityKeyword.size());
00071
00072 size_t numRegions = pvtoTables.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
00084 for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
00085 const auto& pvtoTable = pvtoTables[regionIdx];
00086
00087 const auto& saturatedTable = pvtoTable.getSaturatedTable();
00088 assert(saturatedTable.numRows() > 1);
00089
00090 auto& oilMu = oilMuTable_[regionIdx];
00091 auto& satOilMu = saturatedOilMuTable_[regionIdx];
00092 auto& invOilB = inverseOilBTable_[regionIdx];
00093 auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
00094 auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
00095 std::vector<Scalar> invSatOilBArray;
00096 std::vector<Scalar> satOilMuArray;
00097
00098
00099 for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
00100 Scalar Rs = saturatedTable.get("RS", outerIdx);
00101 Scalar BoSat = saturatedTable.get("BO", outerIdx);
00102 Scalar muoSat = saturatedTable.get("MU", outerIdx);
00103
00104 satOilMuArray.push_back(muoSat);
00105 invSatOilBArray.push_back(1.0/BoSat);
00106
00107 invOilB.appendXPos(Rs);
00108 oilMu.appendXPos(Rs);
00109
00110 assert(invOilB.numX() == outerIdx + 1);
00111 assert(oilMu.numX() == outerIdx + 1);
00112
00113 const auto& underSaturatedTable = pvtoTable.getUnderSaturatedTable(outerIdx);
00114 size_t numRows = underSaturatedTable.numRows();
00115 for (unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
00116 Scalar po = underSaturatedTable.get("P", innerIdx);
00117 Scalar Bo = underSaturatedTable.get("BO", innerIdx);
00118 Scalar muo = underSaturatedTable.get("MU", innerIdx);
00119
00120 invOilB.appendSamplePoint(outerIdx, po, 1.0/Bo);
00121 oilMu.appendSamplePoint(outerIdx, po, muo);
00122 }
00123 }
00124
00125
00126
00127 {
00128 const auto& tmpPressureColumn = saturatedTable.getColumn("P");
00129 const auto& tmpGasSolubilityColumn = saturatedTable.getColumn("RS");
00130
00131 invSatOilB.setXYContainers(tmpPressureColumn, invSatOilBArray);
00132 satOilMu.setXYContainers(tmpPressureColumn, satOilMuArray);
00133 gasDissolutionFac.setXYContainers(tmpPressureColumn, tmpGasSolubilityColumn);
00134 }
00135
00136 updateSaturationPressure_(regionIdx);
00137
00138 for (unsigned xIdx = 0; xIdx < invOilB.numX(); ++xIdx) {
00139
00140 assert(invOilB.numY(xIdx) > 0);
00141
00142
00143
00144 if (invOilB.numY(xIdx) > 1)
00145 continue;
00146
00147
00148
00149
00150 size_t masterTableIdx = xIdx + 1;
00151 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
00152 {
00153 if (pvtoTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
00154 break;
00155 }
00156
00157 if (masterTableIdx >= saturatedTable.numRows())
00158 OPM_THROW(std::runtime_error,
00159 "PVTO tables are invalid: The last table must exhibit at least one "
00160 "entry for undersaturated oil!");
00161
00162
00163 extendPvtoTable_(regionIdx,
00164 xIdx,
00165 pvtoTable.getUnderSaturatedTable(xIdx),
00166 pvtoTable.getUnderSaturatedTable(masterTableIdx));
00167 }
00168 }
00169
00170 vapPar2_ = 0.0;
00171 if (deck.hasKeyword("VAPPARS")) {
00172 const auto& vapParsKeyword = deck.getKeyword("VAPPARS");
00173 vapPar2_ = vapParsKeyword.getRecord(0).getItem("OIL_DENSITY_PROPENSITY").template get<double>(0);
00174 }
00175
00176 initEnd();
00177 }
00178
00179 private:
00180 void extendPvtoTable_(unsigned regionIdx,
00181 unsigned xIdx,
00182 const SimpleTable& curTable,
00183 const SimpleTable& masterTable)
00184 {
00185 std::vector<double> pressuresArray = curTable.getColumn("P").vectorCopy();
00186 std::vector<double> oilBArray = curTable.getColumn("BO").vectorCopy();
00187 std::vector<double> oilMuArray = curTable.getColumn("MU").vectorCopy();
00188
00189 auto& invOilB = inverseOilBTable_[regionIdx];
00190 auto& oilMu = oilMuTable_[regionIdx];
00191
00192 for (unsigned newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
00193 const auto& pressureColumn = masterTable.getColumn("P");
00194 const auto& BOColumn = masterTable.getColumn("BO");
00195 const auto& viscosityColumn = masterTable.getColumn("MU");
00196
00197
00198 Scalar diffPo = pressureColumn[newRowIdx] - pressureColumn[newRowIdx - 1];
00199 Scalar newPo = pressuresArray.back() + diffPo;
00200
00201
00202 Scalar B1 = BOColumn[newRowIdx];
00203 Scalar B2 = BOColumn[newRowIdx - 1];
00204 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
00205
00206
00207
00208 Scalar newBo = oilBArray.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 newMuo = oilMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
00218
00219
00220
00221 pressuresArray.push_back(newPo);
00222 oilBArray.push_back(newBo);
00223 oilMuArray.push_back(newMuo);
00224
00225
00226 invOilB.appendSamplePoint(xIdx, newPo, 1.0/newBo);
00227 oilMu.appendSamplePoint(xIdx, newPo, newMuo);
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 inverseOilBTable_.resize(numRegions);
00239 inverseOilBMuTable_.resize(numRegions);
00240 inverseSaturatedOilBTable_.resize(numRegions);
00241 inverseSaturatedOilBMuTable_.resize(numRegions);
00242 oilMuTable_.resize(numRegions);
00243 saturatedOilMuTable_.resize(numRegions);
00244 saturatedGasDissolutionFactorTable_.resize(numRegions);
00245 saturationPressure_.resize(numRegions);
00246 }
00247
00251 void setReferenceDensities(unsigned regionIdx,
00252 Scalar rhoRefOil,
00253 Scalar rhoRefGas,
00254 Scalar )
00255 {
00256 oilReferenceDensity_[regionIdx] = rhoRefOil;
00257 gasReferenceDensity_[regionIdx] = rhoRefGas;
00258 }
00259
00265 void setSaturatedOilGasDissolutionFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
00266 { saturatedGasDissolutionFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
00267
00277 void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints& samplePoints)
00278 {
00279 Scalar T = 273.15 + 15.56;
00280 auto& invOilB = inverseOilBTable_[regionIdx];
00281
00282 updateSaturationPressure_(regionIdx);
00283
00284
00285 for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
00286 Scalar p1 = std::get<0>(samplePoints[pIdx]);
00287 Scalar p2 = p1 * 2.0;
00288
00289 Scalar Bo1 = std::get<1>(samplePoints[pIdx]);
00290 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
00291 Scalar Bo2 = Bo1/(1.0 + (p2 - p1)*drhoo_dp);
00292
00293 Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
00294
00295 invOilB.appendXPos(Rs);
00296 invOilB.appendSamplePoint(pIdx, p1, 1.0/Bo1);
00297 invOilB.appendSamplePoint(pIdx, p2, 1.0/Bo2);
00298 }
00299 }
00300
00313 void setInverseOilFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction& invBo)
00314 { inverseOilBTable_[regionIdx] = invBo; }
00315
00321 void setOilViscosity(unsigned regionIdx, const TabulatedTwoDFunction& muo)
00322 { oilMuTable_[regionIdx] = muo; }
00323
00331 void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints& samplePoints)
00332 {
00333 Scalar T = 273.15 + 15.56;
00334
00335
00336 saturatedOilMuTable_[regionIdx].setContainerOfTuples(samplePoints);
00337
00338
00339
00340 for (size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
00341 Scalar p1 = std::get<0>(samplePoints[pIdx]);
00342 Scalar p2 = p1 * 2.0;
00343
00344
00345 Scalar mu1 = std::get<1>(samplePoints[pIdx]);
00346 Scalar mu2 = mu1;
00347
00348 Scalar Rs = saturatedGasDissolutionFactor(regionIdx, T, p1);
00349
00350 oilMuTable_[regionIdx].appendXPos(Rs);
00351 oilMuTable_[regionIdx].appendSamplePoint(pIdx, p1, mu1);
00352 oilMuTable_[regionIdx].appendSamplePoint(pIdx, p2, mu2);
00353 }
00354 }
00355
00359 void initEnd()
00360 {
00361
00362 size_t numRegions = oilMuTable_.size();
00363 for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
00364
00365
00366 const auto& oilMu = oilMuTable_[regionIdx];
00367 const auto& satOilMu = saturatedOilMuTable_[regionIdx];
00368 const auto& invOilB = inverseOilBTable_[regionIdx];
00369 assert(oilMu.numX() == invOilB.numX());
00370
00371 auto& invOilBMu = inverseOilBMuTable_[regionIdx];
00372 auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
00373 auto& invSatOilBMu = inverseSaturatedOilBMuTable_[regionIdx];
00374
00375 std::vector<Scalar> satPressuresArray;
00376 std::vector<Scalar> invSatOilBArray;
00377 std::vector<Scalar> invSatOilBMuArray;
00378 for (unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
00379 invOilBMu.appendXPos(oilMu.xAt(rsIdx));
00380
00381 assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
00382
00383 size_t numPressures = oilMu.numY(rsIdx);
00384 for (unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
00385 invOilBMu.appendSamplePoint(rsIdx,
00386 oilMu.yAt(rsIdx, pIdx),
00387 invOilB.valueAt(rsIdx, pIdx)
00388 / oilMu.valueAt(rsIdx, pIdx));
00389
00390
00391
00392
00393 satPressuresArray.push_back(oilMu.yAt(rsIdx, 0));
00394 invSatOilBArray.push_back(invOilB.valueAt(rsIdx, 0));
00395 invSatOilBMuArray.push_back(invSatOilBArray.back()/satOilMu.valueAt(rsIdx));
00396 }
00397
00398 invSatOilB.setXYContainers(satPressuresArray, invSatOilBArray);
00399 invSatOilBMu.setXYContainers(satPressuresArray, invSatOilBMuArray);
00400
00401 updateSaturationPressure_(regionIdx);
00402 }
00403 }
00404
00408 unsigned numRegions() const
00409 { return inverseOilBMuTable_.size(); }
00410
00414 template <class Evaluation>
00415 Evaluation viscosity(unsigned regionIdx,
00416 const Evaluation& ,
00417 const Evaluation& pressure,
00418 const Evaluation& Rs) const
00419 {
00420
00421 const Evaluation& invBo = inverseOilBTable_[regionIdx].eval(Rs, pressure, true);
00422 const Evaluation& invMuoBo = inverseOilBMuTable_[regionIdx].eval(Rs, pressure, true);
00423
00424 return invBo/invMuoBo;
00425 }
00426
00430 template <class Evaluation>
00431 Evaluation saturatedViscosity(unsigned regionIdx,
00432 const Evaluation& ,
00433 const Evaluation& pressure) const
00434 {
00435
00436 const Evaluation& invBo = inverseSaturatedOilBTable_[regionIdx].eval(pressure, true);
00437 const Evaluation& invMuoBo = inverseSaturatedOilBMuTable_[regionIdx].eval(pressure, true);
00438
00439 return invBo/invMuoBo;
00440 }
00441
00445 template <class Evaluation>
00446 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
00447 const Evaluation& ,
00448 const Evaluation& pressure,
00449 const Evaluation& Rs) const
00450 {
00451
00452 return inverseOilBTable_[regionIdx].eval(Rs, pressure, true);
00453 }
00454
00458 template <class Evaluation>
00459 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
00460 const Evaluation& ,
00461 const Evaluation& pressure) const
00462 {
00463
00464 return inverseSaturatedOilBTable_[regionIdx].eval(pressure, true);
00465 }
00466
00470 template <class Evaluation>
00471 Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
00472 const Evaluation& ,
00473 const Evaluation& pressure) const
00474 { return saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure, true); }
00475
00483 template <class Evaluation>
00484 Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
00485 const Evaluation& ,
00486 const Evaluation& pressure,
00487 const Evaluation& oilSaturation,
00488 Scalar maxOilSaturation) const
00489 {
00490 Evaluation tmp =
00491 saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure, true);
00492
00493
00494
00495 maxOilSaturation = std::min(maxOilSaturation, Scalar(1.0));
00496 if (vapPar2_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
00497 static const Scalar eps = 0.001;
00498 const Evaluation& So = Opm::max(oilSaturation, eps);
00499 tmp *= Opm::max(1e-3, Opm::pow(So/maxOilSaturation, vapPar2_));
00500 }
00501
00502 return tmp;
00503 }
00504
00511 template <class Evaluation>
00512 Evaluation saturationPressure(unsigned regionIdx,
00513 const Evaluation& temperature OPM_UNUSED,
00514 const Evaluation& Rs) const
00515 {
00516 typedef Opm::MathToolbox<Evaluation> Toolbox;
00517
00518 const auto& RsTable = saturatedGasDissolutionFactorTable_[regionIdx];
00519 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
00520
00521
00522 Evaluation pSat = saturationPressure_[regionIdx].eval(Rs, true);
00523
00524
00525
00526
00527 bool onProbation = false;
00528 for (int i = 0; i < 20; ++i) {
00529 const Evaluation& f = RsTable.eval(pSat, true) - Rs;
00530 const Evaluation& fPrime = RsTable.evalDerivative(pSat, true);
00531
00532
00533
00534 if (std::abs(Opm::scalarValue(fPrime)) < 1.0e-30) {
00535 return pSat;
00536 }
00537
00538 const Evaluation& delta = f/fPrime;
00539
00540 pSat -= delta;
00541
00542 if (pSat < 0.0) {
00543
00544
00545 if (onProbation)
00546 return 0.0;
00547
00548 onProbation = true;
00549 pSat = 0.0;
00550 }
00551
00552 if (std::abs(Opm::scalarValue(delta)) < std::abs(Opm::scalarValue(pSat))*eps)
00553 return pSat;
00554 }
00555
00556 std::stringstream errlog;
00557 errlog << "Finding saturation pressure did not converge:"
00558 << " pSat = " << pSat
00559 << ", Rs = " << Rs;
00560 OpmLog::debug("Live oil saturation pressure", errlog.str());
00561 OPM_THROW_NOLOG(NumericalProblem, errlog.str());
00562 }
00563
00564 private:
00565 void updateSaturationPressure_(unsigned regionIdx)
00566 {
00567 typedef std::pair<Scalar, Scalar> Pair;
00568 const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
00569
00570
00571
00572 size_t n = gasDissolutionFac.numSamples();
00573 Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/Scalar(n + 1);
00574
00575 SamplingPoints pSatSamplePoints;
00576 Scalar Rs = 0;
00577 for (size_t i=0; i <= n; ++ i) {
00578 Scalar pSat = gasDissolutionFac.xMin() + Scalar(i)*delta;
00579 Rs = saturatedGasDissolutionFactor(regionIdx,
00580 Scalar(1e30),
00581 pSat);
00582
00583 Pair val(Rs, pSat);
00584 pSatSamplePoints.push_back(val);
00585 }
00586
00587
00588 auto x_coord_comparator = [](const Pair& a, const Pair& b) { return a.first == b.first; };
00589 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
00590 pSatSamplePoints.erase(last, pSatSamplePoints.end());
00591
00592 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
00593 }
00594
00595 std::vector<Scalar> gasReferenceDensity_;
00596 std::vector<Scalar> oilReferenceDensity_;
00597 std::vector<TabulatedTwoDFunction> inverseOilBTable_;
00598 std::vector<TabulatedTwoDFunction> oilMuTable_;
00599 std::vector<TabulatedTwoDFunction> inverseOilBMuTable_;
00600 std::vector<TabulatedOneDFunction> saturatedOilMuTable_;
00601 std::vector<TabulatedOneDFunction> inverseSaturatedOilBTable_;
00602 std::vector<TabulatedOneDFunction> inverseSaturatedOilBMuTable_;
00603 std::vector<TabulatedOneDFunction> saturatedGasDissolutionFactorTable_;
00604 std::vector<TabulatedOneDFunction> saturationPressure_;
00605
00606 Scalar vapPar2_;
00607 };
00608
00609 }
00610
00611 #endif