All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Pages
ThermalOilPvtWrapper.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_OIL_PVT_WRAPPER_HPP
21 #define OPM_THERMAL_OIL_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/OilvisctTable.hpp>
28 
29 #include <vector>
30 
31 namespace Opm
32 {
35  class ThermalOilPvtWrapper : 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 
49  int numRegions;
50  auto tables = eclipseState->getTableManager();
51 
52  if (deck->hasKeyword("PVTO"))
53  numRegions = tables->getPvtoTables().size();
54  else if (deck->hasKeyword("PVDO"))
55  numRegions = tables->getPvdoTables().size();
56  else if (deck->hasKeyword("PVCDO"))
57  numRegions = deck->getKeyword("PVCDO").size();
58  else
59  OPM_THROW(std::runtime_error, "Oil phase was not initialized using a known way");
60 
61  // viscosity
62  if (deck->hasKeyword("VISCREF")) {
63  oilvisctTables_ = &tables->getOilvisctTables();
64  const auto& viscrefKeyword = deck->getKeyword("VISCREF");
65 
66  assert(int(oilvisctTables_->size()) == numRegions);
67  assert(int(viscrefKeyword.size()) == numRegions);
68 
69  viscrefPress_.resize(numRegions);
70  viscrefRs_.resize(numRegions);
71  muRef_.resize(numRegions);
72  for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
73  const auto& viscrefRecord = viscrefKeyword.getRecord(regionIdx);
74  viscrefPress_[regionIdx] = viscrefRecord.getItem("REFERENCE_PRESSURE").getSIDouble(0);
75  viscrefRs_[regionIdx] = viscrefRecord.getItem("REFERENCE_RS").getSIDouble(0);
76 
77  // temperature used to calculate the reference viscosity [K]. the
78  // value does not really matter if the underlying PVT object really
79  // is isothermal...
80  double Tref = 273.15 + 20;
81 
82  // compute the reference viscosity using the isothermal PVT object.
83  double tmp1, tmp2;
84  isothermalPvt_->mu(1,
85  &regionIdx,
86  &viscrefPress_[regionIdx],
87  &Tref,
88  &viscrefRs_[regionIdx],
89  &muRef_[regionIdx],
90  &tmp1,
91  &tmp2);
92  }
93  }
94 
95  // quantities required for density. note that we just always use the values
96  // for the first EOS. (since EOS != PVT region.)
97  tref_ = 0.0;
98  if (deck->hasKeyword("THERMEX1")) {
99  oilCompIdx_ = deck->getKeyword("OCOMPIDX").getRecord(0).getItem("OIL_COMPONENT_INDEX").get< int >(0) - 1;
100 
101  // always use the values of the first EOS
102  tref_ = deck->getKeyword("TREF").getRecord(0).getItem("TEMPERATURE").getSIDouble(oilCompIdx_);
103  pref_ = deck->getKeyword("PREF").getRecord(0).getItem("PRESSURE").getSIDouble(oilCompIdx_);
104  cref_ = deck->getKeyword("CREF").getRecord(0).getItem("COMPRESSIBILITY").getSIDouble(oilCompIdx_);
105  thermex1_ = deck->getKeyword("THERMEX1").getRecord(0).getItem("EXPANSION_COEFF").getSIDouble(oilCompIdx_);
106  }
107  }
108 
109  virtual void mu(const int n,
110  const int* pvtRegionIdx,
111  const double* p,
112  const double* T,
113  const double* z,
114  double* output_mu) const
115  {
116  if (oilvisctTables_)
117  // TODO: temperature dependence for viscosity depending on z
118  OPM_THROW(std::runtime_error,
119  "temperature dependent viscosity as a function of z "
120  "is not yet implemented!");
121 
122  // compute the isothermal viscosity
123  isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
124  }
125 
126  virtual void mu(const int n,
127  const int* pvtRegionIdx,
128  const double* p,
129  const double* T,
130  const double* r,
131  double* output_mu,
132  double* output_dmudp,
133  double* output_dmudr) const
134  {
135  // compute the isothermal viscosity and its derivatives
136  isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
137 
138  if (!oilvisctTables_)
139  // isothermal case
140  return;
141 
142  // temperature dependence
143  for (int i = 0; i < n; ++i) {
144  int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
145 
146  // calculate the viscosity of the isothermal keyword for the reference
147  // pressure given by the VISCREF keyword.
148  double muRef = muRef_[regionIdx];
149 
150  // compute the viscosity deviation due to temperature
151  double alpha;
152  {
153  const OilvisctTable& oilvisctTable = oilvisctTables_->getTable<OilvisctTable>(regionIdx);
154  double muOilvisct = oilvisctTable.evaluate("Viscosity", T[i]);
155  alpha = muOilvisct/muRef;
156  }
157 
158  output_mu[i] *= alpha;
159  output_dmudp[i] *= alpha;
160  output_dmudr[i] *= alpha;
161  // TODO (?): derivative of viscosity w.r.t. temperature.
162  }
163  }
164 
165  virtual void mu(const int n,
166  const int* pvtRegionIdx,
167  const double* p,
168  const double* T,
169  const double* r,
170  const PhasePresence* cond,
171  double* output_mu,
172  double* output_dmudp,
173  double* output_dmudr) const
174  {
175  // compute the isothermal viscosity and its derivatives
176  isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
177 
178  if (!oilvisctTables_)
179  // isothermal case
180  return;
181 
182  // temperature dependence
183  for (int i = 0; i < n; ++i) {
184  int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
185 
186  // calculate the viscosity of the isothermal keyword for the reference
187  // pressure given by the VISCREF keyword.
188  double muRef = muRef_[regionIdx];
189 
190  // compute the viscosity deviation due to temperature
191  double alpha;
192  {
193  const OilvisctTable& oilvisctTable = oilvisctTables_->getTable<OilvisctTable>(regionIdx);
194  double muOilvisct = oilvisctTable.evaluate("Viscosity", T[i]);
195  alpha = muOilvisct/muRef;
196  }
197  output_mu[i] *= alpha;
198  output_dmudp[i] *= alpha;
199  output_dmudr[i] *= alpha;
200  // TODO (?): derivative of viscosity w.r.t. temperature.
201  }
202  }
203 
204  virtual void B(const int n,
205  const int* pvtRegionIdx,
206  const double* p,
207  const double* T,
208  const double* z,
209  double* output_B) const
210  {
211  // isothermal case
212  isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
213 
214  if (thermex1_ <= 0.0)
215  // isothermal case
216  return;
217 
218  // deal with the temperature dependence of the oil phase. we use equation
219  // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
220  // using the isothermal keyword instead of using the value for the
221  // components, so the oil compressibility is already dealt with there. Note
222  // that we only do the part for the oil component here, the part for
223  // dissolved gas is ignored so far.
224  double cT1 = thermex1_;
225  double TRef = tref_;
226  for (int i = 0; i < n; ++i) {
227  double alpha = (1 + cT1*(T[i] - TRef));
228  output_B[i] *= alpha;
229  }
230  }
231 
232  virtual void dBdp(const int n,
233  const int* pvtRegionIdx,
234  const double* p,
235  const double* T,
236  const double* z,
237  double* output_B,
238  double* output_dBdp) const
239  {
240  isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
241 
242  if (thermex1_ <= 0.0)
243  // isothermal case
244  return;
245 
246  // deal with the temperature dependence of the oil phase. we use equation
247  // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
248  // using the isothermal keyword instead of using the value for the
249  // components, so the oil compressibility is already dealt with there. Note
250  // that we only do the part for the oil component here, the part for
251  // dissolved gas is ignored so far.
252  double cT1 = thermex1_;
253  double TRef = tref_;
254  for (int i = 0; i < n; ++i) {
255  double alpha = (1 + cT1*(T[i] - TRef));
256  output_B[i] *= alpha;
257  output_dBdp[i] *= alpha;
258  }
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  isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
271 
272  if (thermex1_ <= 0.0)
273  // isothermal case
274  return;
275 
276  // deal with the temperature dependence of the oil phase. we use equation
277  // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
278  // using the isothermal keyword instead of using the value for the
279  // components, so the oil compressibility is already dealt with there. Note
280  // that we only do the part for the oil component here, the part for
281  // dissolved gas is ignored so far.
282  double cT1 = thermex1_;
283  double TRef = tref_;
284  for (int i = 0; i < n; ++i) {
285  double alpha = 1.0/(1 + cT1*(T[i] - TRef));
286  output_b[i] *= alpha;
287  output_dbdp[i] *= alpha;
288  output_dbdr[i] *= alpha;
289  }
290  }
291 
292  virtual void b(const int n,
293  const int* pvtRegionIdx,
294  const double* p,
295  const double* T,
296  const double* r,
297  const PhasePresence* cond,
298  double* output_b,
299  double* output_dbdp,
300  double* output_dbdr) const
301  {
302  isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
303 
304  if (thermex1_ <= 0.0)
305  // isothermal case
306  return;
307 
308  // deal with the temperature dependence of the oil phase. we use equation
309  // (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
310  // using the isothermal keyword instead of using the value for the
311  // components, so the oil compressibility is already dealt with there. Note
312  // that we only do the part for the oil component here, the part for
313  // dissolved gas is ignored so far.
314  double cT1 = thermex1_;
315  double TRef = tref_;
316  for (int i = 0; i < n; ++i) {
317  double alpha = 1.0/(1 + cT1*(T[i] - TRef));
318  output_b[i] *= alpha;
319  output_dbdp[i] *= alpha;
320  output_dbdr[i] *= alpha;
321  }
322  }
323 
324  virtual void rsSat(const int n,
325  const int* pvtRegionIdx,
326  const double* p,
327  double* output_rsSat,
328  double* output_drsSatdp) const
329  {
330  isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
331  }
332 
333  virtual void rvSat(const int n,
334  const int* pvtRegionIdx,
335  const double* p,
336  double* output_rvSat,
337  double* output_drvSatdp) const
338  {
339  isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
340  }
341 
342  virtual void R(const int n,
343  const int* pvtRegionIdx,
344  const double* p,
345  const double* z,
346  double* output_R) const
347  {
348  isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
349  }
350 
351  virtual void dRdp(const int n,
352  const int* pvtRegionIdx,
353  const double* p,
354  const double* z,
355  double* output_R,
356  double* output_dRdp) const
357  {
358  isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
359  }
360 
361  private:
362  int getPvtRegionIndex_(const int* pvtRegionIdx, int cellIdx) const
363  {
364  if (!pvtRegionIdx)
365  return 0;
366  return pvtRegionIdx[cellIdx];
367  }
368 
369  // the PVT propertied for the isothermal case
370  std::shared_ptr<const PvtInterface> isothermalPvt_;
371 
372  // The PVT properties needed for temperature dependence of the viscosity. We need
373  // to store one value per PVT region.
374  std::vector<double> viscrefPress_;
375  std::vector<double> viscrefRs_;
376  std::vector<double> muRef_;
377 
378  const TableContainer* oilvisctTables_;
379 
380  // The PVT properties needed for temperature dependence of the density. This is
381  // specified as one value per EOS in the manual, but we unconditionally use the
382  // expansion coefficient of the first EOS...
383  int oilCompIdx_;
384  double tref_;
385  double pref_;
386  double cref_;
387  double thermex1_;
388  };
389 
390 }
391 
392 #endif
393 
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 oil viscosity
Definition: ThermalOilPvtWrapper.hpp:43
Class which wraps another (i.e., isothermal) PVT object into one which adds temperature dependence of...
Definition: ThermalOilPvtWrapper.hpp:35