OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

adjusted internal energy calculation to include (optional) Joule-Thomson effect

Browse Source
This commit is contained in:
Paul Egberts 2022-04-12 21:07:44 +02:00
parent a721fbba1f
commit a2aa7b9603
11 changed files with 276 additions and 186 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
opm/material/fluidsystems/BlackOilFluidSystem.hpp
View File

@ -1193,37 +1193,6 @@ public:
const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
LhsEval p_st = 1.01325E5; // pressure at standard condition
LhsEval T_st = 20 + 273.15; // pressure at standard condition
LhsEval Joule_Thomson_coefficient;
LhsEval Tref = 399.15;
LhsEval Pref = 20E5;
LhsEval Cp;
LhsEval dns;
LhsEval alpha;
LhsEval JT1;
GasPvtThermal<Scalar> obj;
Scalar JTC;
if (phaseIdx == gasPhaseIdx) {
LhsEval intEn = gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
//std::cout << "Internal Energy: " << intEn << std::endl;
Cp = (intEn- 480.6E3)/T;
//dns = density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
//alpha = rho d(1/rho)/dT
//dns = 9.82;
// alpha = dns * (2.67941015e-03 + 2*6.86181705e-06*(T - Tref))/
// (0.64861 * (gasPvt_->inverseFormationVolumeFactor(regionIdx, Tref, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))));
// JT1 = -(1.0/Cv)*(1.0 - alpha *T)/dns;
//std::cout << "JT1: "<< JT1 << " Cv: "<< Cv << " dns: " << dns << " alpha: " << alpha << std::endl;
//std::cout << " dns: " << dns << " alpha: " << alpha << "\n" << std::endl;
//Joule_Thomson_coefficient = 3.26E-6;
//Joule_Thomson_coefficient = gasPvt_->internalEnergy_JT(regionIdx, Tref, Pref, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
//JTC = Joule_Thomson_coefficient.value(); //3.34718e-06
//std::cout << " Joule_Thomson_coefficient: " << JTC << "\n" << std::endl;
}
switch (phaseIdx) {
case oilPhaseIdx:
return
@ -1232,25 +1201,12 @@ public:
case gasPhaseIdx:
return
//gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
//+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
gasPvt_->internalEnergy_JT(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
//gasPvt_->internalEnergy_JT(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
//gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
//Cv * (T-T_st) - Cv * JT1 * (p - p_st);
//Cp.value() * (T-T_st) - Cp.value() * JTC * (p - p_st);
//gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
//- Cv * Joule_Thomson_coefficient * (p - p_st);
case waterPhaseIdx:
return
waterPvt_->internalEnergy(regionIdx, T, p)
waterPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));

9
opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp
View File

@ -147,15 +147,6 @@ public:
{
return CO2::gasInternalEnergy(temperature, pressure, extrapolate);
}
template <class Evaluation>
Evaluation internalEnergy_JT(unsigned,
const Evaluation&,
const Evaluation&,
const Evaluation&) const
{
throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
}
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.

1
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp
View File

