OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

moving function computeShearMultLog to class PolymerProperties

Browse Source
This commit is contained in:
Kai Bao 2015-06-10 13:36:37 +02:00
parent e44ef196ac
commit 81d9fe7a55
6 changed files with 117 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

95
opm/polymer/PolymerProperties.cpp
View File

@ -20,6 +20,7 @@
#include <config.h> #include <config.h>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/Point2D.hpp>
#include <cmath> #include <cmath>
#include <vector> #include <vector>
#include <opm/core/utility/linearInterpolation.hpp> #include <opm/core/utility/linearInterpolation.hpp>
@ -408,4 +409,98 @@ namespace Opm
dmc_dc = 0.; dmc_dc = 0.;
} }
} }
bool PolymerProperties::computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult) const
{
double refConcentration = plyshlogRefConc();
double refViscMult = viscMult(refConcentration);
std::vector<double> shear_water_vel = shearWaterVelocity();
std::vector<double> shear_vrf = 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 (size_t 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 (size_t 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 (size_t 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]));
size_t 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) {
detail::Point2D lineSegment[2];
lineSegment[0] = detail::Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]};
lineSegment[1] = detail::Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]};
detail::Point2D line[2];
line[0] = detail::Point2D{0, logWaterVelO};
line[1] = detail::Point2D{logWaterVelO, 0};
detail::Point2D intersectionPoint;
bool foundIntersection = detail::Point2D::findIntersection(lineSegment, line, intersectionPoint);
if (foundIntersection) {
shear_mult[i] = std::exp(intersectionPoint.getY());
} else {
std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl;
return false; // failed in finding the solution.
}
} else {
if (std::abs(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;
}
} }

4
opm/polymer/PolymerProperties.hpp
View File

@ -326,6 +326,9 @@ namespace Opm
void computeMcBoth(const double& c, double& mc, void computeMcBoth(const double& c, double& mc,
double& dmc_dc, bool if_with_der) const; double& dmc_dc, bool if_with_der) const;
/// 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) const;
private: private:
double c_max_; double c_max_;
double mix_param_; double mix_param_;
@ -366,6 +369,7 @@ namespace Opm
double& deff_relperm_wat_ds, double& deff_relperm_wat_ds,
double& deff_relperm_wat_dc, double& deff_relperm_wat_dc,
bool if_with_der) const; bool if_with_der) const;
}; };
} // namespace Opm } // namespace Opm

3
opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
View File

@ -285,9 +285,6 @@ namespace Opm {
int nc, int nc,
int nw) const; int nw) const;
/// 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);
/// Computing the water velocity without shear-thinning for the cell faces. /// Computing the water velocity without shear-thinning for the cell faces.
/// The water velocity will be used for shear-thinning calculation. /// The water velocity will be used for shear-thinning calculation.
void computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr, void computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,

101
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
View File

@ -25,7 +25,6 @@
#define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp> #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/polymer/Point2D.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp> #include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp> #include <opm/autodiff/AutoDiffHelpers.hpp>
@ -288,7 +287,7 @@ namespace Opm {
std::vector<double> visc_mult; std::vector<double> visc_mult;
computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult); computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult);
if(!computeShearMultLog(water_vel, visc_mult, shear_mult_faces_)) { if(!polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_)) {
// std::cerr << " failed in calculating the shear-multiplier " << std::endl; // std::cerr << " failed in calculating the shear-multiplier " << std::endl;
OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. "); OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
} }
@ -578,7 +577,7 @@ namespace Opm {
computeWaterShearVelocityWells(state, well_state, aliveWells, water_vel_wells, visc_mult_wells); computeWaterShearVelocityWells(state, well_state, aliveWells, water_vel_wells, visc_mult_wells);
if (!computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_)) { if (!polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_)) {
// std::cout << " failed in calculating the shear factors for wells " << std::endl; // std::cout << " failed in calculating the shear factors for wells " << std::endl;
OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells "); OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells ");
} }
@ -893,102 +892,6 @@ namespace Opm {
extraAddWellEq(state, xw, cq_ps, cmix_s, cqt_is, well_cells); extraAddWellEq(state, xw, cq_ps, cmix_s, cqt_is, well_cells);
} }
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) {
detail::Point2D lineSegment[2];
lineSegment[0] = detail::Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]};
lineSegment[1] = detail::Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]};
detail::Point2D line[2];
line[0] = detail::Point2D{0, logWaterVelO};
line[1] = detail::Point2D{logWaterVelO, 0};
detail::Point2D intersectionPoint;
bool foundIntersection = detail::Point2D::findIntersection(lineSegment, line, intersectionPoint);
if (foundIntersection) {
shear_mult[i] = std::exp(intersectionPoint.getY());
} else {
std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl;
return false; // failed in finding the solution.
}
} else {
if (std::abs(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;
}
template<class Grid> template<class Grid>
void void
BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr, BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,

8
opm/polymer/fullyimplicit/PolymerPropsAd.cpp
View File

@ -317,4 +317,12 @@ namespace Opm {
return krw / rk; return krw / rk;
} }
bool
PolymerPropsAd::computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult) const
{
return polymer_props_.computeShearMultLog(water_vel, visc_mult, shear_mult);
}
}// namespace Opm }// namespace Opm

8
opm/polymer/fullyimplicit/PolymerPropsAd.hpp
View File

@ -130,6 +130,14 @@ namespace Opm {
ADB ADB
effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw) const; effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw) const;
/// \param[in] water_vel Array of the n values of water velocity or shear rate.
/// \param[in] visc_mult Array of the n values of the viscosity multiplier from PLYVISC table.
/// \parma[out] shear_mult Array of the n values of calculated shear multiplier with PLYSHLOG keyword.
/// \return TRUE if the calculation of shear multiplier is sucessful,
/// FALSE if the calculation of shear multplier is failed.
bool computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult) const;
private: private:
const PolymerProperties& polymer_props_; const PolymerProperties& polymer_props_;
}; };