OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Make thp-func constistent with bhp-func and add some damage prevention logic for severe extrapolation

Browse Source
This commit is contained in:
Stein Krogstad
2024-05-28 12:11:44 +02:00
parent c92d80e0e8
commit c5f7d02bf7
6 changed files with 86 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

95
opm/simulators/wells/VFPHelpers.cpp
View File

@@ -363,7 +363,8 @@ template<class Scalar>
Scalar VFPHelpers<Scalar>::
findTHP(const std::vector<Scalar>& bhp_array,
const std::vector<double>& thp_array,
Scalar bhp)
Scalar bhp,
const bool find_largest)
{
int nthp = thp_array.size();
@@ -429,18 +430,60 @@ findTHP(const std::vector<Scalar>& bhp_array,
}
//bhp_array not sorted, raw search.
else {
//Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
//Since the BHP values might not be sorted, first search within
//our interpolation values, and then try to extrapolate.
int i=0;
bool found = false;
for (; i<nthp-1; ++i) {
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
//Here we're into damage prevention territory, and there may be any number of
//solutions (including zero). The well is currently not controlled by THP, and
//since we're doing severe extrapolaton we would also like, if possible, to prevent
//it from switcing to THP. Accordingly, if there are multiple solutions, we return
//the value for the intersection corresponding to the largest (smallest) THP-value
//for producers (injectors).
if (y0 < bhp && bhp <= y1) {
//first check which extrapolations are valid
const bool first_slope_positive = bhp_array[1] >= bhp_array[0];
const bool valid_low = (bhp < bhp_array[0] && first_slope_positive) || (bhp >= bhp_array[0] && !first_slope_positive);
const bool last_slope_positive = bhp_array[nthp-1] >= bhp_array[nthp-2];
const bool valid_high = (bhp > bhp_array[nthp-1] && last_slope_positive) || (bhp <= bhp_array[nthp-1] && !last_slope_positive);
bool found = false;
int i = 0;
if (find_largest){//find intersection corresponding to the largest thp
// high extrap -> table interp -> low extrap
if (valid_high) {
found = true;
break;
i = nthp-2;
} else {
//search backward within table
for (; i>=0; --i) {
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
if (std::min(y0, y1) < bhp && bhp <= std::max(y0, y1)) {
found = true;
break;
}
}
if (!found && valid_low) {
found = true;
i = 0;
}
}
} else {//find intersection corresponding to the smallest thp
//low extrap -> table interp -> high extrap
if (valid_low) {
found = true;
i = 0;
} else {
//search forward within table
for (; i<nthp-1; ++i) {
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
if (std::min(y0, y1) < bhp && bhp <= std::max(y0, y1)) {
found = true;
break;
}
}
if (!found && valid_high) {
found = true;
i = nthp-2;
}
}
}
if (found) {
@@ -449,30 +492,16 @@ findTHP(const std::vector<Scalar>& bhp_array,
const Scalar& y0 = bhp_array[i ];
const Scalar& y1 = bhp_array[i+1];
thp = findX(x0, x1, y0, y1, bhp);
} else {
// no intersection, just return largest/smallest value in table
if (find_largest) {
thp = thp_array[nthp-1];
} else {
thp = thp_array[0];
}
}
else if (bhp <= bhp_array[0]) {
//TODO: LOG extrapolation
const Scalar& x0 = thp_array[0];
const Scalar& x1 = thp_array[1];
const Scalar& y0 = bhp_array[0];
const Scalar& y1 = bhp_array[1];
thp = findX(x0, x1, y0, y1, bhp);
}
//Target bhp greater than all values in array, extrapolate
else if (bhp > bhp_array[nthp-1]) {
//TODO: LOG extrapolation
const Scalar& x0 = thp_array[nthp-2];
const Scalar& x1 = thp_array[nthp-1];
const Scalar& y0 = bhp_array[nthp-2];
const Scalar& y1 = bhp_array[nthp-1];
thp = findX(x0, x1, y0, y1, bhp);
}
else {
OPM_THROW(std::logic_error, "Programmer error: Unable to find THP in THP array");
}
return thp;
}
return thp;
}
template<class Scalar>

3
opm/simulators/wells/VFPHelpers.hpp
View File

@@ -189,7 +189,8 @@ public:
*/
static Scalar findTHP(const std::vector<Scalar>& bhp_array,
const std::vector<double>& thp_array,
Scalar bhp);
Scalar bhp,
const bool find_largest = true);
/**
* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq

2
opm/simulators/wells/VFPInjProperties.cpp
View File

@@ -76,7 +76,7 @@ thp(const int table_id,
bhp_array[i] = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i).value;
}
return VFPHelpers<Scalar>::findTHP(bhp_array, thp_array, bhp_arg);
return VFPHelpers<Scalar>::findTHP(bhp_array, thp_array, bhp_arg, /*find_largest*/ false);
}
template<class Scalar>

