27 #ifndef OPM_GAS_PVT_THERMAL_HPP
28 #define OPM_GAS_PVT_THERMAL_HPP
33 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
34 #include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
35 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
39 template <
class Scalar,
bool enableThermal>
40 class GasPvtMultiplexer;
48 template <
class Scalar>
57 enableThermalDensity_ =
false;
58 enableJouleThomson_ =
false;
59 enableThermalViscosity_ =
false;
60 isothermalPvt_ =
nullptr;
64 const std::vector<TabulatedOneDFunction>& gasvisctCurves,
65 const std::vector<Scalar>& gasdentRefTemp,
66 const std::vector<Scalar>& gasdentCT1,
67 const std::vector<Scalar>& gasdentCT2,
68 const std::vector<Scalar>& gasJTRefPres,
69 const std::vector<Scalar>& gasJTC,
70 const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
72 bool enableJouleThomson,
74 bool enableInternalEnergy)
75 : isothermalPvt_(isothermalPvt)
76 , gasvisctCurves_(gasvisctCurves)
77 , gasdentRefTemp_(gasdentRefTemp)
78 , gasdentCT1_(gasdentCT1)
79 , gasdentCT2_(gasdentCT2)
80 , gasJTRefPres_(gasJTRefPres)
82 , internalEnergyCurves_(internalEnergyCurves)
84 , enableJouleThomson_(enableJouleThomson)
86 , enableInternalEnergy_(enableInternalEnergy)
93 {
delete isothermalPvt_; }
99 void initFromState(
const EclipseState& eclState,
const Schedule& schedule)
105 isothermalPvt_->initFromState(eclState, schedule);
110 const auto& tables = eclState.getTableManager();
112 enableThermalDensity_ = tables.GasDenT().size() > 0;
113 enableJouleThomson_ = tables.GasJT().size() > 0;
114 enableThermalViscosity_ = tables.hasTables(
"GASVISCT");
115 enableInternalEnergy_ = tables.hasTables(
"SPECHEAT");
117 unsigned numRegions = isothermalPvt_->
numRegions();
121 if (enableThermalViscosity_) {
122 const auto& gasvisctTables = tables.getGasvisctTables();
123 auto gasCompIdx = tables.gas_comp_index();
124 std::string gasvisctColumnName =
"Viscosity" + std::to_string(
static_cast<long long>(gasCompIdx));
126 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
127 const auto& T = gasvisctTables[regionIdx].getColumn(
"Temperature").vectorCopy();
128 const auto& mu = gasvisctTables[regionIdx].getColumn(gasvisctColumnName).vectorCopy();
129 gasvisctCurves_[regionIdx].setXYContainers(T, mu);
134 if (enableThermalDensity_) {
135 const auto& gasDenT = tables.GasDenT();
137 assert(gasDenT.size() == numRegions);
138 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
139 const auto& record = gasDenT[regionIdx];
141 gasdentRefTemp_[regionIdx] = record.T0;
142 gasdentCT1_[regionIdx] = record.C1;
143 gasdentCT2_[regionIdx] = record.C2;
148 if (enableJouleThomson_) {
149 const auto& gasJT = tables.GasJT();
151 assert(gasJT.size() == numRegions);
152 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
153 const auto& record = gasJT[regionIdx];
155 gasJTRefPres_[regionIdx] = record.P0;
156 gasJTC_[regionIdx] = record.C1;
159 const auto& densityTable = eclState.getTableManager().getDensityTable();
161 assert(densityTable.size() == numRegions);
162 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
163 rhoRefO_[regionIdx] = densityTable[regionIdx].oil;
167 if (enableInternalEnergy_) {
171 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
172 const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
173 const auto& temperatureColumn = specHeatTable.getColumn(
"TEMPERATURE");
174 const auto& cvGasColumn = specHeatTable.getColumn(
"CV_GAS");
176 std::vector<double> uSamples(temperatureColumn.size());
183 const Scalar hVap = 480.6e3;
185 Scalar u = temperatureColumn[0]*cvGasColumn[0] + hVap;
186 for (
size_t i = 0;; ++i) {
189 if (i >= temperatureColumn.size() - 1)
194 Scalar c_v0 = cvGasColumn[i];
195 Scalar c_v1 = cvGasColumn[i + 1];
196 Scalar T0 = temperatureColumn[i];
197 Scalar T1 = temperatureColumn[i + 1];
198 u += 0.5*(c_v0 + c_v1)*(T1 - T0);
201 internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
212 gasvisctCurves_.resize(numRegions);
213 internalEnergyCurves_.resize(numRegions);
214 gasdentRefTemp_.resize(numRegions);
215 gasdentCT1_.resize(numRegions);
216 gasdentCT2_.resize(numRegions);
217 gasJTRefPres_.resize(numRegions);
218 gasJTC_.resize(numRegions);
219 rhoRefO_.resize(numRegions);
228 size_t numRegions()
const
229 {
return gasvisctCurves_.size(); }
235 {
return enableThermalDensity_; }
241 {
return enableJouleThomson_; }
247 {
return enableThermalViscosity_; }
252 template <
class Evaluation>
254 const Evaluation& temperature,
255 const Evaluation& pressure,
256 const Evaluation& Rv)
const
258 if (!enableInternalEnergy_)
259 throw std::runtime_error(
"Requested the internal energy of gas but it is disabled");
261 if (!enableJouleThomson_) {
265 return internalEnergyCurves_[regionIdx].eval(temperature,
true);
268 Evaluation Tref = gasdentRefTemp_[regionIdx];
269 Evaluation Pref = gasJTRefPres_[regionIdx];
270 Scalar JTC = gasJTC_[regionIdx];
271 Evaluation Rvw = 0.0;
274 constexpr
const Scalar hVap = 480.6e3;
275 Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature,
true) - hVap)/temperature;
276 Evaluation density = invB * (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
278 Evaluation enthalpyPres;
280 enthalpyPres = -Cp * JTC * (pressure -Pref);
282 else if(enableThermalDensity_) {
283 Scalar c1T = gasdentCT1_[regionIdx];
284 Scalar c2T = gasdentCT2_[regionIdx];
286 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
287 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
289 constexpr
const int N = 100;
290 Evaluation deltaP = (pressure - Pref)/N;
291 Evaluation enthalpyPresPrev = 0;
292 for (
size_t i = 0; i < N; ++i) {
293 Evaluation Pnew = Pref + i * deltaP;
295 (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
297 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
298 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
299 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
300 enthalpyPresPrev = enthalpyPres;
304 throw std::runtime_error(
"Requested Joule-thomson calculation but thermal gas density (GASDENT) is not provided");
307 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
309 return enthalpy - pressure/density;
316 template <
class Evaluation>
318 const Evaluation& temperature,
319 const Evaluation& pressure,
320 const Evaluation& Rv,
321 const Evaluation& Rvw)
const
324 return isothermalPvt_->
viscosity(regionIdx, temperature, pressure, Rv, Rvw);
327 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature);
334 template <
class Evaluation>
336 const Evaluation& temperature,
337 const Evaluation& pressure)
const
343 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature,
true);
350 template <
class Evaluation>
352 const Evaluation& temperature,
353 const Evaluation& pressure,
354 const Evaluation& Rv,
355 const Evaluation& )
const
357 const Evaluation& Rvw = 0.0;
371 Scalar TRef = gasdentRefTemp_[regionIdx];
372 Scalar cT1 = gasdentCT1_[regionIdx];
373 Scalar cT2 = gasdentCT2_[regionIdx];
374 const Evaluation& Y = temperature - TRef;
376 return b/(1 + (cT1 + cT2*Y)*Y);
382 template <
class Evaluation>
384 const Evaluation& temperature,
385 const Evaluation& pressure)
const
400 Scalar TRef = gasdentRefTemp_[regionIdx];
401 Scalar cT1 = gasdentCT1_[regionIdx];
402 Scalar cT2 = gasdentCT2_[regionIdx];
403 const Evaluation& Y = temperature - TRef;
405 return b/(1 + (cT1 + cT2*Y)*Y);
411 template <
class Evaluation>
414 const Evaluation& )
const
425 template <
class Evaluation>
427 const Evaluation& temperature,
428 const Evaluation& pressure)
const
438 template <
class Evaluation>
440 const Evaluation& temperature,
441 const Evaluation& pressure,
442 const Evaluation& oilSaturation,
443 const Evaluation& maxOilSaturation)
const
453 template <
class Evaluation>
455 const Evaluation& temperature,
456 const Evaluation& pressure)
const
459 template <
class Evaluation>
460 Evaluation diffusionCoefficient(
const Evaluation& temperature,
461 const Evaluation& pressure,
462 unsigned compIdx)
const
466 const IsothermalPvt* isoThermalPvt()
const
467 {
return isothermalPvt_; }
469 const Scalar gasReferenceDensity(
unsigned regionIdx)
const
472 const std::vector<TabulatedOneDFunction>& gasvisctCurves()
const
473 {
return gasvisctCurves_; }
475 const std::vector<Scalar>& gasdentRefTemp()
const
476 {
return gasdentRefTemp_; }
478 const std::vector<Scalar>& gasdentCT1()
const
479 {
return gasdentCT1_; }
481 const std::vector<Scalar>& gasdentCT2()
const
482 {
return gasdentCT2_; }
484 const std::vector<TabulatedOneDFunction>& internalEnergyCurves()
const
485 {
return internalEnergyCurves_; }
487 bool enableInternalEnergy()
const
488 {
return enableInternalEnergy_; }
490 const std::vector<Scalar>& gasJTRefPres()
const
491 {
return gasJTRefPres_; }
493 const std::vector<Scalar>& gasJTC()
const
497 bool operator==(
const GasPvtThermal<Scalar>& data)
const
499 if (isothermalPvt_ && !data.isothermalPvt_)
501 if (!isothermalPvt_ && data.isothermalPvt_)
504 return (!this->isoThermalPvt() ||
505 (*this->isoThermalPvt() == *data.isoThermalPvt())) &&
506 this->gasvisctCurves() == data.gasvisctCurves() &&
507 this->gasdentRefTemp() == data.gasdentRefTemp() &&
508 this->gasdentCT1() == data.gasdentCT1() &&
509 this->gasdentCT2() == data.gasdentCT2() &&
510 this->gasJTRefPres() == data.gasJTRefPres() &&
511 this->gasJTC() == data.gasJTC() &&
512 this->internalEnergyCurves() == data.internalEnergyCurves() &&
514 this->enableJouleThomson() == data.enableJouleThomson() &&
516 this->enableInternalEnergy() == data.enableInternalEnergy();
519 GasPvtThermal<Scalar>& operator=(
const GasPvtThermal<Scalar>& data)
521 if (data.isothermalPvt_)
522 isothermalPvt_ =
new IsothermalPvt(*data.isothermalPvt_);
524 isothermalPvt_ =
nullptr;
525 gasvisctCurves_ = data.gasvisctCurves_;
526 gasdentRefTemp_ = data.gasdentRefTemp_;
527 gasdentCT1_ = data.gasdentCT1_;
528 gasdentCT2_ = data.gasdentCT2_;
529 gasJTRefPres_ = data.gasJTRefPres_;
530 gasJTC_ = data.gasJTC_;
531 internalEnergyCurves_ = data.internalEnergyCurves_;
532 enableThermalDensity_ = data.enableThermalDensity_;
533 enableJouleThomson_ = data.enableJouleThomson_;
534 enableThermalViscosity_ = data.enableThermalViscosity_;
535 enableInternalEnergy_ = data.enableInternalEnergy_;
541 IsothermalPvt* isothermalPvt_;
545 std::vector<TabulatedOneDFunction> gasvisctCurves_;
547 std::vector<Scalar> gasdentRefTemp_;
548 std::vector<Scalar> gasdentCT1_;
549 std::vector<Scalar> gasdentCT2_;
551 std::vector<Scalar> gasJTRefPres_;
552 std::vector<Scalar> gasJTC_;
554 std::vector<Scalar> rhoRefO_;
557 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
559 bool enableThermalDensity_;
560 bool enableJouleThomson_;
561 bool enableThermalViscosity_;
562 bool enableInternalEnergy_;
Implements a linearly interpolated scalar function that depends on one variable.
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition: GasPvtMultiplexer.hpp:100
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil saturated gas given a set of parameters.
Definition: GasPvtMultiplexer.hpp:272
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: GasPvtMultiplexer.hpp:218
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &Rvw) const
Returns the formation volume factor [-] of the fluid phase.
Definition: GasPvtMultiplexer.hpp:261
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &Rvw) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: GasPvtMultiplexer.hpp:241
Evaluation diffusionCoefficient(const Evaluation &temperature, const Evaluation &pressure, unsigned compIdx) const
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition: GasPvtMultiplexer.hpp:322
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the oil vaporization factor [m^3/m^3] of oil saturated gas.
Definition: GasPvtMultiplexer.hpp:281
const Scalar gasReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition: GasPvtMultiplexer.hpp:224
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &Rv) const
Returns the saturation pressure of the gas phase [Pa] depending on its mass fraction of the oil compo...
Definition: GasPvtMultiplexer.hpp:313
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of oil saturated gas given a set of parameters.
Definition: GasPvtMultiplexer.hpp:252
This class implements temperature dependence of the PVT properties of gas.
Definition: GasPvtThermal.hpp:50
bool enableThermalDensity() const
Returns true iff the density of the gas phase is temperature dependent.
Definition: GasPvtThermal.hpp:234
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &Rvw) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: GasPvtThermal.hpp:317
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &oilSaturation, const Evaluation &maxOilSaturation) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition: GasPvtThermal.hpp:439
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil-saturated gas.
Definition: GasPvtThermal.hpp:383
bool enableJouleThomsony() const
Returns true iff Joule-Thomson effect for the gas phase is active.
Definition: GasPvtThermal.hpp:240
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition: GasPvtThermal.hpp:426
void initEnd()
Finish initializing the thermal part of the gas phase PVT properties.
Definition: GasPvtThermal.hpp:225
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: GasPvtThermal.hpp:210
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &) const
Returns the formation volume factor [-] of the fluid phase.
Definition: GasPvtThermal.hpp:351
bool enableThermalViscosity() const
Returns true iff the viscosity of the gas phase is temperature dependent.
Definition: GasPvtThermal.hpp:246
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: GasPvtThermal.hpp:412
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv) const
Returns the specific internal energy [J/kg] of gas given a set of parameters.
Definition: GasPvtThermal.hpp:253
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of the oil-saturated gas phase given a set of parameters.
Definition: GasPvtThermal.hpp:335
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the saturation pressure of the gas phase [Pa].
Definition: GasPvtThermal.hpp:454
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48