57 enableThermalDensity_ =
false;
58 enableJouleThomson_ =
false;
59 enableThermalViscosity_ =
false;
60 enableInternalEnergy_ =
false;
61 isothermalPvt_ =
nullptr;
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)
89 , pvtwRefPress_(pvtwRefPress)
91 , pvtwCompressibility_(pvtwCompressibility)
92 , pvtwViscosity_(pvtwViscosity)
93 , pvtwViscosibility_(pvtwViscosibility)
94 , watvisctCurves_(watvisctCurves)
95 , internalEnergyCurves_(internalEnergyCurves)
99 , enableInternalEnergy_(enableInternalEnergy)
106 {
delete isothermalPvt_; }
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);
135 void setVapPars(
const Scalar par1,
const Scalar par2)
137 isothermalPvt_->setVapPars(par1, par2);
150 {
return enableThermalDensity_; }
156 {
return enableJouleThomson_; }
162 {
return enableThermalViscosity_; }
164 size_t numRegions()
const
165 {
return pvtwRefPress_.size(); }
170 template <
class Evaluation>
172 const Evaluation& temperature,
173 const Evaluation& pressure,
174 const Evaluation& Rsw,
175 const Evaluation& saltconcentration)
const
177 if (!enableInternalEnergy_)
178 throw std::runtime_error(
"Requested the internal energy of water but it is disabled");
180 if (!enableJouleThomson_) {
184 return internalEnergyCurves_[regionIdx].eval(temperature,
true);
187 Evaluation Tref = watdentRefTemp_[regionIdx];
188 Evaluation Pref = watJTRefPres_[regionIdx];
189 Scalar JTC =watJTC_[regionIdx];
192 Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature,
true)/temperature;
193 Evaluation density = invB * waterReferenceDensity(regionIdx);
195 Evaluation enthalpyPres;
197 enthalpyPres = -Cp * JTC * (pressure -Pref);
199 else if(enableThermalDensity_) {
200 Scalar c1T = watdentCT1_[regionIdx];
201 Scalar c2T = watdentCT2_[regionIdx];
203 Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
204 (1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
206 constexpr const int N = 100;
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;
212 Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
213 Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
214 enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
215 enthalpyPresPrev = enthalpyPres;
219 throw std::runtime_error(
"Requested Joule-thomson calculation but thermal water density (WATDENT) is not provided");
222 Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
224 return enthalpy - pressure/density;
231 template <
class Evaluation>
233 const Evaluation& temperature,
234 const Evaluation& pressure,
235 const Evaluation& Rsw,
236 const Evaluation& saltconcentration)
const
238 const auto& isothermalMu = isothermalPvt_->
viscosity(regionIdx, temperature, pressure, Rsw, saltconcentration);
242 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
243 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
246 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature,
true);
247 return isothermalMu * muWatvisct/muRef;
253 template <
class Evaluation>
255 const Evaluation& temperature,
256 const Evaluation& pressure,
257 const Evaluation& saltconcentration)
const
259 const auto& isothermalMu = isothermalPvt_->
saturatedViscosity(regionIdx, temperature, pressure, saltconcentration);
263 Scalar x = -pvtwViscosibility_[regionIdx]*(viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
264 Scalar muRef = pvtwViscosity_[regionIdx]/(1.0 + x + 0.5*x*x);
267 const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature,
true);
268 return isothermalMu * muWatvisct/muRef;
274 template <
class Evaluation>
276 const Evaluation& temperature,
277 const Evaluation& pressure,
278 const Evaluation& saltconcentration)
const
280 Evaluation Rsw = 0.0;
286 template <
class Evaluation>
288 const Evaluation& temperature,
289 const Evaluation& pressure,
290 const Evaluation& Rsw,
291 const Evaluation& saltconcentration)
const
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;
306 return 1.0/(((1 - X)*(1 + cT1*Y + cT2*Y*Y))*BwRef);
315 template <
class Evaluation>
319 const Evaluation& )
const
325 template <
class Evaluation>
329 const Evaluation& )
const
332 template <
class Evaluation>
333 Evaluation diffusionCoefficient(
const Evaluation& ,
337 throw std::runtime_error(
"Not implemented: The PVT model does not provide a diffusionCoefficient()");
340 const IsothermalPvt* isoThermalPvt()
const
341 {
return isothermalPvt_; }
343 const Scalar waterReferenceDensity(
unsigned regionIdx)
const
346 const std::vector<Scalar>& viscrefPress()
const
347 {
return viscrefPress_; }
349 const std::vector<Scalar>& watdentRefTemp()
const
350 {
return watdentRefTemp_; }
352 const std::vector<Scalar>& watdentCT1()
const
353 {
return watdentCT1_; }
355 const std::vector<Scalar>& watdentCT2()
const
356 {
return watdentCT2_; }
358 const std::vector<Scalar>& pvtwRefPress()
const
359 {
return pvtwRefPress_; }
361 const std::vector<Scalar>& pvtwRefB()
const
362 {
return pvtwRefB_; }
364 const std::vector<Scalar>& pvtwCompressibility()
const
365 {
return pvtwCompressibility_; }
367 const std::vector<Scalar>& pvtwViscosity()
const
368 {
return pvtwViscosity_; }
370 const std::vector<Scalar>& pvtwViscosibility()
const
371 {
return pvtwViscosibility_; }
373 const std::vector<TabulatedOneDFunction>& watvisctCurves()
const
374 {
return watvisctCurves_; }
376 const std::vector<TabulatedOneDFunction> internalEnergyCurves()
const
377 {
return internalEnergyCurves_; }
379 bool enableInternalEnergy()
const
380 {
return enableInternalEnergy_; }
382 const std::vector<Scalar>& watJTRefPres()
const
383 {
return watJTRefPres_; }
385 const std::vector<Scalar>& watJTC()
const
389 bool operator==(
const WaterPvtThermal<Scalar, enableBrine>& data)
const
391 if (isothermalPvt_ && !data.isothermalPvt_)
393 if (!isothermalPvt_ && data.isothermalPvt_)
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() &&
413 this->enableInternalEnergy() == data.enableInternalEnergy();
416 WaterPvtThermal<Scalar, enableBrine>& operator=(
const WaterPvtThermal<Scalar, enableBrine>& data)
418 if (data.isothermalPvt_)
419 isothermalPvt_ =
new IsothermalPvt(*data.isothermalPvt_);
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_;
444 IsothermalPvt* isothermalPvt_;
448 std::vector<Scalar> viscrefPress_;
450 std::vector<Scalar> watdentRefTemp_;
451 std::vector<Scalar> watdentCT1_;
452 std::vector<Scalar> watdentCT2_;
454 std::vector<Scalar> watJTRefPres_;
455 std::vector<Scalar> watJTC_;
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_;
463 std::vector<TabulatedOneDFunction> watvisctCurves_;
466 std::vector<TabulatedOneDFunction> internalEnergyCurves_;
468 bool enableThermalDensity_;
469 bool enableJouleThomson_;
470 bool enableThermalViscosity_;
471 bool enableInternalEnergy_;