diff --git a/opm/core/props/pvt/PvtConstCompr.hpp b/opm/core/props/pvt/PvtConstCompr.hpp index 46e68f8f..542e2093 100644 --- a/opm/core/props/pvt/PvtConstCompr.hpp +++ b/opm/core/props/pvt/PvtConstCompr.hpp @@ -23,6 +23,8 @@ #include #include +#include + #include #include @@ -98,7 +100,9 @@ namespace Opm ref_B_(1, 1.0), comp_(1, 0.0), viscosity_(1, visc), - visc_comp_(1, 0.0) + visc_comp_(1, 0.0), + oilvisctTables_(0), + watvisctTables_(0) { } @@ -119,6 +123,13 @@ namespace Opm int tableIdx = getTableIndex_(pvtRegionIdx, i); double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x); + + if (oilvisctTables_ != 0 || watvisctTables_ != 0) { + // TODO: temperature dependence + OPM_THROW(std::logic_error, + "temperature dependent viscosity as a function of z " + "is not yet implemented!"); + } } } @@ -146,7 +157,7 @@ namespace Opm virtual void mu(const int n, const int* pvtRegionIdx, const double* p, - const double* /*T*/, + const double* T, const double* /*r*/, const PhasePresence* /*cond*/, double* output_mu, @@ -161,6 +172,46 @@ namespace Opm double d = (1.0 + x + 0.5*x*x); output_mu[i] = viscosity_[tableIdx]/d; output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx]; + + if (oilvisctTables_ != 0) { + // this object handles _either_ oil _or_ water + assert(watvisctTables_ == 0); + + // temperature dependence of the oil phase + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double muRef = -visc_comp_[tableIdx]*(pRef - ref_press_[tableIdx]); + + double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } + + if (watvisctTables_ != 0) { + // this object handles _either_ oil _or_ water + assert(oilvisctTables_ == 0); + + // temperature dependence of the water phase + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double muRef = -visc_comp_[tableIdx]*(pRef - ref_press_[tableIdx]); + + double muWatvisct = (*watvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); + double alpha = muWatvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } } std::fill(output_dmudr, output_dmudr + n, 0.0); } @@ -287,6 +338,22 @@ namespace Opm std::fill(output_dRdp, output_dRdp + n, 0.0); } + /// set the tables which specify the temperature dependence of the oil viscosity + void setOilvisctTables(const std::vector& oilvisctTables, + DeckKeywordConstPtr viscrefKeyword) + { + oilvisctTables_ = &oilvisctTables; + viscrefKeyword_ = viscrefKeyword; + } + + /// set the tables which specify the temperature dependence of the water viscosity + void setWatvisctTables(const std::vector& watvisctTables, + DeckKeywordConstPtr viscrefKeyword) + { + watvisctTables_ = &watvisctTables; + viscrefKeyword_ = viscrefKeyword; + } + private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -302,6 +369,10 @@ namespace Opm std::vector comp_; std::vector viscosity_; std::vector visc_comp_; + + const std::vector* oilvisctTables_; + const std::vector* watvisctTables_; + DeckKeywordConstPtr viscrefKeyword_; }; } diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp index 67ed9335..3d711ca6 100644 --- a/opm/core/props/pvt/PvtDead.cpp +++ b/opm/core/props/pvt/PvtDead.cpp @@ -123,13 +123,20 @@ namespace Opm double tempInvB = b_[regionIdx](p[i]); double tempInvBmu = inverseBmu_[regionIdx](p[i]); output_mu[i] = tempInvB / tempInvBmu; + + if (oilvisctTables_ != 0) { + // TODO: temperature dependence + OPM_THROW(std::logic_error, + "temperature dependent viscosity as a function of z " + "is not yet implemented!"); + } } } void PvtDead::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* /*r*/, double* output_mu, double* output_dmudp, @@ -143,6 +150,24 @@ namespace Opm output_mu[i] = tempInvB / tempInvBmu; output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); + + if (oilvisctTables_ != 0) { + // temperature dependence of the oil phase + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double muRef = b_[regionIdx](pRef)/inverseBmu_[regionIdx](pRef); + + double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } + } std::fill(output_dmudr, output_dmudr + n, 0.0); @@ -151,7 +176,7 @@ namespace Opm void PvtDead::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* /*r*/, const PhasePresence* /*cond*/, double* output_mu, @@ -167,6 +192,23 @@ namespace Opm output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); + + if (oilvisctTables_ != 0) { + // temperature dependence of the oil phase + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double muRef = b_[regionIdx](pRef)/inverseBmu_[regionIdx](pRef); + + double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } } std::fill(output_dmudr, output_dmudr + n, 0.0); diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp index b6f04c2c..0bb20ea6 100644 --- a/opm/core/props/pvt/PvtDead.hpp +++ b/opm/core/props/pvt/PvtDead.hpp @@ -44,7 +44,11 @@ namespace Opm class PvtDead : public PvtInterface { public: - PvtDead() {}; + PvtDead() + { + // by default, specify no temperature dependence of the PVT properties + oilvisctTables_ = 0; + } void initFromOil(const std::vector& pvdoTables); void initFromGas(const std::vector& pvdgTables); @@ -150,6 +154,15 @@ namespace Opm const double* z, double* output_R, double* output_dRdp) const; + + /// set the tables which specify the temperature dependence of the oil viscosity + void setOilvisctTables(const std::vector& oilvisctTables, + DeckKeywordConstPtr viscrefKeyword) + { + oilvisctTables_ = &oilvisctTables; + viscrefKeyword_ = viscrefKeyword; + } + private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -163,6 +176,9 @@ namespace Opm std::vector > b_; std::vector > viscosity_; std::vector > inverseBmu_; + + const std::vector* oilvisctTables_; + DeckKeywordConstPtr viscrefKeyword_; }; } diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp index 293caf62..1e0b37e9 100644 --- a/opm/core/props/pvt/PvtDeadSpline.cpp +++ b/opm/core/props/pvt/PvtDeadSpline.cpp @@ -36,7 +36,10 @@ namespace Opm //------------------------------------------------------------------------- PvtDeadSpline::PvtDeadSpline() - {} + { + // by default, specify no temperature dependence of the PVT properties + oilvisctTables_ = 0; + } void PvtDeadSpline::initFromOil(const std::vector& pvdoTables, int numSamples) @@ -120,7 +123,7 @@ namespace Opm void PvtDeadSpline::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* /*r*/, double* output_mu, double* output_dmudp, @@ -131,6 +134,23 @@ namespace Opm int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + + if (oilvisctTables_ != 0) { + // temperature dependence of the oil phase + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double muRef = viscosity_[regionIdx](pRef); + + double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } } std::fill(output_dmudr, output_dmudr + n, 0.0); } @@ -138,7 +158,7 @@ namespace Opm void PvtDeadSpline::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* /*r*/, const PhasePresence* /*cond*/, double* output_mu, @@ -151,6 +171,23 @@ namespace Opm int regionIdx = getTableIndex_(pvtTableIdx, i); output_mu[i] = viscosity_[regionIdx](p[i]); output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); + + if (oilvisctTables_ != 0) { + // temperature dependence of the oil phase + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double muRef = viscosity_[regionIdx](pRef); + + double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } } std::fill(output_dmudr, output_dmudr + n, 0.0); } diff --git a/opm/core/props/pvt/PvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp index 3f60f9f3..366e0caf 100644 --- a/opm/core/props/pvt/PvtDeadSpline.hpp +++ b/opm/core/props/pvt/PvtDeadSpline.hpp @@ -148,6 +148,15 @@ namespace Opm const double* z, double* output_R, double* output_dRdp) const; + + /// set the tables which specify the temperature dependence of the oil viscosity + void setOilvisctTables(const std::vector& oilvisctTables, + DeckKeywordConstPtr viscrefKeyword) + { + oilvisctTables_ = &oilvisctTables; + viscrefKeyword_ = viscrefKeyword; + } + private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -160,6 +169,9 @@ namespace Opm // table per PVT region. std::vector > b_; std::vector > viscosity_; + + const std::vector* oilvisctTables_; + DeckKeywordConstPtr viscrefKeyword_; }; } diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp index ce189c8b..cf60b836 100644 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ b/opm/core/props/pvt/PvtLiveGas.cpp @@ -114,6 +114,9 @@ namespace Opm double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false); output_mu[i] = inverseB / inverseBMu; + + // temperature dependence: Since E100 does not implement temperature + // dependence of gas viscosity, we skip it here as well... } } @@ -162,6 +165,8 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); + // temperature dependence: Since E100 does not implement temperature + // dependence of gas viscosity, we skip it here as well... } } diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp index bf78e48f..ebefdc4f 100644 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ b/opm/core/props/pvt/PvtLiveOil.cpp @@ -36,6 +36,9 @@ namespace Opm //------------------------------------------------------------------------- PvtLiveOil::PvtLiveOil(const std::vector& pvtoTables) { + // by default, specify no temperature dependence of the PVT properties + oilvisctTables_ = 0; + int numTables = pvtoTables.size(); saturated_oil_table_.resize(numTables); undersat_oil_tables_.resize(numTables); @@ -136,7 +139,7 @@ namespace Opm void PvtLiveOil::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* z, double* output_mu) const { @@ -149,13 +152,20 @@ namespace Opm output_mu[i] = inverseB / inverseBMu; } + + if (oilvisctTables_ != 0) { + // TODO: temperature dependence + OPM_THROW(std::logic_error, + "temperature dependent viscosity as a function of z " + "is not yet implemented!"); + } } /// Viscosity and its p and r derivatives as a function of p, T and r. void PvtLiveOil::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* r, double* output_mu, double* output_dmudp, @@ -183,14 +193,32 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); + if (oilvisctTables_ != 0) { + // temperature dependence + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double rsRef = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0); + double muRef = miscible_oil(pRef, rsRef, tableIdx, 3, 0)/miscible_oil(pRef, rsRef, tableIdx, 1, 0); + + double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } } + } /// Viscosity and its p and r derivatives as a function of p, T and r. void PvtLiveOil::mu(const int n, const int* pvtTableIdx, const double* p, - const double* /*T*/, + const double* T, const double* r, const PhasePresence* cond, double* output_mu, @@ -220,6 +248,24 @@ namespace Opm output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) / (inverseBMu * inverseBMu); + + if (oilvisctTables_ != 0) { + // temperature dependence + DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx); + double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0); + double rsRef = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0); + double muRef = miscible_oil(pRef, rsRef, tableIdx, 3, 0)/miscible_oil(pRef, rsRef, tableIdx, 1, 0); + + double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]); + double alpha = muOilvisct/muRef; + + output_mu[i] *= alpha; + output_dmudp[i] *= alpha; + output_dmudr[i] *= alpha; + + // TODO (?): derivative of oil viscosity w.r.t. temperature. + // probably requires a healthy portion of if-spaghetti + } } } @@ -645,5 +691,4 @@ namespace Opm } } } - } // namespace Opm diff --git a/opm/core/props/pvt/PvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp index e94dd7c1..34191875 100644 --- a/opm/core/props/pvt/PvtLiveOil.hpp +++ b/opm/core/props/pvt/PvtLiveOil.hpp @@ -142,6 +142,14 @@ namespace Opm double* output_R, double* output_dRdp) const; + /// set the tables which specify the temperature dependence of the oil viscosity + void setOilvisctTables(const std::vector& oilvisctTables, + DeckKeywordConstPtr viscrefKeyword) + { + oilvisctTables_ = &oilvisctTables; + viscrefKeyword_ = viscrefKeyword; + } + private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { @@ -179,6 +187,9 @@ namespace Opm // store one table per PVT region. std::vector > > saturated_oil_table_; std::vector > > > undersat_oil_tables_; + + const std::vector* oilvisctTables_; + DeckKeywordConstPtr viscrefKeyword_; }; }