@ -155,6 +155,7 @@ public:
*/
template <class Evaluation>
Evaluation internalEnergy(unsigned,
const Evaluation&,
const Evaluation&,
const Evaluation&) const
{

1
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp
View File

@ -165,6 +165,7 @@ public:
*/
template <class Evaluation>
Evaluation internalEnergy(unsigned,
const Evaluation&,
const Evaluation&,
const Evaluation&) const
{

8
opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp
View File

@ -212,14 +212,6 @@ public:
{
throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
}
template <class Evaluation>
Evaluation internalEnergy_JT(unsigned,
const Evaluation&,
const Evaluation&,
const Evaluation&) const
{
throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
}
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.

7
opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp
View File

@ -234,13 +234,6 @@ public:
const Evaluation& Rv) const
{ OPM_GAS_PVT_MULTIPLEXER_CALL(return pvtImpl.internalEnergy(regionIdx, temperature, pressure, Rv)); return 0; }
template <class Evaluation>
Evaluation internalEnergy_JT(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rv) const
{ OPM_GAS_PVT_MULTIPLEXER_CALL(return pvtImpl.internalEnergy_JT(regionIdx, temperature, pressure, Rv)); return 0; }
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
*/

175
opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
View File

@ -61,7 +61,6 @@ public:
{
enableThermalDensity_ = false;
enableThermalViscosity_ = false;
enableJouleThomson_ = false;
isothermalPvt_ = nullptr;
}
@ -70,21 +69,25 @@ public:
const std::vector<Scalar>& gasdentRefTemp,
const std::vector<Scalar>& gasdentCT1,
const std::vector<Scalar>& gasdentCT2,
const std::vector<Scalar>& gasJTRefPres,
const std::vector<Scalar>& gasJT,
const std::vector<Scalar>& gasJTC,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableThermalViscosity,
bool enableInternalEnergy,
bool enableJouleThomson)
bool enableInternalEnergy)
: isothermalPvt_(isothermalPvt)
, gasvisctCurves_(gasvisctCurves)
, gasdentRefTemp_(gasdentRefTemp)
, gasdentCT1_(gasdentCT1)
, gasdentCT2_(gasdentCT2)
, gasJTRefPres_(gasJTRefPres)
, gasJT_(gasJT)
, gasJTC_(gasJTC)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
, enableJouleThomson_(enableJouleThomson)
{ }
GasPvtThermal(const GasPvtThermal& data)
@ -113,7 +116,6 @@ public:
enableThermalDensity_ = tables.GasDenT().size() > 0;
enableThermalViscosity_ = tables.hasTables("GASVISCT");
enableInternalEnergy_ = tables.hasTables("SPECHEAT");
enableJouleThomson_ = tables.hasTables("JOULETHOMSON");
unsigned numRegions = isothermalPvt_->numRegions();
setNumRegions(numRegions);
@ -132,7 +134,7 @@ public:
}
// temperature dependence of gas density
if (tables.GasDenT().size() > 0) {
if (enableThermalDensity_) {
const auto& gasDenT = tables.GasDenT();
assert(gasDenT.size() == numRegions);
@ -143,21 +145,18 @@ public:
gasdentCT1_[regionIdx] = record.C1;
gasdentCT2_[regionIdx] = record.C2;
}
enableThermalDensity_ = true;
}
// temperature dependence of gas density
if (enableJouleThomson_) {
//const auto& joulethomsonTable = tables.joulethomsonTable();
// assert(gasDenT.size() == numRegions);
// for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
// const auto& record = gasDenT[regionIdx];
// Joule Thomson
if (tables.hasTables("JOULETHO")) {
const auto& joulethomsonTables = tables.getJoulethomsonTables();
// gasdentRefTemp_[regionIdx] = record.T0;
// gasdentCT1_[regionIdx] = record.C1;
// gasdentCT2_[regionIdx] = record.C2;
// }
// enableThermalDensity_ = true;
assert(joulethomsonTables.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
gasJTRefPres_[regionIdx] = joulethomsonTables[regionIdx].getColumn("PREF").front();
gasJT_[regionIdx] = joulethomsonTables[regionIdx].getColumn("GAS_JT").front();
gasJTC_[regionIdx] = joulethomsonTables[regionIdx].getColumn("GAS_JTC").front();
}
}
if (enableInternalEnergy_) {
@ -210,8 +209,9 @@ public:
gasdentRefTemp_.resize(numRegions);
gasdentCT1_.resize(numRegions);
gasdentCT2_.resize(numRegions);
joulethomsonRefPres_.resize(numRegions);
joulethomsonGas_.resize(numRegions);
gasJTRefPres_.resize(numRegions);
gasJT_.resize(numRegions);
gasJTC_.resize(numRegions);
}
/*!
@ -239,7 +239,7 @@ public:
* \brief Returns the specific internal energy [J/kg] of gas given a set of parameters.
*/
template <class Evaluation>
Evaluation internalEnergy_JT(unsigned regionIdx,
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rv) const
@ -247,73 +247,57 @@ public:
if (!enableInternalEnergy_)
throw std::runtime_error("Requested the internal energy of gas but it is disabled");
bool JouleThomson = 1 ; //joulethomsonGas_[regionIdx];
bool enableJouleThomson_ = gasJT_[regionIdx];
Evaluation Tref = gasdentRefTemp_[regionIdx];// 399.15;
Evaluation Pref = joulethomsonRefPres_[regionIdx]; //20E5;
Evaluation bref = inverseFormationVolumeFactor(regionIdx, Tref, pressure, Rv);
Evaluation b = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv);
Evaluation intEn = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
Evaluation Cp = (intEn- 480.6E3)/temperature;
Evaluation rho = b * gasReferenceDensity(regionIdx);
Scalar c1T = gasdentCT1_[regionIdx];
Scalar c2T = gasdentCT2_[regionIdx];
Evaluation alpha = (c1T + 2*c2T*(temperature - Tref))/
(1 + c1T*(temperature - Tref) + c2T*(temperature- Tref)*(temperature - Tref));
Evaluation beta = bref.derivative(1)/bref; //derivative wrt to p
//Cp = Cv + alpha*alpha/(beta*rho);
const int N = 1000;
Evaluation deltaP = (pressure - Pref)/N;
Evaluation Hp_old = 0;
Evaluation Hp, JT;
for (size_t i = 0; i < N; ++i) {
Evaluation p_new = Pref + i * deltaP;
Evaluation alpha_new = (c1T + 2*c2T*(temperature- Tref))/(1 + c1T*(temperature - Tref) + c2T*(temperature - Tref)*(temperature - Tref));
// Evaluation rho_new = inverseFormationVolumeFactor(regionIdx, temperature, p_new, Rv) * gasReferenceDensity(regionIdx);
//JT = -(1.0/Cp)*(1.0 - alpha_new.value() * temperature.value())/rho_new.value();
Evaluation rho_new = inverseFormationVolumeFactor(regionIdx, temperature, p_new, Rv) * gasReferenceDensity(regionIdx);
JT = -(1.0/Cp)*(1.0 - alpha_new.value() * temperature.value())/rho_new.value();
Evaluation dHp = -Cp * JT * deltaP;
Hp = Hp_old + dHp;
Hp_old = Hp;
if (!enableJouleThomson_) {
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
std::cout << "P: "<< pressure << " T: " << temperature << " JT: "<< JT << " Hp: " << Hp << "\n"<<std::flush;
else {
Evaluation Tref = gasdentRefTemp_[regionIdx];
Evaluation Pref = gasJTRefPres_[regionIdx];
Scalar JTC = gasJTC_[regionIdx]; // JTC = 0: JTC calculated
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv);
const Scalar hVap = 480.6e3; // [J / kg]
Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true) - hVap)/temperature;
Evaluation density = invB * gasReferenceDensity(regionIdx);
Evaluation a1 = alpha * pressure/rho;
Evaluation a2 = (beta * pressure -alpha * temperature)/rho;
return Cp * (temperature-Tref) + Hp - pressure/rho;
Evaluation enthalpyPres;
if (JTC != 0) {
enthalpyPres = -Cp * JTC * (pressure -Pref);
}
else if(enableThermalDensity_) {
Scalar c1T = gasdentCT1_[regionIdx];
Scalar c2T = gasdentCT2_[regionIdx];
Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
(1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
const int N = 100; // value is experimental
Evaluation deltaP = (pressure - Pref)/N;
Evaluation enthalpyPresPrev = 0;
for (size_t i = 0; i < N; ++i) {
Evaluation Pnew = Pref + i * deltaP;
Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rv) * gasReferenceDensity(regionIdx);
Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
enthalpyPresPrev = enthalpyPres;
}
}
else {
throw std::runtime_error("Requested Joule-thomson calculation but thermal gas density (GASDENT) is not provided");
}
Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
return enthalpy - pressure/density;
}
}
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation&,
const Evaluation&) const
{
if (!enableInternalEnergy_)
throw std::runtime_error("Requested the internal energy of gas but it is disabled");
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
*/
@ -491,14 +475,15 @@ public:
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }
bool enableJouleThomson() const
{ return enableJouleThomson_; }
const std::vector<Scalar>& gasJTRefPres() const
{ return gasJTRefPres_; }
const std::vector<Scalar>& joulethomsonRefPres() const
{ return joulethomsonRefPres_; }
const std::vector<Scalar>& gasJT() const
{ return gasJT_; }
const std::vector<Scalar>& gasJTC() const
{ return gasJTC_; }
const std::vector<Scalar>& joulethomsonGas() const
{ return joulethomsonGas_; }
bool operator==(const GasPvtThermal<Scalar>& data) const
{
@ -513,11 +498,13 @@ public:
this->gasdentRefTemp() == data.gasdentRefTemp() &&
this->gasdentCT1() == data.gasdentCT1() &&
this->gasdentCT2() == data.gasdentCT2() &&
this->gasJTRefPres() == data.gasJTRefPres() &&
this->gasJT() == data.gasJT() &&
this->gasJTC() == data.gasJTC() &&
this->internalEnergyCurves() == data.internalEnergyCurves() &&
this->enableThermalDensity() == data.enableThermalDensity() &&
this->enableThermalViscosity() == data.enableThermalViscosity() &&
this->enableInternalEnergy() == data.enableInternalEnergy();
this->enableJouleThomson() == data.enableJouleThomson();
}
GasPvtThermal<Scalar>& operator=(const GasPvtThermal<Scalar>& data)
@ -530,11 +517,13 @@ public:
gasdentRefTemp_ = data.gasdentRefTemp_;
gasdentCT1_ = data.gasdentCT1_;
gasdentCT2_ = data.gasdentCT2_;
gasJTRefPres_ = data.gasJTRefPres_;
gasJT_ = data.gasJT_;
gasJTC_ = data.gasJTC_;
internalEnergyCurves_ = data.internalEnergyCurves_;
enableThermalDensity_ = data.enableThermalDensity_;
enableThermalViscosity_ = data.enableThermalViscosity_;
enableInternalEnergy_ = data.enableInternalEnergy_;
enableJouleThomson_ = data.enableJouleThomson_;
return *this;
}
@ -550,8 +539,9 @@ private:
std::vector<Scalar> gasdentCT1_;
std::vector<Scalar> gasdentCT2_;
std::vector<Scalar> joulethomsonRefPres_;
std::vector<int> joulethomsonGas_;
std::vector<Scalar> gasJTRefPres_;
std::vector<Scalar> gasJT_;
std::vector<Scalar> gasJTC_;
// piecewise linear curve representing the internal energy of gas
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
@ -559,7 +549,6 @@ private:
bool enableThermalDensity_;
bool enableThermalViscosity_;
bool enableInternalEnergy_;
bool enableJouleThomson_;
};
} // namespace Opm

