ThermalWaterPvtWrapper.hpp
Go to the documentation of this file.
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
31namespace Opm
32{
35 class ThermalWaterPvtWrapper : public PvtInterface
36 {
37 public:
38 ThermalWaterPvtWrapper()
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
Opm::PhasePresence
Definition: BlackoilPhases.hpp:59
Opm::ThermalWaterPvtWrapper
Definition: ThermalWaterPvtWrapper.hpp:36
Opm::ThermalWaterPvtWrapper::initFromDeck
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
Opm::ThermalWaterPvtWrapper::dRdp
virtual void dRdp(const int n, const int *pvtRegionIdx, const double *p, const double *z, double *output_R, double *output_dRdp) const
Definition: ThermalWaterPvtWrapper.hpp:358
Opm::ThermalWaterPvtWrapper::b
virtual void b(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *r, const PhasePresence *cond, double *output_b, double *output_dbdp, double *output_dbdr) const
Definition: ThermalWaterPvtWrapper.hpp:295
Opm::ThermalWaterPvtWrapper::dBdp
virtual void dBdp(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *z, double *output_B, double *output_dBdp) const
Definition: ThermalWaterPvtWrapper.hpp:229
Opm::ThermalWaterPvtWrapper::R
virtual void R(const int n, const int *pvtRegionIdx, const double *p, const double *z, double *output_R) const
Definition: ThermalWaterPvtWrapper.hpp:349
Opm::ThermalWaterPvtWrapper::b
virtual void b(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *r, double *output_b, double *output_dbdp, double *output_dbdr) const
Definition: ThermalWaterPvtWrapper.hpp:261
Opm::ThermalWaterPvtWrapper::mu
virtual void mu(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *r, double *output_mu, double *output_dmudp, double *output_dmudr) const
Definition: ThermalWaterPvtWrapper.hpp:120
Opm::ThermalWaterPvtWrapper::mu
virtual void mu(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *r, const PhasePresence *cond, double *output_mu, double *output_dmudp, double *output_dmudr) const
Definition: ThermalWaterPvtWrapper.hpp:160
Opm::ThermalWaterPvtWrapper::rsSat
virtual void rsSat(const int n, const int *pvtRegionIdx, const double *p, double *output_rsSat, double *output_drsSatdp) const
Definition: ThermalWaterPvtWrapper.hpp:331
Opm::ThermalWaterPvtWrapper::ThermalWaterPvtWrapper
ThermalWaterPvtWrapper()
Definition: ThermalWaterPvtWrapper.hpp:38
Opm::ThermalWaterPvtWrapper::mu
virtual void mu(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *z, double *output_mu) const
Definition: ThermalWaterPvtWrapper.hpp:103
Opm::ThermalWaterPvtWrapper::rvSat
virtual void rvSat(const int n, const int *pvtRegionIdx, const double *p, double *output_rvSat, double *output_drvSatdp) const
Definition: ThermalWaterPvtWrapper.hpp:340
Opm::ThermalWaterPvtWrapper::B
virtual void B(const int n, const int *pvtRegionIdx, const double *p, const double *T, const double *z, double *output_B) const
Definition: ThermalWaterPvtWrapper.hpp:200
Opm
Definition: AnisotropicEikonal.hpp:44