diff --git a/opm/input/eclipse/EclipseState/Tables/GasvisctTable.hpp b/opm/input/eclipse/EclipseState/Tables/GasvisctTable.hpp index f75256741..07dfe5480 100644 --- a/opm/input/eclipse/EclipseState/Tables/GasvisctTable.hpp +++ b/opm/input/eclipse/EclipseState/Tables/GasvisctTable.hpp @@ -24,15 +24,14 @@ namespace Opm { - class Deck; class DeckItem; class GasvisctTable : public SimpleTable { public: - GasvisctTable( const Deck& deck, const DeckItem& deckItem ); + GasvisctTable( const DeckItem& item, const int tableID ); const TableColumn& getTemperatureColumn() const; - const TableColumn& getGasViscosityColumn(size_t compIdx) const; + const TableColumn& getGasViscosityColumn() const; }; } diff --git a/opm/input/eclipse/EclipseState/Tables/TableManager.hpp b/opm/input/eclipse/EclipseState/Tables/TableManager.hpp index a45dbb310..7dda2781f 100644 --- a/opm/input/eclipse/EclipseState/Tables/TableManager.hpp +++ b/opm/input/eclipse/EclipseState/Tables/TableManager.hpp @@ -289,7 +289,6 @@ namespace Opm { void initRTempTables(const Deck& deck); void initDims(const Deck& deck); void initRocktabTables(const Deck& deck); - void initGasvisctTables(const Deck& deck); void initPlymaxTables(const Deck& deck); void initPlyrockTables(const Deck& deck); diff --git a/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp index a30311a5c..94261b329 100644 --- a/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp +++ b/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp @@ -57,11 +57,14 @@ public: enableThermalDensity_ = false; enableJouleThomson_ = false; enableThermalViscosity_ = false; + enableInternalEnergy_ = false; isothermalPvt_ = nullptr; } GasPvtThermal(IsothermalPvt* isothermalPvt, const std::vector& gasvisctCurves, + const std::vector& viscrefPress, + const std::vector& viscRef, const std::vector& gasdentRefTemp, const std::vector& gasdentCT1, const std::vector& gasdentCT2, @@ -74,6 +77,8 @@ public: bool enableInternalEnergy) : isothermalPvt_(isothermalPvt) , gasvisctCurves_(gasvisctCurves) + , viscrefPress_(viscrefPress) + , viscRef_(viscRef) , gasdentRefTemp_(gasdentRefTemp) , gasdentCT1_(gasdentCT1) , gasdentCT2_(gasdentCT2) @@ -105,6 +110,8 @@ public: void setNumRegions(size_t numRegions) { gasvisctCurves_.resize(numRegions); + viscrefPress_.resize(numRegions); + viscRef_.resize(numRegions); internalEnergyCurves_.resize(numRegions); gasdentRefTemp_.resize(numRegions); gasdentCT1_.resize(numRegions); @@ -120,9 +127,6 @@ public: void initEnd() { } - size_t numRegions() const - { return gasvisctCurves_.size(); } - /*! * \brief Returns true iff the density of the gas phase is temperature dependent. */ @@ -141,6 +145,10 @@ public: bool enableThermalViscosity() const { return enableThermalViscosity_; } + size_t numRegions() const + { return viscrefPress_.size(); } + + /*! * \brief Returns the specific internal energy [J/kg] of gas given a set of parameters. */ @@ -216,12 +224,14 @@ public: const Evaluation& Rv, const Evaluation& Rvw) const { + + const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw); if (!enableThermalViscosity()) - return isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw); + return isothermalMu; - // compute the viscosity deviation due to temperature - const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature); - return muGasvisct; + // compute the viscosity deviation due to temperature + const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true); + return muGasvisct/viscRef_[regionIdx]*isothermalMu; } /*! @@ -232,12 +242,13 @@ public: const Evaluation& temperature, const Evaluation& pressure) const { + const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure); if (!enableThermalViscosity()) - return isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure); + return isothermalMu; // compute the viscosity deviation due to temperature const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true); - return muGasvisct; + return muGasvisct/viscRef_[regionIdx]*isothermalMu; } /*! @@ -379,6 +390,12 @@ public: const std::vector& gasvisctCurves() const { return gasvisctCurves_; } + const std::vector& viscrefPress() const + { return viscrefPress_; } + + const std::vector& viscRef() const + { return viscRef_; } + const std::vector& gasdentRefTemp() const { return gasdentRefTemp_; } @@ -409,6 +426,8 @@ public: return false; return this->gasvisctCurves() == data.gasvisctCurves() && + this->viscrefPress() == data.viscrefPress() && + this->viscRef() == data.viscRef() && this->gasdentRefTemp() == data.gasdentRefTemp() && this->gasdentCT1() == data.gasdentCT1() && this->gasdentCT2() == data.gasdentCT2() && @@ -428,6 +447,8 @@ public: else isothermalPvt_ = nullptr; gasvisctCurves_ = data.gasvisctCurves_; + viscrefPress_ = data.viscrefPress_; + viscRef_ = data.viscRef_; gasdentRefTemp_ = data.gasdentRefTemp_; gasdentCT1_ = data.gasdentCT1_; gasdentCT2_ = data.gasdentCT2_; @@ -448,6 +469,8 @@ private: // The PVT properties needed for temperature dependence of the viscosity. We need // to store one value per PVT region. std::vector gasvisctCurves_; + std::vector viscrefPress_; + std::vector viscRef_; std::vector gasdentRefTemp_; std::vector gasdentCT1_; diff --git a/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp index a6843fcff..31bd742c4 100644 --- a/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp @@ -55,6 +55,7 @@ public: WaterPvtThermal() { enableThermalDensity_ = false; + enableJouleThomson_ = false; enableThermalViscosity_ = false; enableInternalEnergy_ = false; isothermalPvt_ = nullptr; diff --git a/src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp b/src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp index 230072dac..71624ccaa 100644 --- a/src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp +++ b/src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp @@ -581,6 +581,7 @@ std::optional make_jfunc(const Deck& deck) { initSimpleTableContainer(deck, "SPECHEAT", m_tabdims.getNumPVTTables()); initSimpleTableContainer(deck, "SPECROCK", m_tabdims.getNumSatTables()); initSimpleTableContainer(deck, "OILVISCT", m_tabdims.getNumPVTTables()); + initSimpleTableContainer(deck, "GASVISCT", m_tabdims.getNumPVTTables()); initSimpleTableContainer(deck, "WATVISCT", m_tabdims.getNumPVTTables()); initSimpleTableContainer(deck, "PLYADS", m_tabdims.getNumSatTables()); @@ -592,7 +593,6 @@ std::optional make_jfunc(const Deck& deck) { initPlyrockTables(deck); initPlymaxTables(deck); - initGasvisctTables(deck); initRTempTables(deck); initRocktabTables(deck); initPlyshlogTables(deck); @@ -618,34 +618,7 @@ std::optional make_jfunc(const Deck& deck) { } else if (hasRTEMPVD) { initSimpleTableContainer(deck, "RTEMPVD", "RTEMPVD", m_eqldims.getNumEquilRegions()); } - } - - - void TableManager::initGasvisctTables(const Deck& deck) { - - const std::string keywordName = "GASVISCT"; - size_t numTables = m_tabdims.getNumPVTTables(); - - if (!deck.hasKeyword(keywordName)) - return; // the table is not featured by the deck... - - auto& container = forceGetTables(keywordName , numTables); - - if (deck.count(keywordName) > 1) { - complainAboutAmbiguousKeyword(deck, keywordName); - return; - } - - const auto& tableKeyword = deck[keywordName].back(); - for (size_t tableIdx = 0; tableIdx < tableKeyword.size(); ++tableIdx) { - const auto& tableRecord = tableKeyword.getRecord( tableIdx ); - const auto& dataItem = tableRecord.getItem( 0 ); - if (dataItem.data_size() > 0) { - std::shared_ptr table = std::make_shared( deck , dataItem ); - container.addTable( tableIdx , table ); - } - } - } + } void TableManager::initPlyshlogTables(const Deck& deck) { diff --git a/src/opm/input/eclipse/EclipseState/Tables/Tables.cpp b/src/opm/input/eclipse/EclipseState/Tables/Tables.cpp index 9b58faf0f..3110f58cd 100644 --- a/src/opm/input/eclipse/EclipseState/Tables/Tables.cpp +++ b/src/opm/input/eclipse/EclipseState/Tables/Tables.cpp @@ -1170,46 +1170,12 @@ WatvisctTable::getWaterViscosityColumn() const return SimpleTable::getColumn(1); } -GasvisctTable::GasvisctTable(const Deck& deck, const DeckItem& deckItem) +GasvisctTable::GasvisctTable(const DeckItem& item, const int tableID) { - int numComponents = deck.get().back().getRecord(0).getItem(0).get(0); - - auto temperatureDimension = deck.getActiveUnitSystem().getDimension("Temperature"); - auto viscosityDimension = deck.getActiveUnitSystem().getDimension("Viscosity"); m_schema.addColumn(ColumnSchema("Temperature", Table::STRICTLY_INCREASING, Table::DEFAULT_NONE)); - for (int compIdx = 0; compIdx < numComponents; ++compIdx) { - std::string columnName = "Viscosity" + std::to_string(compIdx); - m_schema.addColumn(ColumnSchema(columnName, Table::INCREASING, Table::DEFAULT_NONE)); - } - - SimpleTable::addColumns(); - - if (deckItem.data_size() % numColumns() != 0) - throw std::runtime_error("Number of columns in the data file is inconsistent " - "with the expected number for keyword GASVISCT"); - - size_t rows = deckItem.data_size() / m_schema.size(); - for (size_t columnIndex = 0; columnIndex < m_schema.size(); columnIndex++) { - auto& column = getColumn(columnIndex); - for (size_t rowIdx = 0; rowIdx < rows; rowIdx++) { - size_t deckIndex = rowIdx * m_schema.size() + columnIndex; - - if (deckItem.defaultApplied(deckIndex)) - column.addDefault("GASVISCT"); - else { - double rawValue = deckItem.get(deckIndex); - double SIValue; - - if (columnIndex == 0) - SIValue = temperatureDimension.convertRawToSi(rawValue); - else - SIValue = viscosityDimension.convertRawToSi(rawValue); - - column.addValue(SIValue, "GASVISCT"); - } - } - } + m_schema.addColumn(ColumnSchema("Viscosity", Table::RANDOM, Table::DEFAULT_NONE)); + SimpleTable::init("GASVISCT", item, tableID); } const TableColumn& @@ -1219,17 +1185,16 @@ GasvisctTable::getTemperatureColumn() const } const TableColumn& -GasvisctTable::getGasViscosityColumn(size_t compIdx) const +GasvisctTable::getGasViscosityColumn() const { - return SimpleTable::getColumn(1 + compIdx); + return SimpleTable::getColumn(1); } RtempvdTable::RtempvdTable(const DeckItem& item, const int tableID) { m_schema.addColumn(ColumnSchema("Depth", Table::STRICTLY_INCREASING, Table::DEFAULT_NONE)); m_schema.addColumn(ColumnSchema("Temperature", Table::RANDOM, Table::DEFAULT_NONE)); - - SimpleTable::init("GASVISCT", item, tableID); + SimpleTable::init("RTEMPVD", item, tableID); } const TableColumn& diff --git a/src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp b/src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp index bb1484cf0..2881a5d8f 100644 --- a/src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp +++ b/src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp @@ -61,14 +61,45 @@ initFromState(const EclipseState& eclState, const Schedule& schedule) // viscosity if (enableThermalViscosity_) { + if (tables.getViscrefTable().empty()) + OPM_THROW(std::runtime_error, "VISCREF is required when GASVISCT is present"); + const auto& gasvisctTables = tables.getGasvisctTables(); - auto gasCompIdx = tables.gas_comp_index(); - std::string gasvisctColumnName = "Viscosity" + std::to_string(static_cast(gasCompIdx)); + const auto& viscrefTable = tables.getViscrefTable(); + + if (gasvisctTables.size() != numRegions) { + OPM_THROW(std::runtime_error, + fmt::format("Tables sizes mismatch. GASVISCT: {}, NumRegions: {}\n", + gasvisctTables.size(), numRegions)); + } + if (viscrefTable.size() != numRegions) { + OPM_THROW(std::runtime_error, + fmt::format("Tables sizes mismatch. VISCREF: {}, NumRegions: {}\n", + viscrefTable.size(), numRegions)); + } for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - const auto& T = gasvisctTables[regionIdx].getColumn("Temperature").vectorCopy(); - const auto& mu = gasvisctTables[regionIdx].getColumn(gasvisctColumnName).vectorCopy(); - gasvisctCurves_[regionIdx].setXYContainers(T, mu); + const auto& TCol = gasvisctTables[regionIdx].getColumn("Temperature").vectorCopy(); + const auto& muCol = gasvisctTables[regionIdx].getColumn("Viscosity").vectorCopy(); + gasvisctCurves_[regionIdx].setXYContainers(TCol, muCol); + + viscrefPress_[regionIdx] = viscrefTable[regionIdx].reference_pressure; + + // Temperature used to calculate the reference viscosity [K]. + // The value dooes not matter as the underlying PVT object is isothermal. + constexpr const Scalar Tref = 273.15 + 20; + //TODO: For now I just assumed the default references RV and RVW = 0, + // we could add these two parameters to a new item keyword VISCREF + //or create a new keyword for gas. + constexpr const Scalar Rvref = 0.0; + constexpr const Scalar Rvwref = 0.0; + // compute the reference viscosity using the isothermal PVT object. + viscRef_[regionIdx] = + isothermalPvt_->viscosity(regionIdx, + Tref, + viscrefPress_[regionIdx], + Rvref, + Rvwref); } } diff --git a/src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp b/src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp index d59123cad..13437fa07 100644 --- a/src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp +++ b/src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp @@ -52,6 +52,7 @@ initFromState(const EclipseState& eclState, const Schedule& schedule) const auto& tables = eclState.getTableManager(); enableThermalDensity_ = tables.OilDenT().size() > 0; + enableJouleThomson_ = tables.OilJT().size() > 0; enableThermalViscosity_ = tables.hasTables("OILVISCT"); enableInternalEnergy_ = tables.hasTables("SPECHEAT");