Added actual resistance factor. Not tested.

This commit is contained in:
Xavier Raynaud 2012-03-23 17:40:03 +01:00
parent 4c77354d56
commit 1836247dab
3 changed files with 55 additions and 7 deletions

View File

@ -545,6 +545,8 @@ main(int argc, char** argv)
double mix_param = param.getDefault("mix_param", 1.0);
double rock_density = param.getDefault("rock_density", 1000.0);
double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15);
double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability
double c_max_ads = param.getDefault("c_max_ads", 1.);
std::vector<double> c_vals_visc(2, -1e100);
c_vals_visc[0] = 0.0;
c_vals_visc[1] = 7.0;
@ -561,7 +563,8 @@ main(int argc, char** argv)
// polyprop.ads_vals[1] = param.getDefault("c_max_ads", 0.0025);
ads_vals[1] = 0.0015;
ads_vals[2] = 0.0025;
polyprop.set(c_max, mix_param, rock_density, dead_pore_vol, c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
polyprop.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
}
// Initialize polymer inflow function.

View File

@ -41,6 +41,8 @@ namespace Opm
double mix_param,
double rock_density,
double dead_pore_vol,
double res_factor,
double c_max_ads,
const std::vector<double>& c_vals_visc,
const std::vector<double>& visc_mult_vals,
const std::vector<double>& c_vals_ads,
@ -49,6 +51,8 @@ namespace Opm
mix_param_(mix_param),
rock_density_(rock_density),
dead_pore_vol_(dead_pore_vol),
res_factor_(res_factor),
c_max_ads_(c_max_ads),
c_vals_visc_(c_vals_visc),
visc_mult_vals_(visc_mult_vals),
c_vals_ads_(c_vals_ads),
@ -65,6 +69,8 @@ namespace Opm
double mix_param,
double rock_density,
double dead_pore_vol,
double res_factor,
double c_max_ads,
const std::vector<double>& c_vals_visc,
const std::vector<double>& visc_mult_vals,
const std::vector<double>& c_vals_ads,
@ -74,6 +80,8 @@ namespace Opm
mix_param_ = mix_param;
rock_density_ = rock_density;
dead_pore_vol_ = dead_pore_vol;
res_factor_ = res_factor;
c_max_ads_ = c_max_ads;
c_vals_visc_ = c_vals_visc;
visc_mult_vals_ = visc_mult_vals;
c_vals_ads_ = c_vals_ads;
@ -92,7 +100,9 @@ namespace Opm
// We assume NTSFUN=1
const std::vector<double>& plyrock = gridparser.getPLYROCK().plyrock_;
dead_pore_vol_ = plyrock[0];
res_factor_ = plyrock[2];
rock_density_ = plyrock[3];
c_max_ads_ = plyrock[5];
// We assume NTPVT=1
const PLYVISC& plyvisc = gridparser.getPLYVISC();
@ -126,6 +136,16 @@ namespace Opm
return dead_pore_vol_;
}
double resFactor() const
{
return res_factor_;
}
double cMaxAds() const
{
return c_max_ads_;
}
double viscMult(double c) const
{
return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c);
@ -233,6 +253,8 @@ namespace Opm
double mix_param_;
double rock_density_;
double dead_pore_vol_;
double res_factor_;
double c_max_ads_;
std::vector<double> c_vals_visc_;
std::vector<double> visc_mult_vals_;
std::vector<double> c_vals_ads_;

View File

