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

Merge remote-tracking branch 'trine/add-chi-flash-trine' into add-chi-flash-trine

Browse Source
This commit is contained in:
Kai Bao 2022-06-27 10:59:31 +02:00
commit 79d89bd7ea
10 changed files with 350 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
opm/material/components/Component.hpp
View File

@ -103,6 +103,19 @@ public:
static Scalar criticalPressure()
{ throw std::runtime_error("Not implemented: Component::criticalPressure()"); }
/*!
* \brief Returns the acentric factor of the component.
*/
static Scalar acentricFactor()
{ throw std::runtime_error("Not implemented: acentricFactor of the component"); }
/*!
* \brief Returns the critical volume in \f$\mathrm{[m2/kmol]}\f$ of the component.
*/
static Scalar criticalVolume()
{ throw std::runtime_error("Not implemented: criticalVolume of the compoenent"); }
/*!
* \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point.
*/

12
opm/material/components/TabulatedComponent.hpp
View File

@ -233,6 +233,18 @@ public:
static Scalar criticalPressure()
{ return RawComponent::criticalPressure(); }
/*!
* \brief Returns the acentric factor of the component.
*/
static Scalar acentricFactor()
{ throw std::runtime_error("Not implemented: acentricFactor of the component"); }
/*!
* \brief Returns the critical volume in \f$\mathrm{[m2/kmol]}\f$ of the component.
*/
static Scalar criticalVolume()
{ throw std::runtime_error("Not implemented: criticalVolume of the compoenent"); }
/*!
* \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point.
*/

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

@ -1201,12 +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(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
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));

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

