My Project
WaterPvtThermal.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_WATER_PVT_THERMAL_HPP
28 #define OPM_WATER_PVT_THERMAL_HPP
29 
31 
32 #if HAVE_ECL_INPUT
33 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
34 #include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
35 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
36 #endif
37 
38 namespace Opm {
39 template <class Scalar, bool enableThermal, bool enableBrine>
40 class WaterPvtMultiplexer;
41 
48 template <class Scalar, bool enableBrine>
50 {
51 public:
53  using IsothermalPvt = WaterPvtMultiplexer<Scalar, /*enableThermal=*/false, enableBrine>;
54 
56  {
57  enableThermalDensity_ = false;
58  enableThermalViscosity_ = false;
59  enableInternalEnergy_ = false;
60  isothermalPvt_ = nullptr;
61  }
62 
63  WaterPvtThermal(IsothermalPvt* isothermalPvt,
64  const std::vector<Scalar>& viscrefPress,
65  const std::vector<Scalar>& watdentRefTemp,
66  const std::vector<Scalar>& watdentCT1,
67  const std::vector<Scalar>& watdentCT2,
68  const std::vector<Scalar>& watJTRefPres,
69  const std::vector<Scalar>& watJTC,
70  const std::vector<Scalar>& pvtwRefPress,
71  const std::vector<Scalar>& pvtwRefB,
72  const std::vector<Scalar>& pvtwCompressibility,
73  const std::vector<Scalar>& pvtwViscosity,
74  const std::vector<Scalar>& pvtwViscosibility,
75  const std::vector<TabulatedOneDFunction>& watvisctCurves,
76  const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
78  bool enableJouleThomson,
80  bool enableInternalEnergy)
81  : isothermalPvt_(isothermalPvt)
82  , viscrefPress_(viscrefPress)
83  , watdentRefTemp_(watdentRefTemp)
84  , watdentCT1_(watdentCT1)
85  , watdentCT2_(watdentCT2)
86  , watJTRefPres_(watJTRefPres)
87  , watJTC_(watJTC)
88  , pvtwRefPress_(pvtwRefPress)
89  , pvtwRefB_(pvtwRefB)
90  , pvtwCompressibility_(pvtwCompressibility)
91  , pvtwViscosity_(pvtwViscosity)
92  , pvtwViscosibility_(pvtwViscosibility)
93  , watvisctCurves_(watvisctCurves)
94  , internalEnergyCurves_(internalEnergyCurves)
95  , enableThermalDensity_(enableThermalDensity)
96  , enableJouleThomson_(enableJouleThomson)
97  , enableThermalViscosity_(enableThermalViscosity)
98  , enableInternalEnergy_(enableInternalEnergy)
99  { }
100 
101  WaterPvtThermal(const WaterPvtThermal& data)
102  { *this = data; }
103 
104  ~WaterPvtThermal()
105  { delete isothermalPvt_; }
106 
107 #if HAVE_ECL_INPUT
111  void initFromState(const EclipseState& eclState, const Schedule& schedule)
112  {
114  // initialize the isothermal part
116  isothermalPvt_ = new IsothermalPvt;
117  isothermalPvt_->initFromState(eclState, schedule);
118 
120  // initialize the thermal part
122  const auto& tables = eclState.getTableManager();
123 
124  enableThermalDensity_ = tables.WatDenT().size() > 0;
125  enableJouleThomson_ = tables.WatJT().size() > 0;
126  enableThermalViscosity_ = tables.hasTables("WATVISCT");
127  enableInternalEnergy_ = tables.hasTables("SPECHEAT");
128 
129  unsigned numRegions = isothermalPvt_->numRegions();
130  setNumRegions(numRegions);
131 
132  if (enableThermalDensity_) {
133  const auto& watDenT = tables.WatDenT();
134 
135  assert(watDenT.size() == numRegions);
136  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
137  const auto& record = watDenT[regionIdx];
138 
139  watdentRefTemp_[regionIdx] = record.T0;
140  watdentCT1_[regionIdx] = record.C1;
141  watdentCT2_[regionIdx] = record.C2;
142  }
143 
144  const auto& pvtwTables = tables.getPvtwTable();
145 
146  assert(pvtwTables.size() == numRegions);
147 
148  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
149  pvtwRefPress_[regionIdx] = pvtwTables[regionIdx].reference_pressure;
150  pvtwRefB_[regionIdx] = pvtwTables[regionIdx].volume_factor;
151  }
152  }
153 
154  // Joule Thomson
155  if (enableJouleThomson_) {
156  const auto& watJT = tables.WatJT();
157 
158  assert(watJT.size() == numRegions);
159  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
160  const auto& record = watJT[regionIdx];
161 
162  watJTRefPres_[regionIdx] = record.P0;
163  watJTC_[regionIdx] = record.C1;
164  }
165  }
166 
167  if (enableThermalViscosity_) {
168  if (tables.getViscrefTable().empty())
169  throw std::runtime_error("VISCREF is required when WATVISCT is present");
170 
171  const auto& watvisctTables = tables.getWatvisctTables();
172  const auto& viscrefTables = tables.getViscrefTable();
173 
174  const auto& pvtwTables = tables.getPvtwTable();
175 
176  assert(pvtwTables.size() == numRegions);
177  assert(watvisctTables.size() == numRegions);
178  assert(viscrefTables.size() == numRegions);
179 
180  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
181  const auto& T = watvisctTables[regionIdx].getColumn("Temperature").vectorCopy();
182  const auto& mu = watvisctTables[regionIdx].getColumn("Viscosity").vectorCopy();
183  watvisctCurves_[regionIdx].setXYContainers(T, mu);
184 
185  viscrefPress_[regionIdx] = viscrefTables[regionIdx].reference_pressure;
186  }
187 
188  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
189  pvtwViscosity_[regionIdx] = pvtwTables[regionIdx].viscosity;
190  pvtwViscosibility_[regionIdx] = pvtwTables[regionIdx].viscosibility;
191  }
192  }
193 
194  if (enableInternalEnergy_) {
195  // the specific internal energy of liquid water. be aware that ecl only specifies the heat capacity
196  // (via the SPECHEAT keyword) and we need to integrate it ourselfs to get the
197  // internal energy
198  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
199  const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
200  const auto& temperatureColumn = specHeatTable.getColumn("TEMPERATURE");
201  const auto& cvWaterColumn = specHeatTable.getColumn("CV_WATER");
202 
203  std::vector<double> uSamples(temperatureColumn.size());
204 
205  Scalar u = temperatureColumn[0]*cvWaterColumn[0];
206  for (size_t i = 0;; ++i) {
207  uSamples[i] = u;
208 
209  if (i >= temperatureColumn.size() - 1)
210  break;
211 
212  // integrate to the heat capacity from the current sampling point to the next
213  // one. this leads to a quadratic polynomial.
214  Scalar c_v0 = cvWaterColumn[i];
215  Scalar c_v1 = cvWaterColumn[i + 1];
216  Scalar T0 = temperatureColumn[i];
217  Scalar T1 = temperatureColumn[i + 1];
218  u += 0.5*(c_v0 + c_v1)*(T1 - T0);
219  }
220 
221  internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
222  }
223  }
224  }
225 #endif // HAVE_ECL_INPUT
226 
230  void setNumRegions(size_t numRegions)
231  {
232  pvtwRefPress_.resize(numRegions);
233  pvtwRefB_.resize(numRegions);
234  pvtwCompressibility_.resize(numRegions);
235  pvtwViscosity_.resize(numRegions);
236  pvtwViscosibility_.resize(numRegions);
237  viscrefPress_.resize(numRegions);
238  watvisctCurves_.resize(numRegions);
239  watdentRefTemp_.resize(numRegions);
240  watdentCT1_.resize(numRegions);
241  watdentCT2_.resize(numRegions);
242  watJTRefPres_.resize(numRegions);
243  watJTC_.resize(numRegions);
244  internalEnergyCurves_.resize(numRegions);
245  }
246 
250  void initEnd()
251  { }
252 
256  bool enableThermalDensity() const
257  { return enableThermalDensity_; }
258 
262  bool enableJouleThomsony() const
263  { return enableJouleThomson_; }
264 
269  { return enableThermalViscosity_; }
270 
271  size_t numRegions() const
272  { return pvtwRefPress_.size(); }
273 
277  template <class Evaluation>
278  Evaluation internalEnergy(unsigned regionIdx,
279  const Evaluation& temperature,
280  const Evaluation& pressure,
281  const Evaluation& saltconcentration) const
282  {
283  if (!enableInternalEnergy_)
284  throw std::runtime_error("Requested the internal energy of water but it is disabled");
285 
286  if (!enableJouleThomson_) {
287  // compute the specific internal energy for the specified tempature. We use linear
288  // interpolation here despite the fact that the underlying heat capacities are
289  // piecewise linear (which leads to a quadratic function)
290  return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
291  }
292  else {
293  Evaluation Tref = watdentRefTemp_[regionIdx];
294  Evaluation Pref = watJTRefPres_[regionIdx];
295  Scalar JTC =watJTC_[regionIdx]; // if JTC is default then JTC is calculated
296 
297  Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration);
298  Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
299  Evaluation density = invB * waterReferenceDensity(regionIdx);
300 
301  Evaluation enthalpyPres;
302  if (JTC != 0) {
303  enthalpyPres = -Cp * JTC * (pressure -Pref);
304  }
305  else if(enableThermalDensity_) {
306  Scalar c1T = watdentCT1_[regionIdx];
307  Scalar c2T = watdentCT2_[regionIdx];
308 
309  Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
310  (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
311 
312  constexpr const int N = 100; // value is experimental
313  Evaluation deltaP = (pressure - Pref)/N;
314  Evaluation enthalpyPresPrev = 0;
315  for (size_t i = 0; i < N; ++i) {
316  Evaluation Pnew = Pref + i * deltaP;
317  Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, saltconcentration) * waterReferenceDensity(regionIdx);
318  Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
319  Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
320  enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
321  enthalpyPresPrev = enthalpyPres;
322  }
323  }
324  else {
325  throw std::runtime_error("Requested Joule-thomson calculation but thermal water density (WATDENT) is not provided");
326  }
327 
328  Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
329 
330  return enthalpy - pressure/density;
331  }
332  }
333 
337  template <class Evaluation>
338  Evaluation viscosity(unsigned regionIdx,
339  const Evaluation& temperature,
340  const Evaluation& pressure,
341  const Evaluation& saltconcentration) const
342  {
343  const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, saltconcentration);
344  if (!enableThermalViscosity())
345  return isothermalMu;
346 
347  Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
348  Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
349 
350  // compute the viscosity deviation due to temperature
351  const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
352  return isothermalMu * muWatvisct/muRef;
353  }
354 
358  template <class Evaluation>
359  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
360  const Evaluation& temperature,
361  const Evaluation& pressure,
362  const Evaluation& saltconcentration) const
363  {
364  if (!enableThermalDensity())
365  return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration);
366 
367  Scalar BwRef = pvtwRefB_[regionIdx];
368  Scalar TRef = watdentRefTemp_[regionIdx];
369  const Evaluation& X = pvtwCompressibility_[regionIdx]*(pressure - pvtwRefPress_[regionIdx]);
370  Scalar cT1 = watdentCT1_[regionIdx];
371  Scalar cT2 = watdentCT2_[regionIdx];
372  const Evaluation& Y = temperature - TRef;
373 
374  // this is inconsistent with the density calculation of water in the isothermal
375  // case (it misses the quadratic pressure term), but it is the equation given in
376  // the documentation.
377  return 1.0/(((1 - X)*(1 + cT1*Y + cT2*Y*Y))*BwRef);
378  }
379 
380  const IsothermalPvt* isoThermalPvt() const
381  { return isothermalPvt_; }
382 
383  const Scalar waterReferenceDensity(unsigned regionIdx) const
384  { return isothermalPvt_->waterReferenceDensity(regionIdx); }
385 
386  const std::vector<Scalar>& viscrefPress() const
387  { return viscrefPress_; }
388 
389  const std::vector<Scalar>& watdentRefTemp() const
390  { return watdentRefTemp_; }
391 
392  const std::vector<Scalar>& watdentCT1() const
393  { return watdentCT1_; }
394 
395  const std::vector<Scalar>& watdentCT2() const
396  { return watdentCT2_; }
397 
398  const std::vector<Scalar>& pvtwRefPress() const
399  { return pvtwRefPress_; }
400 
401  const std::vector<Scalar>& pvtwRefB() const
402  { return pvtwRefB_; }
403 
404  const std::vector<Scalar>& pvtwCompressibility() const
405  { return pvtwCompressibility_; }
406 
407  const std::vector<Scalar>& pvtwViscosity() const
408  { return pvtwViscosity_; }
409 
410  const std::vector<Scalar>& pvtwViscosibility() const
411  { return pvtwViscosibility_; }
412 
413  const std::vector<TabulatedOneDFunction>& watvisctCurves() const
414  { return watvisctCurves_; }
415 
416  const std::vector<TabulatedOneDFunction> internalEnergyCurves() const
417  { return internalEnergyCurves_; }
418 
419  bool enableInternalEnergy() const
420  { return enableInternalEnergy_; }
421 
422  const std::vector<Scalar>& watJTRefPres() const
423  { return watJTRefPres_; }
424 
425  const std::vector<Scalar>& watJTC() const
426  { return watJTC_; }
427 
428 
429  bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
430  {
431  if (isothermalPvt_ && !data.isothermalPvt_)
432  return false;
433  if (!isothermalPvt_ && data.isothermalPvt_)
434  return false;
435 
436  return (!this->isoThermalPvt() ||
437  (*this->isoThermalPvt() == *data.isoThermalPvt())) &&
438  this->viscrefPress() == data.viscrefPress() &&
439  this->watdentRefTemp() == data.watdentRefTemp() &&
440  this->watdentCT1() == data.watdentCT1() &&
441  this->watdentCT2() == data.watdentCT2() &&
442  this->watdentCT2() == data.watdentCT2() &&
443  this->watJTRefPres() == data.watJTRefPres() &&
444  this->watJT() == data.watJT() &&
445  this->watJTC() == data.watJTC() &&
446  this->pvtwRefPress() == data.pvtwRefPress() &&
447  this->pvtwRefB() == data.pvtwRefB() &&
448  this->pvtwCompressibility() == data.pvtwCompressibility() &&
449  this->pvtwViscosity() == data.pvtwViscosity() &&
450  this->pvtwViscosibility() == data.pvtwViscosibility() &&
451  this->watvisctCurves() == data.watvisctCurves() &&
452  this->internalEnergyCurves() == data.internalEnergyCurves() &&
453  this->enableThermalDensity() == data.enableThermalDensity() &&
454  this->enableJouleThomson() == data.enableJouleThomson() &&
455  this->enableThermalViscosity() == data.enableThermalViscosity() &&
456  this->enableInternalEnergy() == data.enableInternalEnergy();
457  }
458 
459  WaterPvtThermal<Scalar, enableBrine>& operator=(const WaterPvtThermal<Scalar, enableBrine>& data)
460  {
461  if (data.isothermalPvt_)
462  isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
463  else
464  isothermalPvt_ = nullptr;
465  viscrefPress_ = data.viscrefPress_;
466  watdentRefTemp_ = data.watdentRefTemp_;
467  watdentCT1_ = data.watdentCT1_;
468  watdentCT2_ = data.watdentCT2_;
469  watJTRefPres_ = data.watJTRefPres_;
470  watJTC_ = data.watJTC_;
471  pvtwRefPress_ = data.pvtwRefPress_;
472  pvtwRefB_ = data.pvtwRefB_;
473  pvtwCompressibility_ = data.pvtwCompressibility_;
474  pvtwViscosity_ = data.pvtwViscosity_;
475  pvtwViscosibility_ = data.pvtwViscosibility_;
476  watvisctCurves_ = data.watvisctCurves_;
477  internalEnergyCurves_ = data.internalEnergyCurves_;
478  enableThermalDensity_ = data.enableThermalDensity_;
479  enableJouleThomson_ = data.enableJouleThomson_;
480  enableThermalViscosity_ = data.enableThermalViscosity_;
481  enableInternalEnergy_ = data.enableInternalEnergy_;
482 
483  return *this;
484  }
485 
486 private:
487  IsothermalPvt* isothermalPvt_;
488 
489  // The PVT properties needed for temperature dependence. We need to store one
490  // value per PVT region.
491  std::vector<Scalar> viscrefPress_;
492 
493  std::vector<Scalar> watdentRefTemp_;
494  std::vector<Scalar> watdentCT1_;
495  std::vector<Scalar> watdentCT2_;
496 
497  std::vector<Scalar> watJTRefPres_;
498  std::vector<Scalar> watJTC_;
499 
500  std::vector<Scalar> pvtwRefPress_;
501  std::vector<Scalar> pvtwRefB_;
502  std::vector<Scalar> pvtwCompressibility_;
503  std::vector<Scalar> pvtwViscosity_;
504  std::vector<Scalar> pvtwViscosibility_;
505 
506  std::vector<TabulatedOneDFunction> watvisctCurves_;
507 
508  // piecewise linear curve representing the internal energy of water
509  std::vector<TabulatedOneDFunction> internalEnergyCurves_;
510 
511  bool enableThermalDensity_;
512  bool enableJouleThomson_;
513  bool enableThermalViscosity_;
514  bool enableInternalEnergy_;
515 };
516 
517 } // namespace Opm
518 
519 #endif
Implements a linearly interpolated scalar function that depends on one variable.
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
This class represents the Pressure-Volume-Temperature relations of the water phase in the black-oil m...
Definition: WaterPvtMultiplexer.hpp:75
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: WaterPvtMultiplexer.hpp:141
const Scalar waterReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition: WaterPvtMultiplexer.hpp:147
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtMultiplexer.hpp:177
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WaterPvtMultiplexer.hpp:164
This class implements temperature dependence of the PVT properties of water.
Definition: WaterPvtThermal.hpp:50
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtThermal.hpp:359
bool enableThermalViscosity() const
Returns true iff the viscosity of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:268
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WaterPvtThermal.hpp:338
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition: WaterPvtThermal.hpp:250
bool enableThermalDensity() const
Returns true iff the density of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:256
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the specific internal energy [J/kg] of water given a set of parameters.
Definition: WaterPvtThermal.hpp:278
bool enableJouleThomsony() const
Returns true iff Joule-Thomson effect for the water phase is active.
Definition: WaterPvtThermal.hpp:262
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: WaterPvtThermal.hpp:230