Loading...
Searching...
No Matches
WaterPvtThermal.hpp
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_WATER_PVT_THERMAL_HPP
#define OPM_WATER_PVT_THERMAL_HPP
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <cstddef>
namespace Opm {
#if HAVE_ECL_INPUT
class EclipseState;
class Schedule;
#endif
template <class Scalar, bool enableThermal, bool enableBrine>
class WaterPvtMultiplexer;
template <class Scalar, bool enableBrine>
class WaterPvtThermal
{
public:
using TabulatedOneDFunction = Tabulated1DFunction<Scalar>;
using IsothermalPvt = WaterPvtMultiplexer<Scalar, /*enableThermal=*/false, enableBrine>;
WaterPvtThermal() = default;
WaterPvtThermal(IsothermalPvt* isothermalPvt,
const std::vector<Scalar>& viscrefPress,
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,
const std::vector<Scalar>& pvtwViscosity,
const std::vector<Scalar>& pvtwViscosibility,
const std::vector<TabulatedOneDFunction>& watvisctCurves,
const std::vector<TabulatedOneDFunction>& internalEnergyCurves,
bool enableThermalDensity,
bool enableJouleThomson,
bool enableThermalViscosity,
bool enableInternalEnergy)
: isothermalPvt_(isothermalPvt)
, viscrefPress_(viscrefPress)
, watdentRefTemp_(watdentRefTemp)
, watdentCT1_(watdentCT1)
, watdentCT2_(watdentCT2)
, watJTRefPres_(watJTRefPres)
, watJTC_(watJTC)
, pvtwRefPress_(pvtwRefPress)
, pvtwRefB_(pvtwRefB)
, pvtwCompressibility_(pvtwCompressibility)
, pvtwViscosity_(pvtwViscosity)
, pvtwViscosibility_(pvtwViscosibility)
, watvisctCurves_(watvisctCurves)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableJouleThomson_(enableJouleThomson)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
{ }
WaterPvtThermal(const WaterPvtThermal& data)
{ *this = data; }
~WaterPvtThermal()
{ delete isothermalPvt_; }
#if HAVE_ECL_INPUT
void initFromState(const EclipseState& eclState, const Schedule& schedule);
#endif // HAVE_ECL_INPUT
void setNumRegions(std::size_t numRegions);
void setVapPars(const Scalar par1, const Scalar par2)
{
isothermalPvt_->setVapPars(par1, par2);
}
void initEnd()
{ }
bool enableThermalDensity() const
{ return enableThermalDensity_; }
bool enableJouleThomson() const
{ return enableJouleThomson_; }
bool enableThermalViscosity() const
{ return enableThermalViscosity_; }
Scalar hVap(unsigned regionIdx) const
{ return this->hVap_[regionIdx]; }
std::size_t numRegions() const
{ return pvtwRefPress_.size(); }
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rsw,
const Evaluation& saltconcentration) const
{
if (!enableInternalEnergy_) {
throw std::runtime_error("Requested the internal energy of water but it is disabled");
}
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, Rsw, 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));
constexpr const int N = 100; // value is experimental
Evaluation deltaP = (pressure - Pref) / N;
Evaluation enthalpyPresPrev = 0;
for (std::size_t i = 0; i < N; ++i) {
Evaluation Pnew = Pref + i * deltaP;
Evaluation rho = inverseFormationVolumeFactor(regionIdx, temperature,
Pnew, Rsw, 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;
}
}
template <class Evaluation>
Evaluation viscosity(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rsw,
const Evaluation& saltconcentration) const
{
const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature,
pressure, Rsw, saltconcentration);
if (!enableThermalViscosity()) {
return isothermalMu;
}
Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] -
pvtwRefPress_[regionIdx]);
Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);
// compute the viscosity deviation due to temperature
const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
return isothermalMu * muWatvisct / muRef;
}
template <class Evaluation>
Evaluation saturatedViscosity(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& saltconcentration) const
{
const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature,
pressure, saltconcentration);
if (!enableThermalViscosity()) {
return isothermalMu;
}
Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] -
pvtwRefPress_[regionIdx]);
Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);
// compute the viscosity deviation due to temperature
const auto& muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
return isothermalMu * muWatvisct / muRef;
}
template <class Evaluation>
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& saltconcentration) const
{
Evaluation Rsw = 0.0;
return inverseFormationVolumeFactor(regionIdx, temperature, pressure,
Rsw, saltconcentration);
}
template <class Evaluation>
Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rsw,
const Evaluation& saltconcentration) const
{
if (!enableThermalDensity()) {
return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature,
pressure, Rsw, saltconcentration);
}
Scalar BwRef = pvtwRefB_[regionIdx];
Scalar TRef = watdentRefTemp_[regionIdx];
const Evaluation& X = pvtwCompressibility_[regionIdx] * (pressure -
pvtwRefPress_[regionIdx]);
Scalar cT1 = watdentCT1_[regionIdx];
Scalar cT2 = watdentCT2_[regionIdx];
const Evaluation& Y = temperature - TRef;
// this is inconsistent with the density calculation of water in the isothermal
// case (it misses the quadratic pressure term), but it is the equation given in
// the documentation.
return 1.0 / (((1 - X) * (1 + cT1 * Y + cT2 * Y * Y)) * BwRef);
}
template <class FluidState, class LhsEval = typename FluidState::Scalar>
std::pair<LhsEval, LhsEval>
inverseFormationVolumeFactorAndViscosity(const FluidState& fluidState, unsigned regionIdx)
{
auto [b, mu] = isothermalPvt_->inverseFormationVolumeFactorAndViscosity(fluidState, regionIdx);
const LhsEval& pressure = decay<LhsEval>(fluidState.pressure(FluidState::waterPhaseIdx));
const LhsEval& temperature = decay<LhsEval>(fluidState.temperature(FluidState::waterPhaseIdx));
if (enableThermalDensity()) {
Scalar BwRef = pvtwRefB_[regionIdx];
Scalar TRef = watdentRefTemp_[regionIdx];
const LhsEval X = pvtwCompressibility_[regionIdx] * (pressure - pvtwRefPress_[regionIdx]);
Scalar cT1 = watdentCT1_[regionIdx];
Scalar cT2 = watdentCT2_[regionIdx];
const LhsEval Y = temperature - TRef;
// this is inconsistent with the density calculation of water in the isothermal
// case (it misses the quadratic pressure term), but it is the equation given in
// the documentation.
// Note assignment to 'b', not multiplying the isothermal 'b' by a factor!
b = 1.0 / (((1 - X) * (1 + cT1 * Y + cT2 * Y * Y)) * BwRef);
}
if (enableThermalViscosity()) {
Scalar x = -pvtwViscosibility_[regionIdx] * (viscrefPress_[regionIdx] - pvtwRefPress_[regionIdx]);
Scalar muRef = pvtwViscosity_[regionIdx] / (1.0 + x + 0.5 * x * x);
// compute the viscosity deviation due to temperature
const auto muWatvisct = watvisctCurves_[regionIdx].eval(temperature, true);
// Split operations to get same order of operations as in the viscosity() function.
mu *= muWatvisct;
mu /= muRef;
}
return { b, mu };
}
template <class Evaluation>
Evaluation saturationPressure(unsigned /*regionIdx*/,
const Evaluation& /*temperature*/,
const Evaluation& /*Rs*/,
const Evaluation& /*saltconcentration*/) const
{ return 0.0; /* this is dead water, so there isn't any meaningful saturation pressure! */ }
template <class Evaluation>
Evaluation saturatedGasDissolutionFactor(unsigned /*regionIdx*/,
const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
const Evaluation& /*saltconcentration*/) const
{ return 0.0; /* this is dead water! */ }
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
unsigned /*compIdx*/) const
{
throw std::runtime_error("Not implemented: The PVT model does not provide "
"a diffusionCoefficient()");
}
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }
Scalar waterReferenceDensity(unsigned regionIdx) const
{ return isothermalPvt_->waterReferenceDensity(regionIdx); }
const std::vector<Scalar>& viscrefPress() const
{ return viscrefPress_; }
const std::vector<Scalar>& watdentRefTemp() const
{ return watdentRefTemp_; }
const std::vector<Scalar>& watdentCT1() const
{ return watdentCT1_; }
const std::vector<Scalar>& watdentCT2() const
{ return watdentCT2_; }
const std::vector<Scalar>& pvtwRefPress() const
{ return pvtwRefPress_; }
const std::vector<Scalar>& pvtwRefB() const
{ return pvtwRefB_; }
const std::vector<Scalar>& pvtwCompressibility() const
{ return pvtwCompressibility_; }
const std::vector<Scalar>& pvtwViscosity() const
{ return pvtwViscosity_; }
const std::vector<Scalar>& pvtwViscosibility() const
{ return pvtwViscosibility_; }
const std::vector<TabulatedOneDFunction>& watvisctCurves() const
{ return watvisctCurves_; }
const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
{ return internalEnergyCurves_; }
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;
WaterPvtThermal<Scalar, enableBrine>&
operator=(const WaterPvtThermal<Scalar, enableBrine>& data);
private:
IsothermalPvt* isothermalPvt_{nullptr};
// The PVT properties needed for temperature dependence. We need to store one
// value per PVT region.
std::vector<Scalar> viscrefPress_{};
std::vector<Scalar> watdentRefTemp_{};
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_{};
std::vector<Scalar> pvtwViscosity_{};
std::vector<Scalar> pvtwViscosibility_{};
std::vector<TabulatedOneDFunction> watvisctCurves_{};
// piecewise linear curve representing the internal energy of water
std::vector<TabulatedOneDFunction> internalEnergyCurves_{};
std::vector<Scalar> hVap_{};
bool enableThermalDensity_{false};
bool enableJouleThomson_{false};
bool enableThermalViscosity_{false};
bool enableInternalEnergy_{false};
};
} // namespace Opm
#endif
Tabulated1DFunction.hpp
Implements a linearly interpolated scalar function that depends on one variable.
Opm::EclipseState
Definition EclipseState.hpp:62
Opm::Schedule
Definition Schedule.hpp:101
Opm::Tabulated1DFunction
Implements a linearly interpolated scalar function that depends on one variable.
Definition Tabulated1DFunction.hpp:51
Opm::WaterPvtMultiplexer
This class represents the Pressure-Volume-Temperature relations of the water phase in the black-oil m...
Definition WaterPvtMultiplexer.hpp:90
Opm::WaterPvtMultiplexer::waterReferenceDensity
Scalar waterReferenceDensity(unsigned regionIdx) const
Return the reference density which are considered by this PVT-object.
Definition WaterPvtMultiplexer.cpp:111
Opm::WaterPvtThermal::enableThermalViscosity
bool enableThermalViscosity() const
Returns true iff the viscosity of the water phase is temperature dependent.
Definition WaterPvtThermal.hpp:141
Opm::WaterPvtThermal::saturationPressure
Evaluation saturationPressure(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the saturation pressure of the water phase [Pa] depending on its mass fraction of the gas com...
Definition WaterPvtThermal.hpp:349
Opm::WaterPvtThermal::viscosity
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition WaterPvtThermal.hpp:219
Opm::WaterPvtThermal::initEnd
void initEnd()
Finish initializing the thermal part of the water phase PVT properties.
Definition WaterPvtThermal.hpp:123
Opm::WaterPvtThermal::inverseFormationVolumeFactor
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition WaterPvtThermal.hpp:281
Opm::WaterPvtThermal::enableThermalDensity
bool enableThermalDensity() const
Returns true iff the density of the water phase is temperature dependent.
Definition WaterPvtThermal.hpp:129
Opm::WaterPvtThermal::saturatedGasDissolutionFactor
Evaluation saturatedGasDissolutionFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the gas dissolution factor [m^3/m^3] of the water phase.
Definition WaterPvtThermal.hpp:359
Opm::WaterPvtThermal::internalEnergy
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rsw, const Evaluation &saltconcentration) const
Returns the specific internal energy [J/kg] of water given a set of parameters.
Definition WaterPvtThermal.hpp:154
Opm::WaterPvtThermal::saturatedInverseFormationVolumeFactor
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &saltconcentration) const
Returns the formation volume factor [-] of the fluid phase.
Definition WaterPvtThermal.hpp:268
Opm::WaterPvtThermal::enableJouleThomson
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the water phase is active.
Definition WaterPvtThermal.hpp:135
Opm::WaterPvtThermal::saturatedViscosity
Evaluation saturatedViscosity(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:244
Opm::WaterPvtThermal::setNumRegions
void setNumRegions(std::size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition WaterPvtThermal.cpp:185
Opm::WaterPvtThermal::inverseFormationVolumeFactorAndViscosity
std::pair< LhsEval, LhsEval > inverseFormationVolumeFactorAndViscosity(const FluidState &fluidState, unsigned regionIdx)
Returns the formation volume factor [-] and viscosity [Pa s] of the fluid phase.
Definition WaterPvtThermal.hpp:311
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30