My Project
Loading...
Searching...
No Matches
GasPvtThermal.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef OPM_GAS_PVT_THERMAL_HPP
28#define OPM_GAS_PVT_THERMAL_HPP
29
31
32namespace Opm {
33
34#if HAVE_ECL_INPUT
35class EclipseState;
36class Schedule;
37#endif
38
39template <class Scalar, bool enableThermal>
40class GasPvtMultiplexer;
41
48template <class Scalar>
50{
51public:
52 using IsothermalPvt = GasPvtMultiplexer<Scalar, /*enableThermal=*/false>;
54
56 {
57 enableThermalDensity_ = false;
58 enableJouleThomson_ = false;
59 enableThermalViscosity_ = false;
60 enableInternalEnergy_ = false;
61 isothermalPvt_ = nullptr;
62 }
63
64 GasPvtThermal(IsothermalPvt* isothermalPvt,
65 const std::vector<TabulatedOneDFunction>& gasvisctCurves,
66 const std::vector<Scalar>& viscrefPress,
67 const std::vector<Scalar>& viscRef,
68 const std::vector<Scalar>& gasdentRefTemp,
69 const std::vector<Scalar>& gasdentCT1,
70 const std::vector<Scalar>& gasdentCT2,
71 const std::vector<Scalar>& gasJTRefPres,
72 const std::vector<Scalar>& gasJTC,
73 const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
77 bool enableInternalEnergy)
78 : isothermalPvt_(isothermalPvt)
79 , gasvisctCurves_(gasvisctCurves)
80 , viscrefPress_(viscrefPress)
81 , viscRef_(viscRef)
82 , gasdentRefTemp_(gasdentRefTemp)
83 , gasdentCT1_(gasdentCT1)
84 , gasdentCT2_(gasdentCT2)
85 , gasJTRefPres_(gasJTRefPres)
86 , gasJTC_(gasJTC)
87 , internalEnergyCurves_(internalEnergyCurves)
88 , enableThermalDensity_(enableThermalDensity)
89 , enableJouleThomson_(enableJouleThomson)
90 , enableThermalViscosity_(enableThermalViscosity)
91 , enableInternalEnergy_(enableInternalEnergy)
92 { }
93
94 GasPvtThermal(const GasPvtThermal& data)
95 { *this = data; }
96
98 { delete isothermalPvt_; }
99
100#if HAVE_ECL_INPUT
104 void initFromState(const EclipseState& eclState, const Schedule& schedule);
105#endif // HAVE_ECL_INPUT
106
110 void setNumRegions(size_t numRegions)
111 {
112 gasvisctCurves_.resize(numRegions);
113 viscrefPress_.resize(numRegions);
114 viscRef_.resize(numRegions);
115 internalEnergyCurves_.resize(numRegions);
116 gasdentRefTemp_.resize(numRegions);
117 gasdentCT1_.resize(numRegions);
118 gasdentCT2_.resize(numRegions);
119 gasJTRefPres_.resize(numRegions);
120 gasJTC_.resize(numRegions);
121 rhoRefO_.resize(numRegions);
122 }
123
124 void setVapPars(const Scalar par1, const Scalar par2)
125 {
126 isothermalPvt_->setVapPars(par1, par2);
127 }
128
132 void initEnd()
133 { }
134
139 { return enableThermalDensity_; }
140
145 { return enableJouleThomson_; }
146
151 { return enableThermalViscosity_; }
152
153 size_t numRegions() const
154 { return viscrefPress_.size(); }
155
156
160 template <class Evaluation>
161 Evaluation internalEnergy(unsigned regionIdx,
162 const Evaluation& temperature,
163 const Evaluation& pressure,
164 const Evaluation& Rv,
165 const Evaluation& ) const
166 {
167 if (!enableInternalEnergy_)
168 throw std::runtime_error("Requested the internal energy of gas but it is disabled");
169
170 if (!enableJouleThomson_) {
171 // compute the specific internal energy for the specified tempature. We use linear
172 // interpolation here despite the fact that the underlying heat capacities are
173 // piecewise linear (which leads to a quadratic function)
174 return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
175 }
176 else {
177 Evaluation Tref = gasdentRefTemp_[regionIdx];
178 Evaluation Pref = gasJTRefPres_[regionIdx];
179 Scalar JTC = gasJTC_[regionIdx]; // if JTC is default then JTC is calculated
180 Evaluation Rvw = 0.0;
181
182 Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
183 constexpr const Scalar hVap = 480.6e3; // [J / kg]
184 Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true) - hVap)/temperature;
185 Evaluation density = invB * (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
186
187 Evaluation enthalpyPres;
188 if (JTC != 0) {
189 enthalpyPres = -Cp * JTC * (pressure -Pref);
190 }
191 else if(enableThermalDensity_) {
192 Scalar c1T = gasdentCT1_[regionIdx];
193 Scalar c2T = gasdentCT2_[regionIdx];
194
195 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
196 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
197
198 constexpr const int N = 100; // value is experimental
199 Evaluation deltaP = (pressure - Pref)/N;
200 Evaluation enthalpyPresPrev = 0;
201 for (size_t i = 0; i < N; ++i) {
202 Evaluation Pnew = Pref + i * deltaP;
203 Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rv, Rvw) *
204 (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
205 // see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
206 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
207 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
208 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
209 enthalpyPresPrev = enthalpyPres;
210 }
211 }
212 else {
213 throw std::runtime_error("Requested Joule-thomson calculation but thermal gas density (GASDENT) is not provided");
214 }
215
216 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
217
218 return enthalpy - pressure/density;
219 }
220 }
221
225 template <class Evaluation>
226 Evaluation viscosity(unsigned regionIdx,
227 const Evaluation& temperature,
228 const Evaluation& pressure,
229 const Evaluation& Rv,
230 const Evaluation& Rvw) const
231 {
232
233 const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw);
235 return isothermalMu;
236
237 // compute the viscosity deviation due to temperature
238 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
239 return muGasvisct/viscRef_[regionIdx]*isothermalMu;
240 }
241
245 template <class Evaluation>
246 Evaluation saturatedViscosity(unsigned regionIdx,
247 const Evaluation& temperature,
248 const Evaluation& pressure) const
249 {
250 const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
252 return isothermalMu;
253
254 // compute the viscosity deviation due to temperature
255 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true);
256 return muGasvisct/viscRef_[regionIdx]*isothermalMu;
257 }
258
262 template <class Evaluation>
263 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
264 const Evaluation& temperature,
265 const Evaluation& pressure,
266 const Evaluation& Rv,
267 const Evaluation& /*Rvw*/) const
268 {
269 const Evaluation& Rvw = 0.0;
270 const auto& b =
271 isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
272
274 return b;
275
276 // we use the same approach as for the for water here, but with the OPM-specific
277 // GASDENT keyword.
278 //
279 // TODO: Since gas is quite a bit more compressible than water, it might be
280 // necessary to make GASDENT to a table keyword. If the current temperature
281 // is relatively close to the reference temperature, the current approach
282 // should be good enough, though.
283 Scalar TRef = gasdentRefTemp_[regionIdx];
284 Scalar cT1 = gasdentCT1_[regionIdx];
285 Scalar cT2 = gasdentCT2_[regionIdx];
286 const Evaluation& Y = temperature - TRef;
287
288 return b/(1 + (cT1 + cT2*Y)*Y);
289 }
290
294 template <class Evaluation>
295 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
296 const Evaluation& temperature,
297 const Evaluation& pressure) const
298 {
299 const auto& b =
300 isothermalPvt_->saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
301
303 return b;
304
305 // we use the same approach as for the for water here, but with the OPM-specific
306 // GASDENT keyword.
307 //
308 // TODO: Since gas is quite a bit more compressible than water, it might be
309 // necessary to make GASDENT to a table keyword. If the current temperature
310 // is relatively close to the reference temperature, the current approach
311 // should be good enough, though.
312 Scalar TRef = gasdentRefTemp_[regionIdx];
313 Scalar cT1 = gasdentCT1_[regionIdx];
314 Scalar cT2 = gasdentCT2_[regionIdx];
315 const Evaluation& Y = temperature - TRef;
316
317 return b/(1 + (cT1 + cT2*Y)*Y);
318 }
319
323 template <class Evaluation>
324 Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
325 const Evaluation& /*temperature*/,
326 const Evaluation& /*pressure*/) const
327 { return 0.0; }
328
332 template <class Evaluation = Scalar>
333 Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
334 const Evaluation& /*temperature*/,
335 const Evaluation& /*pressure*/,
336 const Evaluation& /*saltConcentration*/) const
337 { return 0.0; }
338
339
340
348 template <class Evaluation>
349 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
350 const Evaluation& temperature,
351 const Evaluation& pressure) const
352 { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure); }
353
361 template <class Evaluation>
362 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
363 const Evaluation& temperature,
364 const Evaluation& pressure,
365 const Evaluation& oilSaturation,
366 const Evaluation& maxOilSaturation) const
367 { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure, oilSaturation, maxOilSaturation); }
368
376 template <class Evaluation>
377 Evaluation saturationPressure(unsigned regionIdx,
378 const Evaluation& temperature,
379 const Evaluation& pressure) const
380 { return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
381
382 template <class Evaluation>
383 Evaluation diffusionCoefficient(const Evaluation& temperature,
384 const Evaluation& pressure,
385 unsigned compIdx) const
386 {
387 return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
388 }
389 const IsothermalPvt* isoThermalPvt() const
390 { return isothermalPvt_; }
391
392 const Scalar gasReferenceDensity(unsigned regionIdx) const
393 { return isothermalPvt_->gasReferenceDensity(regionIdx); }
394
395 const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
396 { return gasvisctCurves_; }
397
398 const std::vector<Scalar>& viscrefPress() const
399 { return viscrefPress_; }
400
401 const std::vector<Scalar>& viscRef() const
402 { return viscRef_; }
403
404 const std::vector<Scalar>& gasdentRefTemp() const
405 { return gasdentRefTemp_; }
406
407 const std::vector<Scalar>& gasdentCT1() const
408 { return gasdentCT1_; }
409
410 const std::vector<Scalar>& gasdentCT2() const
411 { return gasdentCT2_; }
412
413 const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
414 { return internalEnergyCurves_; }
415
416 bool enableInternalEnergy() const
417 { return enableInternalEnergy_; }
418
419 const std::vector<Scalar>& gasJTRefPres() const
420 { return gasJTRefPres_; }
421
422 const std::vector<Scalar>& gasJTC() const
423 { return gasJTC_; }
424
425
426 bool operator==(const GasPvtThermal<Scalar>& data) const
427 {
428 if (isothermalPvt_ && !data.isothermalPvt_)
429 return false;
430 if (!isothermalPvt_ && data.isothermalPvt_)
431 return false;
432
433 return this->gasvisctCurves() == data.gasvisctCurves() &&
434 this->viscrefPress() == data.viscrefPress() &&
435 this->viscRef() == data.viscRef() &&
436 this->gasdentRefTemp() == data.gasdentRefTemp() &&
437 this->gasdentCT1() == data.gasdentCT1() &&
438 this->gasdentCT2() == data.gasdentCT2() &&
439 this->gasJTRefPres() == data.gasJTRefPres() &&
440 this->gasJTC() == data.gasJTC() &&
441 this->internalEnergyCurves() == data.internalEnergyCurves() &&
442 this->enableThermalDensity() == data.enableThermalDensity() &&
443 this->enableJouleThomson() == data.enableJouleThomson() &&
444 this->enableThermalViscosity() == data.enableThermalViscosity() &&
445 this->enableInternalEnergy() == data.enableInternalEnergy();
446 }
447
448 GasPvtThermal<Scalar>& operator=(const GasPvtThermal<Scalar>& data)
449 {
450 if (data.isothermalPvt_)
451 isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
452 else
453 isothermalPvt_ = nullptr;
454 gasvisctCurves_ = data.gasvisctCurves_;
455 viscrefPress_ = data.viscrefPress_;
456 viscRef_ = data.viscRef_;
457 gasdentRefTemp_ = data.gasdentRefTemp_;
458 gasdentCT1_ = data.gasdentCT1_;
459 gasdentCT2_ = data.gasdentCT2_;
460 gasJTRefPres_ = data.gasJTRefPres_;
461 gasJTC_ = data.gasJTC_;
462 internalEnergyCurves_ = data.internalEnergyCurves_;
463 enableThermalDensity_ = data.enableThermalDensity_;
464 enableJouleThomson_ = data.enableJouleThomson_;
465 enableThermalViscosity_ = data.enableThermalViscosity_;
466 enableInternalEnergy_ = data.enableInternalEnergy_;
467
468 return *this;
469 }
470
471private:
472 IsothermalPvt* isothermalPvt_;
473
474 // The PVT properties needed for temperature dependence of the viscosity. We need
475 // to store one value per PVT region.
476 std::vector<TabulatedOneDFunction> gasvisctCurves_;
477 std::vector<Scalar> viscrefPress_;
478 std::vector<Scalar> viscRef_;
479
480 std::vector<Scalar> gasdentRefTemp_;
481 std::vector<Scalar> gasdentCT1_;
482 std::vector<Scalar> gasdentCT2_;
483
484 std::vector<Scalar> gasJTRefPres_;
485 std::vector<Scalar> gasJTC_;
486
487 std::vector<Scalar> rhoRefO_;
488
489 // piecewise linear curve representing the internal energy of gas
490 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
491
492 bool enableThermalDensity_;
493 bool enableJouleThomson_;
494 bool enableThermalViscosity_;
495 bool enableInternalEnergy_;
496};
497
498} // namespace Opm
499
500#endif
Implements a linearly interpolated scalar function that depends on one variable.
Definition EclipseState.hpp:63
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition GasPvtMultiplexer.hpp:109
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:276
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:265
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:245
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:336
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:285
const Scalar gasReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition GasPvtMultiplexer.hpp:227
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:327
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:256
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:138
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:226
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:362
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil-saturated gas.
Definition GasPvtThermal.hpp:295
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:349
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of water saturated gas.
Definition GasPvtThermal.hpp:333
void initEnd()
Finish initializing the thermal part of the gas phase PVT properties.
Definition GasPvtThermal.hpp:132
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition GasPvtThermal.hpp:110
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:263
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the gas phase is active.
Definition GasPvtThermal.hpp:144
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &) const
Returns the specific internal energy [J/kg] of gas given a set of parameters.
Definition GasPvtThermal.hpp:161
bool enableThermalViscosity() const
Returns true iff the viscosity of the gas phase is temperature dependent.
Definition GasPvtThermal.hpp:150
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition GasPvtThermal.hpp:324
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:246
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the saturation pressure of the gas phase [Pa].
Definition GasPvtThermal.hpp:377
Definition Schedule.hpp:88
Implements a linearly interpolated scalar function that depends on one variable.
Definition Tabulated1DFunction.hpp:51
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30