My Project
OilPvtThermal.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_OIL_PVT_THERMAL_HPP
28 #define OPM_OIL_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 OilPvtMultiplexer;
46 
53 template <class Scalar>
55 {
56 public:
58  typedef OilPvtMultiplexer<Scalar, /*enableThermal=*/false> IsothermalPvt;
59 
61  {
62  enableThermalDensity_ = false;
63  enableThermalViscosity_ = false;
64  enableInternalEnergy_ = false;
65  isothermalPvt_ = nullptr;
66  }
67 
68  OilPvtThermal(IsothermalPvt* isothermalPvt,
69  const std::vector<TabulatedOneDFunction>& oilvisctCurves,
70  const std::vector<Scalar>& viscrefPress,
71  const std::vector<Scalar>& viscrefRs,
72  const std::vector<Scalar>& viscRef,
73  const std::vector<Scalar>& oildentRefTemp,
74  const std::vector<Scalar>& oildentCT1,
75  const std::vector<Scalar>& oildentCT2,
76  const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
79  bool enableInternalEnergy)
80  : isothermalPvt_(isothermalPvt)
81  , oilvisctCurves_(oilvisctCurves)
82  , viscrefPress_(viscrefPress)
83  , viscrefRs_(viscrefRs)
84  , viscRef_(viscRef)
85  , oildentRefTemp_(oildentRefTemp)
86  , oildentCT1_(oildentCT1)
87  , oildentCT2_(oildentCT2)
88  , internalEnergyCurves_(internalEnergyCurves)
89  , enableThermalDensity_(enableThermalDensity)
90  , enableThermalViscosity_(enableThermalViscosity)
91  , enableInternalEnergy_(enableInternalEnergy)
92  { }
93 
94  OilPvtThermal(const OilPvtThermal& data)
95  { *this = data; }
96 
97  ~OilPvtThermal()
98  { delete isothermalPvt_; }
99 
100 #if HAVE_ECL_INPUT
104  void initFromState(const EclipseState& eclState, const Schedule& schedule)
105  {
107  // initialize the isothermal part
109  isothermalPvt_ = new IsothermalPvt;
110  isothermalPvt_->initFromState(eclState, schedule);
111 
113  // initialize the thermal part
115  const auto& tables = eclState.getTableManager();
116 
117  enableThermalDensity_ = tables.OilDenT().size() > 0;
118  enableThermalViscosity_ = tables.hasTables("OILVISCT");
119  enableInternalEnergy_ = tables.hasTables("SPECHEAT");
120 
121  unsigned numRegions = isothermalPvt_->numRegions();
122  setNumRegions(numRegions);
123 
124  // viscosity
125  if (enableThermalViscosity_) {
126  if (tables.getViscrefTable().empty())
127  throw std::runtime_error("VISCREF is required when OILVISCT is present");
128 
129  const auto& oilvisctTables = tables.getOilvisctTables();
130  const auto& viscrefTable = tables.getViscrefTable();
131 
132  assert(oilvisctTables.size() == numRegions);
133  assert(viscrefTable.size() == numRegions);
134 
135  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
136  const auto& TCol = oilvisctTables[regionIdx].getColumn("Temperature").vectorCopy();
137  const auto& muCol = oilvisctTables[regionIdx].getColumn("Viscosity").vectorCopy();
138  oilvisctCurves_[regionIdx].setXYContainers(TCol, muCol);
139 
140  viscrefPress_[regionIdx] = viscrefTable[regionIdx].reference_pressure;
141  viscrefRs_[regionIdx] = viscrefTable[regionIdx].reference_rs;
142 
143  // temperature used to calculate the reference viscosity [K]. the
144  // value does not really matter if the underlying PVT object really
145  // is isothermal...
146  Scalar Tref = 273.15 + 20;
147 
148  // compute the reference viscosity using the isothermal PVT object.
149  viscRef_[regionIdx] =
150  isothermalPvt_->viscosity(regionIdx,
151  Tref,
152  viscrefPress_[regionIdx],
153  viscrefRs_[regionIdx]);
154  }
155  }
156 
157  // temperature dependence of oil density
158  const auto& oilDenT = tables.OilDenT();
159  if (oilDenT.size() > 0) {
160  assert(oilDenT.size() == numRegions);
161  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
162  const auto& record = oilDenT[regionIdx];
163 
164  oildentRefTemp_[regionIdx] = record.T0;
165  oildentCT1_[regionIdx] = record.C1;
166  oildentCT2_[regionIdx] = record.C2;
167  }
168  }
169 
170  if (enableInternalEnergy_) {
171  // the specific internal energy of liquid oil. be aware that ecl only specifies the
172  // heat capacity (via the SPECHEAT keyword) and we need to integrate it
173  // ourselfs to get the internal energy
174  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
175  const auto& specheatTable = tables.getSpecheatTables()[regionIdx];
176  const auto& temperatureColumn = specheatTable.getColumn("TEMPERATURE");
177  const auto& cvOilColumn = specheatTable.getColumn("CV_OIL");
178 
179  std::vector<double> uSamples(temperatureColumn.size());
180 
181  Scalar u = temperatureColumn[0]*cvOilColumn[0];
182  for (size_t i = 0;; ++i) {
183  uSamples[i] = u;
184 
185  if (i >= temperatureColumn.size() - 1)
186  break;
187 
188  // integrate to the heat capacity from the current sampling point to the next
189  // one. this leads to a quadratic polynomial.
190  Scalar c_v0 = cvOilColumn[i];
191  Scalar c_v1 = cvOilColumn[i + 1];
192  Scalar T0 = temperatureColumn[i];
193  Scalar T1 = temperatureColumn[i + 1];
194  u += 0.5*(c_v0 + c_v1)*(T1 - T0);
195  }
196 
197  internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
198  }
199  }
200  }
201 #endif // HAVE_ECL_INPUT
202 
206  void setNumRegions(size_t numRegions)
207  {
208  oilvisctCurves_.resize(numRegions);
209  viscrefPress_.resize(numRegions);
210  viscrefRs_.resize(numRegions);
211  viscRef_.resize(numRegions);
212  internalEnergyCurves_.resize(numRegions);
213  }
214 
218  void initEnd()
219  { }
220 
224  bool enableThermalDensity() const
225  { return enableThermalDensity_; }
226 
231  { return enableThermalViscosity_; }
232 
233  size_t numRegions() const
234  { return viscrefRs_.size(); }
235 
239  template <class Evaluation>
240  Evaluation internalEnergy(unsigned regionIdx,
241  const Evaluation& temperature,
242  const Evaluation&,
243  const Evaluation&) const
244  {
245  if (!enableInternalEnergy_)
246  throw std::runtime_error("Requested the internal energy of oil but it is disabled");
247 
248  // compute the specific internal energy for the specified tempature. We use linear
249  // interpolation here despite the fact that the underlying heat capacities are
250  // piecewise linear (which leads to a quadratic function)
251  return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
252  }
253 
257  template <class Evaluation>
258  Evaluation viscosity(unsigned regionIdx,
259  const Evaluation& temperature,
260  const Evaluation& pressure,
261  const Evaluation& Rs) const
262  {
263  const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rs);
264  if (!enableThermalViscosity())
265  return isothermalMu;
266 
267  // compute the viscosity deviation due to temperature
268  const auto& muOilvisct = oilvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
269  return muOilvisct/viscRef_[regionIdx]*isothermalMu;
270  }
271 
275  template <class Evaluation>
276  Evaluation saturatedViscosity(unsigned regionIdx,
277  const Evaluation& temperature,
278  const Evaluation& pressure) const
279  {
280  const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
281  if (!enableThermalViscosity())
282  return isothermalMu;
283 
284  // compute the viscosity deviation due to temperature
285  const auto& muOilvisct = oilvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
286  return muOilvisct/viscRef_[regionIdx]*isothermalMu;
287  }
288 
289 
293  template <class Evaluation>
294  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
295  const Evaluation& temperature,
296  const Evaluation& pressure,
297  const Evaluation& Rs) const
298  {
299  const auto& b =
300  isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rs);
301 
302  if (!enableThermalDensity())
303  return b;
304 
305  // we use the same approach as for the for water here, but with the OPM-specific
306  // OILDENT keyword.
307  Scalar TRef = oildentRefTemp_[regionIdx];
308  Scalar cT1 = oildentCT1_[regionIdx];
309  Scalar cT2 = oildentCT2_[regionIdx];
310  const Evaluation& Y = temperature - TRef;
311 
312  return b/(1 + (cT1 + cT2*Y)*Y);
313  }
314 
318  template <class Evaluation>
319  Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
320  const Evaluation& temperature,
321  const Evaluation& pressure) const
322  {
323  const auto& b =
324  isothermalPvt_->saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
325 
326  if (!enableThermalDensity())
327  return b;
328 
329  // we use the same approach as for the for water here, but with the OPM-specific
330  // OILDENT keyword.
331  Scalar TRef = oildentRefTemp_[regionIdx];
332  Scalar cT1 = oildentCT1_[regionIdx];
333  Scalar cT2 = oildentCT2_[regionIdx];
334  const Evaluation& Y = temperature - TRef;
335 
336  return b/(1 + (cT1 + cT2*Y)*Y);
337  }
338 
346  template <class Evaluation>
347  Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
348  const Evaluation& temperature,
349  const Evaluation& pressure) const
350  { return isothermalPvt_->saturatedGasDissolutionFactor(regionIdx, temperature, pressure); }
351 
359  template <class Evaluation>
360  Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
361  const Evaluation& temperature,
362  const Evaluation& pressure,
363  const Evaluation& oilSaturation,
364  const Evaluation& maxOilSaturation) const
365  { return isothermalPvt_->saturatedGasDissolutionFactor(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 
388  const IsothermalPvt* isoThermalPvt() const
389  { return isothermalPvt_; }
390 
391  const Scalar oilReferenceDensity(unsigned regionIdx) const
392  { return isothermalPvt_->oilReferenceDensity(regionIdx); }
393 
394  const std::vector<TabulatedOneDFunction>& oilvisctCurves() const
395  { return oilvisctCurves_; }
396 
397  const std::vector<Scalar>& viscrefPress() const
398  { return viscrefPress_; }
399 
400  const std::vector<Scalar>& viscrefRs() const
401  { return viscrefRs_; }
402 
403  const std::vector<Scalar>& viscRef() const
404  { return viscRef_; }
405 
406  const std::vector<Scalar>& oildentRefTemp() const
407  { return oildentRefTemp_; }
408 
409  const std::vector<Scalar>& oildentCT1() const
410  { return oildentCT1_; }
411 
412  const std::vector<Scalar>& oildentCT2() const
413  { return oildentCT2_; }
414 
415  const std::vector<TabulatedOneDFunction> internalEnergyCurves() const
416  { return internalEnergyCurves_; }
417 
418  bool enableInternalEnergy() const
419  { return enableInternalEnergy_; }
420 
421  bool operator==(const OilPvtThermal<Scalar>& data) const
422  {
423  if (isothermalPvt_ && !data.isothermalPvt_)
424  return false;
425  if (!isothermalPvt_ && data.isothermalPvt_)
426  return false;
427 
428  return (!this->isoThermalPvt() ||
429  (*this->isoThermalPvt() == *data.isoThermalPvt())) &&
430  this->oilvisctCurves() == data.oilvisctCurves() &&
431  this->viscrefPress() == data.viscrefPress() &&
432  this->viscrefRs() == data.viscrefRs() &&
433  this->viscRef() == data.viscRef() &&
434  this->oildentRefTemp() == data.oildentRefTemp() &&
435  this->oildentCT1() == data.oildentCT1() &&
436  this->oildentCT2() == data.oildentCT2() &&
437  this->internalEnergyCurves() == data.internalEnergyCurves() &&
438  this->enableThermalDensity() == data.enableThermalDensity() &&
439  this->enableThermalViscosity() == data.enableThermalViscosity() &&
440  this->enableInternalEnergy() == data.enableInternalEnergy();
441  }
442 
443  OilPvtThermal<Scalar>& operator=(const OilPvtThermal<Scalar>& data)
444  {
445  if (data.isothermalPvt_)
446  isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
447  else
448  isothermalPvt_ = nullptr;
449  oilvisctCurves_ = data.oilvisctCurves_;
450  viscrefPress_ = data.viscrefPress_;
451  viscrefRs_ = data.viscrefRs_;
452  viscRef_ = data.viscRef_;
453  oildentRefTemp_ = data.oildentRefTemp_;
454  oildentCT1_ = data.oildentCT1_;
455  oildentCT2_ = data.oildentCT2_;
456  internalEnergyCurves_ = data.internalEnergyCurves_;
457  enableThermalDensity_ = data.enableThermalDensity_;
458  enableThermalViscosity_ = data.enableThermalViscosity_;
459  enableInternalEnergy_ = data.enableInternalEnergy_;
460 
461  return *this;
462  }
463 
464 private:
465  IsothermalPvt* isothermalPvt_;
466 
467  // The PVT properties needed for temperature dependence of the viscosity. We need
468  // to store one value per PVT region.
469  std::vector<TabulatedOneDFunction> oilvisctCurves_;
470  std::vector<Scalar> viscrefPress_;
471  std::vector<Scalar> viscrefRs_;
472  std::vector<Scalar> viscRef_;
473 
474  // The PVT properties needed for temperature dependence of the density.
475  std::vector<Scalar> oildentRefTemp_;
476  std::vector<Scalar> oildentCT1_;
477  std::vector<Scalar> oildentCT2_;
478 
479  // piecewise linear curve representing the internal energy of oil
480  std::vector<TabulatedOneDFunction> internalEnergyCurves_;
481 
482  bool enableThermalDensity_;
483  bool enableThermalViscosity_;
484  bool enableInternalEnergy_;
485 };
486 
487 } // namespace Opm
488 
489 #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 oil phase in the black-oil mod...
Definition: OilPvtMultiplexer.hpp:96
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: OilPvtMultiplexer.hpp:177
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rs) const
Returns the formation volume factor [-] of the fluid phase.
Definition: OilPvtMultiplexer.hpp:219
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: OilPvtMultiplexer.hpp:210
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of the fluid phase.
Definition: OilPvtMultiplexer.hpp:229
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the gas dissolution factor [m^3/m^3] of saturated oil.
Definition: OilPvtMultiplexer.hpp:238
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rs) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: OilPvtMultiplexer.hpp:200
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &Rs) const
Returns the saturation pressure [Pa] of oil given the mass fraction of the gas component in the oil p...
Definition: OilPvtMultiplexer.hpp:262
const Scalar oilReferenceDensity(unsigned regionIdx) const
Return the reference density which are considered by this PVT-object.
Definition: OilPvtMultiplexer.hpp:183
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: OilPvtMultiplexer.hpp:271
This class implements temperature dependence of the PVT properties of oil.
Definition: OilPvtThermal.hpp:55
void initEnd()
Finish initializing the thermal part of the oil phase PVT properties.
Definition: OilPvtThermal.hpp:218
bool enableThermalDensity() const
Returns true iff the density of the oil phase is temperature dependent.
Definition: OilPvtThermal.hpp:224
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &, const Evaluation &) const
Returns the specific internal energy [J/kg] of oil given a set of parameters.
Definition: OilPvtThermal.hpp:240
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: OilPvtThermal.hpp:347
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rs) const
Returns the formation volume factor [-] of the fluid phase.
Definition: OilPvtThermal.hpp:294
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: OilPvtThermal.hpp:276
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rs) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: OilPvtThermal.hpp:258
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: OilPvtThermal.hpp:206
bool enableThermalViscosity() const
Returns true iff the viscosity of the oil phase is temperature dependent.
Definition: OilPvtThermal.hpp:230
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the saturation pressure of the oil phase [Pa].
Definition: OilPvtThermal.hpp:375
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &oilSaturation, const Evaluation &maxOilSaturation) const
Returns the gas dissolution factor [m^3/m^3] of the oil phase.
Definition: OilPvtThermal.hpp:360
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of gas-saturated oil phase.
Definition: OilPvtThermal.hpp:319
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47