Loading...
Searching...
No Matches
GasPvtThermal.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_GAS_PVT_THERMAL_HPP
#define OPM_GAS_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>
class GasPvtMultiplexer;
template <class Scalar>
class GasPvtThermal
{
public:
using IsothermalPvt = GasPvtMultiplexer<Scalar, /*enableThermal=*/false>;
using TabulatedOneDFunction = Tabulated1DFunction<Scalar>;
GasPvtThermal() = default;
GasPvtThermal(IsothermalPvt* isothermalPvt,
const std::vector<TabulatedOneDFunction>& gasvisctCurves,
const std::vector<Scalar>& viscrefPress,
const std::vector<Scalar>& viscRef,
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)
, gasvisctCurves_(gasvisctCurves)
, viscrefPress_(viscrefPress)
, viscRef_(viscRef)
, gasdentRefTemp_(gasdentRefTemp)
, gasdentCT1_(gasdentCT1)
, gasdentCT2_(gasdentCT2)
, gasJTRefPres_(gasJTRefPres)
, gasJTC_(gasJTC)
, internalEnergyCurves_(internalEnergyCurves)
, enableThermalDensity_(enableThermalDensity)
, enableJouleThomson_(enableJouleThomson)
, enableThermalViscosity_(enableThermalViscosity)
, enableInternalEnergy_(enableInternalEnergy)
{ }
GasPvtThermal(const GasPvtThermal& data)
{ *this = data; }
~GasPvtThermal()
{ 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_; }
std::size_t numRegions() const
{ return viscrefPress_.size(); }
template <class Evaluation>
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rv,
[[maybe_unused]] const Evaluation& RvW) const
{
if (!enableInternalEnergy_) {
throw std::runtime_error("Requested the internal energy of gas 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 {
// NB should probably be revisited with adding more or restrict to linear internal energy
OpmLog::warning("Experimental code for jouleThomson: simulation will be slower");
Evaluation Tref = gasdentRefTemp_[regionIdx];
Evaluation Pref = gasJTRefPres_[regionIdx];
Scalar JTC = gasJTC_[regionIdx]; // if JTC is default then JTC is calculaited
Evaluation Rvw = 0.0;
Evaluation invB = inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rv, Rvw);
// NB this assumes internalEnergyCurve(0) = 0 derivative should be used add Cp table ??
Evaluation Cp = (internalEnergyCurves_[regionIdx].eval(temperature, /*extrapolate=*/true))/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));
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, 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;
}
}
template <class Evaluation>
Evaluation viscosity(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rv,
const Evaluation& Rvw) const
{
const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw);
if (!enableThermalViscosity())
return isothermalMu;
// compute the viscosity deviation due to temperature
const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
return muGasvisct/viscRef_[regionIdx]*isothermalMu;
}
template <class Evaluation>
Evaluation saturatedViscosity(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{
const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
if (!enableThermalViscosity())
return isothermalMu;
// compute the viscosity deviation due to temperature
const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true);
return muGasvisct/viscRef_[regionIdx]*isothermalMu;
}
template <class Evaluation>
Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rv,
const Evaluation& /*Rvw*/) const
{
const Evaluation& Rvw = 0.0;
const auto& b =
isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature,
pressure, Rv, Rvw);
if (!enableThermalDensity()) {
return b;
}
// we use the same approach as for the for water here, but with the OPM-specific
// GASDENT keyword.
//
// TODO: Since gas is quite a bit more compressible than water, it might be
// necessary to make GASDENT to a table keyword. If the current temperature
// is relatively close to the reference temperature, the current approach
// should be good enough, though.
Scalar TRef = gasdentRefTemp_[regionIdx];
Scalar cT1 = gasdentCT1_[regionIdx];
Scalar cT2 = gasdentCT2_[regionIdx];
const Evaluation& Y = temperature - TRef;
return b / (1 + (cT1 + cT2 * Y) * Y);
}
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& temperature = decay<LhsEval>(fluidState.temperature(FluidState::gasPhaseIdx));
if (enableThermalDensity()) {
// we use the same approach as for the for water here, but with the OPM-specific
// GASDENT keyword.
//
// TODO: Since gas is quite a bit more compressible than water, it might be
// necessary to make GASDENT to a table keyword. If the current temperature
// is relatively close to the reference temperature, the current approach
// should be good enough, though.
Scalar TRef = gasdentRefTemp_[regionIdx];
Scalar cT1 = gasdentCT1_[regionIdx];
Scalar cT2 = gasdentCT2_[regionIdx];
const LhsEval& Y = temperature - TRef;
b /= (1.0 + (cT1 + cT2 * Y) * Y);
}
if (enableThermalViscosity()) {
// compute the viscosity deviation due to temperature
const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
mu *= (muGasvisct / viscRef_[regionIdx]);
}
return { b, mu };
}
template <class Evaluation>
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{
const auto& b =
isothermalPvt_->saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
if (!enableThermalDensity()) {
return b;
}
// we use the same approach as for the for water here, but with the OPM-specific
// GASDENT keyword.
//
// TODO: Since gas is quite a bit more compressible than water, it might be
// necessary to make GASDENT to a table keyword. If the current temperature
// is relatively close to the reference temperature, the current approach
// should be good enough, though.
Scalar TRef = gasdentRefTemp_[regionIdx];
Scalar cT1 = gasdentCT1_[regionIdx];
Scalar cT2 = gasdentCT2_[regionIdx];
const Evaluation& Y = temperature - TRef;
return b / (1 + (cT1 + cT2 * Y) * Y);
}
template <class Evaluation>
Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/) const
{ return 0.0; }
template <class Evaluation = Scalar>
Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
const Evaluation& /*saltConcentration*/) const
{ return 0.0; }
template <class Evaluation>
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{ return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature, pressure); }
template <class Evaluation>
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& oilSaturation,
const Evaluation& maxOilSaturation) const
{
return isothermalPvt_->saturatedOilVaporizationFactor(regionIdx, temperature,
pressure, oilSaturation,
maxOilSaturation);
}
template <class Evaluation>
Evaluation saturationPressure(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{ return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
const Evaluation& pressure,
unsigned compIdx) const
{
return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
}
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }
Scalar gasReferenceDensity(unsigned regionIdx) const
{ return isothermalPvt_->gasReferenceDensity(regionIdx); }
Scalar hVap(unsigned regionIdx) const
{ return this->hVap_[regionIdx]; }
const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
{ return gasvisctCurves_; }
const std::vector<Scalar>& viscrefPress() const
{ return viscrefPress_; }
const std::vector<Scalar>& viscRef() const
{ return viscRef_; }
const std::vector<Scalar>& gasdentRefTemp() const
{ return gasdentRefTemp_; }
const std::vector<Scalar>& gasdentCT1() const
{ return gasdentCT1_; }
const std::vector<Scalar>& gasdentCT2() const
{ return gasdentCT2_; }
const std::vector<TabulatedOneDFunction>& internalEnergyCurves() const
{ return internalEnergyCurves_; }
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;
GasPvtThermal<Scalar>& operator=(const GasPvtThermal<Scalar>& data);
private:
IsothermalPvt* isothermalPvt_{nullptr};
// The PVT properties needed for temperature dependence of the viscosity. We need
// to store one value per PVT region.
std::vector<TabulatedOneDFunction> gasvisctCurves_{};
std::vector<Scalar> viscrefPress_{};
std::vector<Scalar> viscRef_{};
std::vector<Scalar> gasdentRefTemp_{};
std::vector<Scalar> gasdentCT1_{};
std::vector<Scalar> gasdentCT2_{};
std::vector<Scalar> gasJTRefPres_{};
std::vector<Scalar> gasJTC_{};
std::vector<Scalar> rhoRefO_{};
std::vector<Scalar> hVap_{};
// piecewise linear curve representing the internal energy of gas
std::vector<TabulatedOneDFunction> internalEnergyCurves_{};
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::GasPvtMultiplexer
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition GasPvtMultiplexer.hpp:111
Opm::GasPvtMultiplexer::gasReferenceDensity
Scalar gasReferenceDensity(unsigned regionIdx) const
Return the reference density which are considered by this PVT-object.
Definition GasPvtMultiplexer.cpp:57
Opm::GasPvtMultiplexer::diffusionCoefficient
Evaluation diffusionCoefficient(const Evaluation &temperature, const Evaluation &pressure, unsigned compIdx) const
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition GasPvtMultiplexer.hpp:283
Opm::GasPvtThermal::enableThermalDensity
bool enableThermalDensity() const
Returns true iff the density of the gas phase is temperature dependent.
Definition GasPvtThermal.hpp:121
Opm::GasPvtThermal::viscosity
Evaluation viscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &Rvw) const
Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
Definition GasPvtThermal.hpp:212
Opm::GasPvtThermal::saturatedOilVaporizationFactor
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &oilSaturation, const Evaluation &maxOilSaturation) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition GasPvtThermal.hpp:379
Opm::GasPvtThermal::saturatedInverseFormationVolumeFactor
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the formation volume factor [-] of oil-saturated gas.
Definition GasPvtThermal.hpp:313
Opm::GasPvtThermal::saturatedOilVaporizationFactor
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the oil vaporization factor [m^3/m^3] of the gas phase.
Definition GasPvtThermal.hpp:366
Opm::GasPvtThermal::saturatedWaterVaporizationFactor
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of water saturated gas.
Definition GasPvtThermal.hpp:352
Opm::GasPvtThermal::internalEnergy
Evaluation internalEnergy(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &RvW) const
Returns the specific internal energy [J/kg] of gas given a set of parameters.
Definition GasPvtThermal.hpp:143
Opm::GasPvtThermal::setNumRegions
void setNumRegions(std::size_t numRegions)
Set the number of PVT-regions considered by this object.
Definition GasPvtThermal.cpp:193
Opm::GasPvtThermal::initEnd
void initEnd()
Finish initializing the thermal part of the gas phase PVT properties.
Definition GasPvtThermal.hpp:115
Opm::GasPvtThermal::inverseFormationVolumeFactor
Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure, const Evaluation &Rv, const Evaluation &) const
Returns the formation volume factor [-] of the fluid phase.
Definition GasPvtThermal.hpp:248
Opm::GasPvtThermal::enableJouleThomson
bool enableJouleThomson() const
Returns true iff Joule-Thomson effect for the gas phase is active.
Definition GasPvtThermal.hpp:127
Opm::GasPvtThermal::enableThermalViscosity
bool enableThermalViscosity() const
Returns true iff the viscosity of the gas phase is temperature dependent.
Definition GasPvtThermal.hpp:133
Opm::GasPvtThermal::saturatedWaterVaporizationFactor
Evaluation saturatedWaterVaporizationFactor(unsigned, const Evaluation &, const Evaluation &) const
Returns the water vaporization factor [m^3/m^3] of the water phase.
Definition GasPvtThermal.hpp:343
Opm::GasPvtThermal::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 GasPvtThermal.hpp:283
Opm::GasPvtThermal::saturatedViscosity
Evaluation saturatedViscosity(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the dynamic viscosity [Pa s] of the oil-saturated gas phase given a set of parameters.
Definition GasPvtThermal.hpp:231
Opm::GasPvtThermal::saturationPressure
Evaluation saturationPressure(unsigned regionIdx, const Evaluation &temperature, const Evaluation &pressure) const
Returns the saturation pressure of the gas phase [Pa].
Definition GasPvtThermal.hpp:398
Opm::Schedule
Definition Schedule.hpp:101
Opm::Tabulated1DFunction
Implements a linearly interpolated scalar function that depends on one variable.
Definition Tabulated1DFunction.hpp:51
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30