27 #ifndef OPM_LIVE_OIL_PVT_HPP
28 #define OPM_LIVE_OIL_PVT_HPP
37 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
38 #include <opm/parser/eclipse/EclipseState/Schedule/OilVaporizationProperties.hpp>
39 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
40 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
44 #include <opm/common/OpmLog/OpmLog.hpp>
52 template <
class Scalar>
55 typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
66 LiveOilPvt(
const std::vector<Scalar>& gasReferenceDensity,
67 const std::vector<Scalar>& oilReferenceDensity,
68 const std::vector<TabulatedTwoDFunction>& inverseOilBTable,
69 const std::vector<TabulatedTwoDFunction>& oilMuTable,
70 const std::vector<TabulatedTwoDFunction>& inverseOilBMuTable,
71 const std::vector<TabulatedOneDFunction>& saturatedOilMuTable,
72 const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBTable,
73 const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBMuTable,
74 const std::vector<TabulatedOneDFunction>& saturatedGasDissolutionFactorTable,
77 : gasReferenceDensity_(gasReferenceDensity)
78 , oilReferenceDensity_(oilReferenceDensity)
79 , inverseOilBTable_(inverseOilBTable)
80 , oilMuTable_(oilMuTable)
81 , inverseOilBMuTable_(inverseOilBMuTable)
82 , saturatedOilMuTable_(saturatedOilMuTable)
83 , inverseSaturatedOilBTable_(inverseSaturatedOilBTable)
84 , inverseSaturatedOilBMuTable_(inverseSaturatedOilBMuTable)
85 , saturatedGasDissolutionFactorTable_(saturatedGasDissolutionFactorTable)
86 , saturationPressure_(saturationPressure)
94 void initFromState(
const EclipseState& eclState,
const Schedule& schedule)
96 const auto& pvtoTables = eclState.getTableManager().getPvtoTables();
97 const auto& densityTable = eclState.getTableManager().getDensityTable();
99 assert(pvtoTables.size() == densityTable.size());
104 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
105 Scalar rhoRefO = densityTable[regionIdx].oil;
106 Scalar rhoRefG = densityTable[regionIdx].gas;
107 Scalar rhoRefW = densityTable[regionIdx].water;
113 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
114 const auto& pvtoTable = pvtoTables[regionIdx];
116 const auto& saturatedTable = pvtoTable.getSaturatedTable();
117 assert(saturatedTable.numRows() > 1);
119 auto& oilMu = oilMuTable_[regionIdx];
120 auto& satOilMu = saturatedOilMuTable_[regionIdx];
121 auto& invOilB = inverseOilBTable_[regionIdx];
122 auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
123 auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
124 std::vector<Scalar> invSatOilBArray;
125 std::vector<Scalar> satOilMuArray;
128 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
129 Scalar Rs = saturatedTable.get(
"RS", outerIdx);
130 Scalar BoSat = saturatedTable.get(
"BO", outerIdx);
131 Scalar muoSat = saturatedTable.get(
"MU", outerIdx);
133 satOilMuArray.push_back(muoSat);
134 invSatOilBArray.push_back(1.0/BoSat);
136 invOilB.appendXPos(Rs);
137 oilMu.appendXPos(Rs);
139 assert(invOilB.numX() == outerIdx + 1);
140 assert(oilMu.numX() == outerIdx + 1);
142 const auto& underSaturatedTable = pvtoTable.getUnderSaturatedTable(outerIdx);
143 size_t numRows = underSaturatedTable.numRows();
144 for (
unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
145 Scalar po = underSaturatedTable.get(
"P", innerIdx);
146 Scalar Bo = underSaturatedTable.get(
"BO", innerIdx);
147 Scalar muo = underSaturatedTable.get(
"MU", innerIdx);
149 invOilB.appendSamplePoint(outerIdx, po, 1.0/Bo);
150 oilMu.appendSamplePoint(outerIdx, po, muo);
157 const auto& tmpPressureColumn = saturatedTable.getColumn(
"P");
158 const auto& tmpGasSolubilityColumn = saturatedTable.getColumn(
"RS");
160 invSatOilB.setXYContainers(tmpPressureColumn, invSatOilBArray);
161 satOilMu.setXYContainers(tmpPressureColumn, satOilMuArray);
162 gasDissolutionFac.setXYContainers(tmpPressureColumn, tmpGasSolubilityColumn);
165 updateSaturationPressure_(regionIdx);
167 for (
unsigned xIdx = 0; xIdx < invOilB.numX(); ++xIdx) {
169 assert(invOilB.numY(xIdx) > 0);
173 if (invOilB.numY(xIdx) > 1)
179 size_t masterTableIdx = xIdx + 1;
180 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
182 if (pvtoTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
186 if (masterTableIdx >= saturatedTable.numRows())
187 throw std::runtime_error(
"PVTO tables are invalid: The last table must exhibit at least one "
188 "entry for undersaturated oil!");
191 extendPvtoTable_(regionIdx,
193 pvtoTable.getUnderSaturatedTable(xIdx),
194 pvtoTable.getUnderSaturatedTable(masterTableIdx));
199 const auto& oilVap = schedule[0].oilvap();
200 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
201 vapPar2_ = oilVap.vap2();
208 void extendPvtoTable_(
unsigned regionIdx,
210 const SimpleTable& curTable,
211 const SimpleTable& masterTable)
213 std::vector<double> pressuresArray = curTable.getColumn(
"P").vectorCopy();
214 std::vector<double> oilBArray = curTable.getColumn(
"BO").vectorCopy();
215 std::vector<double> oilMuArray = curTable.getColumn(
"MU").vectorCopy();
217 auto& invOilB = inverseOilBTable_[regionIdx];
218 auto& oilMu = oilMuTable_[regionIdx];
220 for (
unsigned newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
221 const auto& pressureColumn = masterTable.getColumn(
"P");
222 const auto& BOColumn = masterTable.getColumn(
"BO");
223 const auto& viscosityColumn = masterTable.getColumn(
"MU");
226 Scalar diffPo = pressureColumn[newRowIdx] - pressureColumn[newRowIdx - 1];
227 Scalar newPo = pressuresArray.back() + diffPo;
230 Scalar B1 = BOColumn[newRowIdx];
231 Scalar B2 = BOColumn[newRowIdx - 1];
232 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
236 Scalar newBo = oilBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
239 Scalar mu1 = viscosityColumn[newRowIdx];
240 Scalar mu2 = viscosityColumn[newRowIdx - 1];
241 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
245 Scalar newMuo = oilMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
249 pressuresArray.push_back(newPo);
250 oilBArray.push_back(newBo);
251 oilMuArray.push_back(newMuo);
254 invOilB.appendSamplePoint(xIdx, newPo, 1.0/newBo);
255 oilMu.appendSamplePoint(xIdx, newPo, newMuo);
268 inverseSaturatedOilBTable_.resize(
numRegions);
269 inverseSaturatedOilBMuTable_.resize(
numRegions);
272 saturatedGasDissolutionFactorTable_.resize(
numRegions);
284 oilReferenceDensity_[regionIdx] = rhoRefOil;
285 gasReferenceDensity_[regionIdx] = rhoRefGas;
294 { saturatedGasDissolutionFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
307 Scalar T = 273.15 + 15.56;
308 auto& invOilB = inverseOilBTable_[regionIdx];
310 updateSaturationPressure_(regionIdx);
313 for (
size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
314 Scalar p1 = std::get<0>(samplePoints[pIdx]);
315 Scalar p2 = p1 * 2.0;
317 Scalar Bo1 = std::get<1>(samplePoints[pIdx]);
318 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
319 Scalar Bo2 = Bo1/(1.0 + (p2 - p1)*drhoo_dp);
323 invOilB.appendXPos(Rs);
324 invOilB.appendSamplePoint(pIdx, p1, 1.0/Bo1);
325 invOilB.appendSamplePoint(pIdx, p2, 1.0/Bo2);
342 { inverseOilBTable_[regionIdx] = invBo; }
350 { oilMuTable_[regionIdx] = muo; }
361 Scalar T = 273.15 + 15.56;
364 saturatedOilMuTable_[regionIdx].setContainerOfTuples(samplePoints);
368 for (
size_t pIdx = 0; pIdx < samplePoints.size(); ++pIdx) {
369 Scalar p1 = std::get<0>(samplePoints[pIdx]);
370 Scalar p2 = p1 * 2.0;
373 Scalar mu1 = std::get<1>(samplePoints[pIdx]);
378 oilMuTable_[regionIdx].appendXPos(Rs);
379 oilMuTable_[regionIdx].appendSamplePoint(pIdx, p1, mu1);
380 oilMuTable_[regionIdx].appendSamplePoint(pIdx, p2, mu2);
391 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
394 const auto& oilMu = oilMuTable_[regionIdx];
395 const auto& satOilMu = saturatedOilMuTable_[regionIdx];
396 const auto& invOilB = inverseOilBTable_[regionIdx];
397 assert(oilMu.numX() == invOilB.numX());
399 auto& invOilBMu = inverseOilBMuTable_[regionIdx];
400 auto& invSatOilB = inverseSaturatedOilBTable_[regionIdx];
401 auto& invSatOilBMu = inverseSaturatedOilBMuTable_[regionIdx];
403 std::vector<Scalar> satPressuresArray;
404 std::vector<Scalar> invSatOilBArray;
405 std::vector<Scalar> invSatOilBMuArray;
406 for (
unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
407 invOilBMu.appendXPos(oilMu.xAt(rsIdx));
409 assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
411 size_t numPressures = oilMu.numY(rsIdx);
412 for (
unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
413 invOilBMu.appendSamplePoint(rsIdx,
414 oilMu.yAt(rsIdx, pIdx),
415 invOilB.valueAt(rsIdx, pIdx)
416 / oilMu.valueAt(rsIdx, pIdx));
421 satPressuresArray.push_back(oilMu.yAt(rsIdx, 0));
422 invSatOilBArray.push_back(invOilB.valueAt(rsIdx, 0));
423 invSatOilBMuArray.push_back(invSatOilBArray.back()/satOilMu.valueAt(rsIdx));
426 invSatOilB.setXYContainers(satPressuresArray, invSatOilBArray);
427 invSatOilBMu.setXYContainers(satPressuresArray, invSatOilBMuArray);
429 updateSaturationPressure_(regionIdx);
437 {
return inverseOilBMuTable_.size(); }
442 template <
class Evaluation>
446 const Evaluation&)
const
448 throw std::runtime_error(
"Requested the enthalpy of oil but the thermal option is not enabled");
454 template <
class Evaluation>
457 const Evaluation& pressure,
458 const Evaluation& Rs)
const
461 const Evaluation& invBo = inverseOilBTable_[regionIdx].eval(Rs, pressure,
true);
462 const Evaluation& invMuoBo = inverseOilBMuTable_[regionIdx].eval(Rs, pressure,
true);
464 return invBo/invMuoBo;
470 template <
class Evaluation>
473 const Evaluation& pressure)
const
476 const Evaluation& invBo = inverseSaturatedOilBTable_[regionIdx].eval(pressure,
true);
477 const Evaluation& invMuoBo = inverseSaturatedOilBMuTable_[regionIdx].eval(pressure,
true);
479 return invBo/invMuoBo;
485 template <
class Evaluation>
488 const Evaluation& pressure,
489 const Evaluation& Rs)
const
492 return inverseOilBTable_[regionIdx].eval(Rs, pressure,
true);
498 template <
class Evaluation>
501 const Evaluation& pressure)
const
504 return inverseSaturatedOilBTable_[regionIdx].eval(pressure,
true);
510 template <
class Evaluation>
513 const Evaluation& pressure)
const
514 {
return saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure,
true); }
523 template <
class Evaluation>
526 const Evaluation& pressure,
527 const Evaluation& oilSaturation,
528 Evaluation maxOilSaturation)
const
531 saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure,
true);
535 maxOilSaturation = min(maxOilSaturation, Scalar(1.0));
536 if (vapPar2_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
537 static const Scalar eps = 0.001;
538 const Evaluation& So = max(oilSaturation, eps);
539 tmp *= max(1e-3, pow(So/maxOilSaturation, vapPar2_));
551 template <
class Evaluation>
554 const Evaluation& Rs)
const
558 const auto& RsTable = saturatedGasDissolutionFactorTable_[regionIdx];
559 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
562 Evaluation pSat = saturationPressure_[regionIdx].eval(Rs,
true);
567 bool onProbation =
false;
568 for (
int i = 0; i < 20; ++i) {
569 const Evaluation& f = RsTable.eval(pSat,
true) - Rs;
570 const Evaluation& fPrime = RsTable.evalDerivative(pSat,
true);
574 if (std::abs(scalarValue(fPrime)) < 1.0e-30) {
578 const Evaluation& delta = f/fPrime;
592 if (std::abs(scalarValue(delta)) < std::abs(scalarValue(pSat))*eps)
596 std::stringstream errlog;
597 errlog <<
"Finding saturation pressure did not converge:"
598 <<
" pSat = " << pSat
601 OpmLog::debug(
"Live oil saturation pressure", errlog.str());
606 template <
class Evaluation>
607 Evaluation diffusionCoefficient(
const Evaluation& ,
611 throw std::runtime_error(
"Not implemented: The PVT model does not provide a diffusionCoefficient()");
613 const Scalar gasReferenceDensity(
unsigned regionIdx)
const
614 {
return gasReferenceDensity_[regionIdx]; }
616 const Scalar oilReferenceDensity(
unsigned regionIdx)
const
617 {
return oilReferenceDensity_[regionIdx]; }
619 const std::vector<TabulatedTwoDFunction>& inverseOilBTable()
const
620 {
return inverseOilBTable_; }
622 const std::vector<TabulatedTwoDFunction>& oilMuTable()
const
623 {
return oilMuTable_; }
625 const std::vector<TabulatedTwoDFunction>& inverseOilBMuTable()
const
626 {
return inverseOilBMuTable_; }
628 const std::vector<TabulatedOneDFunction>& saturatedOilMuTable()
const
629 {
return saturatedOilMuTable_; }
631 const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBTable()
const
632 {
return inverseSaturatedOilBTable_; }
634 const std::vector<TabulatedOneDFunction>& inverseSaturatedOilBMuTable()
const
635 {
return inverseSaturatedOilBMuTable_; }
637 const std::vector<TabulatedOneDFunction>& saturatedGasDissolutionFactorTable()
const
638 {
return saturatedGasDissolutionFactorTable_; }
640 const std::vector<TabulatedOneDFunction>& saturationPressure()
const
641 {
return saturationPressure_; }
643 Scalar vapPar2()
const
646 bool operator==(
const LiveOilPvt<Scalar>& data)
const
648 return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
649 this->oilReferenceDensity_ == data.oilReferenceDensity_ &&
650 this->inverseOilBTable() == data.inverseOilBTable() &&
651 this->oilMuTable() == data.oilMuTable() &&
652 this->inverseOilBMuTable() == data.inverseOilBMuTable() &&
653 this->saturatedOilMuTable() == data.saturatedOilMuTable() &&
654 this->inverseSaturatedOilBTable() == data.inverseSaturatedOilBTable() &&
655 this->inverseSaturatedOilBMuTable() == data.inverseSaturatedOilBMuTable() &&
656 this->saturatedGasDissolutionFactorTable() == data.saturatedGasDissolutionFactorTable() &&
657 this->vapPar2() == data.vapPar2();
661 void updateSaturationPressure_(
unsigned regionIdx)
663 typedef std::pair<Scalar, Scalar> Pair;
664 const auto& gasDissolutionFac = saturatedGasDissolutionFactorTable_[regionIdx];
668 size_t n = gasDissolutionFac.numSamples();
669 Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/Scalar(n + 1);
671 SamplingPoints pSatSamplePoints;
673 for (
size_t i=0; i <= n; ++ i) {
674 Scalar pSat = gasDissolutionFac.xMin() + Scalar(i)*delta;
680 pSatSamplePoints.push_back(val);
684 auto x_coord_comparator = [](
const Pair& a,
const Pair& b) {
return a.first == b.first; };
685 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
686 pSatSamplePoints.erase(last, pSatSamplePoints.end());
688 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
691 std::vector<Scalar> gasReferenceDensity_;
692 std::vector<Scalar> oilReferenceDensity_;
693 std::vector<TabulatedTwoDFunction> inverseOilBTable_;
694 std::vector<TabulatedTwoDFunction> oilMuTable_;
695 std::vector<TabulatedTwoDFunction> inverseOilBMuTable_;
696 std::vector<TabulatedOneDFunction> saturatedOilMuTable_;
697 std::vector<TabulatedOneDFunction> inverseSaturatedOilBTable_;
698 std::vector<TabulatedOneDFunction> inverseSaturatedOilBMuTable_;
699 std::vector<TabulatedOneDFunction> saturatedGasDissolutionFactorTable_;
700 std::vector<TabulatedOneDFunction> saturationPressure_;
A central place for various physical constants occuring in some equations.
This file provides a wrapper around the "final" C++-2011 statement.
Implements a linearly interpolated scalar function that depends on one variable.
This class represents the Pressure-Volume-Temperature relations of the oil phas with dissolved gas.
Definition: LiveOilPvt.hpp:54
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of oil given a set of parameters.
Definition: LiveOilPvt.hpp:443
void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:305
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rs) const
Returns the formation volume factor [-] of the fluid phase.
Definition: LiveOilPvt.hpp:486
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &oilSaturation, Evaluation maxOilSaturation) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: LiveOilPvt.hpp:524
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of the fluid phase.
Definition: LiveOilPvt.hpp:499
void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for gas saturated oil.
Definition: LiveOilPvt.hpp:359
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &, const Evaluation &Rs) const
Returns the saturation pressure of the oil phase [Pa] depending on its mass fraction of the gas compo...
Definition: LiveOilPvt.hpp:552
void setSaturatedOilGasDissolutionFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the gas dissolution factor .
Definition: LiveOilPvt.hpp:293
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rs) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: LiveOilPvt.hpp:455
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: LiveOilPvt.hpp:279
void setOilViscosity(unsigned regionIdx, const TabulatedTwoDFunction &muo)
Initialize the viscosity of the oil phase.
Definition: LiveOilPvt.hpp:349
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: LiveOilPvt.hpp:511
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: LiveOilPvt.hpp:471
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: LiveOilPvt.hpp:436
void initEnd()
Finish initializing the oil phase PVT properties.
Definition: LiveOilPvt.hpp:387
void setInverseOilFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBo)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:341
Definition: Exceptions.hpp:46
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47