From 085ccccb1e95398d1fd4ed9866ec9aa3253ca4fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2022 18:21:11 +0200 Subject: [PATCH] Capture More Dynamic State Per Segment In particular, calculate segment flow velocities, segment holdup fractions, and segment phase viscosities when updating the well state. This is purely for reporting purposes and does not affect the multisegmented model. --- opm/simulators/wells/MultisegmentWellEval.cpp | 133 ++++++++++++++++-- opm/simulators/wells/MultisegmentWellEval.hpp | 2 + opm/simulators/wells/SegmentState.cpp | 64 +++++---- opm/simulators/wells/SegmentState.hpp | 26 +++- 4 files changed, 188 insertions(+), 37 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index a664c77d0..baef79468 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -19,6 +19,7 @@ */ #include + #include #include @@ -33,13 +34,20 @@ #include #include #include +#include #include #include #include #include #include +#include +#include #include +#include +#include +#include +#include namespace Opm { @@ -1397,6 +1405,8 @@ updateWellStateFromPrimaryVariables(WellState& well_state, static constexpr int Oil = BlackoilPhases::Liquid; static constexpr int Water = BlackoilPhases::Aqua; + const auto pvtReg = std::max(this->baseif_.wellEcl().pvt_table_number() - 1, 0); + const PhaseUsage& pu = baseif_.phaseUsage(); assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); const int oil_pos = pu.phase_pos[Oil]; @@ -1409,13 +1419,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state, std::vector fractions(baseif_.numPhases(), 0.0); fractions[oil_pos] = 1.0; - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const int water_pos = pu.phase_pos[Water]; fractions[water_pos] = primary_variables_[seg][WFrac]; fractions[oil_pos] -= fractions[water_pos]; } - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const int gas_pos = pu.phase_pos[Gas]; fractions[gas_pos] = primary_variables_[seg][GFrac]; fractions[oil_pos] -= fractions[gas_pos]; @@ -1445,16 +1455,123 @@ updateWellStateFromPrimaryVariables(WellState& well_state, // update the segment pressure segment_pressure[seg] = primary_variables_[seg][SPres]; + if (seg == 0) { // top segment ws.bhp = segment_pressure[seg]; } + + // Calculate other per-phase dynamic quantities. + + const auto temperature = 0.0; // Ignore thermal effects + const auto saltConc = 0.0; // Ignore salt precipitation + const auto Rvw = 0.0; // Ignore vaporised water. + + auto rsMax = 0.0; + auto rvMax = 0.0; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + // Both oil and gas active. + rsMax = FluidSystem::oilPvt() + .saturatedGasDissolutionFactor(pvtReg, temperature, segment_pressure[seg]); + + rvMax = FluidSystem::gasPvt() + .saturatedOilVaporizationFactor(pvtReg, temperature, segment_pressure[seg]); + } + + // 1) Local condition volume flow rates + const auto& [Rs, Rv] = this->baseif_.rateConverter().inferDissolvedVaporisedRatio + (rsMax, rvMax, segment_rates.begin() + (seg + 0)*this->baseif_.numPhases()); + + { + // Use std::span<> in C++20 and beyond. + const auto rate_start = (seg + 0) * this->baseif_.numPhases(); + const auto* surf_rates = segment_rates.data() + rate_start; + auto* resv_rates = segments.phase_resv_rates.data() + rate_start; + + this->baseif_.rateConverter().calcReservoirVoidageRates + (pvtReg, segment_pressure[seg], + std::max(0.0, Rs), + std::max(0.0, Rv), + temperature, saltConc, surf_rates, resv_rates); + } + + // 2) Local condition holdup fractions. + const auto tot_resv = + std::accumulate(segments.phase_resv_rates.begin() + (seg + 0)*this->baseif_.numPhases(), + segments.phase_resv_rates.begin() + (seg + 1)*this->baseif_.numPhases(), + 0.0); + + std::transform(segments.phase_resv_rates.begin() + (seg + 0)*this->baseif_.numPhases(), + segments.phase_resv_rates.begin() + (seg + 1)*this->baseif_.numPhases(), + segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(), + [tot_resv](const auto qr) { return std::clamp(qr / tot_resv, 0.0, 1.0); }); + + // 3) Local condition flow velocities for segments other than top segment. + if (seg > 0) { + // Possibly poor approximation + // Velocity = Flow rate / cross-sectional area. + // Additionally ignores drift flux. + const auto area = this->baseif_.wellEcl().getSegments() + .getFromSegmentNumber(segments.segment_number()[seg]).crossArea(); + const auto velocity = (area > 0.0) ? tot_resv / area : 0.0; + + std::transform(segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(), + segments.phase_holdup.begin() + (seg + 1)*this->baseif_.numPhases(), + segments.phase_velocity.begin() + (seg + 0)*this->baseif_.numPhases(), + [velocity](const auto hf) { return (hf > 0.0) ? velocity : 0.0; }); + } + + // 4) Local condition phase viscosities. + segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Oil]] = + FluidSystem::oilPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rs); + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Water]] = + FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], saltConc); + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + segments.phase_viscosity[seg*this->baseif_.numPhases() + pu.phase_pos[Gas]] = + FluidSystem::gasPvt().viscosity(pvtReg, temperature, segment_pressure[seg], Rv, Rvw); + } } - WellBhpThpCalculator(baseif_). - updateThp(rho, [this]() { return baseif_.wellEcl().alq_value(); }, - {FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx), - FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx), - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)}, - well_state, deferred_logger); + + // Segment flow velocity in top segment. + { + const auto np = this->baseif_.numPhases(); + auto segVel = [&segments, np](const auto segmentNumber) + { + auto v = 0.0; + const auto* vel = segments.phase_velocity.data() + segmentNumber*np; + for (auto p = 0*np; p < np; ++p) { + if (std::abs(vel[p]) > std::abs(v)) { + v = vel[p]; + } + } + + return v; + }; + + const auto seg = 0; + auto maxVel = 0.0; + for (const auto& inlet : this->segmentSet()[seg].inletSegments()) { + const auto v = segVel(this->segmentNumberToIndex(inlet)); + if (std::abs(v) > std::abs(maxVel)) { + maxVel = v; + } + } + + std::transform(segments.phase_holdup.begin() + (seg + 0)*this->baseif_.numPhases(), + segments.phase_holdup.begin() + (seg + 1)*this->baseif_.numPhases(), + segments.phase_velocity.begin() + (seg + 0)*this->baseif_.numPhases(), + [maxVel](const auto hf) { return (hf > 0.0) ? maxVel : 0.0; }); + } + + WellBhpThpCalculator(this->baseif_) + .updateThp(rho, [this]() { return this->baseif_.wellEcl().alq_value(); }, + {FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx), + FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx), + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)}, + well_state, deferred_logger); } template diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index 8c151a01d..cb6006055 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -35,6 +35,8 @@ #include #include +#include +#include namespace Dune { template class UMFPack; diff --git a/opm/simulators/wells/SegmentState.cpp b/opm/simulators/wells/SegmentState.cpp index 64ad5fea9..e8077cb85 100644 --- a/opm/simulators/wells/SegmentState.cpp +++ b/opm/simulators/wells/SegmentState.cpp @@ -16,45 +16,60 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . */ + #if HAVE_CONFIG_H #include "config.h" #endif // HAVE_CONFIG_H -#include -#include - #include -#include #include +#include +#include +#include +#include +#include + +namespace { + +std::vector make_segment_number(const Opm::WellSegments& segments) +{ + std::vector segment_number; + segment_number.reserve(segments.size()); + + std::transform(segments.begin(), segments.end(), + std::back_inserter(segment_number), + [](const Opm::Segment& segment) + { + return segment.segmentNumber(); + }); + + return segment_number; +} + +} // Anonymous namespace + namespace Opm { -namespace { -std::vector make_segment_number( const WellSegments& segments ) { - std::vector segment_number; - std::transform(segments.begin(), segments.end(), std::back_insert_iterator(segment_number), [](const Segment& segment) { return segment.segmentNumber(); }); - return segment_number; -} - -} - -SegmentState::SegmentState(int num_phases, const WellSegments& segments) : - rates(segments.size() * num_phases), - pressure(segments.size()), - pressure_drop_friction(segments.size()), - pressure_drop_hydrostatic(segments.size()), - pressure_drop_accel(segments.size()), - m_segment_number(make_segment_number(segments)) -{ -} +SegmentState::SegmentState(int num_phases, const WellSegments& segments) + : rates (segments.size() * num_phases) + , phase_resv_rates (segments.size() * num_phases) + , phase_velocity (segments.size() * num_phases) + , phase_holdup (segments.size() * num_phases) + , phase_viscosity (segments.size() * num_phases) + , pressure (segments.size()) + , pressure_drop_friction (segments.size()) + , pressure_drop_hydrostatic(segments.size()) + , pressure_drop_accel (segments.size()) + , m_segment_number (make_segment_number(segments)) +{} double SegmentState::pressure_drop(std::size_t index) const { return this->pressure_drop_friction[index] + this->pressure_drop_hydrostatic[index] + this->pressure_drop_accel[index]; } - bool SegmentState::empty() const { return this->rates.empty(); } @@ -63,7 +78,6 @@ std::size_t SegmentState::size() const { return this->pressure.size(); } - void SegmentState::scale_pressure(const double bhp) { if (this->empty()) throw std::logic_error("Tried to pressure scale empty SegmentState"); @@ -80,4 +94,4 @@ const std::vector& SegmentState::segment_number() const { return this->m_segment_number; } -} +} // namespace Opm diff --git a/opm/simulators/wells/SegmentState.hpp b/opm/simulators/wells/SegmentState.hpp index 7a71d2d27..efcdc88d3 100644 --- a/opm/simulators/wells/SegmentState.hpp +++ b/opm/simulators/wells/SegmentState.hpp @@ -20,35 +20,53 @@ #ifndef OPM_SEGMENTSTATE_HEADER_INCLUDED #define OPM_SEGMENTSTATE_HEADER_INCLUDED +#include #include namespace Opm { - class WellSegments; -class WellConnections; +} // namespace Opm +namespace Opm +{ class SegmentState { public: SegmentState() = default; SegmentState(int num_phases, const WellSegments& segments); + double pressure_drop(std::size_t index) const; bool empty() const; void scale_pressure(double bhp); + const std::vector& segment_number() const; std::size_t size() const; std::vector rates; + + /// Segment condition volume flow rates through segment (per phase) + std::vector phase_resv_rates; + + /// Segment condition flow velocity through segment (per phase) + std::vector phase_velocity; + + /// Segment condition holdup fractions through segment (per phase) + std::vector phase_holdup; + + /// Segment condition phase viscosities. + std::vector phase_viscosity; + std::vector pressure; std::vector pressure_drop_friction; std::vector pressure_drop_hydrostatic; std::vector pressure_drop_accel; + private: std::vector m_segment_number; }; -} +} // namepace Opm -#endif +#endif // OPM_SEGMENTSTATE_HEADER_INCLUDED