adding the function to compute shear-thinning effect based on PLYSHLOG

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -270,7 +270,8 @@ namespace Opm {
         /// return true, if found.
         bool findIntersection (Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point);
 
-
+        /// Computing the shear multiplier based on the water velocity/shear rate with PLYSHLOG keyword
+        bool computeShearMultLog( std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult);
 
     };
 

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -617,6 +617,102 @@ namespace Opm {
         }
     }
 
+    template<class Grid>
+    bool
+    BlackoilPolymerModel<Grid>::computeShearMultLog( std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult){
+
+        double refConcentration = polymer_props_ad_.plyshlogRefConc();
+        double refViscMult = polymer_props_ad_.viscMult(refConcentration);
+
+        std::vector<double> shear_water_vel = polymer_props_ad_.shearWaterVelocity();
+        std::vector<double> shear_vrf = polymer_props_ad_.shearViscosityReductionFactor();
+
+        std::vector<double> logShearWaterVel;
+        std::vector<double> logShearVRF;
+
+        logShearWaterVel.resize(shear_water_vel.size());
+        logShearVRF.resize(shear_water_vel.size());
+
+        // converting the table using the reference condition
+        for(int i = 0; i < shear_vrf.size(); ++i){
+            shear_vrf[i] = (refViscMult * shear_vrf[i] - 1.) / (refViscMult - 1);
+            logShearWaterVel[i] = std::log(shear_water_vel[i]);
+        }
+
+        shear_mult.resize(water_vel.size());
+
+        // the mimum velocity to apply the shear-thinning
+        const double minShearVel = shear_water_vel[0];
+        const double maxShearVel = shear_water_vel.back();
+        const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
+
+        for(int i = 0; i < water_vel.size(); ++i){
+
+            if( visc_mult[i] - 1. < epsilon || std::abs(water_vel[i]) < minShearVel ) {
+                 shear_mult[i] = 1.0;
+                 continue;
+            }
+
+            for(int j = 0; j < shear_vrf.size(); ++j){
+                logShearVRF[j] = (1 + (visc_mult[i] - 1.0) * shear_vrf[j]) / visc_mult[i];
+                logShearVRF[j] = std::log(logShearVRF[j]);
+            }
+
+            // const double logWaterVelO = std::log(water_vel[i]);
+            const double logWaterVelO = std::log(std::abs(water_vel[i]));
+
+            int iIntersection; // finding the intersection on the iIntersectionth table segment
+            bool foundSegment = false;
+
+            for(iIntersection = 0; iIntersection < shear_vrf.size() - 1; ++iIntersection){
+
+                double temp1 = logShearVRF[iIntersection] + logShearWaterVel[iIntersection] - logWaterVelO;
+                double temp2 = logShearVRF[iIntersection + 1] + logShearWaterVel[iIntersection + 1] - logWaterVelO;
+
+                // ignore the cases the temp1 or temp2 is zero first for simplicity.
+                // several more complicated cases remain to be implemented.
+                if( temp1 * temp2 < 0.){
+                    foundSegment = true;
+                    break;
+                }
+            }
+
+            if(foundSegment == true){
+                Point2D lineSegment[2];
+                lineSegment[0] = Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]};
+                lineSegment[1] = Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]};
+
+                Point2D line[2];
+                line[0] = Point2D{0, logWaterVelO};
+                line[1] = Point2D{logWaterVelO, 0};
+
+                Point2D intersectionPoint;
+
+                bool foundIntersection = findIntersection(lineSegment, line, intersectionPoint);
+
+                if(foundIntersection){
+                    shear_mult[i] = std::exp(intersectionPoint.y);
+                }else{
+                    std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl;
+                    return false; // failed in finding the solution.
+                }
+            }else{
+                if (water_vel[i] < maxShearVel) {
+                    std::cout << " the veclocity is " << water_vel[i] << std::endl;
+                    std::cout << " max shear velocity is " << maxShearVel << std::endl;
+                    std::cerr << " something wrong happend in finding segment" << std::endl;
+                    return false;
+                } else {
+                    shear_mult[i] = std::exp(logShearVRF.back());
+                }
+            }
+
+        }
+
+        return true;
+    }
+
+
 } // namespace Opm
 
 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED