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

ECL thermal: assume heat capacities for constant volume instead for constant pressure

Browse Source 
this implies that the internal energy instead of the enthalpy is
specified by the fluid PVT classes.
This commit is contained in:
Andreas Lauser 2018-01-25 19:10:29 +01:00
parent be4245a397
commit 096c22be92
14 changed files with 108 additions and 95 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -971,10 +971,23 @@ public:
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: return oilPvt_->enthalpy(regionIdx, T, p, Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
case gasPhaseIdx: return gasPvt_->enthalpy(regionIdx, T, p, Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
case waterPhaseIdx: return waterPvt_->enthalpy(regionIdx, T, p);
case oilPhaseIdx:
return
oilPvt_->internalEnergy(regionIdx, T, p, Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
case gasPhaseIdx:
return
gasPvt_->internalEnergy(regionIdx, T, p, Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
case waterPhaseIdx:
return
waterPvt_->internalEnergy(regionIdx, T, p)
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
default: OPM_THROW(std::logic_error, "Unhandled phase index " << phaseIdx);
}

2
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp
View File

@ -167,7 +167,7 @@ public:
* \brief Returns the specific enthalpy [J/kg] of oil given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx OPM_UNUSED,
Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rs OPM_UNUSED) const

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

@ -163,7 +163,7 @@ public:
* \brief Returns the specific enthalpy [J/kg] of water given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx OPM_UNUSED,
Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& pressure OPM_UNUSED) const
{

2
opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp
View File

@ -169,7 +169,7 @@ public:
* \brief Returns the specific enthalpy [J/kg] of oil given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx OPM_UNUSED,
Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rs OPM_UNUSED) const

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

@ -195,7 +195,7 @@ public:
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx OPM_UNUSED,
Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rv OPM_UNUSED) const

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

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

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

