PVT properties: allow them to be temperature dependent

Note that this patch does not introduce any real temperature
dependence but only changes the APIs for the viscosity and for the
density related methods. Note that I also don't like the fact that
this requires so many changes to so many files, but with the current
design of the property classes I cannot see a way to avoid this...
This commit is contained in:
Andreas Lauser 2014-11-20 12:15:01 +01:00
parent b6e5dcd706
commit 18c641e0c7
3 changed files with 100 additions and 42 deletions

View File

@ -122,12 +122,15 @@ namespace Opm
*
* \param[in] p Fluid pressure.
*
* \param[in] T Temperature.
*
* \param[in] z Surface volumes of all phases.
*
* \return Phase densities at phase point.
*/
std::vector<double>
operator()(const double p,
const double T,
const std::vector<double>& z) const
{
const int np = props_.numPhases();
@ -136,7 +139,7 @@ namespace Opm
assert (z.size() == std::vector<double>::size_type(np));
double* dAdp = 0;
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
props_.matrix(1, &p, &T, &z[0], &c_[0], &A[0], dAdp);
std::vector<double> rho(np, 0.0);
props_.density(1, &A[0], &c_[0], &rho[0]);
@ -171,11 +174,15 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
virtual double operator()(const double depth,
const double press,
const double temp,
const double sat = 0.0) const = 0;
};
@ -194,6 +201,9 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. In "no mixing
* policy", this is identically zero.
@ -201,6 +211,7 @@ namespace Opm
double
operator()(const double /* depth */,
const double /* press */,
const double /* temp */,
const double /* sat */ = 0.0) const
{
return 0.0;
@ -247,18 +258,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double depth,
const double press,
const double temp,
const double sat_gas = 0.0) const
{
if (sat_gas > 0.0) {
return satRs(press);
return satRs(press, temp);
} else {
return std::min(satRs(press), linearInterpolation(depth_, rs_, depth));
return std::min(satRs(press, temp), linearInterpolation(depth_, rs_, depth));
}
}
@ -270,9 +285,9 @@ namespace Opm
double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press) const
double satRs(const double press, const double temp) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
@ -323,18 +338,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RV
* value.
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double depth,
const double press,
const double temp,
const double sat_oil = 0.0 ) const
{
if (sat_oil > 0.0) {
return satRv(press);
return satRv(press, temp);
} else {
return std::min(satRv(press), linearInterpolation(depth_, rv_, depth));
return std::min(satRv(press, temp), linearInterpolation(depth_, rv_, depth));
}
}
@ -346,9 +365,9 @@ namespace Opm
double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press) const
double satRv(const double press, const double temp) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.
@ -382,15 +401,16 @@ namespace Opm
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
: props_(props), cell_(cell)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
rs_sat_contact_ = satRs(p_contact);
rs_sat_contact_ = satRs(p_contact, T_contact);
}
/**
@ -402,18 +422,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RS
* value.
*
* \param[in] temp Temperature at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double /* depth */,
const double press,
const double temp,
const double sat_gas = 0.0) const
{
if (sat_gas > 0.0) {
return satRs(press);
return satRs(press, temp);
} else {
return std::min(satRs(press), rs_sat_contact_);
return std::min(satRs(press, temp), rs_sat_contact_);
}
}
@ -424,9 +448,9 @@ namespace Opm
double rs_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press) const
double satRs(const double press, const double temp) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
@ -460,15 +484,16 @@ namespace Opm
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
: props_(props), cell_(cell)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
rv_sat_contact_ = satRv(p_contact);
rv_sat_contact_ = satRv(p_contact, T_contact);
}
/**
@ -480,18 +505,22 @@ namespace Opm
* \param[in] press Pressure at which to calculate RV
* value.
*
* \param[in] temp Temperature at which to calculate RV
* value.
*
* \return Dissolved oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double /*depth*/,
const double press,
const double temp,
const double sat_oil = 0.0) const
{
if (sat_oil > 0.0) {
return satRv(press);
return satRv(press, temp);
} else {
return std::min(satRv(press), rv_sat_contact_);
return std::min(satRv(press, temp), rv_sat_contact_);
}
}
@ -502,9 +531,9 @@ namespace Opm
double rv_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press) const
double satRv(const double press, const double temp) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.

View File

@ -171,6 +171,7 @@ namespace Opm
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] oil_pressure Oil pressure for each cell in range.
* \param[in] temperature Temperature for each cell in range.
* \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range.
*/
@ -178,6 +179,7 @@ namespace Opm
std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells,
const std::vector<double> oil_pressure,
const std::vector<double>& temperature,
const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation);
@ -298,7 +300,8 @@ namespace Opm
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press;
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
const double T_contact = 273.15 + 20; // standard temperature for now
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact, T_contact));
}
}
} else {
@ -329,7 +332,8 @@ namespace Opm
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press + rec[i].goc.press;
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
const double T_contact = 273.15 + 20; // standard temperature for now
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact, T_contact));
}
}
} else {
@ -399,6 +403,7 @@ namespace Opm
props.phaseUsage());
PVec press = phasePressures(G, eqreg, cells, grav);
const std::vector<double>& temp = temperature(G, eqreg, cells);
const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, press);
@ -411,8 +416,8 @@ namespace Opm
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
const Vec rs = computeRs(G, cells, press[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_);
}

