27 #ifndef OPM_WATER_PVT_THERMAL_HPP
28 #define OPM_WATER_PVT_THERMAL_HPP
33 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
34 #include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
35 #include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
39 template <
class Scalar,
bool enableThermal,
bool enableBrine>
40 class WaterPvtMultiplexer;
48 template <
class Scalar,
bool enableBrine>
57 enableThermalDensity_ =
false;
58 enableThermalViscosity_ =
false;
59 enableInternalEnergy_ =
false;
60 isothermalPvt_ =
nullptr;
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,
78 bool enableJouleThomson,
80 bool enableInternalEnergy)
81 : isothermalPvt_(isothermalPvt)
82 , viscrefPress_(viscrefPress)
83 , watdentRefTemp_(watdentRefTemp)
84 , watdentCT1_(watdentCT1)
85 , watdentCT2_(watdentCT2)
86 , watJTRefPres_(watJTRefPres)
88 , pvtwRefPress_(pvtwRefPress)
90 , pvtwCompressibility_(pvtwCompressibility)
91 , pvtwViscosity_(pvtwViscosity)
92 , pvtwViscosibility_(pvtwViscosibility)
93 , watvisctCurves_(watvisctCurves)
94 , internalEnergyCurves_(internalEnergyCurves)
96 , enableJouleThomson_(enableJouleThomson)
98 , enableInternalEnergy_(enableInternalEnergy)
105 {
delete isothermalPvt_; }
111 void initFromState(
const EclipseState& eclState,
const Schedule& schedule)
117 isothermalPvt_->initFromState(eclState, schedule);
122 const auto& tables = eclState.getTableManager();
124 enableThermalDensity_ = tables.WatDenT().size() > 0;
125 enableJouleThomson_ = tables.WatJT().size() > 0;
126 enableThermalViscosity_ = tables.hasTables(
"WATVISCT");
127 enableInternalEnergy_ = tables.hasTables(
"SPECHEAT");
129 unsigned numRegions = isothermalPvt_->
numRegions();
132 if (enableThermalDensity_) {
133 const auto& watDenT = tables.WatDenT();
135 assert(watDenT.size() == numRegions);
136 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
137 const auto& record = watDenT[regionIdx];
139 watdentRefTemp_[regionIdx] = record.T0;
140 watdentCT1_[regionIdx] = record.C1;
141 watdentCT2_[regionIdx] = record.C2;
144 const auto& pvtwTables = tables.getPvtwTable();
146 assert(pvtwTables.size() == numRegions);
148 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
149 pvtwRefPress_[regionIdx] = pvtwTables[regionIdx].reference_pressure;
150 pvtwRefB_[regionIdx] = pvtwTables[regionIdx].volume_factor;
155 if (enableJouleThomson_) {
156 const auto& watJT = tables.WatJT();
158 assert(watJT.size() == numRegions);
159 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
160 const auto& record = watJT[regionIdx];
162 watJTRefPres_[regionIdx] = record.P0;
163 watJTC_[regionIdx] = record.C1;
167 if (enableThermalViscosity_) {
168 if (tables.getViscrefTable().empty())
169 throw std::runtime_error(
"VISCREF is required when WATVISCT is present");
171 const auto& watvisctTables = tables.getWatvisctTables();
172 const auto& viscrefTables = tables.getViscrefTable();
174 const auto& pvtwTables = tables.getPvtwTable();
176 assert(pvtwTables.size() == numRegions);
177 assert(watvisctTables.size() == numRegions);
178 assert(viscrefTables.size() == numRegions);
180 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
181 const auto& T = watvisctTables[regionIdx].getColumn(
"Temperature").vectorCopy();
182 const auto& mu = watvisctTables[regionIdx].getColumn(
"Viscosity").vectorCopy();
183 watvisctCurves_[regionIdx].setXYContainers(T, mu);
185 viscrefPress_[regionIdx] = viscrefTables[regionIdx].reference_pressure;
188 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
189 pvtwViscosity_[regionIdx] = pvtwTables[regionIdx].viscosity;
190 pvtwViscosibility_[regionIdx] = pvtwTables[regionIdx].viscosibility;
194 if (enableInternalEnergy_) {
198 for (
unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
199 const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
200 const auto& temperatureColumn = specHeatTable.getColumn(
"TEMPERATURE");
201 const auto& cvWaterColumn = specHeatTable.getColumn(
"CV_WATER");
203 std::vector<double> uSamples(temperatureColumn.size());
205 Scalar u = temperatureColumn[0]*cvWaterColumn[0];
206 for (
size_t i = 0;; ++i) {
209 if (i >= temperatureColumn.size() - 1)
214 Scalar c_v0 = cvWaterColumn[i];
215 Scalar c_v1 = cvWaterColumn[i + 1];
216 Scalar T0 = temperatureColumn[i];
217 Scalar T1 = temperatureColumn[i + 1];
218 u += 0.5*(c_v0 + c_v1)*(T1 - T0);
221 internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
232 pvtwRefPress_.resize(numRegions);
233 pvtwRefB_.resize(numRegions);
234 pvtwCompressibility_.resize(numRegions);
235 pvtwViscosity_.resize(numRegions);
236 pvtwViscosibility_.resize(numRegions);
237 viscrefPress_.resize(numRegions);
238 watvisctCurves_.resize(numRegions);
239 watdentRefTemp_.resize(numRegions);
240 watdentCT1_.resize(numRegions);
241 watdentCT2_.resize(numRegions);
242 watJTRefPres_.resize(numRegions);
243 watJTC_.resize(numRegions);
244 internalEnergyCurves_.resize(numRegions);
257 {
return enableThermalDensity_; }
263 {
return enableJouleThomson_; }
269 {
return enableThermalViscosity_; }
271 size_t numRegions()
const
272 {
return pvtwRefPress_.size(); }
277 template <
class Evaluation>
279 const Evaluation& temperature,
280 const Evaluation& pressure,
281 const Evaluation& saltconcentration)
const
283 if (!enableInternalEnergy_)
284 throw std::runtime_error(
"Requested the internal energy of water but it is disabled");
286 if (!enableJouleThomson_) {
290 return internalEnergyCurves_[regionIdx].eval(temperature,
true);
293 Evaluation Tref = watdentRefTemp_[regionIdx];
294 Evaluation Pref = watJTRefPres_[regionIdx];
295 Scalar JTC =watJTC_[regionIdx];
298 Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature,
true)/temperature;
299 Evaluation density = invB * waterReferenceDensity(regionIdx);
301 Evaluation enthalpyPres;
303 enthalpyPres = -Cp * JTC * (pressure -Pref);
305 else if(enableThermalDensity_) {
306 Scalar c1T = watdentCT1_[regionIdx];
307 Scalar c2T = watdentCT2_[regionIdx];
309 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
310 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
312 constexpr
const int N = 100;
313 Evaluation deltaP = (pressure - Pref)/N;
314 Evaluation enthalpyPresPrev = 0;
315 for (
size_t i = 0; i < N; ++i) {
316 Evaluation Pnew = Pref + i * deltaP;
318 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
319 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
320 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
321 enthalpyPresPrev = enthalpyPres;
325 throw std::runtime_error(
"Requested Joule-thomson calculation but thermal water density (WATDENT) is not provided");
328 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
330 return enthalpy - pressure/density;
337 template <
class Evaluation>
339 const Evaluation& temperature,
340 const Evaluation& pressure,
341 const Evaluation& saltconcentration)
const
343 const auto& isothermalMu = isothermalPvt_->
viscosity(regionIdx, temperature, pressure, saltconcentration);
347 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
348 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
351 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature,
true);
352 return isothermalMu * muWatvisct/muRef;
358 template <
class Evaluation>
360 const Evaluation& temperature,
361 const Evaluation& pressure,
362 const Evaluation& saltconcentration)
const
367 Scalar BwRef = pvtwRefB_[regionIdx];
368 Scalar TRef = watdentRefTemp_[regionIdx];
369 const Evaluation& X = pvtwCompressibility_[regionIdx]*(pressure - pvtwRefPress_[regionIdx]);
370 Scalar cT1 = watdentCT1_[regionIdx];
371 Scalar cT2 = watdentCT2_[regionIdx];
372 const Evaluation& Y = temperature - TRef;
377 return 1.0/(((1 - X)*(1 + cT1*Y + cT2*Y*Y))*BwRef);
380 const IsothermalPvt* isoThermalPvt()
const
381 {
return isothermalPvt_; }
383 const Scalar waterReferenceDensity(
unsigned regionIdx)
const
386 const std::vector<Scalar>& viscrefPress()
const
387 {
return viscrefPress_; }
389 const std::vector<Scalar>& watdentRefTemp()
const
390 {
return watdentRefTemp_; }
392 const std::vector<Scalar>& watdentCT1()
const
393 {
return watdentCT1_; }
395 const std::vector<Scalar>& watdentCT2()
const
396 {
return watdentCT2_; }
398 const std::vector<Scalar>& pvtwRefPress()
const
399 {
return pvtwRefPress_; }
401 const std::vector<Scalar>& pvtwRefB()
const
402 {
return pvtwRefB_; }
404 const std::vector<Scalar>& pvtwCompressibility()
const
405 {
return pvtwCompressibility_; }
407 const std::vector<Scalar>& pvtwViscosity()
const
408 {
return pvtwViscosity_; }
410 const std::vector<Scalar>& pvtwViscosibility()
const
411 {
return pvtwViscosibility_; }
413 const std::vector<TabulatedOneDFunction>& watvisctCurves()
const
414 {
return watvisctCurves_; }
416 const std::vector<TabulatedOneDFunction> internalEnergyCurves()
const
417 {
return internalEnergyCurves_; }
419 bool enableInternalEnergy()
const
420 {
return enableInternalEnergy_; }
422 const std::vector<Scalar>& watJTRefPres()
const
423 {
return watJTRefPres_; }
425 const std::vector<Scalar>& watJTC()
const
429 bool operator==(
const WaterPvtThermal<Scalar, enableBrine>& data)
const
431 if (isothermalPvt_ && !data.isothermalPvt_)
433 if (!isothermalPvt_ && data.isothermalPvt_)
436 return (!this->isoThermalPvt() ||
437 (*this->isoThermalPvt() == *data.isoThermalPvt())) &&
438 this->viscrefPress() == data.viscrefPress() &&
439 this->watdentRefTemp() == data.watdentRefTemp() &&
440 this->watdentCT1() == data.watdentCT1() &&
441 this->watdentCT2() == data.watdentCT2() &&
442 this->watdentCT2() == data.watdentCT2() &&
443 this->watJTRefPres() == data.watJTRefPres() &&
444 this->watJT() == data.watJT() &&
445 this->watJTC() == data.watJTC() &&
446 this->pvtwRefPress() == data.pvtwRefPress() &&
447 this->pvtwRefB() == data.pvtwRefB() &&
448 this->pvtwCompressibility() == data.pvtwCompressibility() &&
449 this->pvtwViscosity() == data.pvtwViscosity() &&
450 this->pvtwViscosibility() == data.pvtwViscosibility() &&
451 this->watvisctCurves() == data.watvisctCurves() &&
452 this->internalEnergyCurves() == data.internalEnergyCurves() &&
454 this->enableJouleThomson() == data.enableJouleThomson() &&
456 this->enableInternalEnergy() == data.enableInternalEnergy();
459 WaterPvtThermal<Scalar, enableBrine>& operator=(
const WaterPvtThermal<Scalar, enableBrine>& data)
461 if (data.isothermalPvt_)
462 isothermalPvt_ =
new IsothermalPvt(*data.isothermalPvt_);
464 isothermalPvt_ =
nullptr;
465 viscrefPress_ = data.viscrefPress_;
466 watdentRefTemp_ = data.watdentRefTemp_;
467 watdentCT1_ = data.watdentCT1_;
468 watdentCT2_ = data.watdentCT2_;
469 watJTRefPres_ = data.watJTRefPres_;
470 watJTC_ = data.watJTC_;
471 pvtwRefPress_ = data.pvtwRefPress_;
472 pvtwRefB_ = data.pvtwRefB_;
473 pvtwCompressibility_ = data.pvtwCompressibility_;
474 pvtwViscosity_ = data.pvtwViscosity_;
475 pvtwViscosibility_ = data.pvtwViscosibility_;
476 watvisctCurves_ = data.watvisctCurves_;
477 internalEnergyCurves_ = data.internalEnergyCurves_;
478 enableThermalDensity_ = data.enableThermalDensity_;
479 enableJouleThomson_ = data.enableJouleThomson_;
480 enableThermalViscosity_ = data.enableThermalViscosity_;
481 enableInternalEnergy_ = data.enableInternalEnergy_;
487 IsothermalPvt* isothermalPvt_;
491 std::vector<Scalar> viscrefPress_;
493 std::vector<Scalar> watdentRefTemp_;
494 std::vector<Scalar> watdentCT1_;
495 std::vector<Scalar> watdentCT2_;
497 std::vector<Scalar> watJTRefPres_;
498 std::vector<Scalar> watJTC_;
500 std::vector<Scalar> pvtwRefPress_;
501 std::vector<Scalar> pvtwRefB_;
502 std::vector<Scalar> pvtwCompressibility_;
503 std::vector<Scalar> pvtwViscosity_;
504 std::vector<Scalar> pvtwViscosibility_;
506 std::vector<TabulatedOneDFunction> watvisctCurves_;
509 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
511 bool enableThermalDensity_;
512 bool enableJouleThomson_;
513 bool enableThermalViscosity_;
514 bool enableInternalEnergy_;
Implements a linearly interpolated scalar function that depends on one variable.
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
This class represents the Pressure-Volume-Temperature relations of the water phase in the black-oil m...
Definition: WaterPvtMultiplexer.hpp:75
unsigned numRegions() const
Return the number of PVT regions which are considered by this PVT-object.
Definition: WaterPvtMultiplexer.hpp:141
const Scalar waterReferenceDensity(unsigned regionIdx)
Return the reference density which are considered by this PVT-object.
Definition: WaterPvtMultiplexer.hpp:147
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtMultiplexer.hpp:177
Evaluation viscosity(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:164
This class implements temperature dependence of the PVT properties of water.
Definition: WaterPvtThermal.hpp:50
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition: WaterPvtThermal.hpp:359
bool enableThermalViscosity() const
Returns true iff the viscosity of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:268
Evaluation viscosity(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:338
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition: WaterPvtThermal.hpp:250
bool enableThermalDensity() const
Returns true iff the density of the water phase is temperature dependent.
Definition: WaterPvtThermal.hpp:256
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the specific internal energy [J/kg] of water given a set of parameters.
Definition: WaterPvtThermal.hpp:278
bool enableJouleThomsony() const
Returns true iff Joule-Thomson effect for the water phase is active.
Definition: WaterPvtThermal.hpp:262
void setNumRegions(size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition: WaterPvtThermal.hpp:230