From 20c36d19c19dacd5fcfd604ff596a6a98c3aeedc Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 23 Jun 2017 14:55:56 +0200 Subject: [PATCH] adding computePropertiesForWellConnectionPressures to StandardWell --- opm/autodiff/StandardWell.hpp | 10 +++ opm/autodiff/StandardWell_impl.hpp | 108 ++++++++++++++++++++++++++++ opm/autodiff/WellInterface.hpp | 5 ++ opm/autodiff/WellInterface_impl.hpp | 1 + 4 files changed, 124 insertions(+) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 1e628b488..2c0f44258 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -178,6 +178,7 @@ namespace Opm using WellInterface::well_efficiency_factor_; using WellInterface::active_; using WellInterface::phase_usage_; + using WellInterface::first_perf_; // densities of the fluid in each perforation std::vector perf_densities_; @@ -209,6 +210,15 @@ namespace Opm // TODO: it is also possible to be moved to the base class. EvalWell getQs(const int phase) const; + + // calculate the properties for the well connections + // to calulate the pressure difference between well connections. + void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) const; }; } diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 53d08ba9e..6a2bfb8a5 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -1416,4 +1416,112 @@ namespace Opm } + + + + template + void + StandardWell:: + computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator, + const WellState& xw, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) const + { + const int nperf = numberOfPerforations(); + // TODO: can make this a member? + const int nw = xw.bhp().size(); + const int numComp = numComponents(); + const PhaseUsage& pu = phase_usage_; + b_perf.resize(nperf*numComp); + surf_dens_perf.resize(nperf*numComp); + const int w = indexOfWell(); + + //rs and rv are only used if both oil and gas is present + if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + } + + // Compute the average pressure in each well block + for (int perf = 0; perf < nperf; ++perf) { + const int cell_idx = wellCells()[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + + // TODO: this is another place to show why WellState need to be a vector of WellState. + // TODO: to check why should be perf - 1 + const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1]; + const double p_avg = (xw.perfPress()[perf] + p_above)/2; + const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); + + if (pu.phase_used[BlackoilPhases::Aqua]) { + b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * numComp] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); + if (oilrate > 0) { + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); + double rv = 0.0; + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + + } else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; + if (pu.phase_used[BlackoilPhases::Vapour]) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); + if (gasrate > 0) { + const double oilrate = std::abs(xw.wellRates()[oilpos_well]); + double rs = 0.0; + if (oilrate > 0) { + rs = gasrate / oilrate; + } + rs = std::min(rs, rsmax_perf[perf]); + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); + } + } + + // Surface density. + for (int p = 0; p < pu.num_phases; ++p) { + surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); + } + + // We use cell values for solvent injector + if (has_solvent) { + b_perf[numComp*perf + solventCompIdx] = intQuants.solventInverseFormationVolumeFactor().value(); + surf_dens_perf[numComp*perf + solventCompIdx] = intQuants.solventRefDensity(); + } + } + } + } diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index 214d836da..0891b984f 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -169,6 +169,11 @@ namespace Opm // number of the perforations for this well int number_of_perforations_; + // record the index of the first perforation + // TODO: it might not be needed if we refactor WellState to be a vector + // of states of individual well. + int first_perf_; + // well index for each perforation std::vector well_index_; diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 52188fa0a..a3a5e0b52 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -65,6 +65,7 @@ namespace Opm const int perf_index_begin = wells->well_connpos[index_well]; const int perf_index_end = wells->well_connpos[index_well + 1]; number_of_perforations_ = perf_index_end - perf_index_begin; + first_perf_ = perf_index_begin; well_cell_.resize(number_of_perforations_); std::copy(wells->well_cells + perf_index_begin,