View File

@ -106,11 +106,13 @@ namespace Opm
template <class Density>
class Water {
public:
Water(const Density& rho,
Water(const double temp,
const Density& rho,
const int np,
const int ix,
const double norm_grav)
: rho_(rho)
: temp_(temp)
, rho_(rho)
, svol_(np, 0)
, ix_(ix)
, g_(norm_grav)
@ -126,6 +128,7 @@ namespace Opm
}
private:
const double temp_;
const Density& rho_;
std::vector<double> svol_;
const int ix_;
@ -134,7 +137,7 @@ namespace Opm
double
density(const double press) const
{
const std::vector<double>& rho = rho_(press, svol_);
const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[ix_];
}
@ -143,13 +146,15 @@ namespace Opm
template <class Density, class RS>
class Oil {
public:
Oil(const Density& rho,
Oil(const double temp,
const Density& rho,
const RS& rs,
const int np,
const int oix,
const int gix,
const double norm_grav)
: rho_(rho)
: temp_(temp)
, rho_(rho)
, rs_(rs)
, svol_(np, 0)
, oix_(oix)
@ -167,6 +172,7 @@ namespace Opm
}
private:
const double temp_;
const Density& rho_;
const RS& rs_;
mutable std::vector<double> svol_;
@ -179,10 +185,10 @@ namespace Opm
const double press) const
{
if (gix_ >= 0) {
svol_[gix_] = rs_(depth, press);
svol_[gix_] = rs_(depth, press, temp_);
}
const std::vector<double>& rho = rho_(press, svol_);
const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[oix_];
}
};
@ -190,13 +196,15 @@ namespace Opm
template <class Density, class RV>
class Gas {
public:
Gas(const Density& rho,
Gas(const double temp,
const Density& rho,
const RV& rv,
const int np,
const int gix,
const int oix,
const double norm_grav)
: rho_(rho)
: temp_(temp)
, rho_(rho)
, rv_(rv)
, svol_(np, 0)
, gix_(gix)
@ -214,6 +222,7 @@ namespace Opm
}
private:
const double temp_;
const Density& rho_;
const RV& rv_;
mutable std::vector<double> svol_;
@ -226,10 +235,10 @@ namespace Opm
const double press) const
{
if (oix_ >= 0) {
svol_[oix_] = rv_(depth, press);
svol_[oix_] = rv_(depth, press, temp_);
}
const std::vector<double>& rho = rho_(press, svol_);
const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[gix_];
}
};
@ -333,7 +342,8 @@ namespace Opm
const PhaseUsage& pu = reg.phaseUsage();
const int wix = PhaseIndex::water(pu);
ODE drho(reg.densityCalculator(), pu.num_phases, wix, grav);
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
const double z0 = reg.zwoc();
const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw
@ -373,7 +383,9 @@ namespace Opm
const int oix = PhaseIndex::oil(pu);
const int gix = PhaseIndex::gas(pu);
ODE drho(reg.densityCalculator(),
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(),
reg.dissolutionCalculator(),
pu.num_phases, oix, gix, grav);
@ -426,7 +438,9 @@ namespace Opm
const int gix = PhaseIndex::gas(pu);
const int oix = PhaseIndex::oil(pu);
ODE drho(reg.densityCalculator(),
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(),
reg.evaporationCalculator(),
pu.num_phases, gix, oix, grav);
@ -573,8 +587,16 @@ namespace Opm
return press;
}
template <class Region,
class CellRange>
std::vector<double>
temperature(const UnstructuredGrid& G,
const Region& reg,
const CellRange& cells)
{
// use the standard temperature for everything for now
return std::vector<double>(cells.size(), 273.15 + 20.0);
}
template <class Region, class CellRange>
std::vector< std::vector<double> >
@ -716,6 +738,7 @@ namespace Opm
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] oil_pressure Oil pressure for each cell in range.
* \param[in] temperature Temperature for each cell in range.
* \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range.
*/
@ -723,6 +746,7 @@ namespace Opm
std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells,
const std::vector<double> oil_pressure,
const std::vector<double>& temperature,
const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation)
{
@ -731,7 +755,7 @@ namespace Opm
int count = 0;
for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
const double depth = grid.cell_centroids[3*(*it) + 2];
rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]);
rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
}
return rs;
}