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 
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, bool enableBrine>
45 class WaterPvtMultiplexer;
46 
53 template <class Scalar, bool enableBrine>
55 {
56 public:
58  typedef WaterPvtMultiplexer<Scalar, /*enableThermal=*/false, enableBrine> IsothermalPvt;
59 
61  {
62  enableThermalDensity_ = false;
63  enableThermalViscosity_ = false;
64  enableInternalEnergy_ = false;
65  isothermalPvt_ = nullptr;
66  }
67 
68  WaterPvtThermal(IsothermalPvt* isothermalPvt,
69  const std::vector<Scalar>& viscrefPress,
70  const std::vector<Scalar>& watdentRefTemp,
71  const std::vector<Scalar>& watdentCT1,
72  const std::vector<Scalar>& watdentCT2,
73  const std::vector<Scalar>& pvtwRefPress,
74  const std::vector<Scalar>& pvtwRefB,
75  const std::vector<Scalar>& pvtwCompressibility,
76  const std::vector<Scalar>& pvtwViscosity,
77  const std::vector<Scalar>& pvtwViscosibility,
78  const std::vector<TabulatedOneDFunction>& watvisctCurves,
79  const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
82  bool enableInternalEnergy)
83  : isothermalPvt_(isothermalPvt)
84  , viscrefPress_(viscrefPress)
85  , watdentRefTemp_(watdentRefTemp)
86  , watdentCT1_(watdentCT1)
87  , watdentCT2_(watdentCT2)
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  , enableThermalViscosity_(enableThermalViscosity)
97  , enableInternalEnergy_(enableInternalEnergy)
98  { }
99 
100  WaterPvtThermal(const WaterPvtThermal& data)
101  { *this = data; }
102 
103  ~WaterPvtThermal()
104  { delete isothermalPvt_; }
105 
106 #if HAVE_ECL_INPUT
110  void initFromState(const EclipseState& eclState, const Schedule& schedule)
111  {
113  // initialize the isothermal part
115  isothermalPvt_ = new IsothermalPvt;
116  isothermalPvt_->initFromState(eclState, schedule);
117 
119  // initialize the thermal part
121  const auto& tables = eclState.getTableManager();
122 
123  enableThermalDensity_ = tables.WatDenT().size() > 0;
124  enableThermalViscosity_ = tables.hasTables("WATVISCT");
125  enableInternalEnergy_ = tables.hasTables("SPECHEAT");
126 
127  unsigned numRegions = isothermalPvt_->numRegions();
128  setNumRegions(numRegions);
129 
130  if (enableThermalDensity_) {
131  const auto& watDenT = tables.WatDenT();
132 
133  assert(watDenT.size() == numRegions);
134  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
135  const auto& record = watDenT[regionIdx];
136 
137  watdentRefTemp_[regionIdx] = record.T0;
138  watdentCT1_[regionIdx] = record.C1;
139  watdentCT2_[regionIdx] = record.C2;
140  }
141 
142  const auto& pvtwTables = tables.getPvtwTable();
143 
144  assert(pvtwTables.size() == numRegions);
145  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
146  pvtwRefPress_[regionIdx] = pvtwTables[regionIdx].reference_pressure;
147  pvtwRefB_[regionIdx] = pvtwTables[regionIdx].volume_factor;
148  }
149  }
150 
151  if (enableThermalViscosity_) {
152  if (tables.getViscrefTable().empty())
153  throw std::runtime_error("VISCREF is required when WATVISCT is present");
154 
155  const auto& watvisctTables = tables.getWatvisctTables();
156  const auto& viscrefTables = tables.getViscrefTable();
157 
158  const auto& pvtwTables = tables.getPvtwTable();
159 
160  assert(pvtwTables.size() == numRegions);
161  assert(watvisctTables.size() == numRegions);
162  assert(viscrefTables.size() == numRegions);
163 
164  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
165  const auto& T = watvisctTables[regionIdx].getColumn("Temperature").vectorCopy();
166  const auto& mu = watvisctTables[regionIdx].getColumn("Viscosity").vectorCopy();
167  watvisctCurves_[regionIdx].setXYContainers(T, mu);
168 
169  viscrefPress_[regionIdx] = viscrefTables[regionIdx].reference_pressure;
170  }
171 
172  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
173  pvtwViscosity_[regionIdx] = pvtwTables[regionIdx].viscosity;
174  pvtwViscosibility_[regionIdx] = pvtwTables[regionIdx].viscosibility;
175  }
176  }
177 
178  if (enableInternalEnergy_) {
179  // the specific internal energy of liquid water. be aware that ecl only specifies the heat capacity
180  // (via the SPECHEAT keyword) and we need to integrate it ourselfs to get the
181  // internal energy
182  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
183  const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
184  const auto& temperatureColumn = specHeatTable.getColumn("TEMPERATURE");
185  const auto& cvWaterColumn = specHeatTable.getColumn("CV_WATER");
186 
187  std::vector<double> uSamples(temperatureColumn.size());
188 
189  Scalar u = temperatureColumn[0]*cvWaterColumn[0];
190  for (size_t i = 0;; ++i) {
191  uSamples[i] = u;
192 
193  if (i >= temperatureColumn.size() - 1)
194  break;
195 
196  // integrate to the heat capacity from the current sampling point to the next
197  // one. this leads to a quadratic polynomial.
198  Scalar c_v0 = cvWaterColumn[i];
199  Scalar c_v1 = cvWaterColumn[i + 1];
200  Scalar T0 = temperatureColumn[i];
201  Scalar T1 = temperatureColumn[i + 1];
202  u += 0.5*(c_v0 + c_v1)*(T1 - T0);
203  }
204 
205  internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
206  }
207  }
208  }
209 #endif // HAVE_ECL_INPUT
210 
214  void setNumRegions(size_t numRegions)
215  {
216  pvtwRefPress_.resize(numRegions);
217  pvtwRefB_.resize(numRegions);
218  pvtwCompressibility_.resize(numRegions);
219  pvtwViscosity_.resize(numRegions);
220  pvtwViscosibility_.resize(numRegions);
221  viscrefPress_.resize(numRegions);
222  watvisctCurves_.resize(numRegions);
223  watdentRefTemp_.resize(numRegions);
224  watdentCT1_.resize(numRegions);
225  watdentCT2_.resize(numRegions);
226  internalEnergyCurves_.resize(numRegions);
227  }
228 
232  void initEnd()
233  { }
234 
238  bool enableThermalDensity() const
239  { return enableThermalDensity_; }
240 
245  { return enableThermalViscosity_; }
246 
247  size_t numRegions() const
248  { return pvtwRefPress_.size(); }
249 
253  template <class Evaluation>
254  Evaluation internalEnergy(unsigned regionIdx,
255  const Evaluation& temperature,
256  const Evaluation&) const
257  {
258  if (!enableInternalEnergy_)
259  throw std::runtime_error("Requested the internal energy of oil but it is disabled");
260 
261  // compute the specific internal energy for the specified tempature. We use linear
262  // interpolation here despite the fact that the underlying heat capacities are
263  // piecewise linear (which leads to a quadratic function)
264  return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
265  }
266 
270  template <class Evaluation>
271  Evaluation viscosity(unsigned regionIdx,
272  const Evaluation& temperature,
273  const Evaluation& pressure,
274  const Evaluation& saltconcentration) const
275  {
276  const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, saltconcentration);
277  if (!enableThermalViscosity())
278  return isothermalMu;
279 
280  Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
281  Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
282 
283  // compute the viscosity deviation due to temperature
284  const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
285  return isothermalMu * muWatvisct/muRef;
286  }
287 
291  template <class Evaluation>
292  Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
293  const Evaluation& temperature,
294  const Evaluation& pressure,
295  const Evaluation& saltconcentration) const
296  {
297  if (!enableThermalDensity())
298  return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration);
299 
300  Scalar BwRef = pvtwRefB_[regionIdx];
301  Scalar TRef = watdentRefTemp_[regionIdx];
302  const Evaluation& X = pvtwCompressibility_[regionIdx]*(pressure - pvtwRefPress_[regionIdx]);
303  Scalar cT1 = watdentCT1_[regionIdx];
304  Scalar cT2 = watdentCT2_[regionIdx];
305  const Evaluation& Y = temperature - TRef;
306 
307  // this is inconsistent with the density calculation of water in the isothermal
308  // case (it misses the quadratic pressure term), but it is the equation given in
309  // the documentation.
310  return 1.0/(((1 - X)*(1 + cT1*Y + cT2*Y*Y))*BwRef);
311  }
312 
313  const IsothermalPvt* isoThermalPvt() const
314  { return isothermalPvt_; }
315 
316  const Scalar waterReferenceDensity(unsigned regionIdx) const
317  { return isothermalPvt_->waterReferenceDensity(regionIdx); }
318 
319  const std::vector<Scalar>& viscrefPress() const
320  { return viscrefPress_; }
321 
322  const std::vector<Scalar>& watdentRefTemp() const
323  { return watdentRefTemp_; }
324 
325  const std::vector<Scalar>& watdentCT1() const
326  { return watdentCT1_; }
327 
328  const std::vector<Scalar>& watdentCT2() const
329  { return watdentCT2_; }
330 
331  const std::vector<Scalar>& pvtwRefPress() const
332  { return pvtwRefPress_; }
333 
334  const std::vector<Scalar>& pvtwRefB() const
335  { return pvtwRefB_; }
336 
337  const std::vector<Scalar>& pvtwCompressibility() const
338  { return pvtwCompressibility_; }
339 
340  const std::vector<Scalar>& pvtwViscosity() const
341  { return pvtwViscosity_; }
342 
343  const std::vector<Scalar>& pvtwViscosibility() const
344  { return pvtwViscosibility_; }
345 
346  const std::vector<TabulatedOneDFunction>& watvisctCurves() const
347  { return watvisctCurves_; }
348 
349  const std::vector<TabulatedOneDFunction> internalEnergyCurves() const
350  { return internalEnergyCurves_; }
351 
352  bool enableInternalEnergy() const
353  { return enableInternalEnergy_; }
354 
355  bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
356  {
357  if (isothermalPvt_ && !data.isothermalPvt_)
358  return false;
359  if (!isothermalPvt_ && data.isothermalPvt_)
360  return false;
361 
362  return (!this->isoThermalPvt() ||
363  (*this->isoThermalPvt() == *data.isoThermalPvt())) &&
364  this->viscrefPress() == data.viscrefPress() &&
365  this->watdentRefTemp() == data.watdentRefTemp() &&
366  this->watdentCT1() == data.watdentCT1() &&
367  this->watdentCT2() == data.watdentCT2() &&
368  this->pvtwRefPress() == data.pvtwRefPress() &&
369  this->pvtwRefB() == data.pvtwRefB() &&
370  this->pvtwCompressibility() == data.pvtwCompressibility() &&
371  this->pvtwViscosity() == data.pvtwViscosity() &&
372  this->pvtwViscosibility() == data.pvtwViscosibility() &&
373  this->watvisctCurves() == data.watvisctCurves() &&
374  this->internalEnergyCurves() == data.internalEnergyCurves() &&
375  this->enableThermalDensity() == data.enableThermalDensity() &&
376  this->enableThermalViscosity() == data.enableThermalViscosity() &&
377  this->enableInternalEnergy() == data.enableInternalEnergy();
378  }
379 
380  WaterPvtThermal<Scalar, enableBrine>& operator=(const WaterPvtThermal<Scalar, enableBrine>& data)
381  {
382  if (data.isothermalPvt_)
383  isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
384  else
385  isothermalPvt_ = nullptr;
386  viscrefPress_ = data.viscrefPress_;
387  watdentRefTemp_ = data.watdentRefTemp_;
388  watdentCT1_ = data.watdentCT1_;
389  watdentCT2_ = data.watdentCT2_;
390  pvtwRefPress_ = data.pvtwRefPress_;
391  pvtwRefB_ = data.pvtwRefB_;
392  pvtwCompressibility_ = data.pvtwCompressibility_;
393  pvtwViscosity_ = data.pvtwViscosity_;
394  pvtwViscosibility_ = data.pvtwViscosibility_;
395  watvisctCurves_ = data.watvisctCurves_;
396  internalEnergyCurves_ = data.internalEnergyCurves_;
397  enableThermalDensity_ = data.enableThermalDensity_;
398  enableThermalViscosity_ = data.enableThermalViscosity_;
399  enableInternalEnergy_ = data.enableInternalEnergy_;
400 
401  return *this;
402  }
403 
404 private:
405  IsothermalPvt* isothermalPvt_;
406 
407  // The PVT properties needed for temperature dependence. We need to store one
408  // value per PVT region.
409  std::vector<Scalar> viscrefPress_;
410 
411  std::vector<Scalar> watdentRefTemp_;
412  std::vector<Scalar> watdentCT1_;
413  std::vector<Scalar> watdentCT2_;
414 
415  std::vector<Scalar> pvtwRefPress_;
416  std::vector<Scalar> pvtwRefB_;
417  std::vector<Scalar> pvtwCompressibility_;
418  std::vector<Scalar> pvtwViscosity_;
419  std::vector<Scalar> pvtwViscosibility_;
420 
421  std::vector<TabulatedOneDFunction> watvisctCurves_;
422 
423  // piecewise linear curve representing the internal energy of water
424  std::vector<TabulatedOneDFunction> internalEnergyCurves_;
425 
426  bool enableThermalDensity_;
427  bool enableThermalViscosity_;
428  bool enableInternalEnergy_;
429 };
430 
431 } // namespace Opm
432 
433 #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...
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
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:176
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:163
This class implements temperature dependence of the PVT properties of water.
Definition: WaterPvtThermal.hpp:55
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:292
bool enableThermalViscosity() const
Returns true iff the viscosity of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:244
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:271
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &) const
Returns the specific internal energy [J/kg] of water given a set of parameters.
Definition: WaterPvtThermal.hpp:254
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition: WaterPvtThermal.hpp:232
bool enableThermalDensity() const
Returns true iff the density of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:238
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: WaterPvtThermal.hpp:214