Merge pull request #4629 from akva2/perf_rate_refactor

Refactor perforation rate calculations
This commit is contained in:
Bård Skaflestad 2023-06-23 13:00:22 +02:00 committed by GitHub
commit 712b00dbdf
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 286 additions and 141 deletions

View File

@ -763,7 +763,7 @@ namespace Opm
// producing perforations // producing perforations
if (drawdown > 0.0) { if (drawdown > 0.0) {
// Do nothing is crossflow is not allowed // Do nothing if crossflow is not allowed
if (!allow_cf && this->isInjector()) { if (!allow_cf && this->isInjector()) {
return; return;
} }

View File

@ -466,6 +466,54 @@ namespace Opm
const double dis_gas_rate, const double dis_gas_rate,
const std::vector<EvalWell>& cq_s, const std::vector<EvalWell>& cq_s,
const IntensiveQuantities& intQuants) const; const IntensiveQuantities& intQuants) const;
template<class Value>
void gasOilPerfRateInj(const std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rv,
const Value& rs,
const Value& pressure,
const Value& rvw,
DeferredLogger& deferred_logger) const;
template<class Value>
void gasOilPerfRateProd(std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rv,
const Value& rs,
const Value& rvw) const;
template<class Value>
void gasWaterPerfRateProd(std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rvw,
const Value& rsw) const;
template<class Value>
void gasWaterPerfRateInj(const std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rvw,
const Value& rsw,
const Value& pressure,
DeferredLogger& deferred_logger) const;
template<class Value>
void disOilVapWatVolumeRatio(Value& volumeRatio,
const Value& rvw,
const Value& rsw,
const Value& pressure,
const std::vector<Value>& cmix_s,
const std::vector<Value>& b_perfcells_dense,
DeferredLogger& deferred_logger) const;
template<class Value>
void gasOilVolumeRatio(Value& volumeRatio,
const Value& rv,
const Value& rs,
const Value& pressure,
const std::vector<Value>& cmix_s,
const std::vector<Value>& b_perfcells_dense,
DeferredLogger& deferred_logger) const;
}; };
} }

View File

