mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #124 from GitPaean/Adding_PLYSHLOG_RELATED
Adding shear-thinning calculation with PLYSHLOG
This commit is contained in:
Atgeirr Flø Rasmussen
commit
eabe1986e6
110
opm/polymer/Point2D.hpp
Normal file
110
opm/polymer/Point2D.hpp
Normal file
@ -0,0 +1,110 @@
|
||||
/*
|
||||
Copyright 2015 Statoil ASA.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_POINT2D_HEADER_INCLUDED
|
||||
#define OPM_POINT2D_HEADER_INCLUDED
|
||||
|
||||
namespace Opm {
|
||||
|
||||
|
||||
|
||||
namespace detail {
|
||||
|
||||
|
||||
class Point2D
|
||||
{
|
||||
public:
|
||||
Point2D(const double xi, const double yi)
|
||||
: x_(xi),
|
||||
y_(yi) {
|
||||
|
||||
}
|
||||
|
||||
Point2D()
|
||||
{
|
||||
}
|
||||
|
||||
const double getX() const
|
||||
{
|
||||
return x_;
|
||||
}
|
||||
|
||||
const double getY() const
|
||||
{
|
||||
return y_;
|
||||
}
|
||||
|
||||
void setX(const double x)
|
||||
{
|
||||
x_ = x;
|
||||
}
|
||||
|
||||
void setY(const double y)
|
||||
{
|
||||
y_ = y;
|
||||
}
|
||||
|
||||
/// Finding the intersection point of a line segment and a line.
|
||||
/// return true, if found.
|
||||
static bool findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point)
|
||||
{
|
||||
|
||||
const double x1 = line_segment1[0].getX();
|
||||
const double y1 = line_segment1[0].getY();
|
||||
const double x2 = line_segment1[1].getX();
|
||||
const double y2 = line_segment1[1].getY();
|
||||
|
||||
const double x3 = line2[0].getX();
|
||||
const double y3 = line2[0].getY();
|
||||
const double x4 = line2[1].getX();
|
||||
const double y4 = line2[1].getY();
|
||||
|
||||
const double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
|
||||
|
||||
if (d == 0.) {
|
||||
return false;
|
||||
}
|
||||
|
||||
const double x = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d;
|
||||
const double y = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d;
|
||||
|
||||
if (x >= std::min(x1,x2) && x <= std::max(x1,x2)) {
|
||||
intersection_point.setX(x);
|
||||
intersection_point.setY(y);
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
double x_;
|
||||
double y_;
|
||||
|
||||
}; // class Point2D
|
||||
|
||||
|
||||
} // namespace detail
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_POINT2D_HEADER_INCLUDED
|
||||
|
||||
@ -20,6 +20,7 @@
|
||||
#include <config.h>
|
||||
|
||||
#include <opm/polymer/PolymerProperties.hpp>
|
||||
#include <opm/polymer/Point2D.hpp>
|
||||
#include <cmath>
|
||||
#include <vector>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
@ -69,6 +70,38 @@ namespace Opm
|
||||
return water_vel_vals_;
|
||||
}
|
||||
|
||||
const std::vector<double>&
|
||||
PolymerProperties::shearViscosityReductionFactor() const
|
||||
{
|
||||
return shear_vrf_vals_;
|
||||
}
|
||||
|
||||
double PolymerProperties:: plyshlogRefConc() const
|
||||
{
|
||||
return plyshlog_ref_conc_;
|
||||
}
|
||||
|
||||
bool PolymerProperties::hasPlyshlogRefSalinity() const
|
||||
{
|
||||
return has_plyshlog_ref_salinity_;
|
||||
}
|
||||
|
||||
bool PolymerProperties::hasPlyshlogRefTemp() const
|
||||
{
|
||||
return has_plyshlog_ref_temp_;
|
||||
}
|
||||
|
||||
double PolymerProperties::plyshlogRefSalinity() const
|
||||
{
|
||||
return plyshlog_ref_salinity_;
|
||||
}
|
||||
|
||||
double PolymerProperties::plyshlogRefTemp() const
|
||||
{
|
||||
return plyshlog_ref_temp_;
|
||||
}
|
||||
|
||||
|
||||
double
|
||||
PolymerProperties::shearVrf(const double velocity) const
|
||||
{
|
||||
@ -376,4 +409,98 @@ namespace Opm
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
@ -153,6 +153,41 @@ namespace Opm
|
||||
|
||||
c_vals_ads_ = plyadsTable.getPolymerConcentrationColumn();
|
||||
ads_vals_ = plyadsTable.getAdsorbedPolymerColumn();
|
||||
|
||||
plyshlog_ = deck->hasKeyword("PLYSHLOG");
|
||||
|
||||
if (plyshlog_) {
|
||||
// Assuming NTPVT == 1 always due to the limitation of the parser
|
||||
const auto& plyshlogTable = eclipseState->getPlyshlogTables()[0];
|
||||
|
||||
water_vel_vals_ = plyshlogTable.getWaterVelocityColumn();
|
||||
shear_vrf_vals_ = plyshlogTable.getShearMultiplierColumn();
|
||||
|
||||
// do the unit version here for the water_vel_vals_
|
||||
Opm::UnitSystem unitSystem = *deck->getActiveUnitSystem();
|
||||
double siFactor = unitSystem.parse("Length/Time")->getSIScaling();
|
||||
|
||||
for (size_t i = 0; i < water_vel_vals_.size(); ++i) {
|
||||
water_vel_vals_[i] *= siFactor;
|
||||
}
|
||||
|
||||
|
||||
plyshlog_ref_conc_ = plyshlogTable.getRefPolymerConcentration();
|
||||
|
||||
if (plyshlogTable.hasRefSalinity()) {
|
||||
has_plyshlog_ref_salinity_ = true;
|
||||
plyshlog_ref_salinity_ = plyshlogTable.getRefSalinity();
|
||||
} else {
|
||||
has_plyshlog_ref_salinity_ = false;
|
||||
}
|
||||
|
||||
if (plyshlogTable.hasRefTemperature()) {
|
||||
has_plyshlog_ref_temp_ = true;
|
||||
plyshlog_ref_temp_ = plyshlogTable.getRefTemperature();
|
||||
} else {
|
||||
has_plyshlog_ref_temp_ = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double cMax() const;
|
||||
@ -168,9 +203,28 @@ namespace Opm
|
||||
double cMaxAds() const;
|
||||
|
||||
int adsIndex() const;
|
||||
|
||||
|
||||
/// the water velocity or water shear rate in PLYSHLOG table
|
||||
const std::vector<double>& shearWaterVelocity() const;
|
||||
|
||||
/// the viscosity reduction factor PLYSHLOG table
|
||||
const std::vector<double>& shearViscosityReductionFactor() const;
|
||||
|
||||
/// the reference polymer concentration in PLYSHLOG
|
||||
double plyshlogRefConc() const;
|
||||
|
||||
/// indicate wheter reference salinity is specified in PLYSHLOG
|
||||
bool hasPlyshlogRefSalinity() const;
|
||||
|
||||
/// indicate whether reference temperature is specified in PLYSHLOG
|
||||
bool hasPlyshlogRefTemp() const;
|
||||
|
||||
/// the reference salinity in PLYSHLOG
|
||||
double plyshlogRefSalinity() const;
|
||||
|
||||
/// the reference temperature in PLYSHLOG
|
||||
double plyshlogRefTemp() const;
|
||||
|
||||
double shearVrf(const double velocity) const;
|
||||
|
||||
double shearVrfWithDer(const double velocity, double& der) const;
|
||||
@ -272,6 +326,9 @@ namespace Opm
|
||||
void computeMcBoth(const double& c, double& mc,
|
||||
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:
|
||||
double c_max_;
|
||||
double mix_param_;
|
||||
@ -279,6 +336,8 @@ namespace Opm
|
||||
double dead_pore_vol_;
|
||||
double res_factor_;
|
||||
double c_max_ads_;
|
||||
|
||||
bool plyshlog_;
|
||||
AdsorptionBehaviour ads_index_;
|
||||
std::vector<double> c_vals_visc_;
|
||||
std::vector<double> visc_mult_vals_;
|
||||
@ -286,6 +345,14 @@ namespace Opm
|
||||
std::vector<double> ads_vals_;
|
||||
std::vector<double> water_vel_vals_;
|
||||
std::vector<double> shear_vrf_vals_;
|
||||
|
||||
double plyshlog_ref_conc_;
|
||||
double plyshlog_ref_salinity_;
|
||||
double plyshlog_ref_temp_;
|
||||
bool has_plyshlog_ref_salinity_;
|
||||
bool has_plyshlog_ref_temp_;
|
||||
|
||||
|
||||
void simpleAdsorptionBoth(double c, double& c_ads,
|
||||
double& dc_ads_dc, bool if_with_der) const;
|
||||
void adsorptionBoth(double c, double cmax,
|
||||
@ -302,6 +369,7 @@ namespace Opm
|
||||
double& deff_relperm_wat_ds,
|
||||
double& deff_relperm_wat_dc,
|
||||
bool if_with_der) const;
|
||||
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
@ -114,6 +114,15 @@ namespace Opm {
|
||||
/// \param[in] iteration current iteration number
|
||||
bool getConvergence(const double dt, const int iteration);
|
||||
|
||||
/// Assemble the residual and Jacobian of the nonlinear system.
|
||||
/// \param[in] reservoir_state reservoir state variables
|
||||
/// \param[in, out] well_state well state variables
|
||||
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
|
||||
void assemble(const ReservoirState& reservoir_state,
|
||||
WellState& well_state,
|
||||
const bool initial_assembly);
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
// --------- Types and enums ---------
|
||||
@ -135,6 +144,11 @@ namespace Opm {
|
||||
std::vector<double> wells_rep_radius_;
|
||||
std::vector<double> wells_perf_length_;
|
||||
|
||||
// shear-thinning factor for cell faces
|
||||
std::vector<double> shear_mult_faces_;
|
||||
// shear-thinning factor for well perforations
|
||||
std::vector<double> shear_mult_wells_;
|
||||
|
||||
// Need to declare Base members we want to use here.
|
||||
using Base::grid_;
|
||||
using Base::fluid_;
|
||||
@ -183,6 +197,11 @@ namespace Opm {
|
||||
using Base::drMaxRel;
|
||||
using Base::maxResidualAllowed;
|
||||
|
||||
using Base::updateWellControls;
|
||||
using Base::computeWellConnectionPressures;
|
||||
using Base::addWellControlEq;
|
||||
using Base::computeRelPerm;
|
||||
|
||||
|
||||
void
|
||||
makeConstantState(SolutionState& state) const;
|
||||
@ -257,7 +276,16 @@ namespace Opm {
|
||||
int nc,
|
||||
int nw) const;
|
||||
|
||||
/// Computing the water velocity without shear-thinning for the cell faces.
|
||||
/// The water velocity will be used for shear-thinning calculation.
|
||||
void computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,
|
||||
const std::vector<ADB>& phasePressure, const SolutionState& state,
|
||||
std::vector<double>& water_vel, std::vector<double>& visc_mult);
|
||||
|
||||
/// Computing the water velocity without shear-thinning for the well perforations based on the water flux rate.
|
||||
/// The water velocity will be used for shear-thinning calculation.
|
||||
void computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw,
|
||||
std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells);
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -266,7 +266,65 @@ namespace Opm {
|
||||
BlackoilPolymerModel<Grid>::
|
||||
assembleMassBalanceEq(const SolutionState& state)
|
||||
{
|
||||
Base::assembleMassBalanceEq(state);
|
||||
// Base::assembleMassBalanceEq(state);
|
||||
|
||||
// Compute b_p and the accumulation term b_p*s_p for each phase,
|
||||
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
|
||||
// These quantities are stored in rq_[phase].accum[1].
|
||||
// The corresponding accumulation terms from the start of
|
||||
// the timestep (b^0_p*s^0_p etc.) were already computed
|
||||
// on the initial call to assemble() and stored in rq_[phase].accum[0].
|
||||
computeAccum(state, 1);
|
||||
|
||||
// Set up the common parts of the mass balance equations
|
||||
// for each active phase.
|
||||
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
|
||||
const std::vector<ADB> kr = computeRelPerm(state);
|
||||
|
||||
|
||||
if (has_plyshlog_) {
|
||||
std::vector<double> water_vel;
|
||||
std::vector<double> visc_mult;
|
||||
|
||||
computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult);
|
||||
if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
|
||||
// std::cerr << " failed in calculating the shear-multiplier " << std::endl;
|
||||
OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
|
||||
}
|
||||
}
|
||||
|
||||
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
|
||||
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
|
||||
|
||||
residual_.material_balance_eq[ phaseIdx ] =
|
||||
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
|
||||
+ ops_.div*rq_[phaseIdx].mflux;
|
||||
}
|
||||
|
||||
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
|
||||
|
||||
// Add the extra (flux) terms to the mass balance equations
|
||||
// From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
|
||||
// The extra terms in the accumulation part of the equation are already handled.
|
||||
if (active_[ Oil ] && active_[ Gas ]) {
|
||||
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
|
||||
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
|
||||
|
||||
const UpwindSelector<double> upwindOil(grid_, ops_,
|
||||
rq_[po].dh.value());
|
||||
const ADB rs_face = upwindOil.select(state.rs);
|
||||
|
||||
const UpwindSelector<double> upwindGas(grid_, ops_,
|
||||
rq_[pg].dh.value());
|
||||
const ADB rv_face = upwindGas.select(state.rv);
|
||||
|
||||
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux);
|
||||
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux);
|
||||
|
||||
// OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
|
||||
|
||||
}
|
||||
|
||||
// Add polymer equation.
|
||||
if (has_polymer_) {
|
||||
residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
|
||||
@ -376,6 +434,14 @@ namespace Opm {
|
||||
rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * (transi * rq_[poly_pos_].dh);
|
||||
// Must recompute water flux since we have to use modified mobilities.
|
||||
rq_[ actph ].mflux = upwind.select(rq_[actph].b * rq_[actph].mob) * (transi * rq_[actph].dh);
|
||||
|
||||
// applying the shear-thinning factors
|
||||
if (has_plyshlog_) {
|
||||
V shear_mult_faces_v = Eigen::Map<V>(shear_mult_faces_.data(), shear_mult_faces_.size());
|
||||
ADB shear_mult_faces_adb = ADB::constant(shear_mult_faces_v);
|
||||
rq_[poly_pos_].mflux = rq_[poly_pos_].mflux / shear_mult_faces_adb;
|
||||
rq_[actph].mflux = rq_[actph].mflux / shear_mult_faces_adb;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -468,6 +534,95 @@ namespace Opm {
|
||||
}
|
||||
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
BlackoilPolymerModel<Grid>::assemble(const ReservoirState& reservoir_state,
|
||||
WellState& well_state,
|
||||
const bool initial_assembly)
|
||||
{
|
||||
using namespace Opm::AutoDiffGrid;
|
||||
|
||||
// Possibly switch well controls and updating well state to
|
||||
// get reasonable initial conditions for the wells
|
||||
updateWellControls(well_state);
|
||||
|
||||
// Create the primary variables.
|
||||
SolutionState state = variableState(reservoir_state, well_state);
|
||||
|
||||
if (initial_assembly) {
|
||||
// Create the (constant, derivativeless) initial state.
|
||||
SolutionState state0 = state;
|
||||
makeConstantState(state0);
|
||||
// Compute initial accumulation contributions
|
||||
// and well connection pressures.
|
||||
computeAccum(state0, 0);
|
||||
computeWellConnectionPressures(state0, well_state);
|
||||
}
|
||||
|
||||
// OPM_AD_DISKVAL(state.pressure);
|
||||
// OPM_AD_DISKVAL(state.saturation[0]);
|
||||
// OPM_AD_DISKVAL(state.saturation[1]);
|
||||
// OPM_AD_DISKVAL(state.saturation[2]);
|
||||
// OPM_AD_DISKVAL(state.rs);
|
||||
// OPM_AD_DISKVAL(state.rv);
|
||||
// OPM_AD_DISKVAL(state.qs);
|
||||
// OPM_AD_DISKVAL(state.bhp);
|
||||
|
||||
// -------- Mass balance equations --------
|
||||
assembleMassBalanceEq(state);
|
||||
|
||||
// -------- Well equations ----------
|
||||
if ( ! wellsActive() ) {
|
||||
return;
|
||||
}
|
||||
|
||||
V aliveWells;
|
||||
|
||||
const int np = wells().number_of_phases;
|
||||
std::vector<ADB> cq_s(np, ADB::null());
|
||||
|
||||
const int nw = wells().number_of_wells;
|
||||
const int nperf = wells().well_connpos[nw];
|
||||
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
|
||||
|
||||
std::vector<ADB> mob_perfcells(np, ADB::null());
|
||||
std::vector<ADB> b_perfcells(np, ADB::null());
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
mob_perfcells[phase] = subset(rq_[phase].mob, well_cells);
|
||||
b_perfcells[phase] = subset(rq_[phase].b, well_cells);
|
||||
}
|
||||
if (param_.solve_welleq_initially_ && initial_assembly) {
|
||||
// solve the well equations as a pre-processing step
|
||||
Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
|
||||
}
|
||||
|
||||
Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
|
||||
|
||||
if (has_plyshlog_) {
|
||||
std::vector<double> water_vel_wells;
|
||||
std::vector<double> visc_mult_wells;
|
||||
|
||||
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
|
||||
computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells);
|
||||
|
||||
if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) {
|
||||
OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells ");
|
||||
}
|
||||
|
||||
// applying the shear-thinning to the water phase
|
||||
V shear_mult_wells_v = Eigen::Map<V>(shear_mult_wells_.data(), shear_mult_wells_.size());
|
||||
ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v);
|
||||
mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
|
||||
}
|
||||
|
||||
Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
|
||||
Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
|
||||
Base::addWellFluxEq(cq_s, state);
|
||||
addWellContributionToMassBalanceEq(cq_s, state, well_state);
|
||||
addWellControlEq(state, well_state, aliveWells);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class Grid>
|
||||
@ -590,6 +745,150 @@ namespace Opm {
|
||||
return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
|
||||
}
|
||||
|
||||
template<class Grid>
|
||||
void
|
||||
BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,
|
||||
const std::vector<ADB>& phasePressure, const SolutionState& state,
|
||||
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];
|
||||
|
||||
const std::vector<PhasePresence> cond = phaseCondition();
|
||||
|
||||
const ADB tr_mult = transMult(state.pressure);
|
||||
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_);
|
||||
rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
|
||||
|
||||
// compute gravity potensial using the face average as in eclipse and MRST
|
||||
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_);
|
||||
const ADB rhoavg = ops_.caver * rho;
|
||||
rq_[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
|
||||
if (use_threshold_pressure_) {
|
||||
applyThresholdPressures(rq_[ phase ].dh);
|
||||
}
|
||||
|
||||
const ADB& b = rq_[ phase ].b;
|
||||
const ADB& mob = rq_[ phase ].mob;
|
||||
const ADB& dh = rq_[ phase ].dh;
|
||||
UpwindSelector<double> upwind(grid_, ops_, dh.value());
|
||||
|
||||
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
|
||||
const ADB mc = computeMc(state);
|
||||
ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
|
||||
cmax,
|
||||
kr[canonicalPhaseIdx]);
|
||||
ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data());
|
||||
rq_[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
|
||||
|
||||
const V& polymer_conc = state.concentration.value();
|
||||
V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
|
||||
V visc_mult_faces = upwind.select(visc_mult_cells);
|
||||
|
||||
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());
|
||||
|
||||
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());
|
||||
|
||||
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) {
|
||||
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;
|
||||
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());
|
||||
|
||||
for (int i = 0; i < nface; ++i) {
|
||||
water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]);
|
||||
}
|
||||
}
|
||||
|
||||
template<class Grid>
|
||||
void
|
||||
BlackoilPolymerModel<Grid>::computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw,
|
||||
std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells)
|
||||
{
|
||||
if( ! wellsActive() ) return ;
|
||||
|
||||
const int nw = wells().number_of_wells;
|
||||
const int nperf = wells().well_connpos[nw];
|
||||
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());
|
||||
|
||||
const V& polymer_conc = state.concentration.value();
|
||||
|
||||
V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
|
||||
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());
|
||||
|
||||
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
|
||||
ADB b_perfcells = subset(rq_[water_pos].b, well_cells);
|
||||
|
||||
const ADB& p_perfcells = subset(state.pressure, well_cells);
|
||||
const V& cdp = well_perforation_pressure_diffs_;
|
||||
const ADB perfpressure = (wops_.w2p * state.bhp) + cdp;
|
||||
// Pressure drawdown (also used to determine direction of flow)
|
||||
const ADB drawdown = p_perfcells - perfpressure;
|
||||
|
||||
// selects injection perforations
|
||||
V selectInjectingPerforations = V::Zero(nperf);
|
||||
for (int c = 0; c < nperf; ++c) {
|
||||
if (drawdown.value()[c] < 0) {
|
||||
selectInjectingPerforations[c] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// for the injection wells
|
||||
for (int i = 0; i < well_cells.size(); ++i) {
|
||||
if (xw.polymerInflow()[well_cells[i]] == 0. && selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
|
||||
visc_mult_wells[i] = 1.;
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
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;
|
||||
b_wells.resize(b_perfcells.size());
|
||||
std::copy(&(b_perfcells.value()[0]), &(b_perfcells.value()[0]) + b_perfcells.size(), b_wells.begin());
|
||||
|
||||
for (int 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]);
|
||||
// TODO: CHECK to make sure this formulation is corectly used. Why muliplied by bW.
|
||||
// Although this formulation works perfectly with the tests compared with other formulations
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
|
||||
|
||||
@ -60,6 +60,64 @@ namespace Opm {
|
||||
return polymer_props_.cMax();
|
||||
}
|
||||
|
||||
const std::vector<double>&
|
||||
PolymerPropsAd::shearWaterVelocity() const
|
||||
{
|
||||
return polymer_props_.shearWaterVelocity();
|
||||
}
|
||||
|
||||
const std::vector<double>&
|
||||
PolymerPropsAd::shearViscosityReductionFactor() const
|
||||
{
|
||||
return polymer_props_.shearViscosityReductionFactor();
|
||||
}
|
||||
|
||||
double
|
||||
PolymerPropsAd::plyshlogRefConc() const
|
||||
{
|
||||
return polymer_props_.plyshlogRefConc();
|
||||
}
|
||||
|
||||
bool
|
||||
PolymerPropsAd::hasPlyshlogRefSalinity() const
|
||||
{
|
||||
return polymer_props_.hasPlyshlogRefSalinity();
|
||||
}
|
||||
|
||||
bool
|
||||
PolymerPropsAd::hasPlyshlogRefTemp() const
|
||||
{
|
||||
return polymer_props_.hasPlyshlogRefTemp();
|
||||
}
|
||||
|
||||
double
|
||||
PolymerPropsAd::plyshlogRefSalinity() const
|
||||
{
|
||||
return polymer_props_.plyshlogRefSalinity();
|
||||
}
|
||||
|
||||
double
|
||||
PolymerPropsAd::plyshlogRefTemp() const
|
||||
{
|
||||
return polymer_props_.plyshlogRefTemp();
|
||||
}
|
||||
|
||||
double
|
||||
PolymerPropsAd::viscMult(double c) const
|
||||
{
|
||||
return polymer_props_.viscMult(c);
|
||||
}
|
||||
|
||||
V
|
||||
PolymerPropsAd::viscMult(const V& c) const
|
||||
{
|
||||
int nc = c.size();
|
||||
V visc_mult(nc);
|
||||
for (int i = 0; i < nc; ++i) {
|
||||
visc_mult[i] = polymer_props_.viscMult(c[i]);
|
||||
}
|
||||
return visc_mult;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -259,4 +317,12 @@ namespace Opm {
|
||||
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
|
||||
|
||||
@ -41,9 +41,36 @@ namespace Opm {
|
||||
/// \return The max concentration injected.
|
||||
double cMax() const;
|
||||
|
||||
/// \ return The water velcoity or shear rate in the PLYSHLOG table
|
||||
const std::vector<double>& shearWaterVelocity() const;
|
||||
|
||||
/// \ return The viscosity reducation factor in the PLYSHLOG table
|
||||
const std::vector<double>& shearViscosityReductionFactor() const;
|
||||
|
||||
/// \ return The reference polymer concentration for PLYSHLOG table
|
||||
double plyshlogRefConc() const;
|
||||
|
||||
/// \ return The flag indicating if reference salinity is specified in PLYSHLOG keyword
|
||||
bool hasPlyshlogRefSalinity() const;
|
||||
|
||||
/// \ return The flag indicating if reference temperature is specified in PLYSHLOG keyword
|
||||
bool hasPlyshlogRefTemp() const;
|
||||
|
||||
/// \ return The reference salinity in PLYSHLOG keyword
|
||||
double plyshlogRefSalinity() const;
|
||||
|
||||
/// \ return The reference temperature in PLYSHLOG keyword
|
||||
double plyshlogRefTemp() const;
|
||||
|
||||
double viscMult(double c) const; // multipler interpolated from PLYVISC table
|
||||
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
|
||||
V viscMult(const V& c) const;
|
||||
/// \param[in] c Array of n polymer concentraion values.
|
||||
/// \return Array of n viscosity multiplier from PLVISC table.
|
||||
|
||||
/// Constructor wrapping a polymer props.
|
||||
PolymerPropsAd(const PolymerProperties& polymer_props);
|
||||
|
||||
@ -103,6 +130,14 @@ namespace Opm {
|
||||
ADB
|
||||
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:
|
||||
const PolymerProperties& polymer_props_;
|
||||
};
|
||||
|
||||
Loading…
Reference in New Issue
Block a user