101
opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp
View File

@ -73,6 +73,9 @@ public:
const std::vector<Scalar>& oildentRefTemp,
const std::vector<Scalar>& oildentCT1,
const std::vector<Scalar>& oildentCT2,
const std::vector<Scalar>& oilJTRefPres,
const std::vector<Scalar>& oilJT,
const std::vector<Scalar>& oilJTC,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableThermalViscosity,
@ -85,6 +88,9 @@ public:
, oildentRefTemp_(oildentRefTemp)
, oildentCT1_(oildentCT1)
, oildentCT2_(oildentCT2)
, oilJTRefPres_(oilJTRefPres)
, oilJT_(oilJT)
, oilJTC_(oilJTC)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableThermalViscosity_(enableThermalViscosity)
@ -167,6 +173,18 @@ public:
}
}
// Joule Thomson
if (tables.hasTables("JOULETHO")) {
const auto& joulethomsonTables = tables.getJoulethomsonTables();
assert(joulethomsonTables.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
oilJTRefPres_[regionIdx] = joulethomsonTables[regionIdx].getColumn("PREF").front();
oilJT_[regionIdx] = joulethomsonTables[regionIdx].getColumn("OIL_JT").front();
oilJTC_[regionIdx] = joulethomsonTables[regionIdx].getColumn("OIL_JTC").front();
}
}
if (enableInternalEnergy_) {
// the specific internal energy of liquid oil. be aware that ecl only specifies the
// heat capacity (via the SPECHEAT keyword) and we need to integrate it
@ -210,6 +228,12 @@ public:
viscrefRs_.resize(numRegions);
viscRef_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
oildentRefTemp_.resize(numRegions);
oildentCT1_.resize(numRegions);
oildentCT2_.resize(numRegions);
oilJTRefPres_.resize(numRegions);
oilJT_.resize(numRegions);
oilJTC_.resize(numRegions);
}
/*!
@ -239,16 +263,60 @@ public:
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation&,
const Evaluation&) const
const Evaluation& pressure,
const Evaluation& Rs) const
{
if (!enableInternalEnergy_)
throw std::runtime_error("Requested the internal energy of oil but it is disabled");
throw std::runtime_error("Requested the internal energy of oil but it is disabled");
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
bool enableJouleThomson_ = oilJT_[regionIdx];
if (!enableJouleThomson_) {
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
else {
Evaluation Tref = oildentRefTemp_[regionIdx];
Evaluation Pref = oilJTRefPres_[regionIdx];
Scalar JTC = oilJTC_[regionIdx]; // JTC = 0: JTC calculated
Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rs);
Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
Evaluation density = invB * oilReferenceDensity(regionIdx);
Evaluation enthalpyPres;
if (JTC != 0) {
enthalpyPres = -Cp * JTC * (pressure -Pref);
}
else if(enableThermalDensity_) {
Scalar c1T = oildentCT1_[regionIdx];
Scalar c2T = oildentCT2_[regionIdx];
Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
(1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
const int N = 100; // value is experimental
Evaluation deltaP = (pressure - Pref)/N;
Evaluation enthalpyPresPrev = 0;
for (size_t i = 0; i < N; ++i) {
Evaluation Pnew = Pref + i * deltaP;
Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, Rs) * oilReferenceDensity(regionIdx);
Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
enthalpyPresPrev = enthalpyPres;
}
}
else {
throw std::runtime_error("Requested Joule-thomson calculation but thermal oil density (OILDENT) is not provided");
}
Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
return enthalpy - pressure/density;
}
}
/*!
@ -418,6 +486,15 @@ public:
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }
const std::vector<Scalar>& oilJTRefPres() const
{ return oilJTRefPres_; }
const std::vector<Scalar>& oilJT() const
{ return oilJT_; }
const std::vector<Scalar>& oilJTC() const
{ return oilJTC_; }
bool operator==(const OilPvtThermal<Scalar>& data) const
{
if (isothermalPvt_ && !data.isothermalPvt_)
@ -434,6 +511,9 @@ public:
this->oildentRefTemp() == data.oildentRefTemp() &&
this->oildentCT1() == data.oildentCT1() &&
this->oildentCT2() == data.oildentCT2() &&
this->oilJTRefPres() == data.oilJTRefPres() &&
this->oilJT() == data.oilJT() &&
this->oilJTC() == data.oilJTC() &&
this->internalEnergyCurves() == data.internalEnergyCurves() &&
this->enableThermalDensity() == data.enableThermalDensity() &&
this->enableThermalViscosity() == data.enableThermalViscosity() &&
@ -453,6 +533,9 @@ public:
oildentRefTemp_ = data.oildentRefTemp_;
oildentCT1_ = data.oildentCT1_;
oildentCT2_ = data.oildentCT2_;
oilJTRefPres_ = data.oilJTRefPres_;
oilJT_ = data.oilJT_;
oilJTC_ = data.oilJTC_;
internalEnergyCurves_ = data.internalEnergyCurves_;
enableThermalDensity_ = data.enableThermalDensity_;
enableThermalViscosity_ = data.enableThermalViscosity_;
@ -476,6 +559,10 @@ private:
std::vector<Scalar> oildentCT1_;
std::vector<Scalar> oildentCT2_;
std::vector<Scalar> oilJTRefPres_;
std::vector<Scalar> oilJT_;
std::vector<Scalar> oilJTC_;
// piecewise linear curve representing the internal energy of oil
std::vector<TabulatedOneDFunction> internalEnergyCurves_;

5
opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp
View File

@ -153,8 +153,9 @@ public:
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{ OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.internalEnergy(regionIdx, temperature, pressure)); return 0; }
const Evaluation& pressure,
const Evaluation& saltconcentration) const
{ OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.internalEnergy(regionIdx, temperature, pressure, saltconcentration)); return 0; }
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.

99
opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp
View File

@ -70,6 +70,9 @@ public:
const std::vector<Scalar>& watdentRefTemp,
const std::vector<Scalar>& watdentCT1,
const std::vector<Scalar>& watdentCT2,
const std::vector<Scalar>& watJTRefPres,
const std::vector<Scalar>& watJT,
const std::vector<Scalar>& watJTC,
const std::vector<Scalar>& pvtwRefPress,
const std::vector<Scalar>& pvtwRefB,
const std::vector<Scalar>& pvtwCompressibility,
@ -85,6 +88,9 @@ public:
, watdentRefTemp_(watdentRefTemp)
, watdentCT1_(watdentCT1)
, watdentCT2_(watdentCT2)
, watJTRefPres_(watJTRefPres)
, watJT_(watJT)
, watJTC_(watJTC)
, pvtwRefPress_(pvtwRefPress)
, pvtwRefB_(pvtwRefB)
, pvtwCompressibility_(pvtwCompressibility)
@ -148,6 +154,18 @@ public:
}
}
// Joule Thomson
if (tables.hasTables("JOULETHO")) {
const auto& joulethomsonTables = tables.getJoulethomsonTables();
assert(joulethomsonTables.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
watJTRefPres_[regionIdx] = joulethomsonTables[regionIdx].getColumn("PREF").front();
watJT_[regionIdx] = joulethomsonTables[regionIdx].getColumn("WAT_JT").front();
watJTC_[regionIdx] = joulethomsonTables[regionIdx].getColumn("WAT_JTC").front();
}
}
if (enableThermalViscosity_) {
if (tables.getViscrefTable().empty())
throw std::runtime_error("VISCREF is required when WATVISCT is present");
@ -223,6 +241,9 @@ public:
watdentRefTemp_.resize(numRegions);
watdentCT1_.resize(numRegions);
watdentCT2_.resize(numRegions);
watJTRefPres_.resize(numRegions);
watJT_.resize(numRegions);
watJTC_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
}
@ -253,15 +274,60 @@ public:
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation&) const
const Evaluation& pressure,
const Evaluation& saltconcentration) const
{
if (!enableInternalEnergy_)
throw std::runtime_error("Requested the internal energy of oil but it is disabled");
throw std::runtime_error("Requested the internal energy of water but it is disabled");
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
bool enableJouleThomson_ = watJT_[regionIdx];
if (!enableJouleThomson_) {
// compute the specific internal energy for the specified tempature. We use linear
// interpolation here despite the fact that the underlying heat capacities are
// piecewise linear (which leads to a quadratic function)
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
else {
Evaluation Tref = watdentRefTemp_[regionIdx];
Evaluation Pref = watJTRefPres_[regionIdx];
Scalar JTC =watJTC_[regionIdx]; // JTC = 0: JTC calculated
Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration);
Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
Evaluation density = invB * waterReferenceDensity(regionIdx);
Evaluation enthalpyPres;
if (JTC != 0) {
enthalpyPres = -Cp * JTC * (pressure -Pref);
}
else if(enableThermalDensity_) {
Scalar c1T = watdentCT1_[regionIdx];
Scalar c2T = watdentCT2_[regionIdx];
Evaluation alpha = (c1T + 2 * c2T * (temperature - Tref)) /
(1 + c1T *(temperature - Tref) + c2T * (temperature - Tref) * (temperature - Tref));
const int N = 100; // value is experimental
Evaluation deltaP = (pressure - Pref)/N;
Evaluation enthalpyPresPrev = 0;
for (size_t i = 0; i < N; ++i) {
Evaluation Pnew = Pref + i * deltaP;
Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature, Pnew, saltconcentration) * waterReferenceDensity(regionIdx);
Evaluation jouleThomsonCoefficient = -(1.0/Cp) * (1.0 - alpha * temperature)/rho;
Evaluation deltaEnthalpyPres = -Cp * jouleThomsonCoefficient * deltaP;
enthalpyPres = enthalpyPresPrev + deltaEnthalpyPres;
enthalpyPresPrev = enthalpyPres;
}
}
else {
throw std::runtime_error("Requested Joule-thomson calculation but thermal water density (WATDENT) is not provided");
}
Evaluation enthalpy = Cp * (temperature - Tref) + enthalpyPres;
return enthalpy - pressure/density;
}
}
/*!
@ -352,6 +418,16 @@ public:
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }
const std::vector<Scalar>& watJTRefPres() const
{ return watJTRefPres_; }
const std::vector<Scalar>& watJT() const
{ return watJT_; }
const std::vector<Scalar>& watJTC() const
{ return watJTC_; }
bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
{
if (isothermalPvt_ && !data.isothermalPvt_)
@ -365,6 +441,10 @@ public:
this->watdentRefTemp() == data.watdentRefTemp() &&
this->watdentCT1() == data.watdentCT1() &&
this->watdentCT2() == data.watdentCT2() &&
this->watdentCT2() == data.watdentCT2() &&
this->watJTRefPres() == data.watJTRefPres() &&
this->watJT() == data.watJT() &&
this->watJTC() == data.watJTC() &&
this->pvtwRefPress() == data.pvtwRefPress() &&
this->pvtwRefB() == data.pvtwRefB() &&
this->pvtwCompressibility() == data.pvtwCompressibility() &&
@ -387,6 +467,9 @@ public:
watdentRefTemp_ = data.watdentRefTemp_;
watdentCT1_ = data.watdentCT1_;
watdentCT2_ = data.watdentCT2_;
watJTRefPres_ = data.watJTRefPres_;
watJT_ = data.watJT_;
watJTC_ = data.watJTC_;
pvtwRefPress_ = data.pvtwRefPress_;
pvtwRefB_ = data.pvtwRefB_;
pvtwCompressibility_ = data.pvtwCompressibility_;
@ -412,6 +495,10 @@ private:
std::vector<Scalar> watdentCT1_;
std::vector<Scalar> watdentCT2_;
std::vector<Scalar> watJTRefPres_;
std::vector<Scalar> watJT_;
std::vector<Scalar> watJTC_;
std::vector<Scalar> pvtwRefPress_;
std::vector<Scalar> pvtwRefB_;
std::vector<Scalar> pvtwCompressibility_;

8
opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp
View File

@ -478,14 +478,6 @@ public:
{
throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
}
template <class Evaluation>
Evaluation internalEnergy_JT(unsigned,
const Evaluation&,
const Evaluation&,
const Evaluation&) const
{
throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
}
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.