27 #ifndef OPM_DRY_HUMID_GAS_PVT_HPP
28 #define OPM_DRY_HUMID_GAS_PVT_HPP
36 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
37 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
42 #include <opm/common/OpmLog/OpmLog.hpp>
51 template <
class Scalar>
54 using SamplingPoints = std::vector<std::pair<Scalar, Scalar>>;
66 const std::vector<Scalar>& waterReferenceDensity,
67 const std::vector<TabulatedTwoDFunction>& inverseGasB,
68 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB,
69 const std::vector<TabulatedTwoDFunction>& gasMu,
70 const std::vector<TabulatedTwoDFunction>& inverseGasBMu,
71 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu,
72 const std::vector<TabulatedOneDFunction>& saturatedWaterVaporizationFactorTable,
75 : gasReferenceDensity_(gasReferenceDensity)
76 , waterReferenceDensity_(waterReferenceDensity)
77 , inverseGasB_(inverseGasB)
78 , inverseSaturatedGasB_(inverseSaturatedGasB)
80 , inverseGasBMu_(inverseGasBMu)
81 , inverseSaturatedGasBMu_(inverseSaturatedGasBMu)
82 , saturatedWaterVaporizationFactorTable_(saturatedWaterVaporizationFactorTable)
83 , saturationPressure_(saturationPressure)
95 void initFromState(
const EclipseState& eclState,
const Schedule&)
97 const auto& pvtgwTables = eclState.getTableManager().getPvtgwTables();
98 const auto& densityTable = eclState.getTableManager().getDensityTable();
100 assert(pvtgwTables.size() == densityTable.size());
105 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
106 Scalar rhoRefO = densityTable[regionIdx].oil;
107 Scalar rhoRefG = densityTable[regionIdx].gas;
108 Scalar rhoRefW = densityTable[regionIdx].water;
113 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
114 const auto& pvtgwTable = pvtgwTables[regionIdx];
116 const auto& saturatedTable = pvtgwTable.getSaturatedTable();
117 assert(saturatedTable.numRows() > 1);
119 auto& gasMu = gasMu_[regionIdx];
120 auto& invGasB = inverseGasB_[regionIdx];
121 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
122 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
123 auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
125 waterVaporizationFac.setXYArrays(saturatedTable.numRows(),
126 saturatedTable.getColumn(
"PG"),
127 saturatedTable.getColumn(
"RW"));
129 std::vector<Scalar> invSatGasBArray;
130 std::vector<Scalar> invSatGasBMuArray;
133 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
134 Scalar pg = saturatedTable.get(
"PG" , outerIdx);
135 Scalar B = saturatedTable.get(
"BG" , outerIdx);
136 Scalar mu = saturatedTable.get(
"MUG" , outerIdx);
138 invGasB.appendXPos(pg);
139 gasMu.appendXPos(pg);
141 invSatGasBArray.push_back(1.0/B);
142 invSatGasBMuArray.push_back(1.0/(mu*B));
144 assert(invGasB.numX() == outerIdx + 1);
145 assert(gasMu.numX() == outerIdx + 1);
147 const auto& underSaturatedTable = pvtgwTable.getUnderSaturatedTable(outerIdx);
148 size_t numRows = underSaturatedTable.numRows();
149 for (
size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
150 Scalar Rw = underSaturatedTable.get(
"RW" , innerIdx);
151 Scalar Bg = underSaturatedTable.get(
"BG" , innerIdx);
152 Scalar mug = underSaturatedTable.get(
"MUG" , innerIdx);
154 invGasB.appendSamplePoint(outerIdx, Rw, 1.0/Bg);
155 gasMu.appendSamplePoint(outerIdx, Rw, mug);
160 std::vector<double> tmpPressure = saturatedTable.getColumn(
"PG").vectorCopy( );
162 invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
163 invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
167 for (
unsigned xIdx = 0; xIdx < invGasB.numX(); ++xIdx) {
169 assert(invGasB.numY(xIdx) > 0);
173 if (invGasB.numY(xIdx) > 1)
179 size_t masterTableIdx = xIdx + 1;
180 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
182 if (pvtgwTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
186 if (masterTableIdx >= saturatedTable.numRows())
187 throw std::runtime_error(
"PVTGW tables are invalid: The last table must exhibit at least one "
188 "entry for undersaturated gas!");
192 extendPvtgwTable_(regionIdx,
194 pvtgwTable.getUnderSaturatedTable(xIdx),
195 pvtgwTable.getUnderSaturatedTable(masterTableIdx));
205 void extendPvtgwTable_(
unsigned regionIdx,
207 const SimpleTable& curTable,
208 const SimpleTable& masterTable)
210 std::vector<double> RwArray = curTable.getColumn(
"RW").vectorCopy();
211 std::vector<double> gasBArray = curTable.getColumn(
"BG").vectorCopy();
212 std::vector<double> gasMuArray = curTable.getColumn(
"MUG").vectorCopy();
214 auto& invGasB = inverseGasB_[regionIdx];
215 auto& gasMu = gasMu_[regionIdx];
217 for (
size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
218 const auto& RWColumn = masterTable.getColumn(
"RW");
219 const auto& BGColumn = masterTable.getColumn(
"BG");
220 const auto& viscosityColumn = masterTable.getColumn(
"MUG");
223 Scalar diffRw = RWColumn[newRowIdx] - RWColumn[newRowIdx - 1];
224 Scalar newRw = RwArray.back() + diffRw;
227 Scalar B1 = BGColumn[newRowIdx];
228 Scalar B2 = BGColumn[newRowIdx - 1];
229 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
233 Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
236 Scalar mu1 = viscosityColumn[newRowIdx];
237 Scalar mu2 = viscosityColumn[newRowIdx - 1];
238 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
242 Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
246 RwArray.push_back(newRw);
247 gasBArray.push_back(newBg);
248 gasMuArray.push_back(newMug);
251 invGasB.appendSamplePoint(xIdx, newRw, 1.0/newBg);
252 gasMu.appendSamplePoint(xIdx, newRw, newMug);
268 saturatedWaterVaporizationFactorTable_.resize(
numRegions);
280 waterReferenceDensity_[regionIdx] = rhoRefWater;
281 gasReferenceDensity_[regionIdx] = rhoRefGas;
290 { saturatedWaterVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
305 { inverseGasB_[regionIdx] = invBg; }
313 { gasMu_[regionIdx] = mug; }
324 auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
326 constexpr
const Scalar RwMin = 0.0;
327 Scalar RwMax = waterVaporizationFac.eval(saturatedWaterVaporizationFactorTable_[regionIdx].xMax(),
true);
329 Scalar poMin = samplePoints.front().first;
330 Scalar poMax = samplePoints.back().first;
332 constexpr
const size_t nRw = 20;
333 size_t nP = samplePoints.size()*2;
340 for (
size_t RwIdx = 0; RwIdx < nRw; ++RwIdx) {
341 Scalar Rw = RwMin + (RwMax - RwMin)*RwIdx/nRw;
343 gasMu_[regionIdx].appendXPos(Rw);
345 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
346 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
347 Scalar mug = mugTable.
eval(pg,
true);
349 gasMu_[regionIdx].appendSamplePoint(RwIdx, pg, mug);
361 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
364 const auto& gasMu = gasMu_[regionIdx];
365 const auto& invGasB = inverseGasB_[regionIdx];
366 assert(gasMu.numX() == invGasB.numX());
368 auto& invGasBMu = inverseGasBMu_[regionIdx];
369 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
370 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
372 std::vector<Scalar> satPressuresArray;
373 std::vector<Scalar> invSatGasBArray;
374 std::vector<Scalar> invSatGasBMuArray;
375 for (
size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
376 invGasBMu.appendXPos(gasMu.xAt(pIdx));
378 assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
380 size_t numRw = gasMu.numY(pIdx);
381 for (
size_t RwIdx = 0; RwIdx < numRw; ++RwIdx)
382 invGasBMu.appendSamplePoint(pIdx,
383 gasMu.yAt(pIdx, RwIdx),
384 invGasB.valueAt(pIdx, RwIdx)
385 / gasMu.valueAt(pIdx, RwIdx));
390 satPressuresArray.push_back(gasMu.xAt(pIdx));
391 invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRw - 1));
392 invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRw - 1));
395 invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
396 invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
398 updateSaturationPressure_(regionIdx);
406 {
return gasReferenceDensity_.size(); }
411 template <
class Evaluation>
415 const Evaluation&)
const
417 throw std::runtime_error(
"Requested the enthalpy of gas but the thermal option is not enabled");
423 template <
class Evaluation>
426 const Evaluation& pressure,
428 const Evaluation& Rvw)
const
430 const Evaluation& invBg = inverseGasB_[regionIdx].eval(pressure, Rvw,
true);
431 const Evaluation& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rvw,
true);
433 return invBg/invMugBg;
439 template <
class Evaluation>
442 const Evaluation& pressure)
const
444 const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure,
true);
445 const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure,
true);
447 return invBg/invMugBg;
453 template <
class Evaluation>
456 const Evaluation& pressure,
458 const Evaluation& Rvw)
const
459 {
return inverseGasB_[regionIdx].eval(pressure, Rvw,
true); }
464 template <
class Evaluation>
467 const Evaluation& pressure)
const
468 {
return inverseSaturatedGasB_[regionIdx].eval(pressure,
true); }
473 template <
class Evaluation>
476 const Evaluation& pressure)
const
478 return saturatedWaterVaporizationFactorTable_[regionIdx].eval(pressure,
true);
484 template <
class Evaluation>
489 const Evaluation& )
const
495 template <
class Evaluation>
498 const Evaluation& )
const
508 template <
class Evaluation>
511 const Evaluation& Rw)
const
515 const auto& RwTable = saturatedWaterVaporizationFactorTable_[regionIdx];
516 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
519 Evaluation pSat = saturationPressure_[regionIdx].eval(Rw,
true);
524 bool onProbation =
false;
525 for (
unsigned i = 0; i < 20; ++i) {
526 const Evaluation& f = RwTable.eval(pSat,
true) - Rw;
527 const Evaluation& fPrime = RwTable.evalDerivative(pSat,
true);
531 if (std::abs(scalarValue(fPrime)) < 1.0e-30) {
535 const Evaluation& delta = f/fPrime;
549 if (std::abs(scalarValue(delta)) < std::abs(scalarValue(pSat))*eps)
553 std::stringstream errlog;
554 errlog <<
"Finding saturation pressure did not converge:"
555 <<
" pSat = " << pSat
558 OpmLog::debug(
"Wet gas saturation pressure", errlog.str());
563 template <
class Evaluation>
564 Evaluation diffusionCoefficient(
const Evaluation& ,
568 throw std::runtime_error(
"Not implemented: The PVT model does not provide a diffusionCoefficient()");
571 const Scalar gasReferenceDensity(
unsigned regionIdx)
const
572 {
return gasReferenceDensity_[regionIdx]; }
574 const Scalar waterReferenceDensity(
unsigned regionIdx)
const
575 {
return waterReferenceDensity_[regionIdx]; }
577 const std::vector<TabulatedTwoDFunction>& inverseGasB()
const {
581 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB()
const {
582 return inverseSaturatedGasB_;
585 const std::vector<TabulatedTwoDFunction>& gasMu()
const {
589 const std::vector<TabulatedTwoDFunction>& inverseGasBMu()
const {
590 return inverseGasBMu_;
593 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu()
const {
594 return inverseSaturatedGasBMu_;
597 const std::vector<TabulatedOneDFunction>& saturatedWaterVaporizationFactorTable()
const {
598 return saturatedWaterVaporizationFactorTable_;
601 const std::vector<TabulatedOneDFunction>& saturationPressure()
const {
602 return saturationPressure_;
605 Scalar vapPar1()
const {
609 bool operator==(
const DryHumidGasPvt<Scalar>& data)
const
611 return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
612 this->waterReferenceDensity_ == data.waterReferenceDensity_ &&
613 this->inverseGasB() == data.inverseGasB() &&
614 this->inverseSaturatedGasB() == data.inverseSaturatedGasB() &&
615 this->gasMu() == data.gasMu() &&
616 this->inverseGasBMu() == data.inverseGasBMu() &&
617 this->inverseSaturatedGasBMu() == data.inverseSaturatedGasBMu() &&
618 this->saturatedWaterVaporizationFactorTable() == data.saturatedWaterVaporizationFactorTable() &&
619 this->saturationPressure() == data.saturationPressure() &&
620 this->vapPar1() == data.vapPar1();
624 void updateSaturationPressure_(
unsigned regionIdx)
626 typedef std::pair<Scalar, Scalar> Pair;
627 const auto& waterVaporizationFac = saturatedWaterVaporizationFactorTable_[regionIdx];
631 size_t n = waterVaporizationFac.numSamples();
632 Scalar delta = (waterVaporizationFac.xMax() - waterVaporizationFac.xMin())/Scalar(n + 1);
634 SamplingPoints pSatSamplePoints;
636 for (
size_t i = 0; i <= n; ++ i) {
637 Scalar pSat = waterVaporizationFac.xMin() + Scalar(i)*delta;
641 pSatSamplePoints.push_back(val);
645 auto x_coord_comparator = [](
const Pair& a,
const Pair& b) {
return a.first == b.first; };
646 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
647 if (std::distance(pSatSamplePoints.begin(), last) > 1)
648 pSatSamplePoints.erase(last, pSatSamplePoints.end());
650 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
653 std::vector<Scalar> gasReferenceDensity_;
654 std::vector<Scalar> waterReferenceDensity_;
655 std::vector<TabulatedTwoDFunction> inverseGasB_;
656 std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
657 std::vector<TabulatedTwoDFunction> gasMu_;
658 std::vector<TabulatedTwoDFunction> inverseGasBMu_;
659 std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
660 std::vector<TabulatedOneDFunction> saturatedWaterVaporizationFactorTable_;
661 std::vector<TabulatedOneDFunction> saturationPressure_;
A central place for various physical constants occuring in some equations.
Implements a linearly interpolated scalar function that depends on one variable.
This class represents the Pressure-Volume-Temperature relations of the gas phase with vaporized water...
Definition: DryHumidGasPvt.hpp:53
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &, const Evaluation &Rvw) const
Returns the formation volume factor [-] of the fluid phase.
Definition: DryHumidGasPvt.hpp:454
Evaluation saturatedOilVaporizationFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the oil vaporization factor [m^3/m^3] of the oil phase.
Definition: DryHumidGasPvt.hpp:485
void setGasViscosity(unsigned regionIdx, const TabulatedTwoDFunction &mug)
Initialize the viscosity of the gas phase.
Definition: DryHumidGasPvt.hpp:312
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &, const Evaluation &Rw) const
Returns the saturation pressure of the gas phase [Pa] depending on its mass fraction of the water com...
Definition: DryHumidGasPvt.hpp:509
void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBg)
Initialize the function for the gas formation volume factor.
Definition: DryHumidGasPvt.hpp:304
void initEnd()
Finish initializing the gas phase PVT properties.
Definition: DryHumidGasPvt.hpp:357
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: DryHumidGasPvt.hpp:405
void setReferenceDensities(unsigned regionIdx, Scalar, Scalar rhoRefGas, Scalar rhoRefWater)
Initialize the reference densities of all fluids for a given PVT region.
Definition: DryHumidGasPvt.hpp:275
Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: DryHumidGasPvt.hpp:474
void setSaturatedGasWaterVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil vaporization factor .
Definition: DryHumidGasPvt.hpp:289
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &, const Evaluation &Rvw) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: DryHumidGasPvt.hpp:424
void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for oil saturated gas.
Definition: DryHumidGasPvt.hpp:322
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of gas given a set of parameters.
Definition: DryHumidGasPvt.hpp:412
Evaluation saturatedOilVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the oil vaporization factor [m^3/m^3] of the oil phase.
Definition: DryHumidGasPvt.hpp:496
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of oil saturated gas at a given pressure.
Definition: DryHumidGasPvt.hpp:440
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of water saturated gas at a given pressure.
Definition: DryHumidGasPvt.hpp:465
Definition: Exceptions.hpp:46
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:259
void setContainerOfTuples(const XYContainer &points, bool sortInputs=true)
Set the sampling points of the piecewise linear function using a STL-compatible container of tuple-li...
Definition: Tabulated1DFunction.hpp:185