diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 154984526..09a019e3d 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -188,7 +188,7 @@ namespace Opm template void computePerfRate(const IntensiveQuantities& int_quants, const std::vector& mob_perfcells, - const double Tw, + const std::vector& Tw, const int seg, const int perf, const Value& segment_pressure, @@ -204,7 +204,7 @@ namespace Opm const Value& rv, const std::vector& b_perfcells, const std::vector& mob_perfcells, - const double Tw, + const std::vector& Tw, const int perf, const Value& segment_pressure, const Value& segment_density, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 881246ae0..838df75f3 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -349,7 +349,7 @@ namespace Opm getMobility(ebosSimulator, perf, mob, deferred_logger); const double trans_mult = ebosSimulator.problem().template wellTransMultiplier(intQuants, cell_idx); const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); const Scalar seg_pressure = segment_pressure[seg]; std::vector cq_s(this->num_components_, 0.); Scalar perf_press = 0.0; @@ -768,7 +768,7 @@ namespace Opm const Value& rv, const std::vector& b_perfcells, const std::vector& mob_perfcells, - const double Tw, + const std::vector& Tw, const int perf, const Value& segment_pressure, const Value& segment_density, @@ -804,7 +804,7 @@ namespace Opm // compute component volumetric rates at standard conditions for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) { - const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown); + const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown); cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; } @@ -828,9 +828,6 @@ namespace Opm total_mob += mob_perfcells[comp_idx]; } - // injection perforations total volume rates - const Value cqt_i = - Tw * (total_mob * drawdown); - // compute volume ratio between connection and at standard conditions Value volume_ratio = 0.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { @@ -871,9 +868,10 @@ namespace Opm } } // injecting connections total volumerates at standard conditions - Value cqt_is = cqt_i / volume_ratio; - for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) { - cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is; + for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) { + const Value cqt_i = - Tw[componentIdx] * (total_mob * drawdown); + Value cqt_is = cqt_i / volume_ratio; + cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; } } // end for injection perforations @@ -907,7 +905,7 @@ namespace Opm MultisegmentWell:: computePerfRate(const IntensiveQuantities& int_quants, const std::vector& mob_perfcells, - const double Tw, + const std::vector& Tw, const int seg, const int perf, const Value& segment_pressure, @@ -1182,11 +1180,11 @@ namespace Opm // the well index associated with the connection const double trans_mult = ebos_simulator.problem().template wellTransMultiplier(int_quantities, cell_idx); const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol); + const std::vector tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol); std::vector ipr_a_perf(this->ipr_a_.size()); std::vector ipr_b_perf(this->ipr_b_.size()); for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { - const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx]; + const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx]; ipr_a_perf[comp_idx] += tw_mob * pressure_diff; ipr_b_perf[comp_idx] += tw_mob; } @@ -1682,7 +1680,7 @@ namespace Opm getMobility(ebosSimulator, perf, mob, deferred_logger); const double trans_mult = ebosSimulator.problem().template wellTransMultiplier(int_quants, cell_idx); const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol); std::vector cq_s(this->num_components_, 0.0); EvalWell perf_press; PerforationRates perfRates; @@ -1995,7 +1993,7 @@ namespace Opm getMobility(ebosSimulator, perf, mob, deferred_logger); const double trans_mult = ebosSimulator.problem().template wellTransMultiplier(int_quants, cell_idx); const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol); std::vector cq_s(this->num_components_, 0.0); Scalar perf_press = 0.0; PerforationRates perf_rates; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 728ecb8d9..a19379b7b 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -286,7 +286,7 @@ namespace Opm void computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob, const Value& bhp, - const double Tw, + const std::vector& Tw, const int perf, const bool allow_cf, std::vector& cq_s, @@ -302,7 +302,7 @@ namespace Opm const Value& rvw, const Value& rsw, std::vector& b_perfcells_dense, - const double Tw, + const std::vector& Tw, const int perf, const bool allow_cf, const Value& skin_pressure, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 4cb54f2ef..320a4e26a 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -119,7 +119,7 @@ namespace Opm computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob, const Value& bhp, - const double Tw, + const std::vector& Tw, const int perf, const bool allow_cf, std::vector& cq_s, @@ -224,7 +224,7 @@ namespace Opm const Value& rvw, const Value& rsw, std::vector& b_perfcells_dense, - const double Tw, + const std::vector& Tw, const int perf, const bool allow_cf, const Value& skin_pressure, @@ -249,7 +249,7 @@ namespace Opm // compute component volumetric rates at standard conditions for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) { - const Value cq_p = - Tw * (mob[componentIdx] * drawdown); + const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown); cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p; } @@ -270,9 +270,6 @@ namespace Opm total_mob_dense += mob[componentIdx]; } - // injection perforations total volume rates - const Value cqt_i = - Tw * (total_mob_dense * drawdown); - // compute volume ratio between connection at standard conditions Value volumeRatio = bhp * 0.0; // initialize it with the correct type ; @@ -306,8 +303,9 @@ namespace Opm } // injecting connections total volumerates at standard conditions - Value cqt_is = cqt_i / volumeRatio; for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) { + const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown); + Value cqt_is = cqt_i / volumeRatio; cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; } @@ -496,7 +494,7 @@ namespace Opm PerforationRates perf_rates; double trans_mult = ebosSimulator.problem().template wellTransMultiplier(intQuants, cell_idx); const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf, cq_s, perf_rates, deferred_logger); @@ -790,12 +788,13 @@ namespace Opm } // the well index associated with the connection - const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template wellTransMultiplier(int_quantities, cell_idx); - + double trans_mult = ebos_simulator.problem().template wellTransMultiplier(int_quantities, cell_idx); + const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); + const std::vector tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol); std::vector ipr_a_perf(this->ipr_a_.size()); std::vector ipr_b_perf(this->ipr_b_.size()); for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { - const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx]; + const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx]; ipr_a_perf[comp_idx] += tw_mob * pressure_diff; ipr_b_perf[comp_idx] += tw_mob; } @@ -1366,7 +1365,7 @@ namespace Opm getMobility(ebosSimulator, perf, mob, deferred_logger); double trans_mult = ebosSimulator.problem().template wellTransMultiplier(intQuants, cell_idx); const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); std::vector cq_s(this->num_components_, 0.); PerforationRates perf_rates; @@ -1669,7 +1668,8 @@ namespace Opm std::vector cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.}); PerforationRates perf_rates; double trans_mult = ebos_simulator.problem().template wellTransMultiplier(int_quant, cell_idx); - const double Tw = this->well_index_[perf] * trans_mult; + const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); + const std::vector Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol); computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s, perf_rates, deferred_logger); // TODO: make area a member @@ -2274,7 +2274,7 @@ namespace Opm std::vector cq_s(this->num_components_, 0.); double trans_mult = ebosSimulator.problem().template wellTransMultiplier(intQuants, cell_idx); const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); PerforationRates perf_rates; computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf, cq_s, perf_rates, deferred_logger); diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index af325b6c7..e13c0cdb5 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -349,7 +349,7 @@ public: return 0; } - double wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const; + std::vector wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const; void updateConnectionDFactor(const Simulator& simulator, SingleWellState& ws) const; @@ -447,6 +447,8 @@ protected: double* connII, DeferredLogger& deferred_logger) const; + double computeConnectionDFactor(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const double total_tw, const SingleWellState& ws) const; + }; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 247e53b01..29fb091e0 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1324,40 +1324,37 @@ namespace Opm } template - double + std::vector WellInterface:: wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const { + std::vector wi(this->num_components_, this->well_index_[perf] * trans_mult); const auto& wdfac = this->well_ecl_.getWDFAC(); if (!wdfac.useDFactor()) { - return this->well_index_[perf] * trans_mult; + return wi; } - + // for gas wells we may want to add a Forchheimer term if the WDFAC or WDFACCOR keyword is used if constexpr (! Indices::gasEnabled) { - return this->well_index_[perf] * trans_mult; + return wi; } - // closed connection are still closed if (this->well_index_[perf] == 0) - return 0.0; + return std::vector(this->num_components_, 0.0); - // for gas wells we may want to add a Forchheimer term if the WDFAC or WDFACCOR keyword is used - const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]]; - // viscosity is evaluated at connection pressure - const double connection_pressure = ws.perf_data.pressure[perf]; - const double mu = FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, getValue(intQuants.fluidState().Rv()), getValue(intQuants.fluidState().Rvw())); - const double phi = getValue(intQuants.porosity()); - //double k = connection.Kh()/h * trans_mult; - double Kh = connection.Kh()* trans_mult; - double Ke = connection.Ke()* trans_mult; - double h = Kh / Ke; - double rw = connection.rw(); - double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx()); - double scaling = 3.141592653589 * Kh; - double d = wdfac.useConnectionDFactor()? connection.dFactor() : wdfac.getDFactor(rho, mu, Ke, phi, rw, h); + double tot_tw = 0.0; + for (const auto& c : this->well_ecl_.getConnections()) { + tot_tw += c.CF(); + } + + double d = computeConnectionDFactor(perf, intQuants, trans_mult, tot_tw, ws); const PhaseUsage& pu = this->phaseUsage(); double Q = std::abs(ws.perf_data.phase_rates[perf*pu.num_phases + pu.phase_pos[Gas]]); - return 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (Q/2 * d / scaling)); + const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]]; + double Kh = connection.Kh()* trans_mult; + double scaling = 3.141592653589 * Kh; + const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (Q/2 * d / scaling)); + return wi; } template @@ -1369,27 +1366,54 @@ namespace Opm if (!wdfac.useDFactor()) { return; } - double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx()); auto& perf_data = ws.perf_data; - + double tot_tw = 0.0; + for (const auto& c : this->well_ecl_.getConnections()) { + tot_tw += c.CF(); + } for (int perf = 0; perf < this->number_of_perforations_; ++perf) { const int cell_idx = this->well_cells_[perf]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); const double trans_mult = simulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); - // viscosity is evaluated at connection pressure - const double connection_pressure = ws.perf_data.pressure[perf]; - const double mu = FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, getValue(intQuants.fluidState().Rv()), getValue(intQuants.fluidState().Rvw())); - const double phi = getValue(intQuants.porosity()); - const auto& connection = this->well_ecl_.getConnections()[perf_data.ecl_index[perf]]; - double Kh = connection.Kh()* trans_mult; - double Ke = connection.Ke()* trans_mult; - double h = Kh / Ke; - double rw = connection.rw(); - double d = wdfac.useConnectionDFactor()? connection.dFactor() : wdfac.getDFactor(rho, mu, Ke, phi, rw, h); - perf_data.connection_d_factor[perf] = d; + perf_data.connection_d_factor[perf] = computeConnectionDFactor(perf, intQuants, trans_mult, tot_tw, ws); } } + template + double + WellInterface:: + computeConnectionDFactor(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const double total_tw, const SingleWellState& ws) const { + const double connection_pressure = ws.perf_data.pressure[perf]; + // viscosity is evaluated at connection pressure + const auto& rv = getValue(intQuants.fluidState().Rv()); + const double psat = FluidSystem::gasPvt().saturationPressure(this->pvtRegionIdx(), ws.temperature, rv); + const double mu = connection_pressure < psat ? + FluidSystem::gasPvt().saturatedViscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure) : + FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, rv, getValue(intQuants.fluidState().Rvw())); + double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx()); + const double phi = getValue(intQuants.porosity()); + const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]]; + double Kh = connection.Kh()* trans_mult; + double Ke = connection.Ke()* trans_mult; + double h = Kh / Ke; + double rw = connection.rw(); + const auto& wdfac = this->well_ecl_.getWDFAC(); + if (wdfac.useConnectionDFactor()) { + double d = connection.dFactor(); + // If a negative d factor is set in COMPDAT individual connection d factors should be used directly. + if (d < 0) + return -d; + // If a positive d factor is set in COMPDAT the connection d factors is treated like a well d factor. + // and thus scaled with the well index + return d * total_tw / connection.CF(); + } else { + + double d = wdfac.getDFactor(rho, mu, Ke, phi, rw, h); + return d * total_tw / connection.CF(); + } + } + + template void WellInterface::