@ -75,7 +75,7 @@ public:
{
setNumRegions(numRegions);
for (size_t i = 0; i < numRegions; ++i) {
gasReferenceDensity_[i] = CO2::gasDensity(T_ref, P_ref, true);
gasReferenceDensity_[i] = CO2::gasDensity(T_ref, P_ref, extrapolate);
}
}
#if HAVE_ECL_INPUT
@ -238,7 +238,7 @@ public:
const Evaluation& pressure,
unsigned /*compIdx*/) const
{
return BinaryCoeffBrineCO2::gasDiffCoeff(temperature, pressure);
return BinaryCoeffBrineCO2::gasDiffCoeff(temperature, pressure, extrapolate);
}
const Scalar gasReferenceDensity(unsigned regionIdx) const

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
{

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

@ -60,8 +60,8 @@ public:
GasPvtThermal()
{
enableThermalDensity_ = false;
enableJouleThomson_ = false;
enableThermalViscosity_ = false;
enableInternalEnergy_ = false;
isothermalPvt_ = nullptr;
}
@ -70,8 +70,11 @@ 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>& gasJTC,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableJouleThomson,
bool enableThermalViscosity,
bool enableInternalEnergy)
: isothermalPvt_(isothermalPvt)
@ -79,8 +82,11 @@ public:
, gasdentRefTemp_(gasdentRefTemp)
, gasdentCT1_(gasdentCT1)
, gasdentCT2_(gasdentCT2)
, gasJTRefPres_(gasJTRefPres)
, gasJTC_(gasJTC)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableJouleThomson_(enableJouleThomson)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
{ }
@ -109,6 +115,7 @@ public:
const auto& tables = eclState.getTableManager();
enableThermalDensity_ = tables.GasDenT().size() > 0;
enableJouleThomson_ = tables.GasJT().size() > 0;
enableThermalViscosity_ = tables.hasTables("GASVISCT");
enableInternalEnergy_ = tables.hasTables("SPECHEAT");
@ -129,7 +136,7 @@ public:
}
// temperature dependence of gas density
if (tables.GasDenT().size() > 0) {
if (enableThermalDensity_) {
const auto& gasDenT = tables.GasDenT();
assert(gasDenT.size() == numRegions);
@ -140,7 +147,26 @@ public:
gasdentCT1_[regionIdx] = record.C1;
gasdentCT2_[regionIdx] = record.C2;
}
enableThermalDensity_ = true;
}
// Joule Thomson
if (enableJouleThomson_) {
const auto& gasJT = tables.GasJT();
assert(gasJT.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& record = gasJT[regionIdx];
gasJTRefPres_[regionIdx] = record.P0;
gasJTC_[regionIdx] = record.C1;
}
const auto& densityTable = eclState.getTableManager().getDensityTable();
assert(densityTable.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
rhoRefO_[regionIdx] = densityTable[regionIdx].oil;
}
}
if (enableInternalEnergy_) {
@ -193,6 +219,9 @@ public:
gasdentRefTemp_.resize(numRegions);
gasdentCT1_.resize(numRegions);
gasdentCT2_.resize(numRegions);
gasJTRefPres_.resize(numRegions);
gasJTC_.resize(numRegions);
rhoRefO_.resize(numRegions);
}
/*!
@ -210,6 +239,12 @@ public:
bool enableThermalDensity() const
{ return enableThermalDensity_; }
/*!
* \brief Returns true iff Joule-Thomson effect for the gas phase is active.
*/
bool enableJouleThomsony() const
{ return enableJouleThomson_; }
/*!
* \brief Returns true iff the viscosity of the gas phase is temperature dependent.
*/
@ -222,16 +257,62 @@ public:
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation&,
const Evaluation&) const
const Evaluation& pressure,
const Evaluation& Rv) 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 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);
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 = gasdentRefTemp_[regionIdx];
Evaluation Pref = gasJTRefPres_[regionIdx];
Scalar JTC = gasJTC_[regionIdx]; // if JTC is default then JTC is calculated
Evaluation Rvw = 0.0;
Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
const Scalar hVap = 480.6e3; // [J / kg]
Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true) - hVap)/temperature;
Evaluation density = invB * (gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
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, Rvw) *
(gasReferenceDensity(regionIdx) + Rv * rhoRefO_[regionIdx]);
// see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
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;
}
}
/*!
@ -411,6 +492,13 @@ public:
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }
const std::vector<Scalar>& gasJTRefPres() const
{ return gasJTRefPres_; }
const std::vector<Scalar>& gasJTC() const
{ return gasJTC_; }
bool operator==(const GasPvtThermal<Scalar>& data) const
{
if (isothermalPvt_ && !data.isothermalPvt_)
@ -424,8 +512,11 @@ public:
this->gasdentRefTemp() == data.gasdentRefTemp() &&
this->gasdentCT1() == data.gasdentCT1() &&
this->gasdentCT2() == data.gasdentCT2() &&
this->gasJTRefPres() == data.gasJTRefPres() &&
this->gasJTC() == data.gasJTC() &&
this->internalEnergyCurves() == data.internalEnergyCurves() &&
this->enableThermalDensity() == data.enableThermalDensity() &&
this->enableJouleThomson() == data.enableJouleThomson() &&
this->enableThermalViscosity() == data.enableThermalViscosity() &&
this->enableInternalEnergy() == data.enableInternalEnergy();
}
@ -440,8 +531,11 @@ public:
gasdentRefTemp_ = data.gasdentRefTemp_;
gasdentCT1_ = data.gasdentCT1_;
gasdentCT2_ = data.gasdentCT2_;
gasJTRefPres_ = data.gasJTRefPres_;
gasJTC_ = data.gasJTC_;
internalEnergyCurves_ = data.internalEnergyCurves_;
enableThermalDensity_ = data.enableThermalDensity_;
enableJouleThomson_ = data.enableJouleThomson_;
enableThermalViscosity_ = data.enableThermalViscosity_;
enableInternalEnergy_ = data.enableInternalEnergy_;
@ -459,10 +553,16 @@ private:
std::vector<Scalar> gasdentCT1_;
std::vector<Scalar> gasdentCT2_;
std::vector<Scalar> gasJTRefPres_;
std::vector<Scalar> gasJTC_;
std::vector<Scalar> rhoRefO_;
// piecewise linear curve representing the internal energy of gas
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
bool enableThermalDensity_;
bool enableJouleThomson_;
bool enableThermalViscosity_;
bool enableInternalEnergy_;
};

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

@ -60,6 +60,7 @@ public:
OilPvtThermal()
{
enableThermalDensity_ = false;
enableJouleThomson_ = false;
enableThermalViscosity_ = false;
enableInternalEnergy_ = false;
isothermalPvt_ = nullptr;
@ -73,8 +74,11 @@ 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>& oilJTC,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableJouleThomson,
bool enableThermalViscosity,
bool enableInternalEnergy)
: isothermalPvt_(isothermalPvt)
@ -85,8 +89,11 @@ public:
, oildentRefTemp_(oildentRefTemp)
, oildentCT1_(oildentCT1)
, oildentCT2_(oildentCT2)
, oilJTRefPres_(oilJTRefPres)
, oilJTC_(oilJTC)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableJouleThomson_(enableJouleThomson)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
{ }
@ -167,6 +174,26 @@ public:
}
}
// Joule Thomson
if (enableJouleThomson_) {
const auto& oilJT = tables.OilJT();
assert(oilJT.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& record = oilJT[regionIdx];
oilJTRefPres_[regionIdx] = record.P0;
oilJTC_[regionIdx] = record.C1;
}
const auto& densityTable = eclState.getTableManager().getDensityTable();
assert(densityTable.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
rhoRefG_[regionIdx] = densityTable[regionIdx].gas;
}
}
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 +237,12 @@ public:
viscrefRs_.resize(numRegions);
viscRef_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
oildentRefTemp_.resize(numRegions);
oildentCT1_.resize(numRegions);
oildentCT2_.resize(numRegions);
oilJTRefPres_.resize(numRegions);
oilJTC_.resize(numRegions);
rhoRefG_.resize(numRegions);
}
/*!
@ -224,6 +257,12 @@ public:
bool enableThermalDensity() const
{ return enableThermalDensity_; }
/*!
* \brief Returns true iff Joule-Thomson effect for the oil phase is active.
*/
bool enableJouleThomsony() const
{ return enableJouleThomson_; }
/*!
* \brief Returns true iff the viscosity of the oil phase is temperature dependent.
*/
@ -239,16 +278,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);
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]; // if JTC is default then JTC is calculated
Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rs);
Evaluation Cp = internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true)/temperature;
Evaluation density = invB * (oilReferenceDensity(regionIdx) + Rs * rhoRefG_[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) + Rs * rhoRefG_[regionIdx]) ;
// see e.g.https://en.wikipedia.org/wiki/Joule-Thomson_effect for a derivation of the Joule-Thomson coeff.
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 +501,12 @@ public:
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }
const std::vector<Scalar>& oilJTRefPres() const
{ return oilJTRefPres_; }
const std::vector<Scalar>& oilJTC() const
{ return oilJTC_; }
bool operator==(const OilPvtThermal<Scalar>& data) const
{
if (isothermalPvt_ && !data.isothermalPvt_)
@ -434,8 +523,11 @@ public:
this->oildentRefTemp() == data.oildentRefTemp() &&
this->oildentCT1() == data.oildentCT1() &&
this->oildentCT2() == data.oildentCT2() &&
this->oilJTRefPres() == data.oilJTRefPres() &&
this->oilJTC() == data.oilJTC() &&
this->internalEnergyCurves() == data.internalEnergyCurves() &&
this->enableThermalDensity() == data.enableThermalDensity() &&
this->enableJouleThomson() == data.enableJouleThomson() &&
this->enableThermalViscosity() == data.enableThermalViscosity() &&
this->enableInternalEnergy() == data.enableInternalEnergy();
}
@ -453,8 +545,11 @@ public:
oildentRefTemp_ = data.oildentRefTemp_;
oildentCT1_ = data.oildentCT1_;
oildentCT2_ = data.oildentCT2_;
oilJTRefPres_ = data.oilJTRefPres_;
oilJTC_ = data.oilJTC_;
internalEnergyCurves_ = data.internalEnergyCurves_;
enableThermalDensity_ = data.enableThermalDensity_;
enableJouleThomson_ = data.enableJouleThomson_;
enableThermalViscosity_ = data.enableThermalViscosity_;
enableInternalEnergy_ = data.enableInternalEnergy_;
@ -476,10 +571,16 @@ private:
std::vector<Scalar> oildentCT1_;
std::vector<Scalar> oildentCT2_;
std::vector<Scalar> oilJTRefPres_;
std::vector<Scalar> oilJTC_;
std::vector<Scalar> rhoRefG_;
// piecewise linear curve representing the internal energy of oil
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
bool enableThermalDensity_;
bool enableJouleThomson_;
bool enableThermalViscosity_;
bool enableInternalEnergy_;
};

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.

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

