My Project
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 isothermalPvt_ = nullptr;
61 }
62
63 GasPvtThermal(IsothermalPvt* isothermalPvt,
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,
74 bool enableInternalEnergy)
75 : isothermalPvt_(isothermalPvt)
76 , gasvisctCurves_(gasvisctCurves)
77 , gasdentRefTemp_(gasdentRefTemp)
78 , gasdentCT1_(gasdentCT1)
79 , gasdentCT2_(gasdentCT2)
80 , gasJTRefPres_(gasJTRefPres)
81 , gasJTC_(gasJTC)
82 , internalEnergyCurves_(internalEnergyCurves)
83 , enableThermalDensity_(enableThermalDensity)
84 , enableJouleThomson_(enableJouleThomson)
85 , enableThermalViscosity_(enableThermalViscosity)
86 , enableInternalEnergy_(enableInternalEnergy)
87 { }
88
89 GasPvtThermal(const GasPvtThermal& data)
90 { *this = data; }
91
93 { delete isothermalPvt_; }
94
95#if HAVE_ECL_INPUT
99 void initFromState(const EclipseState& eclState, const Schedule& schedule);
100#endif // HAVE_ECL_INPUT
101
105 void setNumRegions(size_t numRegions)
106 {
107 gasvisctCurves_.resize(numRegions);
108 internalEnergyCurves_.resize(numRegions);
109 gasdentRefTemp_.resize(numRegions);
110 gasdentCT1_.resize(numRegions);
111 gasdentCT2_.resize(numRegions);
112 gasJTRefPres_.resize(numRegions);
113 gasJTC_.resize(numRegions);
114 rhoRefO_.resize(numRegions);
115 }
116
120 void initEnd()
121 { }
122
123 size_t numRegions() const
124 { return gasvisctCurves_.size(); }
125
130 { return enableThermalDensity_; }
131
136 { return enableJouleThomson_; }
137
142 { return enableThermalViscosity_; }
143
147 template <class Evaluation>
148 Evaluation internalEnergy(unsigned regionIdx,
149 const Evaluation& temperature,
150 const Evaluation& pressure,
151 const Evaluation& Rv,
152 const Evaluation& ) const
153 {
154 if (!enableInternalEnergy_)
155 throw std::runtime_error("Requested the internal energy of gas but it is disabled");
156
157 if (!enableJouleThomson_) {
158 // compute the specific internal energy for the specified tempature. We use linear
159 // interpolation here despite the fact that the underlying heat capacities are
160 // piecewise linear (which leads to a quadratic function)
161 return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
162 }
163 else {
164 Evaluation Tref = gasdentRefTemp_[regionIdx];
165 Evaluation Pref = gasJTRefPres_[regionIdx];
166 Scalar JTC = gasJTC_[regionIdx]; // if JTC is default then JTC is calculated
167 Evaluation Rvw = 0.0;
168
169 Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
170 constexpr const Scalar hVap = 480.6e3; // [J / kg]
171 Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true) - hVap)/temperature;
172 Evaluation density = invB * (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
173
174 Evaluation enthalpyPres;
175 if (JTC != 0) {
176 enthalpyPres = -Cp * JTC * (pressure -Pref);
177 }
178 else if(enableThermalDensity_) {
179 Scalar c1T = gasdentCT1_[regionIdx];
180 Scalar c2T = gasdentCT2_[regionIdx];
181
182 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
183 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
184
185 constexpr const int N = 100; // value is experimental
186 Evaluation deltaP = (pressure - Pref)/N;
187 Evaluation enthalpyPresPrev = 0;
188 for (size_t i = 0; i < N; ++i) {
189 Evaluation Pnew = Pref + i * deltaP;
190 Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rv, Rvw) *
191 (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
192 // see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
193 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
194 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
195 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
196 enthalpyPresPrev = enthalpyPres;
197 }
198 }
199 else {
200 throw std::runtime_error("Requested Joule-thomson calculation but thermal gas density (GASDENT) is not provided");
201 }
202
203 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
204
205 return enthalpy - pressure/density;
206 }
207 }
208
212 template <class Evaluation>
213 Evaluation viscosity(unsigned regionIdx,
214 const Evaluation& temperature,
215 const Evaluation& pressure,
216 const Evaluation& Rv,
217 const Evaluation& Rvw) const
218 {
220 return isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw);
221
222 // compute the viscosity deviation due to temperature
223 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature);
224 return muGasvisct;
225 }
226
230 template <class Evaluation>
231 Evaluation saturatedViscosity(unsigned regionIdx,
232 const Evaluation& temperature,
233 const Evaluation& pressure) const
234 {
236 return isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
237
238 // compute the viscosity deviation due to temperature
239 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true);
240 return muGasvisct;
241 }
242
246 template <class Evaluation>
247 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
248 const Evaluation& temperature,
249 const Evaluation& pressure,
250 const Evaluation& Rv,
251 const Evaluation& /*Rvw*/) const
252 {
253 const Evaluation& Rvw = 0.0;
254 const auto& b =
255 isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
256
258 return b;
259
260 // we use the same approach as for the for water here, but with the OPM-specific
261 // GASDENT keyword.
262 //
263 // TODO: Since gas is quite a bit more compressible than water, it might be
264 // necessary to make GASDENT to a table keyword. If the current temperature
265 // is relatively close to the reference temperature, the current approach
266 // should be good enough, though.
267 Scalar TRef = gasdentRefTemp_[regionIdx];
268 Scalar cT1 = gasdentCT1_[regionIdx];
269 Scalar cT2 = gasdentCT2_[regionIdx];
270 const Evaluation& Y = temperature - TRef;
271
272 return b/(1 + (cT1 + cT2*Y)*Y);
273 }
274
278 template <class Evaluation>
279 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
280 const Evaluation& temperature,
281 const Evaluation& pressure) const
282 {
283 const auto& b =
284 isothermalPvt_->saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
285
287 return b;
288
289 // we use the same approach as for the for water here, but with the OPM-specific
290 // GASDENT keyword.
291 //
292 // TODO: Since gas is quite a bit more compressible than water, it might be
293 // necessary to make GASDENT to a table keyword. If the current temperature
294 // is relatively close to the reference temperature, the current approach
295 // should be good enough, though.
296 Scalar TRef = gasdentRefTemp_[regionIdx];
297 Scalar cT1 = gasdentCT1_[regionIdx];
298 Scalar cT2 = gasdentCT2_[regionIdx];
299 const Evaluation& Y = temperature - TRef;
300
301 return b/(1 + (cT1 + cT2*Y)*Y);
302 }
303
307 template <class Evaluation>
308 Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
309 const Evaluation& /*temperature*/,
310 const Evaluation& /*pressure*/) const
311 { return 0.0; }
312
316 template <class Evaluation = Scalar>
317 Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
318 const Evaluation& /*temperature*/,
319 const Evaluation& /*pressure*/,
320 const Evaluation& /*saltConcentration*/) const
321 { return 0.0; }
322
323
324
332 template <class Evaluation>
333 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
334 const Evaluation& temperature,
335 const Evaluation& pressure) const
336 { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure); }
337
345 template <class Evaluation>
346 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
347 const Evaluation& temperature,
348 const Evaluation& pressure,
349 const Evaluation& oilSaturation,
350 const Evaluation& maxOilSaturation) const
351 { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure, oilSaturation, maxOilSaturation); }
352
360 template <class Evaluation>
361 Evaluation saturationPressure(unsigned regionIdx,
362 const Evaluation& temperature,
363 const Evaluation& pressure) const
364 { return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
365
366 template <class Evaluation>
367 Evaluation diffusionCoefficient(const Evaluation& temperature,
368 const Evaluation& pressure,
369 unsigned compIdx) const
370 {
371 return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
372 }
373 const IsothermalPvt* isoThermalPvt() const
374 { return isothermalPvt_; }
375
376 const Scalar gasReferenceDensity(unsigned regionIdx) const
377 { return isothermalPvt_->gasReferenceDensity(regionIdx); }
378
379 const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
380 { return gasvisctCurves_; }
381
382 const std::vector<Scalar>& gasdentRefTemp() const
383 { return gasdentRefTemp_; }
384
385 const std::vector<Scalar>& gasdentCT1() const
386 { return gasdentCT1_; }
387
388 const std::vector<Scalar>& gasdentCT2() const
389 { return gasdentCT2_; }
390
391 const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
392 { return internalEnergyCurves_; }
393
394 bool enableInternalEnergy() const
395 { return enableInternalEnergy_; }
396
397 const std::vector<Scalar>& gasJTRefPres() const
398 { return gasJTRefPres_; }
399
400 const std::vector<Scalar>& gasJTC() const
401 { return gasJTC_; }
402
403
404 bool operator==(const GasPvtThermal<Scalar>& data) const
405 {
406 if (isothermalPvt_ && !data.isothermalPvt_)
407 return false;
408 if (!isothermalPvt_ && data.isothermalPvt_)
409 return false;
410
411 return this->gasvisctCurves() == data.gasvisctCurves() &&
412 this->gasdentRefTemp() == data.gasdentRefTemp() &&
413 this->gasdentCT1() == data.gasdentCT1() &&
414 this->gasdentCT2() == data.gasdentCT2() &&
415 this->gasJTRefPres() == data.gasJTRefPres() &&
416 this->gasJTC() == data.gasJTC() &&
417 this->internalEnergyCurves() == data.internalEnergyCurves() &&
418 this->enableThermalDensity() == data.enableThermalDensity() &&
419 this->enableJouleThomson() == data.enableJouleThomson() &&
420 this->enableThermalViscosity() == data.enableThermalViscosity() &&
421 this->enableInternalEnergy() == data.enableInternalEnergy();
422 }
423
424 GasPvtThermal<Scalar>& operator=(const GasPvtThermal<Scalar>& data)
425 {
426 if (data.isothermalPvt_)
427 isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
428 else
429 isothermalPvt_ = nullptr;
430 gasvisctCurves_ = data.gasvisctCurves_;
431 gasdentRefTemp_ = data.gasdentRefTemp_;
432 gasdentCT1_ = data.gasdentCT1_;
433 gasdentCT2_ = data.gasdentCT2_;
434 gasJTRefPres_ = data.gasJTRefPres_;
435 gasJTC_ = data.gasJTC_;
436 internalEnergyCurves_ = data.internalEnergyCurves_;
437 enableThermalDensity_ = data.enableThermalDensity_;
438 enableJouleThomson_ = data.enableJouleThomson_;
439 enableThermalViscosity_ = data.enableThermalViscosity_;
440 enableInternalEnergy_ = data.enableInternalEnergy_;
441
442 return *this;
443 }
444
445private:
446 IsothermalPvt* isothermalPvt_;
447
448 // The PVT properties needed for temperature dependence of the viscosity. We need
449 // to store one value per PVT region.
450 std::vector<TabulatedOneDFunction> gasvisctCurves_;
451
452 std::vector<Scalar> gasdentRefTemp_;
453 std::vector<Scalar> gasdentCT1_;
454 std::vector<Scalar> gasdentCT2_;
455
456 std::vector<Scalar> gasJTRefPres_;
457 std::vector<Scalar> gasJTC_;
458
459 std::vector<Scalar> rhoRefO_;
460
461 // piecewise linear curve representing the internal energy of gas
462 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
463
464 bool enableThermalDensity_;
465 bool enableJouleThomson_;
466 bool enableThermalViscosity_;
467 bool enableInternalEnergy_;
468};
469
470} // namespace Opm
471
472#endif
Definition: EclipseState.hpp:55
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition: GasPvtMultiplexer.hpp:102
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:256
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:245
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:225
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:316
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:265
const Scalar gasReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition: GasPvtMultiplexer.hpp:207
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:307
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:236
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:129
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:213
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:346
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil-saturated gas.
Definition: GasPvtThermal.hpp:279
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:333
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:317
void initEnd()
Finish initializing the thermal part of the gas phase PVT properties.
Definition: GasPvtThermal.hpp:120
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: GasPvtThermal.hpp:105
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:247
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the gas phase is active.
Definition: GasPvtThermal.hpp:135
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:148
bool enableThermalViscosity() const
Returns true iff the viscosity of the gas phase is temperature dependent.
Definition: GasPvtThermal.hpp:141
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: GasPvtThermal.hpp:308
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:231
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the saturation pressure of the gas phase [Pa].
Definition: GasPvtThermal.hpp:361
Definition: Schedule.hpp:130
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