All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Pages
ThermalWaterPvtWrapper.hpp
1 /*
2  Copyright 2015 Andreas Lauser
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 3 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 
20 #ifndef OPM_THERMAL_WATER_PVT_WRAPPER_HPP
21 #define OPM_THERMAL_WATER_PVT_WRAPPER_HPP
22 
23 #include <opm/core/props/pvt/PvtInterface.hpp>
24 #include <opm/common/ErrorMacros.hpp>
25 
26 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
27 #include <opm/parser/eclipse/EclipseState/Tables/WatvisctTable.hpp>
28 
29 #include <vector>
30 
31 namespace Opm
32 {
35  class ThermalWaterPvtWrapper : public PvtInterface
36  {
37  public:
39  {}
40 
41 
43  void initFromDeck(std::shared_ptr<const PvtInterface> isothermalPvt,
44  const Opm::Deck& deck,
45  const Opm::EclipseState& eclipseState)
46  {
47  isothermalPvt_ = isothermalPvt;
48  watvisctTables_ = 0;
49 
50  // stuff which we need to get from the PVTW keyword
51  const auto& pvtwKeyword = deck->getKeyword("PVTW");
52  int numRegions = pvtwKeyword.size();
53  pvtwRefPress_.resize(numRegions);
54  pvtwRefB_.resize(numRegions);
55  pvtwCompressibility_.resize(numRegions);
56  pvtwViscosity_.resize(numRegions);
57  pvtwViscosibility_.resize(numRegions);
58  for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
59  const auto& pvtwRecord = pvtwKeyword.getRecord(regionIdx);
60  pvtwRefPress_[regionIdx] = pvtwRecord.getItem("P_REF").getSIDouble(0);
61  pvtwRefB_[regionIdx] = pvtwRecord.getItem("WATER_VOL_FACTOR").getSIDouble(0);
62  pvtwViscosity_[regionIdx] = pvtwRecord.getItem("WATER_VISCOSITY").getSIDouble(0);
63  pvtwViscosibility_[regionIdx] = pvtwRecord.getItem("WATER_VISCOSIBILITY").getSIDouble(0);
64  }
65 
66  // quantities required for the temperature dependence of the viscosity
67  // (basically we expect well-behaved VISCREF and WATVISCT keywords.)
68  if (deck->hasKeyword("VISCREF")) {
69  auto tables = eclipseState->getTableManager();
70  watvisctTables_ = &tables->getWatvisctTables();
71  const auto& viscrefKeyword = deck->getKeyword("VISCREF");
72 
73  assert(int(watvisctTables_->size()) == numRegions);
74  assert(int(viscrefKeyword.size()) == numRegions);
75 
76  viscrefPress_.resize(numRegions);
77  for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
78  const auto& viscrefRecord = viscrefKeyword.getRecord(regionIdx);
79 
80  viscrefPress_[regionIdx] = viscrefRecord.getItem("REFERENCE_PRESSURE").getSIDouble(0);
81  }
82  }
83 
84  // quantities required for the temperature dependence of the density
85  if (deck->hasKeyword("WATDENT")) {
86  const auto& watdentKeyword = deck->getKeyword("WATDENT");
87 
88  assert(int(watdentKeyword.size()) == numRegions);
89 
90  watdentRefTemp_.resize(numRegions);
91  watdentCT1_.resize(numRegions);
92  watdentCT2_.resize(numRegions);
93  for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
94  const auto& watdentRecord = watdentKeyword.getRecord(regionIdx);
95 
96  watdentRefTemp_[regionIdx] = watdentRecord.getItem("REFERENCE_TEMPERATURE").getSIDouble(0);
97  watdentCT1_[regionIdx] = watdentRecord.getItem("EXPANSION_COEFF_LINEAR").getSIDouble(0);
98  watdentCT2_[regionIdx] = watdentRecord.getItem("EXPANSION_COEFF_QUADRATIC").getSIDouble(0);
99  }
100  }
101  }
102 
103  virtual void mu(const int n,
104  const int* pvtRegionIdx,
105  const double* p,
106  const double* T,
107  const double* z,
108  double* output_mu) const
109  {
110  if (watvisctTables_)
111  // TODO: temperature dependence for viscosity depending on z
112  OPM_THROW(std::runtime_error,
113  "temperature dependent viscosity as a function of z "
114  "is not yet implemented!");
115 
116  // compute the isothermal viscosity
117  isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
118  }
119 
120  virtual void mu(const int n,
121  const int* pvtRegionIdx,
122  const double* p,
123  const double* T,
124  const double* r,
125  double* output_mu,
126  double* output_dmudp,
127  double* output_dmudr) const
128  {
129  // compute the isothermal viscosity and its derivatives
130  isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
131 
132  if (!watvisctTables_)
133  // isothermal case
134  return;
135 
136  // temperature dependence
137  for (int i = 0; i < n; ++i) {
138  int tableIdx = getTableIndex_(pvtRegionIdx, i);
139 
140  // calculate the viscosity of the isothermal keyword for the reference
141  // pressure given by the VISCREF keyword.
142  double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]);
143  double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
144 
145  // compute the viscosity deviation due to temperature
146  double alpha;
147  {
148  const WatvisctTable& watVisctTable = watvisctTables_->getTable<WatvisctTable>(tableIdx);
149  double muWatvisct = watVisctTable.evaluate("Viscosity", T[i]);
150  alpha = muWatvisct/muRef;
151  }
152 
153  output_mu[i] *= alpha;
154  output_dmudp[i] *= alpha;
155  output_dmudr[i] *= alpha;
156  // TODO (?): derivative of viscosity w.r.t. temperature.
157  }
158  }
159 
160  virtual void mu(const int n,
161  const int* pvtRegionIdx,
162  const double* p,
163  const double* T,
164  const double* r,
165  const PhasePresence* cond,
166  double* output_mu,
167  double* output_dmudp,
168  double* output_dmudr) const
169  {
170  // compute the isothermal viscosity and its derivatives
171  isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
172 
173  if (!watvisctTables_)
174  // isothermal case
175  return;
176 
177  // temperature dependence
178  for (int i = 0; i < n; ++i) {
179  int tableIdx = getTableIndex_(pvtRegionIdx, i);
180 
181  // calculate the viscosity of the isothermal keyword for the reference
182  // pressure given by the VISCREF keyword.
183  double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]);
184  double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
185 
186  // compute the viscosity deviation due to temperature
187  double alpha;
188  {
189  const WatvisctTable& watVisctTable = watvisctTables_->getTable<WatvisctTable>(tableIdx);
190  double muWatvisct = watVisctTable.evaluate("Viscosity", T[i]);
191  alpha = muWatvisct/muRef;
192  }
193  output_mu[i] *= alpha;
194  output_dmudp[i] *= alpha;
195  output_dmudr[i] *= alpha;
196  // TODO (?): derivative of viscosity w.r.t. temperature.
197  }
198  }
199 
200  virtual void B(const int n,
201  const int* pvtRegionIdx,
202  const double* p,
203  const double* T,
204  const double* z,
205  double* output_B) const
206  {
207  if (watdentRefTemp_.empty()) {
208  // isothermal case
209  isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
210  return;
211  }
212 
213  // This changes how the water density depends on pressure compared to what's
214  // used for the PVTW keyword, but it seems to be what Eclipse does. For
215  // details, see the documentation for the WATDENT keyword in the Eclipse RM.
216  for (int i = 0; i < n; ++i) {
217  int tableIdx = getTableIndex_(pvtRegionIdx, i);
218  double BwRef = pvtwRefB_[tableIdx];
219  double TRef = watdentRefTemp_[tableIdx];
220  double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
221  double cT1 = watdentCT1_[tableIdx];
222  double cT2 = watdentCT2_[tableIdx];
223  double Y = T[i] - TRef;
224  double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
225  output_B[i] = Bw;
226  }
227  }
228 
229  virtual void dBdp(const int n,
230  const int* pvtRegionIdx,
231  const double* p,
232  const double* T,
233  const double* z,
234  double* output_B,
235  double* output_dBdp) const
236  {
237  if (watdentRefTemp_.empty()) {
238  // isothermal case
239  isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
240  return;
241  }
242 
243  // This changes how the water density depends on pressure. This is awkward,
244  // but it seems to be what Eclipse does. See the documentation for the
245  // WATDENT keyword in the Eclipse RM
246  for (int i = 0; i < n; ++i) {
247  int tableIdx = getTableIndex_(pvtRegionIdx, i);
248  double BwRef = pvtwRefB_[tableIdx];
249  double TRef = watdentRefTemp_[tableIdx];
250  double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
251  double cT1 = watdentCT1_[tableIdx];
252  double cT2 = watdentCT2_[tableIdx];
253  double Y = T[i] - TRef;
254  double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
255  output_B[i] = Bw;
256  }
257 
258  std::fill(output_dBdp, output_dBdp + n, 0.0);
259  }
260 
261  virtual void b(const int n,
262  const int* pvtRegionIdx,
263  const double* p,
264  const double* T,
265  const double* r,
266  double* output_b,
267  double* output_dbdp,
268  double* output_dbdr) const
269  {
270  if (watdentRefTemp_.empty()) {
271  // isothermal case
272  isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
273  return;
274  }
275 
276  // This changes how the water density depends on pressure. This is awkward,
277  // but it seems to be what Eclipse does. See the documentation for the
278  // WATDENT keyword in the Eclipse RM
279  for (int i = 0; i < n; ++i) {
280  int tableIdx = getTableIndex_(pvtRegionIdx, i);
281  double BwRef = pvtwRefB_[tableIdx];
282  double TRef = watdentRefTemp_[tableIdx];
283  double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
284  double cT1 = watdentCT1_[tableIdx];
285  double cT2 = watdentCT2_[tableIdx];
286  double Y = T[i] - TRef;
287  double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
288  output_b[i] = 1.0/Bw;
289  }
290 
291  std::fill(output_dbdp, output_dbdp + n, 0.0);
292  std::fill(output_dbdr, output_dbdr + n, 0.0);
293  }
294 
295  virtual void b(const int n,
296  const int* pvtRegionIdx,
297  const double* p,
298  const double* T,
299  const double* r,
300  const PhasePresence* cond,
301  double* output_b,
302  double* output_dbdp,
303  double* output_dbdr) const
304  {
305  if (watdentRefTemp_.empty()) {
306  // isothermal case
307  isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
308  return;
309  }
310 
311  // This changes pressure dependence of the water density, but it seems to be
312  // what Eclipse does. See the documentation for the WATDENT keyword in the
313  // Eclipse RM
314  for (int i = 0; i < n; ++i) {
315  int tableIdx = getTableIndex_(pvtRegionIdx, i);
316  double BwRef = pvtwRefB_[tableIdx];
317  double TRef = watdentRefTemp_[tableIdx];
318  double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
319  double cT1 = watdentCT1_[tableIdx];
320  double cT2 = watdentCT2_[tableIdx];
321  double Y = T[i] - TRef;
322  double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
323  output_b[i] = 1.0/Bw;
324  }
325 
326  std::fill(output_dbdp, output_dbdp + n, 0.0);
327  std::fill(output_dbdr, output_dbdr + n, 0.0);
328 
329  }
330 
331  virtual void rsSat(const int n,
332  const int* pvtRegionIdx,
333  const double* p,
334  double* output_rsSat,
335  double* output_drsSatdp) const
336  {
337  isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
338  }
339 
340  virtual void rvSat(const int n,
341  const int* pvtRegionIdx,
342  const double* p,
343  double* output_rvSat,
344  double* output_drvSatdp) const
345  {
346  isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
347  }
348 
349  virtual void R(const int n,
350  const int* pvtRegionIdx,
351  const double* p,
352  const double* z,
353  double* output_R) const
354  {
355  isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
356  }
357 
358  virtual void dRdp(const int n,
359  const int* pvtRegionIdx,
360  const double* p,
361  const double* z,
362  double* output_R,
363  double* output_dRdp) const
364  {
365  isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
366  }
367 
368  private:
369  int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
370  {
371  if (!pvtTableIdx)
372  return 0;
373  return pvtTableIdx[cellIdx];
374  }
375 
376  // the PVT propertied for the isothermal case
377  std::shared_ptr<const PvtInterface> isothermalPvt_;
378 
379  // The PVT properties needed for temperature dependence. We need to store one
380  // value per PVT region.
381  std::vector<double> viscrefPress_;
382 
383  std::vector<double> watdentRefTemp_;
384  std::vector<double> watdentCT1_;
385  std::vector<double> watdentCT2_;
386 
387  std::vector<double> pvtwRefPress_;
388  std::vector<double> pvtwRefB_;
389  std::vector<double> pvtwCompressibility_;
390  std::vector<double> pvtwViscosity_;
391  std::vector<double> pvtwViscosibility_;
392 
393  const TableContainer* watvisctTables_;
394  };
395 
396 }
397 
398 #endif
399 
Class which wraps another (i.e., isothermal) PVT object into one which adds temperature dependence of...
Definition: ThermalWaterPvtWrapper.hpp:35
void initFromDeck(std::shared_ptr< const PvtInterface > isothermalPvt, const Opm::Deck &deck, const Opm::EclipseState &eclipseState)
set the tables which specify the temperature dependence of the water viscosity
Definition: ThermalWaterPvtWrapper.hpp:43