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
32#include <cstddef>
33
34namespace Opm {
35
36#if HAVE_ECL_INPUT
37class EclipseState;
38class Schedule;
39#endif
40
41template <class Scalar, bool enableThermal, bool enableBrine>
42class WaterPvtMultiplexer;
43
50template <class Scalar, bool enableBrine>
52{
53public:
55 using IsothermalPvt = WaterPvtMultiplexer<Scalar, /*enableThermal=*/false, enableBrine>;
56
57 WaterPvtThermal() = default;
58
59 WaterPvtThermal(IsothermalPvt* isothermalPvt,
60 const std::vector<Scalar>& viscrefPress,
61 const std::vector<Scalar>& watdentRefTemp,
62 const std::vector<Scalar>& watdentCT1,
63 const std::vector<Scalar>& watdentCT2,
64 const std::vector<Scalar>& watJTRefPres,
65 const std::vector<Scalar>& watJTC,
66 const std::vector<Scalar>& pvtwRefPress,
67 const std::vector<Scalar>& pvtwRefB,
68 const std::vector<Scalar>& pvtwCompressibility,
69 const std::vector<Scalar>& pvtwViscosity,
70 const std::vector<Scalar>& pvtwViscosibility,
71 const std::vector<TabulatedOneDFunction>& watvisctCurves,
72 const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
76 bool enableInternalEnergy)
77 : isothermalPvt_(isothermalPvt)
78 , viscrefPress_(viscrefPress)
79 , watdentRefTemp_(watdentRefTemp)
80 , watdentCT1_(watdentCT1)
81 , watdentCT2_(watdentCT2)
82 , watJTRefPres_(watJTRefPres)
83 , watJTC_(watJTC)
84 , pvtwRefPress_(pvtwRefPress)
85 , pvtwRefB_(pvtwRefB)
86 , pvtwCompressibility_(pvtwCompressibility)
87 , pvtwViscosity_(pvtwViscosity)
88 , pvtwViscosibility_(pvtwViscosibility)
89 , watvisctCurves_(watvisctCurves)
90 , internalEnergyCurves_(internalEnergyCurves)
91 , enableThermalDensity_(enableThermalDensity)
92 , enableJouleThomson_(enableJouleThomson)
93 , enableThermalViscosity_(enableThermalViscosity)
94 , enableInternalEnergy_(enableInternalEnergy)
95 { }
96
98 { *this = data; }
99
101 { delete isothermalPvt_; }
102
103#if HAVE_ECL_INPUT
107 void initFromState(const EclipseState& eclState, const Schedule& schedule);
108#endif // HAVE_ECL_INPUT
109
113 void setNumRegions(std::size_t numRegions);
114
115 void setVapPars(const Scalar par1, const Scalar par2)
116 {
117 isothermalPvt_->setVapPars(par1, par2);
118 }
119
123 void initEnd()
124 { }
125
130 { return enableThermalDensity_; }
131
136 { return enableJouleThomson_; }
137
142 { return enableThermalViscosity_; }
143
144 Scalar hVap(unsigned regionIdx) const
145 { return this->hVap_[regionIdx]; }
146
147 std::size_t numRegions() const
148 { return pvtwRefPress_.size(); }
149
153 template <class Evaluation>
154 Evaluation internalEnergy(unsigned regionIdx,
155 const Evaluation& temperature,
156 const Evaluation& pressure,
157 const Evaluation& Rsw,
158 const Evaluation& saltconcentration) const
159 {
160 if (!enableInternalEnergy_) {
161 throw std::runtime_error("Requested the internal energy of water but it is disabled");
162 }
163
164 if (!enableJouleThomson_) {
165 // compute the specific internal energy for the specified tempature. We use linear
166 // interpolation here despite the fact that the underlying heat capacities are
167 // piecewise linear (which leads to a quadratic function)
168 return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
169 }
170 else {
171 Evaluation Tref = watdentRefTemp_[regionIdx];
172 Evaluation Pref = watJTRefPres_[regionIdx];
173 Scalar JTC = watJTC_[regionIdx]; // if JTC is default then JTC is calculated
174
175 Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rsw, saltconcentration);
176 Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
177 Evaluation density = invB * waterReferenceDensity(regionIdx);
178
179 Evaluation enthalpyPres;
180 if (JTC != 0) {
181 enthalpyPres = -Cp * JTC * (pressure - Pref);
182 }
183 else if (enableThermalDensity_) {
184 Scalar c1T = watdentCT1_[regionIdx];
185 Scalar c2T = watdentCT2_[regionIdx];
186
187 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
188 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
189
190 constexpr const int N = 100; // value is experimental
191 Evaluation deltaP = (pressure - Pref) / N;
192 Evaluation enthalpyPresPrev = 0;
193 for (std::size_t i = 0; i < N; ++i) {
194 Evaluation Pnew = Pref + i * deltaP;
195 Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature,
196 Pnew, Rsw, saltconcentration) *
197 waterReferenceDensity(regionIdx);
198 Evaluation jouleThomsonCoefficient = -(1.0 / Cp) * (1.0 - alpha * temperature) / rho;
199 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
200 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
201 enthalpyPresPrev = enthalpyPres;
202 }
203 }
204 else {
205 throw std::runtime_error("Requested Joule-thomson calculation but "
206 "thermal water density (WATDENT) is not provided");
207 }
208
209 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
210
211 return enthalpy - pressure/density;
212 }
213 }
214
218 template <class Evaluation>
219 Evaluation viscosity(unsigned regionIdx,
220 const Evaluation& temperature,
221 const Evaluation& pressure,
222 const Evaluation& Rsw,
223 const Evaluation& saltconcentration) const
224 {
225 const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature,
226 pressure, Rsw, saltconcentration);
227 if (!enableThermalViscosity()) {
228 return isothermalMu;
229 }
230
231 Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] -
232 pvtwRefPress_[regionIdx]);
233 Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);
234
235 // compute the viscosity deviation due to temperature
236 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
237 return isothermalMu * muWatvisct / muRef;
238 }
239
243 template <class Evaluation>
244 Evaluation saturatedViscosity(unsigned regionIdx,
245 const Evaluation& temperature,
246 const Evaluation& pressure,
247 const Evaluation& saltconcentration) const
248 {
249 const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature,
250 pressure, saltconcentration);
251 if (!enableThermalViscosity()) {
252 return isothermalMu;
253 }
254
255 Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] -
256 pvtwRefPress_[regionIdx]);
257 Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);
258
259 // compute the viscosity deviation due to temperature
260 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
261 return isothermalMu * muWatvisct / muRef;
262 }
263
267 template <class Evaluation>
268 Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
269 const Evaluation& temperature,
270 const Evaluation& pressure,
271 const Evaluation& saltconcentration) const
272 {
273 Evaluation Rsw = 0.0;
274 return inverseFormationVolumeFactor(regionIdx, temperature, pressure,
275 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 {
287 if (!enableThermalDensity()) {
288 return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature,
289 pressure, Rsw, saltconcentration);
290 }
291
292 Scalar BwRef = pvtwRefB_[regionIdx];
293 Scalar TRef = watdentRefTemp_[regionIdx];
294 const Evaluation& X = pvtwCompressibility_[regionIdx] * (pressure -
295 pvtwRefPress_[regionIdx]);
296 Scalar cT1 = watdentCT1_[regionIdx];
297 Scalar cT2 = watdentCT2_[regionIdx];
298 const Evaluation& Y = temperature - TRef;
299
300 // this is inconsistent with the density calculation of water in the isothermal
301 // case (it misses the quadratic pressure term), but it is the equation given in
302 // the documentation.
303 return 1.0 / (((1 - X) * (1 + cT1 * Y + cT2 * Y * Y)) * BwRef);
304 }
305
313 template <class Evaluation>
314 Evaluation saturationPressure(unsigned /*regionIdx*/,
315 const Evaluation& /*temperature*/,
316 const Evaluation& /*Rs*/,
317 const Evaluation& /*saltconcentration*/) const
318 { return 0.0; /* this is dead water, so there isn't any meaningful saturation pressure! */ }
319
323 template <class Evaluation>
324 Evaluation saturatedGasDissolutionFactor(unsigned /*regionIdx*/,
325 const Evaluation& /*temperature*/,
326 const Evaluation& /*pressure*/,
327 const Evaluation& /*saltconcentration*/) const
328 { return 0.0; /* this is dead water! */ }
329
330 template <class Evaluation>
331 Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
332 const Evaluation& /*pressure*/,
333 unsigned /*compIdx*/) const
334 {
335 throw std::runtime_error("Not implemented: The PVT model does not provide "
336 "a diffusionCoefficient()");
337 }
338
339 const IsothermalPvt* isoThermalPvt() const
340 { return isothermalPvt_; }
341
342 Scalar waterReferenceDensity(unsigned regionIdx) const
343 { return isothermalPvt_->waterReferenceDensity(regionIdx); }
344
345 const std::vector<Scalar>& viscrefPress() const
346 { return viscrefPress_; }
347
348 const std::vector<Scalar>& watdentRefTemp() const
349 { return watdentRefTemp_; }
350
351 const std::vector<Scalar>& watdentCT1() const
352 { return watdentCT1_; }
353
354 const std::vector<Scalar>& watdentCT2() const
355 { return watdentCT2_; }
356
357 const std::vector<Scalar>& pvtwRefPress() const
358 { return pvtwRefPress_; }
359
360 const std::vector<Scalar>& pvtwRefB() const
361 { return pvtwRefB_; }
362
363 const std::vector<Scalar>& pvtwCompressibility() const
364 { return pvtwCompressibility_; }
365
366 const std::vector<Scalar>& pvtwViscosity() const
367 { return pvtwViscosity_; }
368
369 const std::vector<Scalar>& pvtwViscosibility() const
370 { return pvtwViscosibility_; }
371
372 const std::vector<TabulatedOneDFunction>& watvisctCurves() const
373 { return watvisctCurves_; }
374
375 const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
376 { return internalEnergyCurves_; }
377
378 bool enableInternalEnergy() const
379 { return enableInternalEnergy_; }
380
381 const std::vector<Scalar>& watJTRefPres() const
382 { return watJTRefPres_; }
383
384 const std::vector<Scalar>& watJTC() const
385 { return watJTC_; }
386
387 bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const;
388
389 WaterPvtThermal<Scalar, enableBrine>&
390 operator=(const WaterPvtThermal<Scalar, enableBrine>& data);
391
392private:
393 IsothermalPvt* isothermalPvt_{nullptr};
394
395 // The PVT properties needed for temperature dependence. We need to store one
396 // value per PVT region.
397 std::vector<Scalar> viscrefPress_{};
398
399 std::vector<Scalar> watdentRefTemp_{};
400 std::vector<Scalar> watdentCT1_{};
401 std::vector<Scalar> watdentCT2_{};
402
403 std::vector<Scalar> watJTRefPres_{};
404 std::vector<Scalar> watJTC_{};
405
406 std::vector<Scalar> pvtwRefPress_{};
407 std::vector<Scalar> pvtwRefB_{};
408 std::vector<Scalar> pvtwCompressibility_{};
409 std::vector<Scalar> pvtwViscosity_{};
410 std::vector<Scalar> pvtwViscosibility_{};
411
412 std::vector<TabulatedOneDFunction> watvisctCurves_{};
413
414 // piecewise linear curve representing the internal energy of water
415 std::vector<TabulatedOneDFunction> internalEnergyCurves_{};
416 std::vector<Scalar> hVap_{};
417
418 bool enableThermalDensity_{false};
419 bool enableJouleThomson_{false};
420 bool enableThermalViscosity_{false};
421 bool enableInternalEnergy_{false};
422};
423
424} // namespace Opm
425
426#endif
Implements a linearly interpolated scalar function that depends on one variable.
Definition EclipseState.hpp:63
Definition Schedule.hpp:101
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:90
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: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 WaterPvtMultiplexer.hpp:168
Scalar waterReferenceDensity(unsigned regionIdx) const
Return the reference density which are considered by this PVT-object.
Definition WaterPvtMultiplexer.cpp:111
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:180
This class implements temperature dependence of the PVT properties of water.
Definition WaterPvtThermal.hpp:52
bool enableThermalViscosity() const
Returns true iff the viscosity of the water phase is temperature dependent.
Definition WaterPvtThermal.hpp:141
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:314
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:219
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition WaterPvtThermal.hpp:123
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:129
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:324
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:154
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:268
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the water phase is active.
Definition WaterPvtThermal.hpp:135
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:244
void setNumRegions(std::size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition WaterPvtThermal.cpp:185
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30