30
opm/simulators/wells/VFPProdProperties.cpp
View File

@@ -38,25 +38,19 @@ thp(const int table_id,
const Scalar liquid,
const Scalar vapour,
const Scalar bhp_arg,
const Scalar alq) const
const Scalar alq,
const Scalar explicit_wfr,
const Scalar explicit_gfr,
const bool use_expvfp) const
{
const VFPProdTable& table = detail::getTable(m_tables, table_id);
// Find interpolation variables.
Scalar flo = 0.0;
Scalar wfr = 0.0;
Scalar gfr = 0.0;
if (aqua == 0.0 && liquid == 0.0 && vapour == 0.0) {
// All zero, likely at initial state.
// Set FLO variable to minimum to avoid extrapolation.
// The water and gas fractions are kept at 0.0.
flo = table.getFloAxis().front();
} else {
// The usual case.
// Recall that production rate is negative in Opm, so switch the sign.
flo = -detail::getFlo(table, aqua, liquid, vapour);
wfr = detail::getWFR(table, aqua, liquid, vapour);
gfr = detail::getGFR(table, aqua, liquid, vapour);
Scalar flo = detail::getFlo(table, aqua, liquid, vapour);
Scalar wfr = detail::getWFR(table, aqua, liquid, vapour);
Scalar gfr = detail::getGFR(table, aqua, liquid, vapour);
if (use_expvfp || -flo < table.getFloAxis().front()) {
wfr = explicit_wfr;
gfr = explicit_gfr;
}
const std::vector<double> thp_array = table.getTHPAxis();
@@ -67,7 +61,7 @@ thp(const int table_id,
* by interpolating for every value of thp. This might be somewhat
* expensive, but let us assome that nthp is small.
*/
auto flo_i = VFPHelpers<Scalar>::findInterpData( flo, table.getFloAxis());
auto flo_i = VFPHelpers<Scalar>::findInterpData(-flo, table.getFloAxis());
auto wfr_i = VFPHelpers<Scalar>::findInterpData( wfr, table.getWFRAxis());
auto gfr_i = VFPHelpers<Scalar>::findInterpData( gfr, table.getGFRAxis());
auto alq_i = VFPHelpers<Scalar>::findInterpData( alq, table.getALQAxis());
@@ -77,7 +71,7 @@ thp(const int table_id,
bhp_array[i] = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
}
return VFPHelpers<Scalar>::findTHP(bhp_array, thp_array, bhp_arg);
return VFPHelpers<Scalar>::findTHP(bhp_array, thp_array, bhp_arg, /*find_largest*/ true);
}
template<class Scalar>

5
opm/simulators/wells/VFPProdProperties.hpp
View File

@@ -109,7 +109,10 @@ public:
const Scalar liquid,
const Scalar vapour,
const Scalar bhp,
const Scalar alq) const;
const Scalar alq,
const Scalar explicit_wrf,
const Scalar explicit_gfr,
const bool use_expvfp) const;
/**
* Returns the table associated with the ID, or throws an exception if

7
opm/simulators/wells/WellBhpThpCalculator.cpp
View File

@@ -143,11 +143,14 @@ calculateThpFromBhp(const std::vector<Scalar>& rates,
else if (well_.isProducer()) {
const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
const bool use_vfpexp = well_.useVfpExplicit();
const double wfr = well_.vfpProperties()->getExplicitWFR(table_id, well_.indexOfWell());
const double gfr = well_.vfpProperties()->getExplicitGFR(table_id, well_.indexOfWell());
auto thp_func =
[this, table_id, aqua, liquid, vapour, dp, &alq]
[this, table_id, aqua, liquid, vapour, dp, &alq, wfr, gfr, use_vfpexp]
(const Scalar bhp_value, const Scalar pressure_loss) {
return this->well_.vfpProperties()->getProd()->thp(
table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value());
table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value(), wfr, gfr, use_vfpexp);
};
thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger);
}