diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 6a5eda68f..31f2e7e8b 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/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& water_vel, std::vector& visc_mult, std::vector& shear_mult); }; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index dea3eeee8..98a1e4249 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -617,6 +617,102 @@ namespace Opm { } } + template + bool + BlackoilPolymerModel::computeShearMultLog( std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult){ + + double refConcentration = polymer_props_ad_.plyshlogRefConc(); + double refViscMult = polymer_props_ad_.viscMult(refConcentration); + + std::vector shear_water_vel = polymer_props_ad_.shearWaterVelocity(); + std::vector shear_vrf = polymer_props_ad_.shearViscosityReductionFactor(); + + std::vector logShearWaterVel; + std::vector 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::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