diff --git a/examples/flow_polymer.cpp b/examples/flow_polymer.cpp index 20e36fd47..0c6f912b5 100644 --- a/examples/flow_polymer.cpp +++ b/examples/flow_polymer.cpp @@ -281,6 +281,7 @@ try deck->hasKeyword("VAPOIL"), polymer, deck->hasKeyword("PLYSHLOG"), + deck->hasKeyword("SHRATE"), eclipseState, outputWriter, deck, diff --git a/opm/polymer/PolymerProperties.cpp b/opm/polymer/PolymerProperties.cpp index 87b8e8b11..191c78a93 100644 --- a/opm/polymer/PolymerProperties.cpp +++ b/opm/polymer/PolymerProperties.cpp @@ -101,6 +101,22 @@ namespace Opm return plyshlog_ref_temp_; } + bool PolymerProperties::hasPlyshlog() const + { + return has_plyshlog_; + } + + bool PolymerProperties::hasShrate() const + { + return has_shrate_; + } + + + double PolymerProperties::shrate() const + { + return shrate_; + } + double PolymerProperties::shearVrf(const double velocity) const diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index 7feff246c..ba1190c4b 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -25,6 +25,7 @@ #include #include +#include namespace Opm @@ -154,9 +155,10 @@ namespace Opm c_vals_ads_ = plyadsTable.getPolymerConcentrationColumn(); ads_vals_ = plyadsTable.getAdsorbedPolymerColumn(); - plyshlog_ = deck->hasKeyword("PLYSHLOG"); + has_plyshlog_ = deck->hasKeyword("PLYSHLOG"); + has_shrate_ = deck->hasKeyword("SHRATE"); - if (plyshlog_) { + if (has_plyshlog_) { // Assuming NTPVT == 1 always due to the limitation of the parser const auto& plyshlogTable = eclipseState->getPlyshlogTables()[0]; @@ -165,7 +167,21 @@ namespace Opm // do the unit version here for the water_vel_vals_ Opm::UnitSystem unitSystem = *deck->getActiveUnitSystem(); - double siFactor = unitSystem.parse("Length/Time")->getSIScaling(); + double siFactor; + if (has_shrate_) { + siFactor = unitSystem.parse("1/Time")->getSIScaling(); + DeckKeywordConstPtr shrateKeyword = deck->getKeyword("SHRATE"); + std::vector shrate = shrateKeyword->getSIDoubleData(); + if (shrate.size() == 1) { + shrate_ = shrate[0]; + } else if (shrate.size() == 0) { + shrate_ = 4.8; // default value + } else { + OPM_THROW(std::logic_error, "Only NTPVT == 1 is allowed for SHRATE keyword now !\n"); + } + } else { + siFactor = unitSystem.parse("Length/Time")->getSIScaling(); + } for (size_t i = 0; i < water_vel_vals_.size(); ++i) { water_vel_vals_[i] *= siFactor; @@ -204,6 +220,9 @@ namespace Opm int adsIndex() const; + /// indicate whehter PLYSHLOG is specified + bool hasPlyshlog() const; + /// the water velocity or water shear rate in PLYSHLOG table const std::vector& shearWaterVelocity() const; @@ -225,6 +244,12 @@ namespace Opm /// the reference temperature in PLYSHLOG double plyshlogRefTemp() const; + /// indicate whether SHRATE keyword is specified + bool hasShrate() const; + + /// the value of SHRATE + double shrate() const; + double shearVrf(const double velocity) const; double shearVrfWithDer(const double velocity, double& der) const; @@ -337,7 +362,12 @@ namespace Opm double res_factor_; double c_max_ads_; - bool plyshlog_; + bool has_plyshlog_; + bool has_shrate_; + // Assuming NTPVT == 1 always due to the limitation of the parser + // only one SHRATE value + // TODO: to be extended later when parser is improved. + double shrate_; AdsorptionBehaviour ads_index_; std::vector c_vals_visc_; std::vector visc_mult_vals_; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 5521800d6..da17671db 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -56,17 +56,21 @@ namespace Opm { /// Construct the model. It will retain references to the /// arguments of this functions, and they are expected to /// remain in scope for the lifetime of the solver. - /// \param[in] param parameters - /// \param[in] grid grid data structure - /// \param[in] fluid fluid properties - /// \param[in] geo rock properties - /// \param[in] rock_comp_props if non-null, rock compressibility properties - /// \param[in] wells well structure - /// \param[in] linsolver linear solver - /// \param[in] has_disgas turn on dissolved gas - /// \param[in] has_vapoil turn on vaporized oil feature - /// \param[in] has_polymer turn on polymer feature - /// \param[in] has_plyshlog true when PLYSHLOG keyword available + /// \param[in] param parameters + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] geo rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] wells well structure + /// \param[in] linsolver linear solver + /// \param[in] has_disgas turn on dissolved gas + /// \param[in] has_vapoil turn on vaporized oil feature + /// \param[in] has_polymer turn on polymer feature + /// \param[in] has_plyshlog true when PLYSHLOG keyword available + /// \param[in] has_shrate true when PLYSHLOG keyword available + /// \param[in] wells_rep_radius representative radius of well perforations during shear effects calculation + /// \param[in] wells_perf_length perforation length for well perforations + /// \param[in] wells_bore_diameter wellbore diameters for well performations /// \param[in] terminal_output request output to cout/cerr BlackoilPolymerModel(const typename Base::ModelParameters& param, const Grid& grid, @@ -80,8 +84,10 @@ namespace Opm { const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, + const bool has_shrate, const std::vector& wells_rep_radius, const std::vector& wells_perf_length, + const std::vector& wells_bore_diameter, const bool terminal_output); /// Called once before each time step. @@ -136,6 +142,7 @@ namespace Opm { const PolymerPropsAd& polymer_props_ad_; const bool has_polymer_; const bool has_plyshlog_; + const bool has_shrate_; const int poly_pos_; V cmax_; @@ -143,6 +150,8 @@ namespace Opm { // to be used in shear-thinning computation. std::vector wells_rep_radius_; std::vector wells_perf_length_; + // wellbore diameters + std::vector wells_bore_diameter_; // shear-thinning factor for cell faces std::vector shear_mult_faces_; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 0a650e0c7..ccd906bc0 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -86,17 +86,21 @@ namespace Opm { const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, + const bool has_shrate, const std::vector& wells_rep_radius, const std::vector& wells_perf_length, + const std::vector& wells_bore_diameter, const bool terminal_output) : Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver, has_disgas, has_vapoil, terminal_output), polymer_props_ad_(polymer_props_ad), has_polymer_(has_polymer), has_plyshlog_(has_plyshlog), + has_shrate_(has_shrate), poly_pos_(detail::polymerPos(fluid.phaseUsage())), wells_rep_radius_(wells_rep_radius), - wells_perf_length_(wells_perf_length) + wells_perf_length_(wells_perf_length), + wells_bore_diameter_(wells_bore_diameter) { if (has_polymer_) { if (!active_[Water]) { @@ -752,8 +756,6 @@ namespace Opm { std::vector& water_vel, std::vector& visc_mult) { - std::vector b_faces; - const int phase = fluid_.phaseUsage().phase_pos[Water]; // water position const int canonicalPhaseIdx = canph_[phase]; @@ -791,36 +793,73 @@ namespace Opm { size_t nface = visc_mult_faces.size(); visc_mult.resize(nface); - std::copy(&(visc_mult_faces[0]), &(visc_mult_faces[0]) + nface, visc_mult.begin()); + std::copy(visc_mult_faces.data(), visc_mult_faces.data() + nface, visc_mult.begin()); rq_[ phase ].mflux = upwind.select(b * mob) * (transi * dh); const auto& b_faces_adb = upwind.select(b); - b_faces.resize(b_faces_adb.size()); - std::copy(&(b_faces_adb.value()[0]), &(b_faces_adb.value()[0]) + b_faces_adb.size(), b_faces.begin()); + std::vector b_faces(b_faces_adb.value().data(), b_faces_adb.value().data() + b_faces_adb.size()); const auto& internal_faces = ops_.internal_faces; std::vector internal_face_areas; internal_face_areas.resize(internal_faces.size()); - for (int i = 0; i < internal_faces.size(); ++i) { + for (size_t i = 0; i < internal_faces.size(); ++i) { internal_face_areas[i] = grid_.face_areas[internal_faces[i]]; } const ADB phi = Opm::AutoDiffBlock::constant(Eigen::Map(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1)); const ADB phiavg_adb = ops_.caver * phi; - std::vector phiavg; - phiavg.resize(phiavg_adb.size()); - std::copy(&(phiavg_adb.value()[0]), &(phiavg_adb.value()[0]) + phiavg_adb.size(), phiavg.begin()); + std::vector phiavg(phiavg_adb.value().data(), phiavg_adb.value().data() + phiavg_adb.size()); water_vel.resize(nface); - std::copy(&(rq_[0].mflux.value()[0]), &(rq_[0].mflux.value()[0]) + nface, water_vel.begin()); + std::copy(rq_[0].mflux.value().data(), rq_[0].mflux.value().data() + nface, water_vel.begin()); for (size_t i = 0; i < nface; ++i) { water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]); } + + // for SHRATE keyword treatment + if (has_shrate_) { + + // get the upwind water saturation + const Opm::PhaseUsage pu = fluid_.phaseUsage(); + const ADB& sw = state.saturation[pu.phase_pos[ Water ]]; + const ADB& sw_upwind_adb = upwind.select(sw); + std::vector sw_upwind(sw_upwind_adb.value().data(), sw_upwind_adb.value().data() + sw_upwind_adb.size()); + + // get the absolute permeability for the faces + std::vector perm; + perm.resize(transi.size()); + + for (size_t i = 0; i < transi.size(); ++i) { + perm[i] = transi[i] / internal_faces[i]; + } + + // get the upwind krw_eff + const ADB& krw_adb = upwind.select(krw_eff); + std::vector krw_upwind(krw_adb.value().data(), krw_adb.value().data() + krw_adb.size()); + + const double& shrate_const = polymer_props_ad_.shrate(); + + const double epsilon = std::numeric_limits::epsilon(); + // std::cout << "espilon is " << epsilon << std::endl; + // std::cin.ignore(); + + for (size_t i = 0; i < water_vel.size(); ++i) { + // assuming only when upwinding water saturation is not zero + // there will be non-zero water velocity + if (std::abs(water_vel[i]) < epsilon) { + continue; + } + + water_vel[i] *= shrate_const * std::sqrt(phiavg[i] / (perm[i] * sw_upwind[i] * krw_upwind[i])); + + } + } + } template @@ -835,7 +874,7 @@ namespace Opm { const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); water_vel_wells.resize(cq_sw.size()); - std::copy(&(cq_sw.value()[0]), &(cq_sw.value()[0]) + cq_sw.size(), water_vel_wells.begin()); + std::copy(cq_sw.value().data(), cq_sw.value().data() + cq_sw.size(), water_vel_wells.begin()); const V& polymer_conc = state.concentration.value(); @@ -843,7 +882,7 @@ namespace Opm { V visc_mult_wells_v = subset(visc_mult_cells, well_cells); visc_mult_wells.resize(visc_mult_wells_v.size()); - std::copy(&(visc_mult_wells_v[0]), &(visc_mult_wells_v[0]) + visc_mult_wells_v.size(), visc_mult_wells.begin()); + std::copy(visc_mult_wells_v.data(), visc_mult_wells_v.data() + visc_mult_wells_v.size(), visc_mult_wells.begin()); const int water_pos = fluid_.phaseUsage().phase_pos[Water]; ADB b_perfcells = subset(rq_[water_pos].b, well_cells); @@ -872,13 +911,9 @@ namespace Opm { const ADB phi = Opm::AutoDiffBlock::constant(Eigen::Map(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1)); const ADB phi_wells_adb = subset(phi, well_cells); - std::vector phi_wells; - phi_wells.resize(phi_wells_adb.size()); - std::copy(&(phi_wells_adb.value()[0]), &(phi_wells_adb.value()[0]) + phi_wells_adb.size(), phi_wells.begin()); + std::vector phi_wells(phi_wells_adb.value().data(), phi_wells_adb.value().data() + phi_wells_adb.size()); - std::vector b_wells; - b_wells.resize(b_perfcells.size()); - std::copy(&(b_perfcells.value()[0]), &(b_perfcells.value()[0]) + b_perfcells.size(), b_wells.begin()); + std::vector b_wells(b_perfcells.value().data(), b_perfcells.value().data() + b_perfcells.size()); for (size_t i = 0; i < water_vel_wells.size(); ++i) { water_vel_wells[i] = b_wells[i] * water_vel_wells[i] / (phi_wells[i] * 2. * M_PI * wells_rep_radius_[i] * wells_perf_length_[i]); @@ -886,6 +921,14 @@ namespace Opm { // Although this formulation works perfectly with the tests compared with other formulations } + // for SHRATE treatment + if (has_shrate_) { + const double& shrate_const = polymer_props_ad_.shrate(); + for (size_t i = 0; i < water_vel_wells.size(); ++i) { + water_vel_wells[i] = shrate_const * water_vel_wells[i] / wells_bore_diameter_[i]; + } + } + return; } diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp index 37dc56ba4..cad97caa2 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -102,6 +102,16 @@ namespace Opm { return polymer_props_.plyshlogRefTemp(); } + double PolymerPropsAd::shrate() const + { + if (polymer_props_.hasShrate()) { + return polymer_props_.shrate(); + } else { + OPM_THROW(std::logic_error, "the SHRATE keyword is not specified while requested \n"); + } + } + + double PolymerPropsAd::viscMult(double c) const { diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index 7bf90ea3c..cc54ce282 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -62,6 +62,9 @@ namespace Opm { /// \ return The reference temperature in PLYSHLOG keyword double plyshlogRefTemp() const; + /// \ return the value of SHRATE + double shrate() const; + double viscMult(double c) const; // multipler interpolated from PLYVISC table typedef AutoDiffBlock ADB; diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp index 600dd9768..49e7be303 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -117,6 +117,7 @@ namespace Opm const bool vapoil, const bool polymer, const bool plyshlog, + const bool shrate, std::shared_ptr eclipse_state, BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, @@ -136,10 +137,13 @@ namespace Opm bool has_polymer_; // flag for PLYSHLOG keyword bool has_plyshlog_; + // flag for SHRATE keyword + bool has_shrate_; DeckConstPtr deck_; std::vector wells_rep_radius_; std::vector wells_perf_length_; + std::vector wells_bore_diameter_; // generate the mapping from Cartesian grid cells to global compressed cells, // copied from opm-core, to be used in function computeRepRadiusPerfLength() @@ -147,13 +151,15 @@ namespace Opm setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed); // calculate the representative radius and length for for well peforations + // and store the wellbore diameters // it will be used in the shear-thinning calcluation only. void computeRepRadiusPerfLength(const Opm::EclipseStateConstPtr eclipseState, const size_t timeStep, const GridT& grid, std::vector& wells_rep_radius, - std::vector& wells_perf_length); + std::vector& wells_perf_length, + std::vector& wells_bore_diameter); }; diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 8387819b6..d3e01b984 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -35,6 +35,7 @@ namespace Opm const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, + const bool has_shrate, std::shared_ptr eclipse_state, BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, @@ -54,6 +55,7 @@ namespace Opm , polymer_props_(polymer_props) , has_polymer_(has_polymer) , has_plyshlog_(has_plyshlog) + , has_shrate_(has_shrate) , deck_(deck) { } @@ -78,8 +80,10 @@ namespace Opm BaseType::has_vapoil_, has_polymer_, has_plyshlog_, + has_shrate_, wells_rep_radius_, wells_perf_length_, + wells_bore_diameter_, BaseType::terminal_output_)); if (!BaseType::threshold_pressures_by_face_.empty()) { @@ -117,7 +121,7 @@ namespace Opm polymer_inflow_c); well_state.polymerInflow() = polymer_inflow_c; - computeRepRadiusPerfLength(BaseType::eclipse_state_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_); + computeRepRadiusPerfLength(BaseType::eclipse_state_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_, wells_bore_diameter_); } @@ -146,7 +150,8 @@ namespace Opm const size_t timeStep, const GridT& grid, std::vector& wells_rep_radius, - std::vector& wells_perf_length) + std::vector& wells_perf_length, + std::vector& wells_bore_diameter) { // TODO, the function does not work for parallel running @@ -167,9 +172,11 @@ namespace Opm wells_rep_radius.clear(); wells_perf_length.clear(); + wells_bore_diameter.clear(); wells_rep_radius.reserve(n_perf); wells_perf_length.reserve(n_perf); + wells_bore_diameter.reserve(n_perf); std::map cartesian_to_compressed; @@ -240,6 +247,7 @@ namespace Opm double repR = std::sqrt(re * radius); wells_rep_radius.push_back(repR); wells_perf_length.push_back(perf_length); + wells_bore_diameter.push_back(2. * radius); } } else { if (completion->getState() != WellCompletion::SHUT) {