My Project
Loading...
Searching...
No Matches
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 enableJouleThomson_ = false;
59 enableThermalViscosity_ = false;
60 enableInternalEnergy_ = false;
61 isothermalPvt_ = nullptr;
62 }
63
64 WaterPvtThermal(IsothermalPvt* isothermalPvt,
65 const std::vector<Scalar>& viscrefPress,
66 const std::vector<Scalar>& watdentRefTemp,
67 const std::vector<Scalar>& watdentCT1,
68 const std::vector<Scalar>& watdentCT2,
69 const std::vector<Scalar>& watJTRefPres,
70 const std::vector<Scalar>& watJTC,
71 const std::vector<Scalar>& pvtwRefPress,
72 const std::vector<Scalar>& pvtwRefB,
73 const std::vector<Scalar>& pvtwCompressibility,
74 const std::vector<Scalar>& pvtwViscosity,
75 const std::vector<Scalar>& pvtwViscosibility,
76 const std::vector<TabulatedOneDFunction>& watvisctCurves,
77 const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
81 bool enableInternalEnergy)
82 : isothermalPvt_(isothermalPvt)
83 , viscrefPress_(viscrefPress)
84 , watdentRefTemp_(watdentRefTemp)
85 , watdentCT1_(watdentCT1)
86 , watdentCT2_(watdentCT2)
87 , watJTRefPres_(watJTRefPres)
88 , watJTC_(watJTC)
89 , pvtwRefPress_(pvtwRefPress)
90 , pvtwRefB_(pvtwRefB)
91 , pvtwCompressibility_(pvtwCompressibility)
92 , pvtwViscosity_(pvtwViscosity)
93 , pvtwViscosibility_(pvtwViscosibility)
94 , watvisctCurves_(watvisctCurves)
95 , internalEnergyCurves_(internalEnergyCurves)
96 , enableThermalDensity_(enableThermalDensity)
97 , enableJouleThomson_(enableJouleThomson)
98 , enableThermalViscosity_(enableThermalViscosity)
99 , enableInternalEnergy_(enableInternalEnergy)
100 { }
101
103 { *this = data; }
104
106 { delete isothermalPvt_; }
107
108#if HAVE_ECL_INPUT
112 void initFromState(const EclipseState& eclState, const Schedule& schedule);
113#endif // HAVE_ECL_INPUT
114
118 void setNumRegions(size_t numRegions)
119 {
120 pvtwRefPress_.resize(numRegions);
121 pvtwRefB_.resize(numRegions);
122 pvtwCompressibility_.resize(numRegions);
123 pvtwViscosity_.resize(numRegions);
124 pvtwViscosibility_.resize(numRegions);
125 viscrefPress_.resize(numRegions);
126 watvisctCurves_.resize(numRegions);
127 watdentRefTemp_.resize(numRegions);
128 watdentCT1_.resize(numRegions);
129 watdentCT2_.resize(numRegions);
130 watJTRefPres_.resize(numRegions);
131 watJTC_.resize(numRegions);
132 internalEnergyCurves_.resize(numRegions);
133 }
134
135 void setVapPars(const Scalar par1, const Scalar par2)
136 {
137 isothermalPvt_->setVapPars(par1, par2);
138 }
139
143 void initEnd()
144 { }
145
150 { return enableThermalDensity_; }
151
156 { return enableJouleThomson_; }
157
162 { return enableThermalViscosity_; }
163
164 size_t numRegions() const
165 { return pvtwRefPress_.size(); }
166
170 template <class Evaluation>
171 Evaluation internalEnergy(unsigned regionIdx,
172 const Evaluation& temperature,
173 const Evaluation& pressure,
174 const Evaluation& Rsw,
175 const Evaluation& saltconcentration) const
176 {
177 if (!enableInternalEnergy_)
178 throw std::runtime_error("Requested the internal energy of water but it is disabled");
179
180 if (!enableJouleThomson_) {
181 // compute the specific internal energy for the specified tempature. We use linear
182 // interpolation here despite the fact that the underlying heat capacities are
183 // piecewise linear (which leads to a quadratic function)
184 return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
185 }
186 else {
187 Evaluation Tref = watdentRefTemp_[regionIdx];
188 Evaluation Pref = watJTRefPres_[regionIdx];
189 Scalar JTC =watJTC_[regionIdx]; // if JTC is default then JTC is calculated
190
191 Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
192 Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
193 Evaluation density = invB * waterReferenceDensity(regionIdx);
194
195 Evaluation enthalpyPres;
196 if (JTC != 0) {
197 enthalpyPres = -Cp * JTC * (pressure -Pref);
198 }
199 else if(enableThermalDensity_) {
200 Scalar c1T = watdentCT1_[regionIdx];
201 Scalar c2T = watdentCT2_[regionIdx];
202
203 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
204 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
205
206 constexpr const int N = 100; // value is experimental
207 Evaluation deltaP = (pressure - Pref)/N;
208 Evaluation enthalpyPresPrev = 0;
209 for (size_t i = 0; i < N; ++i) {
210 Evaluation Pnew = Pref + i * deltaP;
211 Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rsw, saltconcentration) * waterReferenceDensity(regionIdx);
212 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
213 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
214 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
215 enthalpyPresPrev = enthalpyPres;
216 }
217 }
218 else {
219 throw std::runtime_error("Requested Joule-thomson calculation but thermal water density (WATDENT) is not provided");
220 }
221
222 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
223
224 return enthalpy - pressure/density;
225 }
226 }
227
231 template <class Evaluation>
232 Evaluation viscosity(unsigned regionIdx,
233 const Evaluation& temperature,
234 const Evaluation& pressure,
235 const Evaluation& Rsw,
236 const Evaluation& saltconcentration) const
237 {
238 const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rsw, saltconcentration);
240 return isothermalMu;
241
242 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
243 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
244
245 // compute the viscosity deviation due to temperature
246 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
247 return isothermalMu * muWatvisct/muRef;
248 }
249
253 template <class Evaluation>
254 Evaluation saturatedViscosity(unsigned regionIdx,
255 const Evaluation& temperature,
256 const Evaluation& pressure,
257 const Evaluation& saltconcentration) const
258 {
259 const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure, saltconcentration);
261 return isothermalMu;
262
263 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
264 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
265
266 // compute the viscosity deviation due to temperature
267 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
268 return isothermalMu * muWatvisct/muRef;
269 }
270
274 template <class Evaluation>
275 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
276 const Evaluation& temperature,
277 const Evaluation& pressure,
278 const Evaluation& saltconcentration) const
279 {
280 Evaluation Rsw = 0.0;
281 return inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
282 }
286 template <class Evaluation>
287 Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
288 const Evaluation& temperature,
289 const Evaluation& pressure,
290 const Evaluation& Rsw,
291 const Evaluation& saltconcentration) const
292 {
294 return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
295
296 Scalar BwRef = pvtwRefB_[regionIdx];
297 Scalar TRef = watdentRefTemp_[regionIdx];
298 const Evaluation& X = pvtwCompressibility_[regionIdx]*(pressure - pvtwRefPress_[regionIdx]);
299 Scalar cT1 = watdentCT1_[regionIdx];
300 Scalar cT2 = watdentCT2_[regionIdx];
301 const Evaluation& Y = temperature - TRef;
302
303 // this is inconsistent with the density calculation of water in the isothermal
304 // case (it misses the quadratic pressure term), but it is the equation given in
305 // the documentation.
306 return 1.0/(((1 - X)*(1 + cT1*Y + cT2*Y*Y))*BwRef);
307 }
308
315 template <class Evaluation>
316 Evaluation saturationPressure(unsigned /*regionIdx*/,
317 const Evaluation& /*temperature*/,
318 const Evaluation& /*Rs*/,
319 const Evaluation& /*saltconcentration*/) const
320 { return 0.0; /* this is dead water, so there isn't any meaningful saturation pressure! */ }
321
325 template <class Evaluation>
326 Evaluation saturatedGasDissolutionFactor(unsigned /*regionIdx*/,
327 const Evaluation& /*temperature*/,
328 const Evaluation& /*pressure*/,
329 const Evaluation& /*saltconcentration*/) const
330 { return 0.0; /* this is dead water! */ }
331
332 template <class Evaluation>
333 Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
334 const Evaluation& /*pressure*/,
335 unsigned /*compIdx*/) const
336 {
337 throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
338 }
339
340 const IsothermalPvt* isoThermalPvt() const
341 { return isothermalPvt_; }
342
343 const Scalar waterReferenceDensity(unsigned regionIdx) const
344 { return isothermalPvt_->waterReferenceDensity(regionIdx); }
345
346 const std::vector<Scalar>& viscrefPress() const
347 { return viscrefPress_; }
348
349 const std::vector<Scalar>& watdentRefTemp() const
350 { return watdentRefTemp_; }
351
352 const std::vector<Scalar>& watdentCT1() const
353 { return watdentCT1_; }
354
355 const std::vector<Scalar>& watdentCT2() const
356 { return watdentCT2_; }
357
358 const std::vector<Scalar>& pvtwRefPress() const
359 { return pvtwRefPress_; }
360
361 const std::vector<Scalar>& pvtwRefB() const
362 { return pvtwRefB_; }
363
364 const std::vector<Scalar>& pvtwCompressibility() const
365 { return pvtwCompressibility_; }
366
367 const std::vector<Scalar>& pvtwViscosity() const
368 { return pvtwViscosity_; }
369
370 const std::vector<Scalar>& pvtwViscosibility() const
371 { return pvtwViscosibility_; }
372
373 const std::vector<TabulatedOneDFunction>& watvisctCurves() const
374 { return watvisctCurves_; }
375
376 const std::vector<TabulatedOneDFunction> internalEnergyCurves() const
377 { return internalEnergyCurves_; }
378
379 bool enableInternalEnergy() const
380 { return enableInternalEnergy_; }
381
382 const std::vector<Scalar>& watJTRefPres() const
383 { return watJTRefPres_; }
384
385 const std::vector<Scalar>& watJTC() const
386 { return watJTC_; }
387
388
389 bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
390 {
391 if (isothermalPvt_ && !data.isothermalPvt_)
392 return false;
393 if (!isothermalPvt_ && data.isothermalPvt_)
394 return false;
395
396 return this->viscrefPress() == data.viscrefPress() &&
397 this->watdentRefTemp() == data.watdentRefTemp() &&
398 this->watdentCT1() == data.watdentCT1() &&
399 this->watdentCT2() == data.watdentCT2() &&
400 this->watdentCT2() == data.watdentCT2() &&
401 this->watJTRefPres() == data.watJTRefPres() &&
402 this->watJTC() == data.watJTC() &&
403 this->pvtwRefPress() == data.pvtwRefPress() &&
404 this->pvtwRefB() == data.pvtwRefB() &&
405 this->pvtwCompressibility() == data.pvtwCompressibility() &&
406 this->pvtwViscosity() == data.pvtwViscosity() &&
407 this->pvtwViscosibility() == data.pvtwViscosibility() &&
408 this->watvisctCurves() == data.watvisctCurves() &&
409 this->internalEnergyCurves() == data.internalEnergyCurves() &&
410 this->enableThermalDensity() == data.enableThermalDensity() &&
411 this->enableJouleThomson() == data.enableJouleThomson() &&
412 this->enableThermalViscosity() == data.enableThermalViscosity() &&
413 this->enableInternalEnergy() == data.enableInternalEnergy();
414 }
415
416 WaterPvtThermal<Scalar, enableBrine>& operator=(const WaterPvtThermal<Scalar, enableBrine>& data)
417 {
418 if (data.isothermalPvt_)
419 isothermalPvt_ = new IsothermalPvt(*data.isothermalPvt_);
420 else
421 isothermalPvt_ = nullptr;
422 viscrefPress_ = data.viscrefPress_;
423 watdentRefTemp_ = data.watdentRefTemp_;
424 watdentCT1_ = data.watdentCT1_;
425 watdentCT2_ = data.watdentCT2_;
426 watJTRefPres_ = data.watJTRefPres_;
427 watJTC_ = data.watJTC_;
428 pvtwRefPress_ = data.pvtwRefPress_;
429 pvtwRefB_ = data.pvtwRefB_;
430 pvtwCompressibility_ = data.pvtwCompressibility_;
431 pvtwViscosity_ = data.pvtwViscosity_;
432 pvtwViscosibility_ = data.pvtwViscosibility_;
433 watvisctCurves_ = data.watvisctCurves_;
434 internalEnergyCurves_ = data.internalEnergyCurves_;
435 enableThermalDensity_ = data.enableThermalDensity_;
436 enableJouleThomson_ = data.enableJouleThomson_;
437 enableThermalViscosity_ = data.enableThermalViscosity_;
438 enableInternalEnergy_ = data.enableInternalEnergy_;
439
440 return *this;
441 }
442
443private:
444 IsothermalPvt* isothermalPvt_;
445
446 // The PVT properties needed for temperature dependence. We need to store one
447 // value per PVT region.
448 std::vector<Scalar> viscrefPress_;
449
450 std::vector<Scalar> watdentRefTemp_;
451 std::vector<Scalar> watdentCT1_;
452 std::vector<Scalar> watdentCT2_;
453
454 std::vector<Scalar> watJTRefPres_;
455 std::vector<Scalar> watJTC_;
456
457 std::vector<Scalar> pvtwRefPress_;
458 std::vector<Scalar> pvtwRefB_;
459 std::vector<Scalar> pvtwCompressibility_;
460 std::vector<Scalar> pvtwViscosity_;
461 std::vector<Scalar> pvtwViscosibility_;
462
463 std::vector<TabulatedOneDFunction> watvisctCurves_;
464
465 // piecewise linear curve representing the internal energy of water
466 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
467
468 bool enableThermalDensity_;
469 bool enableJouleThomson_;
470 bool enableThermalViscosity_;
471 bool enableInternalEnergy_;
472};
473
474} // namespace Opm
475
476#endif
Implements a linearly interpolated scalar function that depends on one variable.
Definition EclipseState.hpp:63
Definition Schedule.hpp:88
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:89
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:179
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:193
const Scalar waterReferenceDensity(unsigned regionIdx) const
Return the reference density which are considered by this PVT-object.
Definition WaterPvtMultiplexer.hpp:161
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:206
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:161
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:316
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:232
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition WaterPvtThermal.hpp:143
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:287
bool enableThermalDensity() const
Returns true iff the density of the water phase is temperature dependent.
Definition WaterPvtThermal.hpp:149
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:326
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:171
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:275
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the water phase is active.
Definition WaterPvtThermal.hpp:155
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:254
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition WaterPvtThermal.hpp:118
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30