@ -37,6 +37,26 @@
#include <functional> #include <functional>
#include <numeric> #include <numeric>
namespace {
template<class dValue, class Value>
auto dValueError(const dValue& d,
const std::string& name,
const std::string& methodName,
const Value& Rs,
const Value& Rv,
const Value& pressure)
{
return fmt::format("Problematic d value {} obtained for well {}"
" during {} calculations with rs {}"
", rv {} and pressure {}."
" Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
" for this connection.", d, name, methodName, Rs, Rv, pressure);
}
}
namespace Opm namespace Opm
{ {
@ -232,44 +252,10 @@ namespace Opm
} }
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const Value cq_sOil = cq_s[oilCompIdx];
const Value cq_sGas = cq_s[gasCompIdx];
const Value dis_gas = rs * cq_sOil;
const Value vap_oil = rv * cq_sGas;
cq_s[gasCompIdx] += dis_gas;
cq_s[oilCompIdx] += vap_oil;
// recording the perforation solution gas rate and solution oil rates
if (this->isProducer()) {
perf_rates.dis_gas = getValue(dis_gas);
perf_rates.vap_oil = getValue(vap_oil);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const Value vap_wat = rvw * cq_sGas;
cq_s[waterCompIdx] += vap_wat;
if (this->isProducer())
perf_rates.vap_wat = getValue(vap_wat);
}
} else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { } else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const Value cq_sWat = cq_s[waterCompIdx];
const Value cq_sGas = cq_s[gasCompIdx];
const Value vap_wat = rvw * cq_sGas;
const Value dis_gas_wat = rsw * cq_sWat;
cq_s[waterCompIdx] += vap_wat;
cq_s[gasCompIdx] += dis_gas_wat;
if (this->isProducer()) {
perf_rates.vap_wat = getValue(vap_wat);
perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
} }
}
} else { } else {
// Do nothing if crossflow is not allowed // Do nothing if crossflow is not allowed
if (!allow_cf && this->isProducer()) { if (!allow_cf && this->isProducer()) {
@ -289,26 +275,8 @@ namespace Opm
Value volumeRatio = bhp * 0.0; // initialize it with the correct type Value volumeRatio = bhp * 0.0; // initialize it with the correct type
; ;
if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) { if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); cmix_s, b_perfcells_dense, deferred_logger);
// Incorporate RSW/RVW factors if both water and gas active
const Value d = 1.0 - rvw * rsw;
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << this->name()
<< " during computePerfRate calculations with rsw " << rsw
<< ", rvw " << rvw << " and pressure " << pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
}
const Value tmp_wat = d > 0.0? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d : cmix_s[waterCompIdx];
volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d : cmix_s[gasCompIdx];
volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
} else { } else {
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
@ -321,26 +289,8 @@ namespace Opm
} }
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); cmix_s, b_perfcells_dense, deferred_logger);
// Incorporate RS/RV factors if both oil and gas active
const Value d = 1.0 - rv * rs;
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << this->name()
<< " during computePerfRate calculations with rs " << rs
<< ", rv " << rv << " and pressure " << pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
}
const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
} }
else { else {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
@ -362,68 +312,13 @@ namespace Opm
// calculating the perforation solution gas rate and solution oil rates // calculating the perforation solution gas rate and solution oil rates
if (this->isProducer()) { if (this->isProducer()) {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); gasOilPerfRateInj(cq_s, perf_rates,
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); rv, rs, pressure, rvw, deferred_logger);
// TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
// s means standard condition, r means reservoir condition
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * b_g + rs * q_or * b_o
// d = 1.0 - rs * rv
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - getValue(rv) * getValue(rs);
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << this->name()
<< " during computePerfRate calculations with rs " << rs
<< ", rv " << rv << " and pressure " << pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
} else {
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
// q_ws = q_wr * b_w + rvw * q_gr * b_g
// q_wr = 1 / b_w * (q_ws - rvw * q_gr * b_g) = 1 / b_w * (q_ws - rvw * 1 / d (q_gs - rs * q_os))
// vaporized water in gas
// rvw * q_gr * b_g = q_ws -q_wr *b_w = rvw * (q_gs -rs *q_os) / d
perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
}
} }
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
//no oil //no oil
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); pressure, deferred_logger);
const double dw = 1.0 - getValue(rvw) * getValue(rsw);
if (dw <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << dw << " obtained for well " << this->name()
<< " during computePerfRate calculations with rsw " << rsw
<< ", rvw " << rvw << " and pressure " << pressure
<< " obtaining d " << dw
<< " Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
} else {
// vaporized water into gas
// rvw * q_gr * b_g = rvw * (q_gs - rsw * q_ws) / dw
perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
// dissolved gas in water
// rsw * q_wr * b_w = rsw * (q_ws - rvw * q_gs) / dw
perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
}
} }
} }
} }
@ -2471,4 +2366,206 @@ namespace Opm
return {Base::restrictEval(cq_s_zfrac_effective), cq_s_zfrac_effective}; return {Base::restrictEval(cq_s_zfrac_effective), cq_s_zfrac_effective};
} }
template <typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
gasOilPerfRateInj(const std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rv,
const Value& rs,
const Value& pressure,
const Value& rvw,
DeferredLogger& deferred_logger) const
{
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
// s means standard condition, r means reservoir condition
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * b_g + rs * q_or * b_o
// d = 1.0 - rs * rv
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - getValue(rv) * getValue(rs);
if (d <= 0.0) {
deferred_logger.debug(dValueError(d, this->name(),
"gasOilPerfRateInj",
rs, rv, pressure));
} else {
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_rates.dis_gas = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
// q_ws = q_wr * b_w + rvw * q_gr * b_g
// q_wr = 1 / b_w * (q_ws - rvw * q_gr * b_g) = 1 / b_w * (q_ws - rvw * 1 / d (q_gs - rs * q_os))
// vaporized water in gas
// rvw * q_gr * b_g = q_ws -q_wr *b_w = rvw * (q_gs -rs *q_os) / d
perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
}
}
template <typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
gasOilPerfRateProd(std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rv,
const Value& rs,
const Value& rvw) const
{
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const Value cq_sOil = cq_s[oilCompIdx];
const Value cq_sGas = cq_s[gasCompIdx];
const Value dis_gas = rs * cq_sOil;
const Value vap_oil = rv * cq_sGas;
cq_s[gasCompIdx] += dis_gas;
cq_s[oilCompIdx] += vap_oil;
// recording the perforation solution gas rate and solution oil rates
if (this->isProducer()) {
perf_rates.dis_gas = getValue(dis_gas);
perf_rates.vap_oil = getValue(vap_oil);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const Value vap_wat = rvw * cq_sGas;
cq_s[waterCompIdx] += vap_wat;
if (this->isProducer())
perf_rates.vap_wat = getValue(vap_wat);
}
}
template <typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
gasWaterPerfRateProd(std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rvw,
const Value& rsw) const
{
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const Value cq_sWat = cq_s[waterCompIdx];
const Value cq_sGas = cq_s[gasCompIdx];
const Value vap_wat = rvw * cq_sGas;
const Value dis_gas_wat = rsw * cq_sWat;
cq_s[waterCompIdx] += vap_wat;
cq_s[gasCompIdx] += dis_gas_wat;
if (this->isProducer()) {
perf_rates.vap_wat = getValue(vap_wat);
perf_rates.dis_gas_in_water = getValue(dis_gas_wat);
}
}
template <typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
gasWaterPerfRateInj(const std::vector<Value>& cq_s,
PerforationRates& perf_rates,
const Value& rvw,
const Value& rsw,
const Value& pressure,
DeferredLogger& deferred_logger) const
{
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const double dw = 1.0 - getValue(rvw) * getValue(rsw);
if (dw <= 0.0) {
deferred_logger.debug(dValueError(dw, this->name(),
"gasWaterPerfRateInj",
rsw, rvw, pressure));
} else {
// vaporized water into gas
// rvw * q_gr * b_g = rvw * (q_gs - rsw * q_ws) / dw
perf_rates.vap_wat = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rsw) * getValue(cq_s[waterCompIdx])) / dw;
// dissolved gas in water
// rsw * q_wr * b_w = rsw * (q_ws - rvw * q_gs) / dw
perf_rates.dis_gas_in_water = getValue(rsw) * (getValue(cq_s[waterCompIdx]) - getValue(rvw) * getValue(cq_s[gasCompIdx])) / dw;
}
}
template <typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
disOilVapWatVolumeRatio(Value& volumeRatio,
const Value& rvw,
const Value& rsw,
const Value& pressure,
const std::vector<Value>& cmix_s,
const std::vector<Value>& b_perfcells_dense,
DeferredLogger& deferred_logger) const
{
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// Incorporate RSW/RVW factors if both water and gas active
const Value d = 1.0 - rvw * rsw;
if (d <= 0.0) {
deferred_logger.debug(dValueError(d, this->name(),
"disOilVapWatVolumeRatio",
rsw, rvw, pressure));
}
const Value tmp_wat = d > 0.0 ? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d
: cmix_s[waterCompIdx];
volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx];
const Value tmp_gas = d > 0.0 ? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d
: cmix_s[waterCompIdx];
volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
}
template <typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
gasOilVolumeRatio(Value& volumeRatio,
const Value& rv,
const Value& rs,
const Value& pressure,
const std::vector<Value>& cmix_s,
const std::vector<Value>& b_perfcells_dense,
DeferredLogger& deferred_logger) const
{
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// Incorporate RS/RV factors if both oil and gas active
const Value d = 1.0 - rv * rs;
if (d <= 0.0) {
deferred_logger.debug(dValueError(d, this->name(),
"gasOilVolumeRatio",
rs, rv, pressure));
}
const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
}
} // namespace Opm } // namespace Opm