Merge pull request #133 from GitPaean/Adding_SHRATE_RELATED

Adding shear-thinning effects based on keyword SHRATE

examples/flow_polymer.cpp

@@ -281,6 +281,7 @@ try
                   deck->hasKeyword("VAPOIL"),
                   polymer,
                   deck->hasKeyword("PLYSHLOG"),
+                  deck->hasKeyword("SHRATE"),
                   eclipseState,
                   outputWriter,
                   deck,

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

opm/polymer/PolymerProperties.hpp

@@ -25,6 +25,7 @@
 
 #include <cmath>
 #include <vector>
+#include <opm/core/utility/ErrorMacros.hpp>
 
 
 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<double> 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<double>& 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<double> c_vals_visc_;
         std::vector<double> visc_mult_vals_;

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -67,6 +67,10 @@ namespace Opm {
         /// \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<double>&              wells_rep_radius,
                              const std::vector<double>&              wells_perf_length,
+                             const std::vector<double>&              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<double> wells_rep_radius_;
         std::vector<double> wells_perf_length_;
+        // wellbore diameters
+        std::vector<double> wells_bore_diameter_;
 
         // shear-thinning factor for cell faces
         std::vector<double> shear_mult_faces_;

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<double>&              wells_rep_radius,
                                                      const std::vector<double>&              wells_perf_length,
+                                                     const std::vector<double>&              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<double>& water_vel, std::vector<double>& visc_mult)
     {
 
-        std::vector<double> 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::vector<double> b_faces(b_faces_adb.value().data(), b_faces_adb.value().data() + b_faces_adb.size());
-        std::copy(&(b_faces_adb.value()[0]), &(b_faces_adb.value()[0]) + b_faces_adb.size(), b_faces.begin());
 
         const auto& internal_faces = ops_.internal_faces;
 
         std::vector<double> 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<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
         const ADB phiavg_adb  = ops_.caver * phi;
 
-        std::vector<double> phiavg;
+        std::vector<double> phiavg(phiavg_adb.value().data(), phiavg_adb.value().data() + phiavg_adb.size());
-        phiavg.resize(phiavg_adb.size());
-        std::copy(&(phiavg_adb.value()[0]), &(phiavg_adb.value()[0]) + phiavg_adb.size(), phiavg.begin());
 
         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<double> 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<double> 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<double> 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<double>::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<class Grid>
@@ -835,7 +874,7 @@ namespace Opm {
         const std::vector<int> 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<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
         const ADB phi_wells_adb = subset(phi, well_cells);
 
-        std::vector<double> phi_wells;
+        std::vector<double> phi_wells(phi_wells_adb.value().data(), phi_wells_adb.value().data() + phi_wells_adb.size());
-        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<double> b_wells;
+        std::vector<double> b_wells(b_perfcells.value().data(), b_perfcells.value().data() + b_perfcells.size());
-        b_wells.resize(b_perfcells.size());
-        std::copy(&(b_perfcells.value()[0]), &(b_perfcells.value()[0]) + b_perfcells.size(), b_wells.begin());
 
         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;
     }
 

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
     {

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<double> ADB;

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<EclipseState> 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<double> wells_rep_radius_;
         std::vector<double> wells_perf_length_;
+        std::vector<double> 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<int,int>& 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<double>&            wells_rep_radius,
-                                   std::vector<double>&            wells_perf_length);
+                                   std::vector<double>&            wells_perf_length,
+                                   std::vector<double>&            wells_bore_diameter);
 
     };
 

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<EclipseState> 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<double>&            wells_rep_radius,
-                               std::vector<double>&            wells_perf_length)
+                               std::vector<double>&            wells_perf_length,
+                               std::vector<double>&            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<int,int> 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) {