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
32namespace Opm {
33
34#if HAVE_ECL_INPUT
35class EclipseState;
36class Schedule;
37#endif
38
39template <class Scalar, bool enableThermal, bool enableBrine>
40class WaterPvtMultiplexer;
41
48template <class Scalar, bool enableBrine>
50{
51public:
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,
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
102 { *this = data; }
103
105 { delete isothermalPvt_; }
106
107#if HAVE_ECL_INPUT
111 void initFromState(const EclipseState& eclState, const Schedule& schedule);
112#endif // HAVE_ECL_INPUT
113
117 void setNumRegions(size_t numRegions)
118 {
119 pvtwRefPress_.resize(numRegions);
120 pvtwRefB_.resize(numRegions);
121 pvtwCompressibility_.resize(numRegions);
122 pvtwViscosity_.resize(numRegions);
123 pvtwViscosibility_.resize(numRegions);
124 viscrefPress_.resize(numRegions);
125 watvisctCurves_.resize(numRegions);
126 watdentRefTemp_.resize(numRegions);
127 watdentCT1_.resize(numRegions);
128 watdentCT2_.resize(numRegions);
129 watJTRefPres_.resize(numRegions);
130 watJTC_.resize(numRegions);
131 internalEnergyCurves_.resize(numRegions);
132 }
133
137 void initEnd()
138 { }
139
144 { return enableThermalDensity_; }
145
150 { return enableJouleThomson_; }
151
156 { return enableThermalViscosity_; }
157
158 size_t numRegions() const
159 { return pvtwRefPress_.size(); }
160
164 template <class Evaluation>
165 Evaluation internalEnergy(unsigned regionIdx,
166 const Evaluation& temperature,
167 const Evaluation& pressure,
168 const Evaluation& Rsw,
169 const Evaluation& saltconcentration) const
170 {
171 if (!enableInternalEnergy_)
172 throw std::runtime_error("Requested the internal energy of water but it is disabled");
173
174 if (!enableJouleThomson_) {
175 // compute the specific internal energy for the specified tempature. We use linear
176 // interpolation here despite the fact that the underlying heat capacities are
177 // piecewise linear (which leads to a quadratic function)
178 return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
179 }
180 else {
181 Evaluation Tref = watdentRefTemp_[regionIdx];
182 Evaluation Pref = watJTRefPres_[regionIdx];
183 Scalar JTC =watJTC_[regionIdx]; // if JTC is default then JTC is calculated
184
185 Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
186 Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
187 Evaluation density = invB * waterReferenceDensity(regionIdx);
188
189 Evaluation enthalpyPres;
190 if (JTC != 0) {
191 enthalpyPres = -Cp * JTC * (pressure -Pref);
192 }
193 else if(enableThermalDensity_) {
194 Scalar c1T = watdentCT1_[regionIdx];
195 Scalar c2T = watdentCT2_[regionIdx];
196
197 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
198 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
199
200 constexpr const int N = 100; // value is experimental
201 Evaluation deltaP = (pressure - Pref)/N;
202 Evaluation enthalpyPresPrev = 0;
203 for (size_t i = 0; i < N; ++i) {
204 Evaluation Pnew = Pref + i * deltaP;
205 Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rsw, saltconcentration) * waterReferenceDensity(regionIdx);
206 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
207 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
208 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
209 enthalpyPresPrev = enthalpyPres;
210 }
211 }
212 else {
213 throw std::runtime_error("Requested Joule-thomson calculation but thermal water density (WATDENT) is not provided");
214 }
215
216 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
217
218 return enthalpy - pressure/density;
219 }
220 }
221
225 template <class Evaluation>
226 Evaluation viscosity(unsigned regionIdx,
227 const Evaluation& temperature,
228 const Evaluation& pressure,
229 const Evaluation& Rsw,
230 const Evaluation& saltconcentration) const
231 {
232 const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rsw, saltconcentration);
234 return isothermalMu;
235
236 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
237 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
238
239 // compute the viscosity deviation due to temperature
240 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
241 return isothermalMu * muWatvisct/muRef;
242 }
243
247 template <class Evaluation>
248 Evaluation saturatedViscosity(unsigned regionIdx,
249 const Evaluation& temperature,
250 const Evaluation& pressure,
251 const Evaluation& saltconcentration) const
252 {
253 const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure, saltconcentration);
255 return isothermalMu;
256
257 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
258 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
259
260 // compute the viscosity deviation due to temperature
261 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
262 return isothermalMu * muWatvisct/muRef;
263 }
264
268 template <class Evaluation>
269 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
270 const Evaluation& temperature,
271 const Evaluation& pressure,
272 const Evaluation& saltconcentration) const
273 {
274 Evaluation Rsw = 0.0;
275 return inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
276 }
280 template <class Evaluation>
281 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
282 const Evaluation& temperature,
283 const Evaluation& pressure,
284 const Evaluation& Rsw,
285 const Evaluation& saltconcentration) const
286 {
288 return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
289
290 Scalar BwRef = pvtwRefB_[regionIdx];
291 Scalar TRef = watdentRefTemp_[regionIdx];
292 const Evaluation& X = pvtwCompressibility_[regionIdx]*(pressure - pvtwRefPress_[regionIdx]);
293 Scalar cT1 = watdentCT1_[regionIdx];
294 Scalar cT2 = watdentCT2_[regionIdx];
295 const Evaluation& Y = temperature - TRef;
296
297 // this is inconsistent with the density calculation of water in the isothermal
298 // case (it misses the quadratic pressure term), but it is the equation given in
299 // the documentation.
300 return 1.0/(((1 - X)*(1 + cT1*Y + cT2*Y*Y))*BwRef);
301 }
302
309 template <class Evaluation>
310 Evaluation saturationPressure(unsigned /*regionIdx*/,
311 const Evaluation& /*temperature*/,
312 const Evaluation& /*Rs*/,
313 const Evaluation& /*saltconcentration*/) const
314 { return 0.0; /* this is dead water, so there isn't any meaningful saturation pressure! */ }
315
319 template <class Evaluation>
320 Evaluation saturatedGasDissolutionFactor(unsigned /*regionIdx*/,
321 const Evaluation& /*temperature*/,
322 const Evaluation& /*pressure*/,
323 const Evaluation& /*saltconcentration*/) const
324 { return 0.0; /* this is dead water! */ }
325
326 template <class Evaluation>
327 Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
328 const Evaluation& /*pressure*/,
329 unsigned /*compIdx*/) const
330 {
331 throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
332 }
333
334 const IsothermalPvt* isoThermalPvt() const
335 { return isothermalPvt_; }
336
337 const Scalar waterReferenceDensity(unsigned regionIdx) const
338 { return isothermalPvt_->waterReferenceDensity(regionIdx); }
339
340 const std::vector<Scalar>& viscrefPress() const
341 { return viscrefPress_; }
342
343 const std::vector<Scalar>& watdentRefTemp() const
344 { return watdentRefTemp_; }
345
346 const std::vector<Scalar>& watdentCT1() const
347 { return watdentCT1_; }
348
349 const std::vector<Scalar>& watdentCT2() const
350 { return watdentCT2_; }
351
352 const std::vector<Scalar>& pvtwRefPress() const
353 { return pvtwRefPress_; }
354
355 const std::vector<Scalar>& pvtwRefB() const
356 { return pvtwRefB_; }
357
358 const std::vector<Scalar>& pvtwCompressibility() const
359 { return pvtwCompressibility_; }
360
361 const std::vector<Scalar>& pvtwViscosity() const
362 { return pvtwViscosity_; }
363
364 const std::vector<Scalar>& pvtwViscosibility() const
365 { return pvtwViscosibility_; }
366
367 const std::vector<TabulatedOneDFunction>& watvisctCurves() const
368 { return watvisctCurves_; }
369
370 const std::vector<TabulatedOneDFunction> internalEnergyCurves() const
371 { return internalEnergyCurves_; }
372
373 bool enableInternalEnergy() const
374 { return enableInternalEnergy_; }
375
376 const std::vector<Scalar>& watJTRefPres() const
377 { return watJTRefPres_; }
378
379 const std::vector<Scalar>& watJTC() const
380 { return watJTC_; }
381
382
383 bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
384 {
385 if (isothermalPvt_ && !data.isothermalPvt_)
386 return false;
387 if (!isothermalPvt_ && data.isothermalPvt_)
388 return false;
389
390 return this->viscrefPress() == data.viscrefPress() &&
391 this->watdentRefTemp() == data.watdentRefTemp() &&
392 this->watdentCT1() == data.watdentCT1() &&
393 this->watdentCT2() == data.watdentCT2() &&
394 this->watdentCT2() == data.watdentCT2() &&
395 this->watJTRefPres() == data.watJTRefPres() &&
396 this->watJTC() == data.watJTC() &&
397 this->pvtwRefPress() == data.pvtwRefPress() &&
398 this->pvtwRefB() == data.pvtwRefB() &&
399 this->pvtwCompressibility() == data.pvtwCompressibility() &&
400 this->pvtwViscosity() == data.pvtwViscosity() &&
401 this->pvtwViscosibility() == data.pvtwViscosibility() &&
402 this->watvisctCurves() == data.watvisctCurves() &&
403 this->internalEnergyCurves() == data.internalEnergyCurves() &&
404 this->enableThermalDensity() == data.enableThermalDensity() &&
405 this->enableJouleThomson() == data.enableJouleThomson() &&
406 this->enableThermalViscosity() == data.enableThermalViscosity() &&
407 this->enableInternalEnergy() == data.enableInternalEnergy();
408 }
409
410 WaterPvtThermal<Scalar, enableBrine>& operator=(const WaterPvtThermal<Scalar, enableBrine>& data)
411 {
412 if (data.isothermalPvt_)
413 isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
414 else
415 isothermalPvt_ = nullptr;
416 viscrefPress_ = data.viscrefPress_;
417 watdentRefTemp_ = data.watdentRefTemp_;
418 watdentCT1_ = data.watdentCT1_;
419 watdentCT2_ = data.watdentCT2_;
420 watJTRefPres_ = data.watJTRefPres_;
421 watJTC_ = data.watJTC_;
422 pvtwRefPress_ = data.pvtwRefPress_;
423 pvtwRefB_ = data.pvtwRefB_;
424 pvtwCompressibility_ = data.pvtwCompressibility_;
425 pvtwViscosity_ = data.pvtwViscosity_;
426 pvtwViscosibility_ = data.pvtwViscosibility_;
427 watvisctCurves_ = data.watvisctCurves_;
428 internalEnergyCurves_ = data.internalEnergyCurves_;
429 enableThermalDensity_ = data.enableThermalDensity_;
430 enableJouleThomson_ = data.enableJouleThomson_;
431 enableThermalViscosity_ = data.enableThermalViscosity_;
432 enableInternalEnergy_ = data.enableInternalEnergy_;
433
434 return *this;
435 }
436
437private:
438 IsothermalPvt* isothermalPvt_;
439
440 // The PVT properties needed for temperature dependence. We need to store one
441 // value per PVT region.
442 std::vector<Scalar> viscrefPress_;
443
444 std::vector<Scalar> watdentRefTemp_;
445 std::vector<Scalar> watdentCT1_;
446 std::vector<Scalar> watdentCT2_;
447
448 std::vector<Scalar> watJTRefPres_;
449 std::vector<Scalar> watJTC_;
450
451 std::vector<Scalar> pvtwRefPress_;
452 std::vector<Scalar> pvtwRefB_;
453 std::vector<Scalar> pvtwCompressibility_;
454 std::vector<Scalar> pvtwViscosity_;
455 std::vector<Scalar> pvtwViscosibility_;
456
457 std::vector<TabulatedOneDFunction> watvisctCurves_;
458
459 // piecewise linear curve representing the internal energy of water
460 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
461
462 bool enableThermalDensity_;
463 bool enableJouleThomson_;
464 bool enableThermalViscosity_;
465 bool enableInternalEnergy_;
466};
467
468} // namespace Opm
469
470#endif
Definition: EclipseState.hpp:55
Definition: Schedule.hpp:130
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:51
This class represents the Pressure-Volume-Temperature relations of the water phase in the black-oil m...
Definition: WaterPvtMultiplexer.hpp:82
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WaterPvtMultiplexer.hpp:163
Evaluation saturatedViscosity(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:177
const Scalar waterReferenceDensity(unsigned regionIdx) const
Return the reference density which are considered by this PVT-object.
Definition: WaterPvtMultiplexer.hpp:145
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtMultiplexer.hpp:190
This class implements temperature dependence of the PVT properties of water.
Definition: WaterPvtThermal.hpp:50
bool enableThermalViscosity() const
Returns true iff the viscosity of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:155
Evaluation saturationPressure(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the saturation pressure of the water phase [Pa] depending on its mass fraction of the gas com...
Definition: WaterPvtThermal.hpp:310
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition: WaterPvtThermal.hpp:226
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition: WaterPvtThermal.hpp:137
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtThermal.hpp:281
bool enableThermalDensity() const
Returns true iff the density of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:143
Evaluation saturatedGasDissolutionFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the gas dissolution factor [m^3/m^3] of the water phase.
Definition: WaterPvtThermal.hpp:320
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the specific internal energy [J/kg] of water given a set of parameters.
Definition: WaterPvtThermal.hpp:165
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtThermal.hpp:269
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the water phase is active.
Definition: WaterPvtThermal.hpp:149
Evaluation saturatedViscosity(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:248
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: WaterPvtThermal.hpp:117
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition: Exceptions.hpp:30