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 
36 
37 #if HAVE_ECL_INPUT
38 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
39 #include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
40 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
41 #endif
42 
43 namespace Opm {
44 template <class Scalar, bool enableThermal>
45 class GasPvtMultiplexer;
46 
53 template <class Scalar>
55 {
56 public:
57  typedef GasPvtMultiplexer<Scalar, /*enableThermal=*/false> IsothermalPvt;
59 
61  {
62  enableThermalDensity_ = false;
63  enableThermalViscosity_ = false;
64  enableInternalEnergy_ = false;
65  isothermalPvt_ = nullptr;
66  }
67 
68  GasPvtThermal(IsothermalPvt* isothermalPvt,
69  const std::vector<TabulatedOneDFunction>& gasvisctCurves,
70  const std::vector<Scalar>& gasdentRefTemp,
71  const std::vector<Scalar>& gasdentCT1,
72  const std::vector<Scalar>& gasdentCT2,
73  const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
76  bool enableInternalEnergy)
77  : isothermalPvt_(isothermalPvt)
78  , gasvisctCurves_(gasvisctCurves)
79  , gasdentRefTemp_(gasdentRefTemp)
80  , gasdentCT1_(gasdentCT1)
81  , gasdentCT2_(gasdentCT2)
82  , internalEnergyCurves_(internalEnergyCurves)
83  , enableThermalDensity_(enableThermalDensity)
84  , enableThermalViscosity_(enableThermalViscosity)
85  , enableInternalEnergy_(enableInternalEnergy)
86  { }
87 
88  GasPvtThermal(const GasPvtThermal& data)
89  { *this = data; }
90 
91  ~GasPvtThermal()
92  { delete isothermalPvt_; }
93 
94 #if HAVE_ECL_INPUT
98  void initFromState(const EclipseState& eclState, const Schedule& schedule)
99  {
101  // initialize the isothermal part
103  isothermalPvt_ = new IsothermalPvt;
104  isothermalPvt_->initFromState(eclState, schedule);
105 
107  // initialize the thermal part
109  const auto& tables = eclState.getTableManager();
110 
111  enableThermalDensity_ = tables.GasDenT().size() > 0;
112  enableThermalViscosity_ = tables.hasTables("GASVISCT");
113  enableInternalEnergy_ = tables.hasTables("SPECHEAT");
114 
115  unsigned numRegions = isothermalPvt_->numRegions();
116  setNumRegions(numRegions);
117 
118  // viscosity
119  if (enableThermalViscosity_) {
120  const auto& gasvisctTables = tables.getGasvisctTables();
121  auto gasCompIdx = tables.gas_comp_index();
122  std::string gasvisctColumnName = "Viscosity" + std::to_string(static_cast<long long>(gasCompIdx));
123 
124  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
125  const auto& T = gasvisctTables[regionIdx].getColumn("Temperature").vectorCopy();
126  const auto& mu = gasvisctTables[regionIdx].getColumn(gasvisctColumnName).vectorCopy();
127  gasvisctCurves_[regionIdx].setXYContainers(T, mu);
128  }
129  }
130 
131  // temperature dependence of gas density
132  if (tables.GasDenT().size() > 0) {
133  const auto& gasDenT = tables.GasDenT();
134 
135  assert(gasDenT.size() == numRegions);
136  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
137  const auto& record = gasDenT[regionIdx];
138 
139  gasdentRefTemp_[regionIdx] = record.T0;
140  gasdentCT1_[regionIdx] = record.C1;
141  gasdentCT2_[regionIdx] = record.C2;
142  }
143  enableThermalDensity_ = true;
144  }
145 
146  if (enableInternalEnergy_) {
147  // the specific internal energy of gas. be aware that ecl only specifies the heat capacity
148  // (via the SPECHEAT keyword) and we need to integrate it ourselfs to get the
149  // internal energy
150  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
151  const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
152  const auto& temperatureColumn = specHeatTable.getColumn("TEMPERATURE");
153  const auto& cvGasColumn = specHeatTable.getColumn("CV_GAS");
154 
155  std::vector<double> uSamples(temperatureColumn.size());
156 
157  // the specific enthalpy of vaporization. since ECL does not seem to
158  // feature a proper way to specify this quantity, we use the value for
159  // methane. A proper model would also need to consider the enthalpy
160  // change due to dissolution, i.e. the enthalpies of the gas and oil
161  // phases should depend on the phase composition
162  const Scalar hVap = 480.6e3; // [J / kg]
163 
164  Scalar u = temperatureColumn[0]*cvGasColumn[0] + hVap;
165  for (size_t i = 0;; ++i) {
166  uSamples[i] = u;
167 
168  if (i >= temperatureColumn.size() - 1)
169  break;
170 
171  // integrate to the heat capacity from the current sampling point to the next
172  // one. this leads to a quadratic polynomial.
173  Scalar c_v0 = cvGasColumn[i];
174  Scalar c_v1 = cvGasColumn[i + 1];
175  Scalar T0 = temperatureColumn[i];
176  Scalar T1 = temperatureColumn[i + 1];
177  u += 0.5*(c_v0 + c_v1)*(T1 - T0);
178  }
179 
180  internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
181  }
182  }
183  }
184 #endif // HAVE_ECL_INPUT
185 
189  void setNumRegions(size_t numRegions)
190  {
191  gasvisctCurves_.resize(numRegions);
192  internalEnergyCurves_.resize(numRegions);
193  gasdentRefTemp_.resize(numRegions);
194  gasdentCT1_.resize(numRegions);
195  gasdentCT2_.resize(numRegions);
196  }
197 
201  void initEnd()
202  { }
203 
204  size_t numRegions() const
205  { return gasvisctCurves_.size(); }
206 
210  bool enableThermalDensity() const
211  { return enableThermalDensity_; }
212 
217  { return enableThermalViscosity_; }
218 
222  template <class Evaluation>
223  Evaluation internalEnergy(unsigned regionIdx,
224  const Evaluation& temperature,
225  const Evaluation&,
226  const Evaluation&) const
227  {
228  if (!enableInternalEnergy_)
229  throw std::runtime_error("Requested the internal energy of oil but it is disabled");
230 
231  // compute the specific internal energy for the specified tempature. We use linear
232  // interpolation here despite the fact that the underlying heat capacities are
233  // piecewise linear (which leads to a quadratic function)
234  return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
235  }
236 
240  template <class Evaluation>
241  Evaluation viscosity(unsigned regionIdx,
242  const Evaluation& temperature,
243  const Evaluation& pressure,
244  const Evaluation& Rv) const
245  {
246  if (!enableThermalViscosity())
247  return isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv);
248 
249  // compute the viscosity deviation due to temperature
250  const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature);
251  return muGasvisct;
252  }
253 
257  template <class Evaluation>
258  Evaluation saturatedViscosity(unsigned regionIdx,
259  const Evaluation& temperature,
260  const Evaluation& pressure) const
261  {
262  if (!enableThermalViscosity())
263  return isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
264 
265  // compute the viscosity deviation due to temperature
266  const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true);
267  return muGasvisct;
268  }
269 
273  template <class Evaluation>
274  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
275  const Evaluation& temperature,
276  const Evaluation& pressure,
277  const Evaluation& Rv) const
278  {
279  const auto& b =
280  isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv);
281 
282  if (!enableThermalDensity())
283  return b;
284 
285  // we use the same approach as for the for water here, but with the OPM-specific
286  // GASDENT keyword.
287  //
288  // TODO: Since gas is quite a bit more compressible than water, it might be
289  // necessary to make GASDENT to a table keyword. If the current temperature
290  // is relatively close to the reference temperature, the current approach
291  // should be good enough, though.
292  Scalar TRef = gasdentRefTemp_[regionIdx];
293  Scalar cT1 = gasdentCT1_[regionIdx];
294  Scalar cT2 = gasdentCT2_[regionIdx];
295  const Evaluation& Y = temperature - TRef;
296 
297  return b/(1 + (cT1 + cT2*Y)*Y);
298  }
299 
303  template <class Evaluation>
304  Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
305  const Evaluation& temperature,
306  const Evaluation& pressure) const
307  {
308  const auto& b =
309  isothermalPvt_->saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
310 
311  if (!enableThermalDensity())
312  return b;
313 
314  // we use the same approach as for the for water here, but with the OPM-specific
315  // GASDENT keyword.
316  //
317  // TODO: Since gas is quite a bit more compressible than water, it might be
318  // necessary to make GASDENT to a table keyword. If the current temperature
319  // is relatively close to the reference temperature, the current approach
320  // should be good enough, though.
321  Scalar TRef = gasdentRefTemp_[regionIdx];
322  Scalar cT1 = gasdentCT1_[regionIdx];
323  Scalar cT2 = gasdentCT2_[regionIdx];
324  const Evaluation& Y = temperature - TRef;
325 
326  return b/(1 + (cT1 + cT2*Y)*Y);
327  }
328 
332  template <class Evaluation>
333  Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
334  const Evaluation& /*temperature*/,
335  const Evaluation& /*pressure*/) const
336  { return 0.0; }
337 
338 
346  template <class Evaluation>
347  Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
348  const Evaluation& temperature,
349  const Evaluation& pressure) const
350  { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure); }
351 
359  template <class Evaluation>
360  Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
361  const Evaluation& temperature,
362  const Evaluation& pressure,
363  const Evaluation& oilSaturation,
364  const Evaluation& maxOilSaturation) const
365  { return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure, oilSaturation, maxOilSaturation); }
366 
374  template <class Evaluation>
375  Evaluation saturationPressure(unsigned regionIdx,
376  const Evaluation& temperature,
377  const Evaluation& pressure) const
378  { return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
379 
380  template <class Evaluation>
381  Evaluation diffusionCoefficient(const Evaluation& temperature,
382  const Evaluation& pressure,
383  unsigned compIdx) const
384  {
385  return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
386  }
387  const IsothermalPvt* isoThermalPvt() const
388  { return isothermalPvt_; }
389 
390  const Scalar gasReferenceDensity(unsigned regionIdx) const
391  { return isothermalPvt_->gasReferenceDensity(regionIdx); }
392 
393  const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
394  { return gasvisctCurves_; }
395 
396  const std::vector<Scalar>& gasdentRefTemp() const
397  { return gasdentRefTemp_; }
398 
399  const std::vector<Scalar>& gasdentCT1() const
400  { return gasdentCT1_; }
401 
402  const std::vector<Scalar>& gasdentCT2() const
403  { return gasdentCT2_; }
404 
405  const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
406  { return internalEnergyCurves_; }
407 
408  bool enableInternalEnergy() const
409  { return enableInternalEnergy_; }
410 
411  bool operator==(const GasPvtThermal<Scalar>& data) const
412  {
413  if (isothermalPvt_ && !data.isothermalPvt_)
414  return false;
415  if (!isothermalPvt_ && data.isothermalPvt_)
416  return false;
417 
418  return (!this->isoThermalPvt() ||
419  (*this->isoThermalPvt() == *data.isoThermalPvt())) &&
420  this->gasvisctCurves() == data.gasvisctCurves() &&
421  this->gasdentRefTemp() == data.gasdentRefTemp() &&
422  this->gasdentCT1() == data.gasdentCT1() &&
423  this->gasdentCT2() == data.gasdentCT2() &&
424  this->internalEnergyCurves() == data.internalEnergyCurves() &&
425  this->enableThermalDensity() == data.enableThermalDensity() &&
426  this->enableThermalViscosity() == data.enableThermalViscosity() &&
427  this->enableInternalEnergy() == data.enableInternalEnergy();
428  }
429 
430  GasPvtThermal<Scalar>& operator=(const GasPvtThermal<Scalar>& data)
431  {
432  if (data.isothermalPvt_)
433  isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
434  else
435  isothermalPvt_ = nullptr;
436  gasvisctCurves_ = data.gasvisctCurves_;
437  gasdentRefTemp_ = data.gasdentRefTemp_;
438  gasdentCT1_ = data.gasdentCT1_;
439  gasdentCT2_ = data.gasdentCT2_;
440  internalEnergyCurves_ = data.internalEnergyCurves_;
441  enableThermalDensity_ = data.enableThermalDensity_;
442  enableThermalViscosity_ = data.enableThermalViscosity_;
443  enableInternalEnergy_ = data.enableInternalEnergy_;
444 
445  return *this;
446  }
447 
448 private:
449  IsothermalPvt* isothermalPvt_;
450 
451  // The PVT properties needed for temperature dependence of the viscosity. We need
452  // to store one value per PVT region.
453  std::vector<TabulatedOneDFunction> gasvisctCurves_;
454 
455  std::vector<Scalar> gasdentRefTemp_;
456  std::vector<Scalar> gasdentCT1_;
457  std::vector<Scalar> gasdentCT2_;
458 
459  // piecewise linear curve representing the internal energy of gas
460  std::vector<TabulatedOneDFunction> internalEnergyCurves_;
461 
462  bool enableThermalDensity_;
463  bool enableThermalViscosity_;
464  bool enableInternalEnergy_;
465 };
466 
467 } // namespace Opm
468 
469 #endif
A central place for various physical constants occuring in some equations.
This file provides a wrapper around the "final" C++-2011 statement.
Class implementing cubic splines.
Implements a linearly interpolated scalar function that depends on one variable.
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition: GasPvtMultiplexer.hpp:93
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv) const
Returns the formation volume factor [-] of the fluid phase.
Definition: GasPvtMultiplexer.hpp:242
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:252
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: GasPvtMultiplexer.hpp:200
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:302
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:261
const Scalar gasReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition: GasPvtMultiplexer.hpp:206
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: GasPvtMultiplexer.hpp:223
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:293
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:233
This class implements temperature dependence of the PVT properties of gas.
Definition: GasPvtThermal.hpp:55
bool enableThermalDensity() const
Returns true iff the density of the gas phase is temperature dependent.
Definition: GasPvtThermal.hpp:210
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:360
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil-saturated gas.
Definition: GasPvtThermal.hpp:304
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &, const Evaluation &) const
Returns the specific internal energy [J/kg] of gas given a set of parameters.
Definition: GasPvtThermal.hpp:223
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:347
void initEnd()
Finish initializing the thermal part of the gas phase PVT properties.
Definition: GasPvtThermal.hpp:201
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: GasPvtThermal.hpp:189
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: GasPvtThermal.hpp:241
bool enableThermalViscosity() const
Returns true iff the viscosity of the gas phase is temperature dependent.
Definition: GasPvtThermal.hpp:216
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition: GasPvtThermal.hpp:333
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv) const
Returns the formation volume factor [-] of the fluid phase.
Definition: GasPvtThermal.hpp:274
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:258
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the saturation pressure of the gas phase [Pa].
Definition: GasPvtThermal.hpp:375
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47