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

Merge pull request #765 from andlaus/implement_temperature_dependent_density

Browse Source 
Implement temperature dependent density
This commit is contained in:
Atgeirr Flø Rasmussen committed 2015-03-26 15:10:13 +01:00
2 parents b199d2b61b a386621495
commit ca9c198368
12 files changed
+1097 -252
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/core/props/pvt/PvtConstCompr.hpp
+2 -75
View File
@@ -23,8 +23,6 @@
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
#include <algorithm>
@@ -47,8 +45,6 @@ namespace Opm
{
public:
PvtConstCompr()
: oilvisctTables_(0),
watvisctTables_(0)
{}
void initFromWater(Opm::DeckKeywordConstPtr pvtwKeyword)
@@ -102,9 +98,7 @@ namespace Opm
ref_B_(1, 1.0),
comp_(1, 0.0),
viscosity_(1, visc),
visc_comp_(1, 0.0),
oilvisctTables_(0),
watvisctTables_(0)
visc_comp_(1, 0.0)
{
}
@@ -125,13 +119,6 @@ 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!");
}
}
}
@@ -159,7 +146,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,
@@ -174,46 +161,6 @@ 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);
}
@@ -340,22 +287,6 @@ 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<Opm::OilvisctTable>& oilvisctTables,
DeckKeywordConstPtr viscrefKeyword)
{
oilvisctTables_ = &oilvisctTables;
viscrefKeyword_ = viscrefKeyword;
}
/// set the tables which specify the temperature dependence of the water viscosity
void setWatvisctTables(const std::vector<Opm::WatvisctTable>& watvisctTables,
DeckKeywordConstPtr viscrefKeyword)
{
watvisctTables_ = &watvisctTables;
viscrefKeyword_ = viscrefKeyword;
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
@@ -371,10 +302,6 @@ namespace Opm
std::vector<double> comp_;
std::vector<double> viscosity_;
std::vector<double> visc_comp_;
const std::vector<Opm::OilvisctTable>* oilvisctTables_;
const std::vector<Opm::WatvisctTable>* watvisctTables_;
DeckKeywordConstPtr viscrefKeyword_;
};
}
opm/core/props/pvt/PvtDead.cpp
+2 -44
View File
@@ -123,20 +123,13 @@ 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,
@@ -150,24 +143,6 @@ 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);
@@ -176,7 +151,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,
@@ -192,23 +167,6 @@ 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);
opm/core/props/pvt/PvtDead.hpp
+1 -17
View File
@@ -44,11 +44,7 @@ namespace Opm
class PvtDead : public PvtInterface
{
public:
PvtDead()
{
// by default, specify no temperature dependence of the PVT properties
oilvisctTables_ = 0;
}
PvtDead() {};
void initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables);
void initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables);
@@ -154,15 +150,6 @@ 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<Opm::OilvisctTable>& oilvisctTables,
DeckKeywordConstPtr viscrefKeyword)
{
oilvisctTables_ = &oilvisctTables;
viscrefKeyword_ = viscrefKeyword;
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
@@ -176,9 +163,6 @@ namespace Opm
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
std::vector<NonuniformTableLinear<double> > inverseBmu_;
const std::vector<Opm::OilvisctTable>* oilvisctTables_;
DeckKeywordConstPtr viscrefKeyword_;
};
}
opm/core/props/pvt/PvtDeadSpline.cpp
+3 -40
View File
@@ -36,10 +36,7 @@ namespace Opm
//-------------------------------------------------------------------------
PvtDeadSpline::PvtDeadSpline()
{
// by default, specify no temperature dependence of the PVT properties
oilvisctTables_ = 0;
}
{}
void PvtDeadSpline::initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables,
int numSamples)
@@ -123,7 +120,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,
@@ -134,23 +131,6 @@ 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);
}
@@ -158,7 +138,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,
@@ -171,23 +151,6 @@ 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);
}
opm/core/props/pvt/PvtDeadSpline.hpp
-12
View File
@@ -148,15 +148,6 @@ 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<Opm::OilvisctTable>& oilvisctTables,
DeckKeywordConstPtr viscrefKeyword)
{
oilvisctTables_ = &oilvisctTables;
viscrefKeyword_ = viscrefKeyword;
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
@@ -169,9 +160,6 @@ namespace Opm
// table per PVT region.
std::vector<UniformTableLinear<double> > b_;
std::vector<UniformTableLinear<double> > viscosity_;
const std::vector<Opm::OilvisctTable>* oilvisctTables_;
DeckKeywordConstPtr viscrefKeyword_;
};
}
opm/core/props/pvt/PvtLiveGas.cpp
-5
View File
@@ -114,9 +114,6 @@ 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...
}
}
@@ -165,8 +162,6 @@ 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...
}
}
opm/core/props/pvt/PvtLiveOil.cpp
+3 -48
View File
@@ -36,9 +36,6 @@ namespace Opm
//-------------------------------------------------------------------------
PvtLiveOil::PvtLiveOil(const std::vector<Opm::PvtoTable>& 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);
@@ -152,20 +149,13 @@ 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,
@@ -193,32 +183,14 @@ 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,
@@ -248,24 +220,6 @@ 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
}
}
}
@@ -691,4 +645,5 @@ namespace Opm
}
}
}
} // namespace Opm
opm/core/props/pvt/PvtLiveOil.hpp
-11
View File
@@ -142,14 +142,6 @@ 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<Opm::OilvisctTable>& oilvisctTables,
DeckKeywordConstPtr viscrefKeyword)
{
oilvisctTables_ = &oilvisctTables;
viscrefKeyword_ = viscrefKeyword;
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
@@ -187,9 +179,6 @@ namespace Opm
// store one table per PVT region.
std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
const std::vector<Opm::OilvisctTable>* oilvisctTables_;
DeckKeywordConstPtr viscrefKeyword_;
};
}
opm/core/props/pvt/ThermalGasPvtWrapper.hpp
+310
View File
@@ -0,0 +1,310 @@
/*
Copyright 2015 Andreas Lauser
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 3 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/>.
*/
#ifndef OPM_THERMAL_GAS_PVT_WRAPPER_HPP
#define OPM_THERMAL_GAS_PVT_WRAPPER_HPP
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
/// Class which wraps another (i.e., isothermal) PVT object into one which adds
/// temperature dependence of gas
class ThermalGasPvtWrapper : public PvtInterface
{
public:
ThermalGasPvtWrapper()
{}
/// extract the quantities needed specify the temperature dependence of the gas
/// viscosity and density from the deck
void initFromDeck(std::shared_ptr<const PvtInterface> isothermalPvt,
Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState)
{
isothermalPvt_ = isothermalPvt;
int numRegions;
if (deck->hasKeyword("PVTG"))
numRegions = eclipseState->getPvtgTables().size();
else if (deck->hasKeyword("PVDG"))
numRegions = eclipseState->getPvdgTables().size();
else
OPM_THROW(std::runtime_error, "Gas phase was not initialized using a known way");
// viscosity
if (deck->hasKeyword("GASVISCT")) {
gasvisctTables_ = &eclipseState->getGasvisctTables();
assert(gasvisctTables_->size() == numRegions);
gasCompIdx_ = deck->getKeyword("GCOMPIDX")->getRecord(0)->getItem("GAS_COMPONENT_INDEX")->getInt(0) - 1;
gasvisctColumnName_ = "Viscosity"+std::to_string(static_cast<long long>(gasCompIdx_));
}
// density
if (deck->hasKeyword("TREF")) {
tref_ = deck->getKeyword("TREF")->getRecord(0)->getItem("TEMPERATURE")->getSIDouble(0);
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const
{
if (gasvisctTables_)
// TODO: temperature dependence for viscosity depending on z
OPM_THROW(std::runtime_error,
"temperature dependent viscosity as a function of z "
"is not yet implemented!");
// compute the isothermal viscosity
isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
if (gasvisctTables_ != 0) {
for (int i = 0; i < n; ++i) {
// temperature dependence of the gas phase. this assumes that the gas
// component index has been set properly, and it also looses the
// pressure dependence of gas. (This does not make much sense, but it
// seems to be what the documentation for the GASVISCT keyword in the
// RM says.)
int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
double muGasvisct = (*gasvisctTables_)[regionIdx].evaluate(gasvisctColumnName_, T[i]);
output_mu[i] = muGasvisct;
output_dmudp[i] = 0.0;
output_dmudr[i] = 0.0;
// TODO (?): derivative of gas viscosity w.r.t. temperature.
}
}
else {
// compute the isothermal viscosity and its derivatives
isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
if (gasvisctTables_ != 0) {
for (int i = 0; i < n; ++i) {
// temperature dependence of the gas phase. this assumes that the gas
// component index has been set properly, and it also looses the
// pressure dependence of gas. (This does not make much sense, but it
// seems to be what the documentation for the GASVISCT keyword in the
// RM says.)
int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
double muGasvisct = (*gasvisctTables_)[regionIdx].evaluate(gasvisctColumnName_, T[i]);
output_mu[i] = muGasvisct;
output_dmudp[i] = 0.0;
output_dmudr[i] = 0.0;
// TODO (?): derivative of gas viscosity w.r.t. temperature.
}
}
else {
// compute the isothermal viscosity and its derivatives
isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
}
}
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const
{
// isothermal case
isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
if (tref_ > 0.0) {
// the Eclipse TD/RM do not explicitly specify the relation of the gas
// density and the temperature, but equation (69.49) (for Eclipse 2011.1)
// implies that the temperature dependence of the gas phase is rho(T, p) =
// rho(tref_, p)/tref_*T ...
for (int i = 0; i < n; ++i) {
double alpha = tref_/T[i];
output_B[i] *= alpha;
}
}
}
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const
{
isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
if (tref_ > 0.0) {
// the Eclipse TD/RM do not explicitly specify the relation of the gas
// density and the temperature, but equation (69.49) (for Eclipse 2011.1)
// implies that the temperature dependence of the gas phase is rho(T, p) =
// rho(tref_, p)/tref_*T ...
for (int i = 0; i < n; ++i) {
double alpha = tref_/T[i];
output_B[i] *= alpha;
output_dBdp[i] *= alpha;
}
}
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
if (tref_ > 0.0) {
// the Eclipse TD/RM do not explicitly specify the relation of the gas
// density and the temperature, but equation (69.49) (for Eclipse 2011.1)
// implies that the temperature dependence of the gas phase is rho(T, p) =
// rho(tref_, p)/tref_*T ...
for (int i = 0; i < n; ++i) {
double alpha = T[i]/tref_;
output_b[i] *= alpha;
output_dbdp[i] *= alpha;
output_dbdr[i] *= alpha;
}
}
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
if (tref_ > 0.0) {
// the Eclipse TD/RM do not explicitly specify the relation of the gas
// density and the temperature, but equation (69.49) (for Eclipse 2011.1)
// implies that the temperature dependence of the gas phase is rho(T, p) =
// rho(tref_, p)/tref_*T ...
for (int i = 0; i < n; ++i) {
double alpha = T[i]/tref_;
output_b[i] *= alpha;
output_dbdp[i] *= alpha;
output_dbdr[i] *= alpha;
}
}
}
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const
{
isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
}
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
}
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const
{
isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
}
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
}
private:
int getPvtRegionIndex_(const int* pvtRegionIdx, int cellIdx) const
{
if (!pvtRegionIdx)
return 0;
return pvtRegionIdx[cellIdx];
}
// the PVT propertied for the isothermal case
std::shared_ptr<const PvtInterface> isothermalPvt_;
// The PVT properties needed for temperature dependence of the viscosity. We need
// to store one value per PVT region.
const std::vector<Opm::GasvisctTable>* gasvisctTables_;
std::string gasvisctColumnName_;
int gasCompIdx_;
// The PVT properties needed for temperature dependence of the density.
double tref_;
};
}
#endif
opm/core/props/pvt/ThermalOilPvtWrapper.hpp
+383
View File
@@ -0,0 +1,383 @@
/*
Copyright 2015 Andreas Lauser
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 3 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/>.
*/
#ifndef OPM_THERMAL_OIL_PVT_WRAPPER_HPP
#define OPM_THERMAL_OIL_PVT_WRAPPER_HPP
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
/// Class which wraps another (i.e., isothermal) PVT object into one which adds
/// temperature dependence of oil
class ThermalOilPvtWrapper : public PvtInterface
{
public:
ThermalOilPvtWrapper()
{}
/// set the tables which specify the temperature dependence of the oil viscosity
void initFromDeck(std::shared_ptr<const PvtInterface> isothermalPvt,
Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState)
{
isothermalPvt_ = isothermalPvt;
int numRegions;
if (deck->hasKeyword("PVTO"))
numRegions = eclipseState->getPvtoTables().size();
else if (deck->hasKeyword("PVDO"))
numRegions = eclipseState->getPvdoTables().size();
else if (deck->hasKeyword("PVCDO"))
numRegions = deck->getKeyword("PVCDO")->size();
else
OPM_THROW(std::runtime_error, "Oil phase was not initialized using a known way");
// viscosity
if (deck->hasKeyword("VISCREF")) {
oilvisctTables_ = &eclipseState->getOilvisctTables();
Opm::DeckKeywordConstPtr viscrefKeyword = deck->getKeyword("VISCREF");
assert(oilvisctTables_->size() == numRegions);
assert(viscrefKeyword->size() == numRegions);
viscrefPress_.resize(numRegions);
viscrefRs_.resize(numRegions);
muRef_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
DeckRecordConstPtr viscrefRecord = viscrefKeyword->getRecord(regionIdx);
viscrefPress_[regionIdx] = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
viscrefRs_[regionIdx] = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0);
// temperature used to calculate the reference viscosity [K]. the
// value does not really matter if the underlying PVT object really
// is isothermal...
double Tref = 273.15 + 20;
// compute the reference viscosity using the isothermal PVT object.
double tmp1, tmp2;
isothermalPvt_->mu(1,
&regionIdx,
&viscrefPress_[regionIdx],
&Tref,
&viscrefRs_[regionIdx],
&muRef_[regionIdx],
&tmp1,
&tmp2);
}
}
// quantities required for density. note that we just always use the values
// for the first EOS. (since EOS != PVT region.)
tref_ = 0.0;
if (deck->hasKeyword("THERMEX1")) {
oilCompIdx_ = deck->getKeyword("OCOMPIDX")->getRecord(0)->getItem("OIL_COMPONENT_INDEX")->getInt(0) - 1;
// always use the values of the first EOS
tref_ = deck->getKeyword("TREF")->getRecord(0)->getItem("TEMPERATURE")->getSIDouble(oilCompIdx_);
pref_ = deck->getKeyword("PREF")->getRecord(0)->getItem("PRESSURE")->getSIDouble(oilCompIdx_);
cref_ = deck->getKeyword("CREF")->getRecord(0)->getItem("COMPRESSIBILITY")->getSIDouble(oilCompIdx_);
thermex1_ = deck->getKeyword("THERMEX1")->getRecord(0)->getItem("EXPANSION_COEFF")->getSIDouble(oilCompIdx_);
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const
{
if (oilvisctTables_)
// TODO: temperature dependence for viscosity depending on z
OPM_THROW(std::runtime_error,
"temperature dependent viscosity as a function of z "
"is not yet implemented!");
// compute the isothermal viscosity
isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// compute the isothermal viscosity and its derivatives
isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
if (!oilvisctTables_)
// isothermal case
return;
// temperature dependence
for (int i = 0; i < n; ++i) {
int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
// calculate the viscosity of the isothermal keyword for the reference
// pressure given by the VISCREF keyword.
double muRef = muRef_[regionIdx];
// compute the viscosity deviation due to temperature
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 viscosity w.r.t. temperature.
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// compute the isothermal viscosity and its derivatives
isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
if (!oilvisctTables_)
// isothermal case
return;
// temperature dependence
for (int i = 0; i < n; ++i) {
int regionIdx = getPvtRegionIndex_(pvtRegionIdx, i);
// calculate the viscosity of the isothermal keyword for the reference
// pressure given by the VISCREF keyword.
double muRef = muRef_[regionIdx];
// compute the viscosity deviation due to temperature
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 viscosity w.r.t. temperature.
}
}
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const
{
// isothermal case
isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
if (thermex1_ <= 0.0)
// isothermal case
return;
// deal with the temperature dependence of the oil phase. we use equation
// (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
// using the isothermal keyword instead of using the value for the
// components, so the oil compressibility is already dealt with there. Note
// that we only do the part for the oil component here, the part for
// dissolved gas is ignored so far.
double cT1 = thermex1_;
double TRef = tref_;
for (int i = 0; i < n; ++i) {
double alpha = (1 + cT1*(T[i] - TRef));
output_B[i] *= alpha;
}
}
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const
{
isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
if (thermex1_ <= 0.0)
// isothermal case
return;
// deal with the temperature dependence of the oil phase. we use equation
// (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
// using the isothermal keyword instead of using the value for the
// components, so the oil compressibility is already dealt with there. Note
// that we only do the part for the oil component here, the part for
// dissolved gas is ignored so far.
double cT1 = thermex1_;
double TRef = tref_;
for (int i = 0; i < n; ++i) {
double alpha = (1 + cT1*(T[i] - TRef));
output_B[i] *= alpha;
output_dBdp[i] *= alpha;
}
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
if (thermex1_ <= 0.0)
// isothermal case
return;
// deal with the temperature dependence of the oil phase. we use equation
// (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
// using the isothermal keyword instead of using the value for the
// components, so the oil compressibility is already dealt with there. Note
// that we only do the part for the oil component here, the part for
// dissolved gas is ignored so far.
double cT1 = thermex1_;
double TRef = tref_;
for (int i = 0; i < n; ++i) {
double alpha = 1.0/(1 + cT1*(T[i] - TRef));
output_b[i] *= alpha;
output_dbdp[i] *= alpha;
output_dbdr[i] *= alpha;
}
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
if (thermex1_ <= 0.0)
// isothermal case
return;
// deal with the temperature dependence of the oil phase. we use equation
// (3.208) from the Eclipse 2011.1 Reference Manual, but we calculate rho_ref
// using the isothermal keyword instead of using the value for the
// components, so the oil compressibility is already dealt with there. Note
// that we only do the part for the oil component here, the part for
// dissolved gas is ignored so far.
double cT1 = thermex1_;
double TRef = tref_;
for (int i = 0; i < n; ++i) {
double alpha = 1.0/(1 + cT1*(T[i] - TRef));
output_b[i] *= alpha;
output_dbdp[i] *= alpha;
output_dbdr[i] *= alpha;
}
}
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const
{
isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
}
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
}
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const
{
isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
}
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
}
private:
int getPvtRegionIndex_(const int* pvtRegionIdx, int cellIdx) const
{
if (!pvtRegionIdx)
return 0;
return pvtRegionIdx[cellIdx];
}
// the PVT propertied for the isothermal case
std::shared_ptr<const PvtInterface> isothermalPvt_;
// The PVT properties needed for temperature dependence of the viscosity. We need
// to store one value per PVT region.
std::vector<double> viscrefPress_;
std::vector<double> viscrefRs_;
std::vector<double> muRef_;
const std::vector<Opm::OilvisctTable>* oilvisctTables_;
// The PVT properties needed for temperature dependence of the density. This is
// specified as one value per EOS in the manual, but we unconditionally use the
// expansion coefficient of the first EOS...
int oilCompIdx_;
double tref_;
double pref_;
double cref_;
double thermex1_;
};
}
#endif
opm/core/props/pvt/ThermalWaterPvtWrapper.hpp
+390
View File
@@ -0,0 +1,390 @@
/*
Copyright 2015 Andreas Lauser
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 3 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/>.
*/
#ifndef OPM_THERMAL_WATER_PVT_WRAPPER_HPP
#define OPM_THERMAL_WATER_PVT_WRAPPER_HPP
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
/// Class which wraps another (i.e., isothermal) PVT object into one which adds
/// temperature dependence of water
class ThermalWaterPvtWrapper : public PvtInterface
{
public:
ThermalWaterPvtWrapper()
{}
/// set the tables which specify the temperature dependence of the water viscosity
void initFromDeck(std::shared_ptr<const PvtInterface> isothermalPvt,
Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState)
{
isothermalPvt_ = isothermalPvt;
watvisctTables_ = 0;
// stuff which we need to get from the PVTW keyword
Opm::DeckKeywordConstPtr pvtwKeyword = deck->getKeyword("PVTW");
int numRegions = pvtwKeyword->size();
pvtwRefPress_.resize(numRegions);
pvtwRefB_.resize(numRegions);
pvtwCompressibility_.resize(numRegions);
pvtwViscosity_.resize(numRegions);
pvtwViscosibility_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
Opm::DeckRecordConstPtr pvtwRecord = pvtwKeyword->getRecord(regionIdx);
pvtwRefPress_[regionIdx] = pvtwRecord->getItem("P_REF")->getSIDouble(0);
pvtwRefB_[regionIdx] = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0);
pvtwViscosity_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0);
pvtwViscosibility_[regionIdx] = pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0);
}
// quantities required for the temperature dependence of the viscosity
// (basically we expect well-behaved VISCREF and WATVISCT keywords.)
if (deck->hasKeyword("VISCREF")) {
watvisctTables_ = &eclipseState->getWatvisctTables();
Opm::DeckKeywordConstPtr viscrefKeyword = deck->getKeyword("VISCREF");
assert(watvisctTables_->size() == numRegions);
assert(viscrefKeyword->size() == numRegions);
viscrefPress_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
Opm::DeckRecordConstPtr viscrefRecord = viscrefKeyword->getRecord(regionIdx);
viscrefPress_[regionIdx] = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
}
}
// quantities required for the temperature dependence of the density
if (deck->hasKeyword("WATDENT")) {
DeckKeywordConstPtr watdentKeyword = deck->getKeyword("WATDENT");
assert(watdentKeyword->size() == numRegions);
watdentRefTemp_.resize(numRegions);
watdentCT1_.resize(numRegions);
watdentCT2_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
Opm::DeckRecordConstPtr watdentRecord = watdentKeyword->getRecord(regionIdx);
watdentRefTemp_[regionIdx] = watdentRecord->getItem("REFERENCE_TEMPERATURE")->getSIDouble(0);
watdentCT1_[regionIdx] = watdentRecord->getItem("EXPANSION_COEFF_LINEAR")->getSIDouble(0);
watdentCT2_[regionIdx] = watdentRecord->getItem("EXPANSION_COEFF_QUADRATIC")->getSIDouble(0);
}
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const
{
if (watvisctTables_)
// TODO: temperature dependence for viscosity depending on z
OPM_THROW(std::runtime_error,
"temperature dependent viscosity as a function of z "
"is not yet implemented!");
// compute the isothermal viscosity
isothermalPvt_->mu(n, pvtRegionIdx, p, T, z, output_mu);
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// compute the isothermal viscosity and its derivatives
isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, output_mu, output_dmudp, output_dmudr);
if (!watvisctTables_)
// isothermal case
return;
// temperature dependence
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
// calculate the viscosity of the isothermal keyword for the reference
// pressure given by the VISCREF keyword.
double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]);
double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
// compute the viscosity deviation due to temperature
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 viscosity w.r.t. temperature.
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// compute the isothermal viscosity and its derivatives
isothermalPvt_->mu(n, pvtRegionIdx, p, T, r, cond, output_mu, output_dmudp, output_dmudr);
if (!watvisctTables_)
// isothermal case
return;
// temperature dependence
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
// calculate the viscosity of the isothermal keyword for the reference
// pressure given by the VISCREF keyword.
double x = -pvtwViscosibility_[tableIdx]*(viscrefPress_[tableIdx] - pvtwRefPress_[tableIdx]);
double muRef = pvtwViscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
// compute the viscosity deviation due to temperature
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 viscosity w.r.t. temperature.
}
}
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const
{
if (watdentRefTemp_.empty()) {
// isothermal case
isothermalPvt_->B(n, pvtRegionIdx, p, T, z, output_B);
return;
}
// This changes how the water density depends on pressure compared to what's
// used for the PVTW keyword, but it seems to be what Eclipse does. For
// details, see the documentation for the WATDENT keyword in the Eclipse RM.
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double BwRef = pvtwRefB_[tableIdx];
double TRef = watdentRefTemp_[tableIdx];
double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
double cT1 = watdentCT1_[tableIdx];
double cT2 = watdentCT2_[tableIdx];
double Y = T[i] - TRef;
double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
output_B[i] = Bw;
}
}
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const
{
if (watdentRefTemp_.empty()) {
// isothermal case
isothermalPvt_->dBdp(n, pvtRegionIdx, p, T, z, output_B, output_dBdp);
return;
}
// This changes how the water density depends on pressure. This is awkward,
// but it seems to be what Eclipse does. See the documentation for the
// WATDENT keyword in the Eclipse RM
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double BwRef = pvtwRefB_[tableIdx];
double TRef = watdentRefTemp_[tableIdx];
double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
double cT1 = watdentCT1_[tableIdx];
double cT2 = watdentCT2_[tableIdx];
double Y = T[i] - TRef;
double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
output_B[i] = Bw;
}
std::fill(output_dBdp, output_dBdp + n, 0.0);
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
if (watdentRefTemp_.empty()) {
// isothermal case
isothermalPvt_->b(n, pvtRegionIdx, p, T, r, output_b, output_dbdp, output_dbdr);
return;
}
// This changes how the water density depends on pressure. This is awkward,
// but it seems to be what Eclipse does. See the documentation for the
// WATDENT keyword in the Eclipse RM
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double BwRef = pvtwRefB_[tableIdx];
double TRef = watdentRefTemp_[tableIdx];
double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
double cT1 = watdentCT1_[tableIdx];
double cT2 = watdentCT2_[tableIdx];
double Y = T[i] - TRef;
double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
output_b[i] = 1.0/Bw;
}
std::fill(output_dbdp, output_dbdp + n, 0.0);
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
if (watdentRefTemp_.empty()) {
// isothermal case
isothermalPvt_->b(n, pvtRegionIdx, p, T, r, cond, output_b, output_dbdp, output_dbdr);
return;
}
// This changes pressure dependence of the water density, but it seems to be
// what Eclipse does. See the documentation for the WATDENT keyword in the
// Eclipse RM
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double BwRef = pvtwRefB_[tableIdx];
double TRef = watdentRefTemp_[tableIdx];
double X = pvtwCompressibility_[tableIdx]*(p[i] - pvtwRefPress_[tableIdx]);
double cT1 = watdentCT1_[tableIdx];
double cT2 = watdentCT2_[tableIdx];
double Y = T[i] - TRef;
double Bw = BwRef*(1 - X)*(1 + cT1*Y + cT2*Y*Y);
output_b[i] = 1.0/Bw;
}
std::fill(output_dbdp, output_dbdp + n, 0.0);
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const
{
isothermalPvt_->rsSat(n, pvtRegionIdx, p, output_rsSat, output_drsSatdp);
}
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
isothermalPvt_->rvSat(n, pvtRegionIdx, p, output_rvSat, output_drvSatdp);
}
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const
{
isothermalPvt_->R(n, pvtRegionIdx, p, z, output_R);
}
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
isothermalPvt_->dRdp(n, pvtRegionIdx, p, z, output_R, output_dRdp);
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// the PVT propertied for the isothermal case
std::shared_ptr<const PvtInterface> isothermalPvt_;
// The PVT properties needed for temperature dependence. We need to store one
// value per PVT region.
std::vector<double> viscrefPress_;
std::vector<double> watdentRefTemp_;
std::vector<double> watdentCT1_;
std::vector<double> watdentCT2_;
std::vector<double> pvtwRefPress_;
std::vector<double> pvtwRefB_;
std::vector<double> pvtwCompressibility_;
std::vector<double> pvtwViscosity_;
std::vector<double> pvtwViscosibility_;
const std::vector<Opm::WatvisctTable>* watvisctTables_;
};
}
#endif