diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 977c6c059..e20e5dfd3 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -187,6 +187,7 @@ namespace Opm using Base::wellHasTHPConstraints; using Base::mostStrictBhpFromBhpLimits; using Base::scalingFactor; + using Base::scaleProductivityIndex; // protected member variables from the Base class using Base::current_step_; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index bc52550ea..950e937fc 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -461,6 +461,10 @@ namespace Opm well_state.wellVaporizedOilRates()[index_of_well_] = 0.; well_state.wellDissolvedGasRates()[index_of_well_] = 0.; + for (int p = 0; p < np; ++p) { + well_state.productivityIndex()[np*index_of_well_ + p] = 0.; + } + for (int perf = 0; perf < number_of_perforations_; ++perf) { const int cell_idx = well_cells_[perf]; @@ -601,8 +605,19 @@ namespace Opm // Store the perforation pressure for later usage. well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf]; + + // Compute Productivity index + for (int p = 0; p < np; ++p) { + const unsigned int compIdx = flowPhaseToEbosCompIdx(p); + const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value(); + double productivity_index = cq_s[compIdx].value() / drawdown; + scaleProductivityIndex(perf, productivity_index); + well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index; + } + } + // add vol * dF/dt + Q to the well equations; for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) { EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt; diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index d613b2126..4fcf5e4e4 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -333,6 +333,8 @@ namespace Opm void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector& B_avg, bool terminal_output); + void scaleProductivityIndex(const int perfIdx, double& productivity_index) const; + }; } diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 129cafc00..e111c63a7 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -103,6 +103,7 @@ namespace Opm saturation_table_number_.begin() ); } well_efficiency_factor_ = 1.0; + } @@ -1048,4 +1049,37 @@ namespace Opm } } } + + template + void + WellInterface::scaleProductivityIndex(const int perfIdx, double& productivity_index) const + { + + const auto& connection = well_ecl_->getConnections(current_step_)[perfIdx]; + + if (well_ecl_->getDrainageRadius(current_step_) < 0) { + OpmLog::warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero"); + productivity_index = 0.0; + return; + } + + if (connection.r0() > well_ecl_->getDrainageRadius(current_step_)) { + OpmLog::info("PRODUCTIVITY_INDEX_INFO", "The effective radius is larger then the well drainage radius for well " + name() + + " They are set to equal in the well productivity index calculations"); + return; + } + + // For zero drainage radius the productivity index is just the transmissibility times the mobility + if (well_ecl_->getDrainageRadius(current_step_) == 0) { + return; + } + + // Scale the productivity index to account for the drainage radius. + // Assumes steady radial flow only valied for horizontal wells + productivity_index *= + (std::log(connection.r0() / connection.rw()) + connection.skinFactor()) / + (std::log(well_ecl_->getDrainageRadius(current_step_) / connection.rw()) + connection.skinFactor()); + } + + } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index fb80ca5bb..8f27b5d00 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -128,6 +128,8 @@ namespace Opm perfRateSolvent_.clear(); perfRateSolvent_.resize(nperf, 0.0); + productivity_index_.resize(nw * np, 0.0); + // intialize wells that have been there before // order may change so the mapping is based on the well name if(prevState && !prevState->wellMap().empty()) { @@ -262,6 +264,8 @@ namespace Opm well_dissolved_gas_rates_.resize(nw, 0.0); well_vaporized_oil_rates_.resize(nw, 0.0); + productivity_index_.resize(nw * np, 0.0); + // Ensure that we start out with zero rates by default. perfphaserates_.clear(); perfphaserates_.resize(nperf * np, 0.0); @@ -480,6 +484,18 @@ namespace Opm well.rates.set( rt::reservoir_gas, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Gas]] ); } + if ( pu.phase_used[Water] ) { + well.rates.set( rt::productivity_index_water, this->productivity_index_[well_rate_index + pu.phase_pos[Water]] ); + } + + if ( pu.phase_used[Oil] ) { + well.rates.set( rt::productivity_index_oil, this->productivity_index_[well_rate_index + pu.phase_pos[Oil]] ); + } + + if ( pu.phase_used[Gas] ) { + well.rates.set( rt::productivity_index_gas, this->productivity_index_[well_rate_index + pu.phase_pos[Gas]] ); + } + well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] ); well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] ); @@ -761,6 +777,14 @@ namespace Opm return top_segment_index_[w]; } + std::vector& productivityIndex() { + return productivity_index_; + } + + const std::vector& productivityIndex() const { + return productivity_index_; + } + private: std::vector perfphaserates_; std::vector current_controls_; @@ -793,6 +817,9 @@ namespace Opm std::vector top_segment_index_; int nseg_; // total number of the segments + // Productivity Index + std::vector productivity_index_; + }; } // namespace Opm