25 #ifndef OPM_LIVE_OIL_PVT_HPP
26 #define OPM_LIVE_OIL_PVT_HPP
36 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
40 template <
class Scalar>
41 class GasPvtMultiplexer;
47 template <
class Scalar>
55 typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
62 void initFromDeck(DeckConstPtr deck, EclipseStateConstPtr eclState)
64 const auto& pvtoTables = eclState->getTableManager()->getPvtoTables();
65 DeckKeywordConstPtr densityKeyword = deck->getKeyword(
"DENSITY");
67 assert(pvtoTables.size() == densityKeyword->size());
69 size_t numRegions = pvtoTables.size();
72 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
73 Scalar rhoRefO = densityKeyword->getRecord(regionIdx)->getItem(
"OIL")->getSIDouble(0);
74 Scalar rhoRefG = densityKeyword->getRecord(regionIdx)->getItem(
"GAS")->getSIDouble(0);
75 Scalar rhoRefW = densityKeyword->getRecord(regionIdx)->getItem(
"WATER")->getSIDouble(0);
81 Scalar T = 273.15 + 15.56;
89 const auto& pvtoTable = pvtoTables[regionIdx];
91 const auto saturatedTable = pvtoTable.getOuterTable();
92 assert(saturatedTable->numRows() > 1);
94 auto& oilMu = oilMuTable_[regionIdx];
95 auto& invOilB = inverseOilBTable_[regionIdx];
96 auto& gasDissolutionFac = gasDissolutionFactorTable_[regionIdx];
98 gasDissolutionFac.setXYArrays(saturatedTable->numRows(),
99 saturatedTable->getPressureColumn(),
100 saturatedTable->getGasSolubilityColumn());
103 for (
unsigned outerIdx = 0; outerIdx < saturatedTable->numRows(); ++ outerIdx) {
104 Scalar Rs = saturatedTable->getGasSolubilityColumn()[outerIdx];
106 invOilB.appendXPos(Rs);
107 oilMu.appendXPos(Rs);
109 assert(invOilB.numX() == outerIdx + 1);
110 assert(oilMu.numX() == outerIdx + 1);
112 const auto underSaturatedTable = pvtoTable.getInnerTable(outerIdx);
113 size_t numRows = underSaturatedTable->numRows();
114 for (
unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
115 Scalar po = underSaturatedTable->getPressureColumn()[innerIdx];
116 Scalar Bo = underSaturatedTable->getOilFormationFactorColumn()[innerIdx];
117 Scalar muo = underSaturatedTable->getOilViscosityColumn()[innerIdx];
119 invOilB.appendSamplePoint(outerIdx, po, 1.0/Bo);
120 oilMu.appendSamplePoint(outerIdx, po, muo);
125 for (
unsigned xIdx = 0; xIdx < invOilB.numX(); ++xIdx) {
127 assert(invOilB.numY(xIdx) > 0);
131 if (invOilB.numY(xIdx) > 1)
137 size_t masterTableIdx = xIdx + 1;
138 for (; masterTableIdx < pvtoTable.getOuterTable()->numRows(); ++masterTableIdx)
140 if (pvtoTable.getInnerTable(masterTableIdx)->numRows() > 1)
144 if (masterTableIdx >= pvtoTable.getOuterTable()->numRows())
145 OPM_THROW(std::runtime_error,
146 "PVTO tables are invalid: The last table must exhibit at least one "
147 "entry for undersaturated oil!");
153 const auto masterTable = pvtoTable.getInnerTable(masterTableIdx);
154 const auto curTable = pvtoTable.getInnerTable(xIdx);
155 for (
unsigned newRowIdx = 1; newRowIdx < masterTable->numRows(); ++ newRowIdx) {
157 masterTable->getPressureColumn()[newRowIdx]
158 / masterTable->getPressureColumn()[0];
161 masterTable->getOilFormationFactorColumn()[newRowIdx]
162 / masterTable->getOilFormationFactorColumn()[0];
165 masterTable->getOilViscosityColumn()[newRowIdx]
166 / masterTable->getOilViscosityColumn()[0];
168 Scalar newPo = curTable->getPressureColumn()[0]*alphaPo;
169 Scalar newBo = curTable->getOilFormationFactorColumn()[0]*alphaBo;
170 Scalar newMuo = curTable->getOilViscosityColumn()[0]*alphaMuo;
172 invOilB.appendSamplePoint(xIdx, newPo, 1.0/newBo);
173 oilMu.appendSamplePoint(xIdx, newPo, newMuo);
178 #endif // HAVE_OPM_PARSER
182 oilMolarMass_.resize(numRegions);
183 gasMolarMass_.resize(numRegions);
184 oilReferenceDensity_.resize(numRegions);
185 gasReferenceDensity_.resize(numRegions);
186 inverseOilBTable_.resize(numRegions);
187 inverseOilBMuTable_.resize(numRegions);
188 oilMuTable_.resize(numRegions);
189 gasDissolutionFactorTable_.resize(numRegions);
190 saturationPressureSpline_.resize(numRegions);
201 oilReferenceDensity_[regionIdx] = rhoRefOil;
202 gasReferenceDensity_[regionIdx] = rhoRefGas;
213 oilMolarMass_[regionIdx] = MOil;
214 gasMolarMass_[regionIdx] = MGas;
224 { gasDissolutionFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
237 auto& invOilB = inverseOilBTable_[regionIdx];
239 auto &RsTable = gasDissolutionFactorTable_[regionIdx];
241 Scalar T = 273.15 + 15.56;
244 Scalar RsMax = RsTable.eval(gasDissolutionFactorTable_[regionIdx].xMax(),
true);
246 Scalar poMin = samplePoints.front().first;
247 Scalar poMax = samplePoints.back().first;
250 size_t nP = samplePoints.size()*2;
252 Scalar rhogRef = gasReferenceDensity_[regionIdx];
253 Scalar rhooRef = oilReferenceDensity_[regionIdx];
255 Spline oilFormationVolumeFactorSpline;
258 updateSaturationPressureSpline_(regionIdx);
262 for (
size_t RsIdx = 0; RsIdx < nRs; ++RsIdx) {
263 Scalar Rs = RsMin + (RsMax - RsMin)*RsIdx/nRs;
264 Scalar XoG = Rs/(rhooRef/rhogRef + Rs);
266 invOilB.appendXPos(Rs);
268 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
269 Scalar po = poMin + (poMax - poMin)*pIdx/nP;
272 Scalar BoSat = oilFormationVolumeFactorSpline.
eval(poSat,
true);
273 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
274 Scalar rhoo = oilReferenceDensity_[regionIdx]/BoSat*(1 + drhoo_dp*(po - poSat));
276 Scalar Bo = oilReferenceDensity_[regionIdx]/rhoo;
278 invOilB.appendSamplePoint(RsIdx, po, 1.0/Bo);
296 { inverseOilBTable_[regionIdx] = invBo; }
304 { oilMuTable_[regionIdx] = muo; }
315 auto& gasDissolutionFac = gasDissolutionFactorTable_[regionIdx];
318 Scalar RsMax = gasDissolutionFac.eval(gasDissolutionFactorTable_[regionIdx].xMax(),
true);
320 Scalar poMin = samplePoints.front().first;
321 Scalar poMax = samplePoints.back().first;
324 size_t nP = samplePoints.size()*2;
331 for (
size_t RsIdx = 0; RsIdx < nRs; ++RsIdx) {
332 Scalar Rs = RsMin + (RsMax - RsMin)*RsIdx/nRs;
334 oilMuTable_[regionIdx].appendXPos(Rs);
336 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
337 Scalar po = poMin + (poMax - poMin)*pIdx/nP;
338 Scalar muo = muoSpline.
eval(po,
true);
340 oilMuTable_[regionIdx].appendSamplePoint(RsIdx, po, muo);
353 size_t numRegions = oilMuTable_.size();
354 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
357 const auto& oilMu = oilMuTable_[regionIdx];
358 const auto& invOilB = inverseOilBTable_[regionIdx];
359 assert(oilMu.numX() == invOilB.numX());
361 auto& invOilBMu = inverseOilBMuTable_[regionIdx];
363 for (
unsigned rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) {
364 invOilBMu.appendXPos(oilMu.xAt(rsIdx));
366 assert(oilMu.numY(rsIdx) == invOilB.numY(rsIdx));
368 size_t numPressures = oilMu.numY(rsIdx);
369 for (
unsigned pIdx = 0; pIdx < numPressures; ++pIdx)
370 invOilBMu.appendSamplePoint(rsIdx,
371 oilMu.yAt(rsIdx, pIdx),
372 invOilB.valueAt(rsIdx, pIdx)*
373 1/oilMu.valueAt(rsIdx, pIdx));
376 updateSaturationPressureSpline_(regionIdx);
383 template <
class Evaluation>
386 const Evaluation& pressure,
387 const Evaluation& XoG)
const
389 const Evaluation& Rs =
390 XoG/(1 - XoG)*(oilReferenceDensity_[regionIdx]/gasReferenceDensity_[regionIdx]);
393 const Evaluation& invBo = inverseOilBTable_[regionIdx].eval(Rs, pressure,
true);
394 const Evaluation& invMuoBo = inverseOilBMuTable_[regionIdx].eval(Rs, pressure,
true);
396 return invBo/invMuoBo;
402 template <
class Evaluation>
404 const Evaluation& temperature,
405 const Evaluation& pressure,
406 const Evaluation& XoG)
const
408 Scalar rhooRef = oilReferenceDensity_[regionIdx];
409 Scalar rhogRef = gasReferenceDensity_[regionIdx];
416 Evaluation rhoo = rhooRef/Bo;
421 const Evaluation Rs = XoG/(1 - XoG) * rhooRef/rhogRef;
422 rhoo += rhogRef*Rs/Bo;
430 template <
class Evaluation>
433 const Evaluation& pressure,
434 const Evaluation& XoG)
const
436 const Evaluation& Rs =
437 XoG/(1-XoG)*(oilReferenceDensity_[regionIdx]/gasReferenceDensity_[regionIdx]);
442 return 1.0 / inverseOilBTable_[regionIdx].eval(Rs, pressure,
true);
449 template <
class Evaluation>
452 const Evaluation& pressure)
const
457 return 20e3/pressure;
460 template <
class Evaluation>
462 const Evaluation& temperature,
463 const Evaluation& pressure)
const
472 template <
class Evaluation>
474 const Evaluation& temperature,
475 const Evaluation& pressure)
const
489 return phi_gG / x_oGSat;
495 template <
class Evaluation>
498 const Evaluation& pressure)
const
499 {
return gasDissolutionFactorTable_[regionIdx].eval(pressure,
true); }
507 template <
class Evaluation>
509 const Evaluation& temperature,
510 const Evaluation& XoG)
const
515 Evaluation pSat = saturationPressureSpline_[regionIdx].eval(XoG,
true);
516 Evaluation eps = pSat*1e-11;
521 for (
int i = 0; i < 20; ++i) {
525 const Evaluation& delta = f/fPrime;
528 Scalar absDelta =
std::abs(Toolbox::value(delta));
529 if (absDelta < Toolbox::value(pSat) * 1e-10 || absDelta < 1e-4)
533 OPM_THROW(NumericalProblem,
"Could find the oil saturation pressure for X_o^G = " << XoG);
536 template <
class Evaluation>
538 const Evaluation& temperature,
539 const Evaluation& pressure)
const
541 Scalar rho_gRef = gasReferenceDensity_[regionIdx];
542 Scalar rho_oRef = oilReferenceDensity_[regionIdx];
551 return rho_oG/(rho_oRef + rho_oG);
554 template <
class Evaluation>
556 const Evaluation& temperature,
557 const Evaluation& pressure)
const
564 Scalar MG = gasMolarMass_[regionIdx];
565 Scalar MO = oilMolarMass_[regionIdx];
567 Evaluation avgMolarMass = MO/(1 + XoG*(MO/MG - 1));
568 return XoG*avgMolarMass/MG;
572 void updateSaturationPressureSpline_(
unsigned regionIdx)
574 auto& gasDissolutionFac = gasDissolutionFactorTable_[regionIdx];
578 size_t n = gasDissolutionFac.numSamples()*5;
579 Scalar delta = (gasDissolutionFac.xMax() - gasDissolutionFac.xMin())/(n + 1);
581 SamplingPoints pSatSamplePoints;
583 for (
size_t i=0; i <= n; ++ i) {
584 Scalar pSat = gasDissolutionFac.xMin() + i*delta;
589 std::pair<Scalar, Scalar> val(XoG, pSat);
590 pSatSamplePoints.push_back(val);
592 saturationPressureSpline_[regionIdx].setContainerOfTuples(pSatSamplePoints,
596 const GasPvtMultiplexer *gasPvt_;
598 std::vector<Scalar> gasMolarMass_;
599 std::vector<Scalar> oilMolarMass_;
600 std::vector<Scalar> gasReferenceDensity_;
601 std::vector<Scalar> oilReferenceDensity_;
602 std::vector<TabulatedTwoDFunction> inverseOilBTable_;
603 std::vector<TabulatedTwoDFunction> oilMuTable_;
604 std::vector<TabulatedTwoDFunction> inverseOilBMuTable_;
605 std::vector<TabulatedOneDFunction> gasDissolutionFactorTable_;
606 std::vector<Spline> saturationPressureSpline_;
Definition: Spline.hpp:104
bool CheckDefined(const T &value OPM_UNUSED)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition: Valgrind.hpp:74
This class represents the Pressure-Volume-Temperature relations of the oil phas with dissolved gas...
Definition: LiveOilPvt.hpp:48
Definition: Air_Mesitylene.hpp:31
void setSaturatedOilViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for gas saturated oil.
Definition: LiveOilPvt.hpp:313
Evaluation fugacityCoefficientGas(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Definition: LiveOilPvt.hpp:473
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: LiveOilPvt.hpp:196
void setContainerOfTuples(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=false)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:474
Evaluation fugacityCoefficientGas(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the fugacity coefficient [Pa] of the gas component in the gas phase given a set of parameters...
Definition: GasPvtMultiplexer.hpp:175
void setOilViscosity(unsigned regionIdx, const TabulatedTwoDFunction &muo)
Initialize the viscosity of the oil phase.
Definition: LiveOilPvt.hpp:303
Evaluation formationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &XoG) const
Returns the formation volume factor [-] of the fluid phase.
Definition: LiveOilPvt.hpp:431
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition: ConstantCompressibilityOilPvt.hpp:39
Evaluation fugacityCoefficientOil(unsigned, const Evaluation &, const Evaluation &pressure) const
Returns the fugacity coefficient [Pa] of a component in the fluid phase given a set of parameters...
Definition: LiveOilPvt.hpp:450
Class implementing cubic splines.
Evaluation saturatedOilGasMassFraction(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Definition: LiveOilPvt.hpp:537
void setMolarMasses(unsigned regionIdx, Scalar MOil, Scalar MGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: LiveOilPvt.hpp:208
Implements a linearly interpolated scalar function that depends on one variable.
Evaluation density(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &XoG) const
Returns the density [kg/m^3] of the fluid phase given a set of parameters.
Definition: LiveOilPvt.hpp:403
Evaluation fugacityCoefficientWater(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Definition: LiveOilPvt.hpp:461
A central place for various physical constants occuring in some equations.
void setSaturatedOilGasDissolutionFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the gas dissolution factor .
Definition: LiveOilPvt.hpp:223
This file provides a wrapper around the "final" C++-2011 statement.
void initEnd(const GasPvtMultiplexer *gasPvt)
Finish initializing the oil phase PVT properties.
Definition: LiveOilPvt.hpp:348
void setSaturatedOilFormationVolumeFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil formation volume factor.
Definition: LiveOilPvt.hpp:235
Evaluation< Scalar, VarSetTag, numVars > abs(const Evaluation< Scalar, VarSetTag, numVars > &)
Definition: Math.hpp:41
Class implementing cubic splines.
Definition: Spline.hpp:89
Evaluation oilSaturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &XoG) const
Returns the saturation pressure of the oil phase [Pa] depending on its mass fraction of the gas compo...
Definition: LiveOilPvt.hpp:508
Scalar eval(Scalar x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Spline.hpp:809
void setNumRegions(size_t numRegions)
Definition: LiveOilPvt.hpp:180
Evaluation gasDissolutionFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: LiveOilPvt.hpp:496
void setInverseOilFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBo)
Initialize the spline for the oil formation volume factor.
Definition: LiveOilPvt.hpp:295
Evaluation saturatedOilGasMoleFraction(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Definition: LiveOilPvt.hpp:555
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:44
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &XoG) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: LiveOilPvt.hpp:384
A central place for various physical constants occuring in some equations.
Definition: Constants.hpp:39