@ -70,6 +70,8 @@ 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>& watJTC,
const std::vector<Scalar>& pvtwRefPress,
const std::vector<Scalar>& pvtwRefB,
const std::vector<Scalar>& pvtwCompressibility,
@ -78,6 +80,7 @@ public:
const std::vector<TabulatedOneDFunction>& watvisctCurves,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableJouleThomson,
bool enableThermalViscosity,
bool enableInternalEnergy)
: isothermalPvt_(isothermalPvt)
@ -85,6 +88,8 @@ public:
, watdentRefTemp_(watdentRefTemp)
, watdentCT1_(watdentCT1)
, watdentCT2_(watdentCT2)
, watJTRefPres_(watJTRefPres)
, watJTC_(watJTC)
, pvtwRefPress_(pvtwRefPress)
, pvtwRefB_(pvtwRefB)
, pvtwCompressibility_(pvtwCompressibility)
@ -93,6 +98,7 @@ public:
, watvisctCurves_(watvisctCurves)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableJouleThomson_(enableJouleThomson)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
{ }
@ -121,6 +127,7 @@ public:
const auto& tables = eclState.getTableManager();
enableThermalDensity_ = tables.WatDenT().size() > 0;
enableJouleThomson_ = tables.WatJT().size() > 0;
enableThermalViscosity_ = tables.hasTables("WATVISCT");
enableInternalEnergy_ = tables.hasTables("SPECHEAT");
@ -142,12 +149,26 @@ public:
const auto& pvtwTables = tables.getPvtwTable();
assert(pvtwTables.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
pvtwRefPress_[regionIdx] = pvtwTables[regionIdx].reference_pressure;
pvtwRefB_[regionIdx] = pvtwTables[regionIdx].volume_factor;
}
}
// Joule Thomson
if (enableJouleThomson_) {
const auto& watJT = tables.WatJT();
assert(watJT.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& record = watJT[regionIdx];
watJTRefPres_[regionIdx] = record.P0;
watJTC_[regionIdx] = record.C1;
}
}
if (enableThermalViscosity_) {
if (tables.getViscrefTable().empty())
throw std::runtime_error("VISCREF is required when WATVISCT is present");
@ -223,6 +244,8 @@ public:
watdentRefTemp_.resize(numRegions);
watdentCT1_.resize(numRegions);
watdentCT2_.resize(numRegions);
watJTRefPres_.resize(numRegions);
watJTC_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
}
@ -238,6 +261,12 @@ public:
bool enableThermalDensity() const
{ return enableThermalDensity_; }
/*!
* \brief Returns true iff Joule-Thomson effect for the water phase is active.
*/
bool enableJouleThomsony() const
{ return enableJouleThomson_; }
/*!
* \brief Returns true iff the viscosity of the water phase is temperature dependent.
*/
@ -253,15 +282,58 @@ 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);
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]; // if JTC is default then JTC is 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 +424,13 @@ public:
bool enableInternalEnergy() const
{ return enableInternalEnergy_; }
const std::vector<Scalar>& watJTRefPres() const
{ return watJTRefPres_; }
const std::vector<Scalar>& watJTC() const
{ return watJTC_; }
bool operator==(const WaterPvtThermal<Scalar, enableBrine>& data) const
{
if (isothermalPvt_ && !data.isothermalPvt_)
@ -365,6 +444,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() &&
@ -373,6 +456,7 @@ public:
this->watvisctCurves() == data.watvisctCurves() &&
this->internalEnergyCurves() == data.internalEnergyCurves() &&
this->enableThermalDensity() == data.enableThermalDensity() &&
this->enableJouleThomson() == data.enableJouleThomson() &&
this->enableThermalViscosity() == data.enableThermalViscosity() &&
this->enableInternalEnergy() == data.enableInternalEnergy();
}
@ -387,6 +471,8 @@ public:
watdentRefTemp_ = data.watdentRefTemp_;
watdentCT1_ = data.watdentCT1_;
watdentCT2_ = data.watdentCT2_;
watJTRefPres_ = data.watJTRefPres_;
watJTC_ = data.watJTC_;
pvtwRefPress_ = data.pvtwRefPress_;
pvtwRefB_ = data.pvtwRefB_;
pvtwCompressibility_ = data.pvtwCompressibility_;
@ -395,6 +481,7 @@ public:
watvisctCurves_ = data.watvisctCurves_;
internalEnergyCurves_ = data.internalEnergyCurves_;
enableThermalDensity_ = data.enableThermalDensity_;
enableJouleThomson_ = data.enableJouleThomson_;
enableThermalViscosity_ = data.enableThermalViscosity_;
enableInternalEnergy_ = data.enableInternalEnergy_;
@ -412,6 +499,9 @@ private:
std::vector<Scalar> watdentCT1_;
std::vector<Scalar> watdentCT2_;
std::vector<Scalar> watJTRefPres_;
std::vector<Scalar> watJTC_;
std::vector<Scalar> pvtwRefPress_;
std::vector<Scalar> pvtwRefB_;
std::vector<Scalar> pvtwCompressibility_;
@ -424,6 +514,7 @@ private:
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
bool enableThermalDensity_;
bool enableJouleThomson_;
bool enableThermalViscosity_;
bool enableInternalEnergy_;
};