@ -62,7 +62,7 @@ public:
{
enableThermalDensity_ = false;
enableThermalViscosity_ = false;
enableEnthalpy_ = false;
enableInternalEnergy_ = false;
}
~GasPvtThermal()
@ -88,7 +88,7 @@ public:
enableThermalDensity_ = deck.hasKeyword("GASDENT");
enableThermalViscosity_ = deck.hasKeyword("GASVISCT");
enableEnthalpy_ = deck.hasKeyword("SPECHEAT");
enableInternalEnergy_ = deck.hasKeyword("SPECHEAT");
unsigned numRegions = isothermalPvt_->numRegions();
setNumRegions(numRegions);
@ -121,15 +121,15 @@ public:
}
if (deck.hasKeyword("SPECHEAT")) {
// the specific enthalpy of gas. be aware that ecl only specifies the heat capacity
// the specific internal energy of gas. be aware that ecl only specifies the heat capacity
// (via the SPECHEAT keyword) and we need to integrate it ourselfs to get the
// enthalpy
// internal energy
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
const auto& temperatureColumn = specHeatTable.getColumn("TEMPERATURE");
const auto& cpGasColumn = specHeatTable.getColumn("CP_GAS");
const auto& cvGasColumn = specHeatTable.getColumn("CV_GAS");
std::vector<double> hSamples(temperatureColumn.size());
std::vector<double> uSamples(temperatureColumn.size());
// the specific enthalpy of vaporization. since ECL does not seem to
// feature a proper way to specify this quantity, we use the value for
@ -138,25 +138,25 @@ public:
// phases should depend on the phase composition
const Scalar hVap = 480.6e3; // [J / kg]
Scalar h = temperatureColumn[0]*cpGasColumn[0] + hVap;
Scalar u = temperatureColumn[0]*cvGasColumn[0] + hVap;
for (size_t i = 0;; ++i) {
hSamples[i] = h;
uSamples[i] = u;
if (i >= temperatureColumn.size() - 1)
break;
// integrate to the heat capacity from the current sampling point to the next
// one. this leads to a quadratic polynomial.
Scalar h0 = cpGasColumn[i];
Scalar h1 = cpGasColumn[i + 1];
Scalar h0 = cvGasColumn[i];
Scalar h1 = cvGasColumn[i + 1];
Scalar T0 = temperatureColumn[i];
Scalar T1 = temperatureColumn[i + 1];
Scalar m = (h1 - h0)/(T1 - T0);
Scalar deltaH = 0.5*m*(T1*T1 - T0*T0) + h0*(T1 - T0);
h += deltaH;
Scalar deltaU = 0.5*m*(T1*T1 - T0*T0) + h0*(T1 - T0);
u += deltaU;
}
enthalpyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), hSamples);
internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
}
}
}
@ -168,7 +168,7 @@ public:
void setNumRegions(size_t numRegions)
{
gasvisctCurves_.resize(numRegions);
enthalpyCurves_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
gasdentRefTemp_.resize(numRegions);
gasdentCT1_.resize(numRegions);
gasdentCT2_.resize(numRegions);
@ -196,22 +196,22 @@ public:
{ return enableThermalViscosity_; }
/*!
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
* \brief Returns the specific internal energy [J/kg] of gas given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rv OPM_UNUSED) const
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rv OPM_UNUSED) const
{
if (!enableEnthalpy_)
if (!enableInternalEnergy_)
OPM_THROW(std::runtime_error,
"Requested the enthalpy of oil but it is disabled");
"Requested the internal energy of oil but it is disabled");
// compute the specific enthalpy for the specified tempature. We use linear
// 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 enthalpyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
/*!
@ -358,12 +358,12 @@ private:
std::vector<Scalar> gasdentCT1_;
std::vector<Scalar> gasdentCT2_;
// piecewise linear curve representing the enthalpy of gas
std::vector<TabulatedOneDFunction> enthalpyCurves_;
// piecewise linear curve representing the internal energy of gas
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
bool enableThermalDensity_;
bool enableThermalViscosity_;
bool enableEnthalpy_;
bool enableInternalEnergy_;
};
} // namespace Opm

2
opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp
View File

@ -412,7 +412,7 @@ public:
* \brief Returns the specific enthalpy [J/kg] of oil given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx OPM_UNUSED,
Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rs OPM_UNUSED) const

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

@ -154,11 +154,11 @@ public:
* \brief Returns the specific enthalpy [J/kg] oil given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx,
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rs) const
{ OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.enthalpy(regionIdx, temperature, pressure, Rs)); return 0; }
{ OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.internalEnergy(regionIdx, temperature, pressure, Rs)); return 0; }
/*!
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.

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

@ -62,7 +62,7 @@ public:
{
enableThermalDensity_ = false;
enableThermalViscosity_ = false;
enableEnthalpy_ = false;
enableInternalEnergy_ = false;
}
~OilPvtThermal()
@ -88,7 +88,7 @@ public:
enableThermalDensity_ = deck.hasKeyword("OILDENT");
enableThermalViscosity_ = deck.hasKeyword("VISCREF");
enableEnthalpy_ = deck.hasKeyword("SPECHEAT");
enableInternalEnergy_ = deck.hasKeyword("SPECHEAT");
unsigned numRegions = isothermalPvt_->numRegions();
setNumRegions(numRegions);
@ -139,35 +139,35 @@ public:
}
if (deck.hasKeyword("SPECHEAT")) {
// the specific enthalpy of liquid oil. be aware that ecl only specifies the
// 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
// ourselfs to get the enthalpy
// ourselfs to get the internal energy
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& specheatTable = tables.getSpecheatTables()[regionIdx];
const auto& temperatureColumn = specheatTable.getColumn("TEMPERATURE");
const auto& cpOilColumn = specheatTable.getColumn("CP_OIL");
const auto& cvOilColumn = specheatTable.getColumn("CV_OIL");
std::vector<double> hSamples(temperatureColumn.size());
std::vector<double> uSamples(temperatureColumn.size());
Scalar h = temperatureColumn[0]*cpOilColumn[0];
Scalar u = temperatureColumn[0]*cvOilColumn[0];
for (size_t i = 0;; ++i) {
hSamples[i] = h;
uSamples[i] = u;
if (i >= temperatureColumn.size() - 1)
break;
// integrate to the heat capacity from the current sampling point to the next
// one. this leads to a quadratic polynomial.
Scalar h0 = cpOilColumn[i];
Scalar h1 = cpOilColumn[i + 1];
Scalar h0 = cvOilColumn[i];
Scalar h1 = cvOilColumn[i + 1];
Scalar T0 = temperatureColumn[i];
Scalar T1 = temperatureColumn[i + 1];
Scalar m = (h1 - h0)/(T1 - T0);
Scalar deltaH = 0.5*m*(T1*T1 - T0*T0) + h0*(T1 - T0);
h += deltaH;
Scalar deltaU = 0.5*m*(T1*T1 - T0*T0) + h0*(T1 - T0);
u += deltaU;
}
enthalpyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), hSamples);
internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
}
}
}
@ -182,7 +182,7 @@ public:
viscrefPress_.resize(numRegions);
viscrefRs_.resize(numRegions);
viscRef_.resize(numRegions);
enthalpyCurves_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
}
/*!
@ -207,22 +207,22 @@ public:
{ return viscrefRs_.size(); }
/*!
* \brief Returns the specific enthalpy [J/kg] of oil given a set of parameters.
* \brief Returns the specific internal energy [J/kg] of oil given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rs OPM_UNUSED) const
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rs OPM_UNUSED) const
{
if (!enableEnthalpy_)
if (!enableInternalEnergy_)
OPM_THROW(std::runtime_error,
"Requested the enthalpy of oil but it is disabled");
"Requested the internal energy of oil but it is disabled");
// compute the specific enthalpy for the specified tempature. We use linear
// 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 enthalpyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
/*!
@ -366,12 +366,12 @@ private:
std::vector<Scalar> oildentCT1_;
std::vector<Scalar> oildentCT2_;
// piecewise linear curve representing the enthalpy of oil
std::vector<TabulatedOneDFunction> enthalpyCurves_;
// piecewise linear curve representing the internal energy of oil
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
bool enableThermalDensity_;
bool enableThermalViscosity_;
bool enableEnthalpy_;
bool enableInternalEnergy_;
};
} // namespace Opm

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

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

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

@ -62,7 +62,7 @@ public:
{
enableThermalDensity_ = false;
enableThermalViscosity_ = false;
enableEnthalpy_ = false;
enableInternalEnergy_ = false;
}
~WaterPvtThermal()
@ -88,7 +88,7 @@ public:
enableThermalDensity_ = deck.hasKeyword("WATDENT");
enableThermalViscosity_ = deck.hasKeyword("VISCREF");
enableEnthalpy_ = deck.hasKeyword("SPECHEAT");
enableInternalEnergy_ = deck.hasKeyword("SPECHEAT");
unsigned numRegions = isothermalPvt_->numRegions();
setNumRegions(numRegions);
@ -125,35 +125,35 @@ public:
}
if (deck.hasKeyword("SPECHEAT")) {
// the specific enthalpy of liquid water. be aware that ecl only specifies the heat capacity
// the specific internal energy of liquid water. be aware that ecl only specifies the heat capacity
// (via the SPECHEAT keyword) and we need to integrate it ourselfs to get the
// enthalpy
// internal energy
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& specHeatTable = tables.getSpecheatTables()[regionIdx];
const auto& temperatureColumn = specHeatTable.getColumn("TEMPERATURE");
const auto& cpWaterColumn = specHeatTable.getColumn("CP_WATER");
const auto& cvWaterColumn = specHeatTable.getColumn("CV_WATER");
std::vector<double> hSamples(temperatureColumn.size());
std::vector<double> uSamples(temperatureColumn.size());
Scalar h = temperatureColumn[0]*cpWaterColumn[0];
Scalar u = temperatureColumn[0]*cvWaterColumn[0];
for (size_t i = 0;; ++i) {
hSamples[i] = h;
uSamples[i] = u;
if (i >= temperatureColumn.size() - 1)
break;
// integrate to the heat capacity from the current sampling point to the next
// one. this leads to a quadratic polynomial.
Scalar h0 = cpWaterColumn[i];
Scalar h1 = cpWaterColumn[i + 1];
Scalar h0 = cvWaterColumn[i];
Scalar h1 = cvWaterColumn[i + 1];
Scalar T0 = temperatureColumn[i];
Scalar T1 = temperatureColumn[i + 1];
Scalar m = (h1 - h0)/(T1 - T0);
Scalar deltaH = 0.5*m*(T1*T1 - T0*T0) + h0*(T1 - T0);
h += deltaH;
Scalar deltaU = 0.5*m*(T1*T1 - T0*T0) + h0*(T1 - T0);
u += deltaU;
}
enthalpyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), hSamples);
internalEnergyCurves_[regionIdx].setXYContainers(temperatureColumn.vectorCopy(), uSamples);
}
}
}
@ -174,7 +174,7 @@ public:
watdentRefTemp_.resize(numRegions);
watdentCT1_.resize(numRegions);
watdentCT2_.resize(numRegions);
enthalpyCurves_.resize(numRegions);
internalEnergyCurves_.resize(numRegions);
}
/*!
@ -198,22 +198,22 @@ public:
size_t numRegions() const
{ return pvtwRefPress_.size(); }
/*!
* \brief Returns the specific enthalpy [J/kg] of water given a set of parameters.
/*!
* \brief Returns the specific internal energy [J/kg] of water given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure OPM_UNUSED) const
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure OPM_UNUSED) const
{
if (!enableEnthalpy_)
if (!enableInternalEnergy_)
OPM_THROW(std::runtime_error,
"Requested the enthalpy of oil but it is disabled");
"Requested the internal energy of oil but it is disabled");
// compute the specific enthalpy for the specified tempature. We use linear
// 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 enthalpyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
return internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
}
/*!
@ -279,12 +279,12 @@ private:
std::vector<TabulatedOneDFunction> watvisctCurves_;
// piecewise linear curve representing the enthalpy of water
std::vector<TabulatedOneDFunction> enthalpyCurves_;
// piecewise linear curve representing the internal energy of water
std::vector<TabulatedOneDFunction> internalEnergyCurves_;
bool enableThermalDensity_;
bool enableThermalViscosity_;
bool enableEnthalpy_;
bool enableInternalEnergy_;
};
} // namespace Opm

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

@ -442,7 +442,7 @@ public:
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
*/
template <class Evaluation>
Evaluation enthalpy(unsigned regionIdx OPM_UNUSED,
Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED,
const Evaluation& temperature OPM_UNUSED,
const Evaluation& pressure OPM_UNUSED,
const Evaluation& Rv OPM_UNUSED) const

4
opm/material/thermal/EclThermalLawManager.hpp
View File

@ -207,8 +207,8 @@ private:
auto& specrockParams = multiplexerParams.template getRealParams<SolidEnergyLawParams::specrockApproach>();
const auto& temperatureColumn = specrockTable.getColumn("TEMPERATURE");
const auto& cpRockColumn = specrockTable.getColumn("CP_ROCK");
specrockParams.setHeatCapacities(temperatureColumn, cpRockColumn);
const auto& cvRockColumn = specrockTable.getColumn("CV_ROCK");
specrockParams.setHeatCapacities(temperatureColumn, cvRockColumn);
specrockParams.finalize();
multiplexerParams.finalize();