@ -38,6 +38,8 @@ public:
double porosity;
double dtpv; // dt/pv(i)
double dps;
double res_factor;
double c_max_ads;
double rhor;
double ads0;
int gradient_method;
@ -197,6 +199,7 @@ namespace Opm
const TransportModelPolymer& tm_;
const int cell_;
const double s0_;
const double cmax0_;
const double influx_; // sum_j min(v_ij, 0)*f(s_j)
const double outflux_; // sum_j max(v_ij, 0)
const double dtpv_; // dt/pv(i)
@ -204,6 +207,7 @@ namespace Opm
explicit ResidualS(const TransportModelPolymer& tmodel,
const int cell,
const double s0,
const double cmax0,
const double influx,
const double outflux,
const double dtpv,
@ -211,6 +215,7 @@ namespace Opm
: tm_(tmodel),
cell_(cell),
s0_(s0),
cmax0_(cmax0),
influx_(influx),
outflux_(outflux),
dtpv_(dtpv),
@ -301,7 +306,7 @@ namespace Opm
double operator()(double c) const
{
double dps = tm.polyprops_.deadPoreVol();
ResidualS res_s(tm, cell, s0, influx, outflux, dtpv, c);
ResidualS res_s(tm, cell, s0, cmax0, influx, outflux, dtpv, c);
int iters_used;
// Solve for s first.
s = modifiedRegulaFalsi(res_s, std::max(tm.smin_[2*cell], dps), tm.smax_[2*cell],
@ -341,7 +346,8 @@ namespace Opm
dps = tm.polyprops_.deadPoreVol();
rhor = tm.polyprops_.rockDensity();
ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
res_factor = tm.polyprops_.resFactor();
c_max_ads = tm.polyprops_.cMaxAds();
double dflux = -tm.source_[cell];
bool src_is_inflow = dflux < 0.0;
influx = src_is_inflow ? dflux : 0.0;
@ -843,6 +849,10 @@ namespace Opm
{
double c_max_limit = polyprops_.cMax();
double cbar = c/c_max_limit;
double c_ads = polyprops_.adsorbtion(std::max(c, cmax_[cell]));
double c_max_ads = polyprops_.cMaxAds();
double res_factor = polyprops_.resFactor();
double res_k = 1 + (res_factor - 1)*c_ads/c_max_ads;
double mu_w = visc_[0];
double mu_m = polyprops_.viscMult(c)*mu_w;
double mu_p = polyprops_.viscMult(polyprops_.cMax())*mu_w;
@ -855,7 +865,7 @@ namespace Opm
double sat[2] = { s, 1.0 - s };
double mob[2];
props_.relperm(1, sat, &cell, mob, 0);
mob[0] *= inv_visc_eff[0];
mob[0] *= inv_visc_eff[0]/res_k;
mob[1] *= inv_visc_eff[1];
return mob[0]/(mob[0] + mob[1]);
}
@ -864,6 +874,19 @@ namespace Opm
{
double c_max_limit = polyprops_.cMax();
double cbar = c/c_max_limit;
double c_ads;
double c_ads_dc;
double c_max = cmax_[cell];
if (c < c_max) {
c_ads = polyprops_.adsorbtion(c_max);
c_ads_dc = 0;
} else {
c_ads = polyprops_.adsorbtionWithDer(c, &c_ads_dc);
}
double c_max_ads = polyprops_.cMaxAds();
double res_factor = polyprops_.resFactor();
double res_k = 1 + (res_factor - 1)*c_ads/c_max_ads;
double res_k_dc = (res_factor - 1)*c_ads_dc/c_max_ads;
double mu_w = visc_[0];
double mu_m_dc; // derivative of mu_m with respect to c
double mu_m = polyprops_.viscMultWithDer(c, &mu_m_dc)*mu_w;
@ -886,11 +909,11 @@ namespace Opm
double perm[2];
double perm_ds[4];
props_.relperm(1, sat, &cell, perm, perm_ds);
mob[0] = perm[0]/visc_eff[0];
mob[0] = perm[0]/visc_eff[0]/res_k;
mob[1] = perm[1]/visc_eff[1];
mob_ds[0] = perm_ds[0]/visc_eff[0];
mob_ds[0] = perm_ds[0]/visc_eff[0]/res_k;
mob_ds[1] = perm_ds[1]/visc_eff[1];
mob_dc[0] = - perm[0]*mu_w_eff_dc/(mu_w_eff*mu_w_eff);
mob_dc[0] = - perm[0]*(mu_w_eff_dc/(mu_w_eff*mu_w_eff*res_k) + res_k_dc/(res_k*res_k*mu_w_eff));
mob_dc[1] = 0.;
der[0] = (mob_ds[0]*mob[1] - mob_ds[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1]));
der[1] = (mob_dc[0]*mob[1] - mob_dc[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1]));