From eeb1b7e36c9ae3c30ba5d262fa623dc0eea2e10b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 20 Apr 2021 15:42:27 +0200 Subject: [PATCH] Initialize empty producers using the mobility ratio and the transmissbility ratio --- opm/simulators/wells/StandardWell.hpp | 5 +- opm/simulators/wells/StandardWell_impl.hpp | 75 ++++++++++++++++------ 2 files changed, 57 insertions(+), 23 deletions(-) diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 75e67680d..8af52ede7 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -458,7 +458,8 @@ namespace Opm void computeConnectionPressureDelta(); - void computeWellConnectionDensitesPressures(const WellState& well_state, + void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator, + const WellState& well_state, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, @@ -468,7 +469,7 @@ namespace Opm void computeAccumWell(); void computeWellConnectionPressures(const Simulator& ebosSimulator, - const WellState& well_state); + const WellState& well_state); void computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index dc40c5c37..ea4ba0ca2 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1854,20 +1854,25 @@ namespace Opm } } else { assert(this->isProducer()); - // Using the preferred phase to decide the mix initialization. - switch (this->wellEcl().getPreferredPhase()) { - case Phase::OIL: - mix[FluidSystem::oilCompIdx] = 1.0; - break; - case Phase::GAS: - mix[FluidSystem::gasCompIdx] = 1.0; - break; - case Phase::WATER: - mix[FluidSystem::waterCompIdx] = 1.0; - break; - default: - // No others supported. - break; + // For the frist perforation without flow we use the preferred phase to decide the mix initialization. + if (perf == 0) { // + switch (this->wellEcl().getPreferredPhase()) { + case Phase::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case Phase::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case Phase::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + default: + // No others supported. + break; + } + // For the rest of the perforation without flow we use mix from the above perforation. + } else { + mix = x; } } @@ -2088,12 +2093,11 @@ namespace Opm - - template void StandardWell:: - computeWellConnectionDensitesPressures(const WellState& well_state, + computeWellConnectionDensitesPressures(const Simulator& ebosSimulator, + const WellState& well_state, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, @@ -2103,7 +2107,6 @@ namespace Opm const int nperf = number_of_perforations_; const int np = number_of_phases_; std::vector perfRates(b_perf.size(),0.0); - for (int perf = 0; perf < nperf; ++perf) { for (int comp = 0; comp < np; ++comp) { perfRates[perf * num_components_ + comp] = well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(comp)]; @@ -2113,6 +2116,38 @@ namespace Opm } } + // for producers where all perforations have zero rate we + // approximate the perforation mixture using the mobility ratio + // and weight the perforations using the well transmissibility. + bool all_zero = std::all_of(perfRates.begin(), perfRates.end(), [](double val) { return val == 0.0; }); + if ( all_zero && this->isProducer() ) { + double total_tw = 0; + for (int perf = 0; perf < nperf; ++perf) { + total_tw += well_index_[perf]; + } + for (int perf = 0; perf < nperf; ++perf) { + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + const double well_tw_fraction = well_index_[perf] / total_tw; + double total_mobility = 0.0; + for (int p = 0; p < np; ++p) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(p); + total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value(); + } + if(has_solvent) { + total_mobility += intQuants.solventInverseFormationVolumeFactor().value() * intQuants.solventMobility().value(); + } + for (int p = 0; p < np; ++p) { + int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(p); + perfRates[perf * num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility; + } + if(has_solvent) { + perfRates[perf * num_components_ + contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility; + } + } + } + computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); computeConnectionPressureDelta(); @@ -2137,7 +2172,7 @@ namespace Opm std::vector rvmax_perf; std::vector surf_dens_perf; computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - computeWellConnectionDensitesPressures(well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); } @@ -2179,8 +2214,6 @@ namespace Opm - - template void StandardWell::