27 #ifndef OPM_WET_GAS_PVT_HPP
28 #define OPM_WET_GAS_PVT_HPP
37 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
38 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
39 #include <opm/parser/eclipse/EclipseState/Schedule/OilVaporizationProperties.hpp>
43 #include <opm/common/OpmLog/OpmLog.hpp>
52 template <
class Scalar>
55 typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
66 WetGasPvt(
const std::vector<Scalar>& gasReferenceDensity,
67 const std::vector<Scalar>& oilReferenceDensity,
68 const std::vector<TabulatedTwoDFunction>& inverseGasB,
69 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB,
70 const std::vector<TabulatedTwoDFunction>& gasMu,
71 const std::vector<TabulatedTwoDFunction>& inverseGasBMu,
72 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu,
73 const std::vector<TabulatedOneDFunction>& saturatedOilVaporizationFactorTable,
76 : gasReferenceDensity_(gasReferenceDensity)
77 , oilReferenceDensity_(oilReferenceDensity)
78 , inverseGasB_(inverseGasB)
79 , inverseSaturatedGasB_(inverseSaturatedGasB)
81 , inverseGasBMu_(inverseGasBMu)
82 , inverseSaturatedGasBMu_(inverseSaturatedGasBMu)
83 , saturatedOilVaporizationFactorTable_(saturatedOilVaporizationFactorTable)
84 , saturationPressure_(saturationPressure)
96 void initFromState(
const EclipseState& eclState,
const Schedule& schedule)
98 const auto& pvtgTables = eclState.getTableManager().getPvtgTables();
99 const auto& densityTable = eclState.getTableManager().getDensityTable();
101 assert(pvtgTables.size() == densityTable.size());
106 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
107 Scalar rhoRefO = densityTable[regionIdx].oil;
108 Scalar rhoRefG = densityTable[regionIdx].gas;
109 Scalar rhoRefW = densityTable[regionIdx].water;
114 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
115 const auto& pvtgTable = pvtgTables[regionIdx];
117 const auto& saturatedTable = pvtgTable.getSaturatedTable();
118 assert(saturatedTable.numRows() > 1);
120 auto& gasMu = gasMu_[regionIdx];
121 auto& invGasB = inverseGasB_[regionIdx];
122 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
123 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
124 auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
126 oilVaporizationFac.setXYArrays(saturatedTable.numRows(),
127 saturatedTable.getColumn(
"PG"),
128 saturatedTable.getColumn(
"RV"));
130 std::vector<Scalar> invSatGasBArray;
131 std::vector<Scalar> invSatGasBMuArray;
134 for (
unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
135 Scalar pg = saturatedTable.get(
"PG" , outerIdx);
136 Scalar B = saturatedTable.get(
"BG" , outerIdx);
137 Scalar mu = saturatedTable.get(
"MUG" , outerIdx);
139 invGasB.appendXPos(pg);
140 gasMu.appendXPos(pg);
142 invSatGasBArray.push_back(1.0/B);
143 invSatGasBMuArray.push_back(1.0/(mu*B));
145 assert(invGasB.numX() == outerIdx + 1);
146 assert(gasMu.numX() == outerIdx + 1);
148 const auto& underSaturatedTable = pvtgTable.getUnderSaturatedTable(outerIdx);
149 size_t numRows = underSaturatedTable.numRows();
150 for (
size_t innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
151 Scalar Rv = underSaturatedTable.get(
"RV" , innerIdx);
152 Scalar Bg = underSaturatedTable.get(
"BG" , innerIdx);
153 Scalar mug = underSaturatedTable.get(
"MUG" , innerIdx);
155 invGasB.appendSamplePoint(outerIdx, Rv, 1.0/Bg);
156 gasMu.appendSamplePoint(outerIdx, Rv, mug);
161 std::vector<double> tmpPressure = saturatedTable.getColumn(
"PG").vectorCopy( );
163 invSatGasB.setXYContainers(tmpPressure, invSatGasBArray);
164 invSatGasBMu.setXYContainers(tmpPressure, invSatGasBMuArray);
168 for (
unsigned xIdx = 0; xIdx < invGasB.numX(); ++xIdx) {
170 assert(invGasB.numY(xIdx) > 0);
174 if (invGasB.numY(xIdx) > 1)
180 size_t masterTableIdx = xIdx + 1;
181 for (; masterTableIdx < saturatedTable.numRows(); ++masterTableIdx)
183 if (pvtgTable.getUnderSaturatedTable(masterTableIdx).numRows() > 1)
187 if (masterTableIdx >= saturatedTable.numRows())
188 throw std::runtime_error(
"PVTG tables are invalid: The last table must exhibit at least one "
189 "entry for undersaturated gas!");
193 extendPvtgTable_(regionIdx,
195 pvtgTable.getUnderSaturatedTable(xIdx),
196 pvtgTable.getUnderSaturatedTable(masterTableIdx));
201 const auto& oilVap = schedule[0].oilvap();
202 if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) {
203 vapPar1_ = oilVap.vap1();
210 void extendPvtgTable_(
unsigned regionIdx,
212 const SimpleTable& curTable,
213 const SimpleTable& masterTable)
215 std::vector<double> RvArray = curTable.getColumn(
"RV").vectorCopy();
216 std::vector<double> gasBArray = curTable.getColumn(
"BG").vectorCopy();
217 std::vector<double> gasMuArray = curTable.getColumn(
"MUG").vectorCopy();
219 auto& invGasB = inverseGasB_[regionIdx];
220 auto& gasMu = gasMu_[regionIdx];
222 for (
size_t newRowIdx = 1; newRowIdx < masterTable.numRows(); ++ newRowIdx) {
223 const auto& RVColumn = masterTable.getColumn(
"RV");
224 const auto& BGColumn = masterTable.getColumn(
"BG");
225 const auto& viscosityColumn = masterTable.getColumn(
"MUG");
228 Scalar diffRv = RVColumn[newRowIdx] - RVColumn[newRowIdx - 1];
229 Scalar newRv = RvArray.back() + diffRv;
232 Scalar B1 = BGColumn[newRowIdx];
233 Scalar B2 = BGColumn[newRowIdx - 1];
234 Scalar x = (B1 - B2)/( (B1 + B2)/2.0 );
238 Scalar newBg = gasBArray.back()*(1.0 + x/2.0)/(1.0 - x/2.0);
241 Scalar mu1 = viscosityColumn[newRowIdx];
242 Scalar mu2 = viscosityColumn[newRowIdx - 1];
243 Scalar xMu = (mu1 - mu2)/( (mu1 + mu2)/2.0 );
247 Scalar newMug = gasMuArray.back()*(1.0 + xMu/2)/(1.0 - xMu/2.0);
251 RvArray.push_back(newRv);
252 gasBArray.push_back(newBg);
253 gasMuArray.push_back(newMug);
256 invGasB.appendSamplePoint(xIdx, newRv, 1.0/newBg);
257 gasMu.appendSamplePoint(xIdx, newRv, newMug);
273 saturatedOilVaporizationFactorTable_.resize(
numRegions);
285 oilReferenceDensity_[regionIdx] = rhoRefOil;
286 gasReferenceDensity_[regionIdx] = rhoRefGas;
295 { saturatedOilVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); }
308 auto& invGasB = inverseGasB_[regionIdx];
310 const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
312 Scalar T = 273.15 + 15.56;
315 Scalar RvMax = RvTable.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(),
true);
317 Scalar poMin = samplePoints.front().first;
318 Scalar poMax = samplePoints.back().first;
321 size_t nP = samplePoints.size()*2;
323 Scalar rhooRef = oilReferenceDensity_[regionIdx];
328 updateSaturationPressure_(regionIdx);
334 for (
size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
335 Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
337 invGasB.appendXPos(Rv);
339 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
340 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
342 Scalar poSat = saturationPressure(regionIdx, T, Rv);
343 Scalar BgSat = gasFormationVolumeFactor.
eval(poSat,
true);
344 Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76);
345 Scalar rhoo = rhooRef/BgSat*(1 + drhoo_dp*(pg - poSat));
347 Scalar Bg = rhooRef/rhoo;
349 invGasB.appendSamplePoint(RvIdx, pg, 1.0/Bg);
367 { inverseGasB_[regionIdx] = invBg; }
375 { gasMu_[regionIdx] = mug; }
386 auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
389 Scalar RvMax = oilVaporizationFac.eval(saturatedOilVaporizationFactorTable_[regionIdx].xMax(),
true);
391 Scalar poMin = samplePoints.front().first;
392 Scalar poMax = samplePoints.back().first;
395 size_t nP = samplePoints.size()*2;
402 for (
size_t RvIdx = 0; RvIdx < nRv; ++RvIdx) {
403 Scalar Rv = RvMin + (RvMax - RvMin)*RvIdx/nRv;
405 gasMu_[regionIdx].appendXPos(Rv);
407 for (
size_t pIdx = 0; pIdx < nP; ++pIdx) {
408 Scalar pg = poMin + (poMax - poMin)*pIdx/nP;
409 Scalar mug = mugTable.
eval(pg,
true);
411 gasMu_[regionIdx].appendSamplePoint(RvIdx, pg, mug);
423 for (
unsigned regionIdx = 0; regionIdx <
numRegions; ++ regionIdx) {
426 const auto& gasMu = gasMu_[regionIdx];
427 const auto& invGasB = inverseGasB_[regionIdx];
428 assert(gasMu.numX() == invGasB.numX());
430 auto& invGasBMu = inverseGasBMu_[regionIdx];
431 auto& invSatGasB = inverseSaturatedGasB_[regionIdx];
432 auto& invSatGasBMu = inverseSaturatedGasBMu_[regionIdx];
434 std::vector<Scalar> satPressuresArray;
435 std::vector<Scalar> invSatGasBArray;
436 std::vector<Scalar> invSatGasBMuArray;
437 for (
size_t pIdx = 0; pIdx < gasMu.numX(); ++pIdx) {
438 invGasBMu.appendXPos(gasMu.xAt(pIdx));
440 assert(gasMu.numY(pIdx) == invGasB.numY(pIdx));
442 size_t numRv = gasMu.numY(pIdx);
443 for (
size_t rvIdx = 0; rvIdx < numRv; ++rvIdx)
444 invGasBMu.appendSamplePoint(pIdx,
445 gasMu.yAt(pIdx, rvIdx),
446 invGasB.valueAt(pIdx, rvIdx)
447 / gasMu.valueAt(pIdx, rvIdx));
452 satPressuresArray.push_back(gasMu.xAt(pIdx));
453 invSatGasBArray.push_back(invGasB.valueAt(pIdx, numRv - 1));
454 invSatGasBMuArray.push_back(invGasBMu.valueAt(pIdx, numRv - 1));
457 invSatGasB.setXYContainers(satPressuresArray, invSatGasBArray);
458 invSatGasBMu.setXYContainers(satPressuresArray, invSatGasBMuArray);
460 updateSaturationPressure_(regionIdx);
468 {
return gasReferenceDensity_.size(); }
473 template <
class Evaluation>
477 const Evaluation&)
const
479 throw std::runtime_error(
"Requested the enthalpy of gas but the thermal option is not enabled");
485 template <
class Evaluation>
488 const Evaluation& pressure,
489 const Evaluation& Rv)
const
491 const Evaluation& invBg = inverseGasB_[regionIdx].eval(pressure, Rv,
true);
492 const Evaluation& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rv,
true);
494 return invBg/invMugBg;
500 template <
class Evaluation>
503 const Evaluation& pressure)
const
505 const Evaluation& invBg = inverseSaturatedGasB_[regionIdx].eval(pressure,
true);
506 const Evaluation& invMugBg = inverseSaturatedGasBMu_[regionIdx].eval(pressure,
true);
508 return invBg/invMugBg;
514 template <
class Evaluation>
517 const Evaluation& pressure,
518 const Evaluation& Rv)
const
519 {
return inverseGasB_[regionIdx].eval(pressure, Rv,
true); }
524 template <
class Evaluation>
527 const Evaluation& pressure)
const
528 {
return inverseSaturatedGasB_[regionIdx].eval(pressure,
true); }
533 template <
class Evaluation>
536 const Evaluation& pressure)
const
538 return saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure,
true);
548 template <
class Evaluation>
551 const Evaluation& pressure,
552 const Evaluation& oilSaturation,
553 Evaluation maxOilSaturation)
const
556 saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure,
true);
560 maxOilSaturation = min(maxOilSaturation, Scalar(1.0));
561 if (vapPar1_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
562 static const Scalar eps = 0.001;
563 const Evaluation& So = max(oilSaturation, eps);
564 tmp *= max(1e-3, pow(So/maxOilSaturation, vapPar1_));
580 template <
class Evaluation>
583 const Evaluation& Rv)
const
587 const auto& RvTable = saturatedOilVaporizationFactorTable_[regionIdx];
588 const Scalar eps = std::numeric_limits<typename Toolbox::Scalar>::epsilon()*1e6;
591 Evaluation pSat = saturationPressure_[regionIdx].eval(Rv,
true);
596 bool onProbation =
false;
597 for (
unsigned i = 0; i < 20; ++i) {
598 const Evaluation& f = RvTable.eval(pSat,
true) - Rv;
599 const Evaluation& fPrime = RvTable.evalDerivative(pSat,
true);
603 if (std::abs(scalarValue(fPrime)) < 1.0e-30) {
607 const Evaluation& delta = f/fPrime;
621 if (std::abs(scalarValue(delta)) < std::abs(scalarValue(pSat))*eps)
625 std::stringstream errlog;
626 errlog <<
"Finding saturation pressure did not converge:"
627 <<
" pSat = " << pSat
630 OpmLog::debug(
"Wet gas saturation pressure", errlog.str());
635 template <
class Evaluation>
636 Evaluation diffusionCoefficient(
const Evaluation& ,
640 throw std::runtime_error(
"Not implemented: The PVT model does not provide a diffusionCoefficient()");
643 const Scalar gasReferenceDensity(
unsigned regionIdx)
const
644 {
return gasReferenceDensity_[regionIdx]; }
646 const Scalar oilReferenceDensity(
unsigned regionIdx)
const
647 {
return oilReferenceDensity_[regionIdx]; }
649 const std::vector<TabulatedTwoDFunction>& inverseGasB()
const {
653 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasB()
const {
654 return inverseSaturatedGasB_;
657 const std::vector<TabulatedTwoDFunction>& gasMu()
const {
661 const std::vector<TabulatedTwoDFunction>& inverseGasBMu()
const {
662 return inverseGasBMu_;
665 const std::vector<TabulatedOneDFunction>& inverseSaturatedGasBMu()
const {
666 return inverseSaturatedGasBMu_;
669 const std::vector<TabulatedOneDFunction>& saturatedOilVaporizationFactorTable()
const {
670 return saturatedOilVaporizationFactorTable_;
673 const std::vector<TabulatedOneDFunction>& saturationPressure()
const {
674 return saturationPressure_;
677 Scalar vapPar1()
const {
681 bool operator==(
const WetGasPvt<Scalar>& data)
const
683 return this->gasReferenceDensity_ == data.gasReferenceDensity_ &&
684 this->oilReferenceDensity_ == data.oilReferenceDensity_ &&
685 this->inverseGasB() == data.inverseGasB() &&
686 this->inverseSaturatedGasB() == data.inverseSaturatedGasB() &&
687 this->gasMu() == data.gasMu() &&
688 this->inverseGasBMu() == data.inverseGasBMu() &&
689 this->inverseSaturatedGasBMu() == data.inverseSaturatedGasBMu() &&
690 this->saturatedOilVaporizationFactorTable() == data.saturatedOilVaporizationFactorTable() &&
691 this->saturationPressure() == data.saturationPressure() &&
692 this->vapPar1() == data.vapPar1();
696 void updateSaturationPressure_(
unsigned regionIdx)
698 typedef std::pair<Scalar, Scalar> Pair;
699 const auto& oilVaporizationFac = saturatedOilVaporizationFactorTable_[regionIdx];
703 size_t n = oilVaporizationFac.numSamples();
704 Scalar delta = (oilVaporizationFac.xMax() - oilVaporizationFac.xMin())/Scalar(n + 1);
706 SamplingPoints pSatSamplePoints;
708 for (
size_t i = 0; i <= n; ++ i) {
709 Scalar pSat = oilVaporizationFac.xMin() + Scalar(i)*delta;
713 pSatSamplePoints.push_back(val);
717 auto x_coord_comparator = [](
const Pair& a,
const Pair& b) {
return a.first == b.first; };
718 auto last = std::unique(pSatSamplePoints.begin(), pSatSamplePoints.end(), x_coord_comparator);
719 pSatSamplePoints.erase(last, pSatSamplePoints.end());
721 saturationPressure_[regionIdx].setContainerOfTuples(pSatSamplePoints);
724 std::vector<Scalar> gasReferenceDensity_;
725 std::vector<Scalar> oilReferenceDensity_;
726 std::vector<TabulatedTwoDFunction> inverseGasB_;
727 std::vector<TabulatedOneDFunction> inverseSaturatedGasB_;
728 std::vector<TabulatedTwoDFunction> gasMu_;
729 std::vector<TabulatedTwoDFunction> inverseGasBMu_;
730 std::vector<TabulatedOneDFunction> inverseSaturatedGasBMu_;
731 std::vector<TabulatedOneDFunction> saturatedOilVaporizationFactorTable_;
732 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.
Definition: Exceptions.hpp:46
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:258
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:184
This class represents the Pressure-Volume-Temperature relations of the gas phas with vaporized oil.
Definition: WetGasPvt.hpp:54
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rv) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WetGasPvt.hpp:515
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil saturated gas at a given pressure.
Definition: WetGasPvt.hpp:525
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition: WetGasPvt.hpp:534
void initEnd()
Finish initializing the gas phase PVT properties.
Definition: WetGasPvt.hpp:419
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &, const Evaluation &Rv) const
Returns the saturation pressure of the gas phase [Pa] depending on its mass fraction of the oil compo...
Definition: WetGasPvt.hpp:581
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &oilSaturation, Evaluation maxOilSaturation) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition: WetGasPvt.hpp:549
Evaluation viscosity(unsigned regionIdx, const Evaluation &, const Evaluation &pressure, const Evaluation &Rv) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WetGasPvt.hpp:486
void setGasViscosity(unsigned regionIdx, const TabulatedTwoDFunction &mug)
Initialize the viscosity of the gas phase.
Definition: WetGasPvt.hpp:374
Evaluation internalEnergy(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the specific enthalpy [J/kg] of gas given a set of parameters.
Definition: WetGasPvt.hpp:474
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: WetGasPvt.hpp:467
void setSaturatedGasViscosity(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the phase viscosity for oil saturated gas.
Definition: WetGasPvt.hpp:384
void setSaturatedGasFormationVolumeFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the gas formation volume factor.
Definition: WetGasPvt.hpp:306
void setSaturatedGasOilVaporizationFactor(unsigned regionIdx, const SamplingPoints &samplePoints)
Initialize the function for the oil vaporization factor .
Definition: WetGasPvt.hpp:294
void setInverseGasFormationVolumeFactor(unsigned regionIdx, const TabulatedTwoDFunction &invBg)
Initialize the function for the gas formation volume factor.
Definition: WetGasPvt.hpp:366
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: WetGasPvt.hpp:501
void setReferenceDensities(unsigned regionIdx, Scalar rhoRefOil, Scalar rhoRefGas, Scalar)
Initialize the reference densities of all fluids for a given PVT region.
Definition: WetGasPvt.hpp:280