Added functions for b and mu with isSaturated as input.

Functions for volume factor and viscosity that explicitly take a boolean
variable indicating whether the fluid is saturated or not is added to
the SinglePvtInterface.

Corresponding changes are done in the dependent PVT files.

The new functionality is tested in test_blackoilfluid
This commit is contained in:
Tor Harald Sandve 2013-11-28 10:29:24 +01:00
parent 78301dbc6d
commit a5a4d7ba39
11 changed files with 389 additions and 2 deletions

View File

@ -108,6 +108,30 @@ namespace Opm
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void mu(const int n,
const double* p,
const double* /*r*/,
const bool* /*isSat*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
if (visc_comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = -visc_comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_/d;
output_dmudp[i] = (viscosity_/(d*d))*(1+x) * visc_comp_;
}
} else {
std::fill(output_mu, output_mu + n, viscosity_);
std::fill(output_dmudp, output_dmudp + n, 0.0);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void B(const int n,
const double* p,
const double* /*z*/,
@ -172,6 +196,34 @@ namespace Opm
}
virtual void b(const int n,
const double* p,
const double* /*r*/,
const bool* /*isSat*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
if (comp_) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(p[i] - ref_press_);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_B;
output_b[i] = d/ref_B_;
output_dbdp[i] = (1 + x) * comp_/ref_B_;
}
} else {
std::fill(output_b, output_b + n, 1/ref_B_);
std::fill(output_dbdp, output_dbdp + n, 0.0);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void rbub(const int n,
const double* /*p*/,
double* output_rbub,

View File

@ -96,6 +96,23 @@ namespace Opm
}
void SinglePvtDead::mu(const int n,
const double* p,
const double* /*r*/,
const bool* /*isSat*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
output_dmudp[i] = viscosity_.derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void SinglePvtDead::B(const int n,
const double* p,
const double* /*z*/,
@ -140,6 +157,25 @@ namespace Opm
}
void SinglePvtDead::b(const int n,
const double* p,
const double* /*r*/,
const bool* /*isSat*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = b_(p[i]);
output_dbdp[i] = b_.derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void SinglePvtDead::rbub(const int n,
const double* /*p*/,
double* output_rbub,

View File

@ -49,6 +49,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void mu(const int n,
const double* p,
const double* r,
@ -56,6 +57,16 @@ namespace Opm
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const double* p,
@ -70,6 +81,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void b(const int n,
const double* p,
const double* r,
@ -77,6 +89,15 @@ namespace Opm
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Gas resolution and its derivatives at bublepoint as a function of p.

View File

@ -98,6 +98,24 @@ namespace Opm
}
void SinglePvtDeadSpline::mu(const int n,
const double* p,
const double* /*r*/,
const bool* /*isSat*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = viscosity_(p[i]);
output_dmudp[i] = viscosity_.derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void SinglePvtDeadSpline::B(const int n,
const double* p,
const double* /*z*/,
@ -141,6 +159,25 @@ namespace Opm
}
void SinglePvtDeadSpline::b(const int n,
const double* p,
const double* /*r*/,
const bool* /*isSat*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = b_(p[i]);
output_dbdp[i] = b_.derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void SinglePvtDeadSpline::rbub(const int n,
const double* /*p*/,
double* output_rbub,

View File

@ -50,6 +50,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void mu(const int n,
const double* p,
const double* r,
@ -57,6 +58,16 @@ namespace Opm
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const double* p,
@ -71,6 +82,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void b(const int n,
const double* p,
const double* r,
@ -78,6 +90,16 @@ namespace Opm
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Gas resolution and its derivatives at bublepoint as a function of p.
virtual void rbub(const int n,
const double* p,

View File

@ -56,6 +56,7 @@ namespace Opm
double* output_mu) const = 0;
/// Viscosity as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void mu(const int n,
const double* p,
const double* r,
@ -63,6 +64,16 @@ namespace Opm
double* output_dmudp,
double* output_dmudr) const = 0;
/// Viscosity as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const = 0;
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const double* p,
@ -77,6 +88,7 @@ namespace Opm
double* output_dBdp) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void b(const int n,
const double* p,
const double* r,
@ -84,6 +96,16 @@ namespace Opm
double* output_dbdp,
double* output_dpdr) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dpdr) const = 0;
/// Gas resolution at bublepoint as a function of pressure
virtual void rbub(const int n,
const double* p,

View File

@ -109,6 +109,18 @@ namespace Opm
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveGas::mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// Formation volume factor as a function of p and z.
void SinglePvtLiveGas::B(const int n,
@ -149,6 +161,19 @@ namespace Opm
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void SinglePvtLiveGas::b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// Gas resolution and its derivatives at bublepoint as a function of p.
void SinglePvtLiveGas::rbub(const int n,
const double* p,

View File

@ -47,6 +47,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void mu(const int n,
const double* p,
const double* r,
@ -54,6 +55,16 @@ namespace Opm
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const double* p,
@ -68,6 +79,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void b(const int n,
const double* p,
const double* r,
@ -75,6 +87,16 @@ namespace Opm
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Gas resolution and its derivatives at bublepoint as a function of p.

View File

@ -140,6 +140,24 @@ namespace Opm
}
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveOil::mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[i] = miscible_oil(p[i], r[i],isSat[i], 2, 0);
output_dmudp[i] = miscible_oil(p[i], r[i],isSat[i], 2, 1);
output_dmudr[i] = miscible_oil(p[i], r[i],isSat[i], 2, 2);
}
}
/// Formation volume factor as a function of p and z.
void SinglePvtLiveOil::B(const int n,
@ -185,6 +203,24 @@ namespace Opm
}
}
void SinglePvtLiveOil::b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_b[i] = miscible_oil(p[i], r[i], isSat[i],1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], isSat[i], 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], isSat[i],1, 2);
}
}
void SinglePvtLiveOil::rbub(const int n,
const double* p,
double* output_rbub,
@ -423,4 +459,81 @@ namespace Opm
}
}
double SinglePvtLiveOil::miscible_oil(const double press,
const double r,
const bool isSat,
const int item,
const int deriv) const
{
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], r);
double w = (r - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (isSat) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], r);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[3][is+1]-saturated_oil_table_[3][is]);
return val;
}
} else {
if (isSat) { // Saturated case
return linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], r);
double w = (r - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

View File

@ -48,6 +48,7 @@ namespace Opm
double* output_mu) const;
/// Viscosity and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void mu(const int n,
const double* p,
const double* r,
@ -55,6 +56,16 @@ namespace Opm
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void mu(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p and z.
virtual void B(const int n,
const double* p,
@ -69,6 +80,7 @@ namespace Opm
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// The fluid is considered saturated if r >= rbub(p).
virtual void b(const int n,
const double* p,
const double* r,
@ -76,6 +88,16 @@ namespace Opm
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
/// Whether the fluid is saturated or not is given explicitly by isSat.
virtual void b(const int n,
const double* p,
const double* r,
const bool* isSat,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Gas resolution and its derivatives at bublepoint as a function of p.
virtual void rbub(const int n,
const double* p,
@ -112,6 +134,12 @@ namespace Opm
const int item,
const int deriv = 0) const;
double miscible_oil(const double press,
const double r,
const bool isSat,
const int item,
const int deriv = 0) const;
// PVT properties of live oil (with dissolved gas)
std::vector<std::vector<double> > saturated_oil_table_;
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;

View File

@ -97,6 +97,7 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
std::vector<double> p(n);
std::vector<double> r(n);
bool isSat[n];
std::vector<double> z(n * np);
std::vector<double> mu(n);
@ -122,6 +123,11 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
r[1] = 200;
r[2] = 200 + h_r;
isSat[0] = true;
isSat[1] = true;
isSat[2] = true;
// undersaturated
p[3] = p[0];
p[4] = p[1];
@ -131,6 +137,9 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
r[4] = 50;
r[5] = 50 +h_r;
isSat[3] = false;
isSat[4] = false;
isSat[5] = false;
// Corresponing z factors, used to compare with the [p,z] interface
for (int i = 0; i < n; ++i) {
@ -141,7 +150,7 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
// test mu
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->mu(n, &p[0], &r[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &p[0], &r[0], &isSat[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &p[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/h_p;
dmudr_diff = (mu_new[2]-mu_new[0])/h_r;
@ -175,7 +184,7 @@ BOOST_AUTO_TEST_CASE(test_blackoilfluid)
double dbdr_diff_u;
for (int phase = 1; phase < 2; ++phase) {
props_[phase]->b(n, &p[0], &r[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->b(n, &p[0], &r[0], &isSat[0], &b[0], &dbdp[0], &dbdr[0]);
//props_[phase]->B(n, p, z, B);
props_[phase]->dBdp(n, &p[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/h_p;