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
32#include <cstddef>
33
34namespace Opm {
35
36#if HAVE_ECL_INPUT
37class EclipseState;
38class Schedule;
39#endif
40
41template <class Scalar, bool enableThermal>
42class GasPvtMultiplexer;
43
50template <class Scalar>
52{
53public:
54 using IsothermalPvt = GasPvtMultiplexer<Scalar, /*enableThermal=*/false>;
56
57 GasPvtThermal() = default;
58
59 GasPvtThermal(IsothermalPvt* isothermalPvt,
60 const std::vector<TabulatedOneDFunction>& gasvisctCurves,
61 const std::vector<Scalar>& viscrefPress,
62 const std::vector<Scalar>& viscRef,
63 const std::vector<Scalar>& gasdentRefTemp,
64 const std::vector<Scalar>& gasdentCT1,
65 const std::vector<Scalar>& gasdentCT2,
66 const std::vector<Scalar>& gasJTRefPres,
67 const std::vector<Scalar>& gasJTC,
68 const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
72 bool enableInternalEnergy)
73 : isothermalPvt_(isothermalPvt)
74 , gasvisctCurves_(gasvisctCurves)
75 , viscrefPress_(viscrefPress)
76 , viscRef_(viscRef)
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(std::size_t numRegions);
106
107 void setVapPars(const Scalar par1, const Scalar par2)
108 {
109 isothermalPvt_->setVapPars(par1, par2);
110 }
111
115 void initEnd()
116 { }
117
122 { return enableThermalDensity_; }
123
128 { return enableJouleThomson_; }
129
134 { return enableThermalViscosity_; }
135
136 std::size_t numRegions() const
137 { return viscrefPress_.size(); }
138
142 template <class Evaluation>
143 Evaluation internalEnergy(unsigned regionIdx,
144 const Evaluation& temperature,
145 const Evaluation& pressure,
146 const Evaluation& Rv,
147 [[maybe_unused]] const Evaluation& RvW) const
148 {
149 if (!enableInternalEnergy_) {
150 throw std::runtime_error("Requested the internal energy of gas but it is disabled");
151 }
152
153 if (!enableJouleThomson_) {
154 // compute the specific internal energy for the specified tempature. We use linear
155 // interpolation here despite the fact that the underlying heat capacities are
156 // piecewise linear (which leads to a quadratic function)
157 return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
158 }
159 else {
160 // NB should probably be revisited with adding more or restrict to linear internal energy
161 OpmLog::warning("Experimental code for jouleThomson: simulation will be slower");
162 Evaluation Tref = gasdentRefTemp_[regionIdx];
163 Evaluation Pref = gasJTRefPres_[regionIdx];
164 Scalar JTC = gasJTC_[regionIdx]; // if JTC is default then JTC is calculaited
165 Evaluation Rvw = 0.0;
166
167 Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
168 // NB this assumes internalEnergyCurve(0) = 0 derivative should be used add Cp table ??
169 Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true))/temperature;
170 Evaluation density = invB * (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
171
172 Evaluation enthalpyPres;
173 if (JTC != 0) {
174 enthalpyPres = -Cp * JTC * (pressure -Pref);
175 }
176 else if (enableThermalDensity_) {
177 Scalar c1T = gasdentCT1_[regionIdx];
178 Scalar c2T = gasdentCT2_[regionIdx];
179
180 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
181 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
182
183 constexpr const int N = 100; // value is experimental
184 Evaluation deltaP = (pressure - Pref)/N;
185 Evaluation enthalpyPresPrev = 0;
186 for (std::size_t i = 0; i < N; ++i) {
187 Evaluation Pnew = Pref + i * deltaP;
188 Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rv, Rvw) *
189 (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
190 // see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
191 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
192 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
193 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
194 enthalpyPresPrev = enthalpyPres;
195 }
196 }
197 else {
198 throw std::runtime_error("Requested Joule-thomson calculation but thermal "
199 "gas density (GASDENT) is not provided");
200 }
201
202 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
203
204 return enthalpy - pressure/density;
205 }
206 }
207
211 template <class Evaluation>
212 Evaluation viscosity(unsigned regionIdx,
213 const Evaluation& temperature,
214 const Evaluation& pressure,
215 const Evaluation& Rv,
216 const Evaluation& Rvw) const
217 {
218 const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw);
220 return isothermalMu;
221
222 // compute the viscosity deviation due to temperature
223 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
224 return muGasvisct/viscRef_[regionIdx]*isothermalMu;
225 }
226
230 template <class Evaluation>
231 Evaluation saturatedViscosity(unsigned regionIdx,
232 const Evaluation& temperature,
233 const Evaluation& pressure) const
234 {
235 const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
237 return isothermalMu;
238
239 // compute the viscosity deviation due to temperature
240 const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true);
241 return muGasvisct/viscRef_[regionIdx]*isothermalMu;
242 }
243
247 template <class Evaluation>
248 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
249 const Evaluation& temperature,
250 const Evaluation& pressure,
251 const Evaluation& Rv,
252 const Evaluation& /*Rvw*/) const
253 {
254 const Evaluation& Rvw = 0.0;
255 const auto& b =
256 isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature,
257 pressure, Rv, Rvw);
258
259 if (!enableThermalDensity()) {
260 return b;
261 }
262
263 // we use the same approach as for the for water here, but with the OPM-specific
264 // GASDENT keyword.
265 //
266 // TODO: Since gas is quite a bit more compressible than water, it might be
267 // necessary to make GASDENT to a table keyword. If the current temperature
268 // is relatively close to the reference temperature, the current approach
269 // should be good enough, though.
270 Scalar TRef = gasdentRefTemp_[regionIdx];
271 Scalar cT1 = gasdentCT1_[regionIdx];
272 Scalar cT2 = gasdentCT2_[regionIdx];
273 const Evaluation& Y = temperature - TRef;
274
275 return b / (1 + (cT1 + cT2 * Y) * Y);
276 }
277
281 template <class Evaluation>
282 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
283 const Evaluation& temperature,
284 const Evaluation& pressure) const
285 {
286 const auto& b =
287 isothermalPvt_->saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
288
289 if (!enableThermalDensity()) {
290 return b;
291 }
292
293 // we use the same approach as for the for water here, but with the OPM-specific
294 // GASDENT keyword.
295 //
296 // TODO: Since gas is quite a bit more compressible than water, it might be
297 // necessary to make GASDENT to a table keyword. If the current temperature
298 // is relatively close to the reference temperature, the current approach
299 // should be good enough, though.
300 Scalar TRef = gasdentRefTemp_[regionIdx];
301 Scalar cT1 = gasdentCT1_[regionIdx];
302 Scalar cT2 = gasdentCT2_[regionIdx];
303 const Evaluation& Y = temperature - TRef;
304
305 return b / (1 + (cT1 + cT2 * Y) * Y);
306 }
307
311 template <class Evaluation>
312 Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
313 const Evaluation& /*temperature*/,
314 const Evaluation& /*pressure*/) const
315 { return 0.0; }
316
320 template <class Evaluation = Scalar>
321 Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
322 const Evaluation& /*temperature*/,
323 const Evaluation& /*pressure*/,
324 const Evaluation& /*saltConcentration*/) const
325 { return 0.0; }
326
334 template <class Evaluation>
335 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
336 const Evaluation& temperature,
337 const Evaluation& pressure) const
338 { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure); }
339
347 template <class Evaluation>
348 Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
349 const Evaluation& temperature,
350 const Evaluation& pressure,
351 const Evaluation& oilSaturation,
352 const Evaluation& maxOilSaturation) const
353 {
354 return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature,
355 pressure, oilSaturation,
356 maxOilSaturation);
357 }
358
366 template <class Evaluation>
367 Evaluation saturationPressure(unsigned regionIdx,
368 const Evaluation& temperature,
369 const Evaluation& pressure) const
370 { return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
371
372 template <class Evaluation>
373 Evaluation diffusionCoefficient(const Evaluation& temperature,
374 const Evaluation& pressure,
375 unsigned compIdx) const
376 {
377 return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
378 }
379 const IsothermalPvt* isoThermalPvt() const
380 { return isothermalPvt_; }
381
382 Scalar gasReferenceDensity(unsigned regionIdx) const
383 { return isothermalPvt_->gasReferenceDensity(regionIdx); }
384
385 Scalar hVap(unsigned regionIdx) const
386 { return this->hVap_[regionIdx]; }
387
388 const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
389 { return gasvisctCurves_; }
390
391 const std::vector<Scalar>& viscrefPress() const
392 { return viscrefPress_; }
393
394 const std::vector<Scalar>& viscRef() const
395 { return viscRef_; }
396
397 const std::vector<Scalar>& gasdentRefTemp() const
398 { return gasdentRefTemp_; }
399
400 const std::vector<Scalar>& gasdentCT1() const
401 { return gasdentCT1_; }
402
403 const std::vector<Scalar>& gasdentCT2() const
404 { return gasdentCT2_; }
405
406 const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
407 { return internalEnergyCurves_; }
408
409 bool enableInternalEnergy() const
410 { return enableInternalEnergy_; }
411
412 const std::vector<Scalar>& gasJTRefPres() const
413 { return gasJTRefPres_; }
414
415 const std::vector<Scalar>& gasJTC() const
416 { return gasJTC_; }
417
418 bool operator==(const GasPvtThermal<Scalar>& data) const;
419
420 GasPvtThermal<Scalar>& operator=(const GasPvtThermal<Scalar>& data);
421
422private:
423 IsothermalPvt* isothermalPvt_{nullptr};
424
425 // The PVT properties needed for temperature dependence of the viscosity. We need
426 // to store one value per PVT region.
427 std::vector<TabulatedOneDFunction> gasvisctCurves_{};
428 std::vector<Scalar> viscrefPress_{};
429 std::vector<Scalar> viscRef_{};
430
431 std::vector<Scalar> gasdentRefTemp_{};
432 std::vector<Scalar> gasdentCT1_{};
433 std::vector<Scalar> gasdentCT2_{};
434
435 std::vector<Scalar> gasJTRefPres_{};
436 std::vector<Scalar> gasJTC_{};
437
438 std::vector<Scalar> rhoRefO_{};
439 std::vector<Scalar> hVap_{};
440
441 // piecewise linear curve representing the internal energy of gas
442 std::vector<TabulatedOneDFunction> internalEnergyCurves_{};
443
444 bool enableThermalDensity_{false};
445 bool enableJouleThomson_{false};
446 bool enableThermalViscosity_{false};
447 bool enableInternalEnergy_{false};
448};
449
450} // namespace Opm
451
452#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:111
Scalar gasReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition GasPvtMultiplexer.cpp:57
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:209
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:198
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:178
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:269
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:218
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:260
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:189
This class implements temperature dependence of the PVT properties of gas.
Definition GasPvtThermal.hpp:52
bool enableThermalDensity() const
Returns true iff the density of the gas phase is temperature dependent.
Definition GasPvtThermal.hpp:121
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:212
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:348
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil-saturated gas.
Definition GasPvtThermal.hpp:282
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:335
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:321
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &RvW) const
Returns the specific internal energy [J/kg] of gas given a set of parameters.
Definition GasPvtThermal.hpp:143
void setNumRegions(std::size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition GasPvtThermal.cpp:192
void initEnd()
Finish initializing the thermal part of the gas phase PVT properties.
Definition GasPvtThermal.hpp:115
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:248
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the gas phase is active.
Definition GasPvtThermal.hpp:127
bool enableThermalViscosity() const
Returns true iff the viscosity of the gas phase is temperature dependent.
Definition GasPvtThermal.hpp:133
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition GasPvtThermal.hpp:312
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:367
Definition Schedule.hpp:101
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