Added support for adsorbtion index (desorbtion allowed or not).

This commit is contained in:
Xavier Raynaud 2012-03-26 16:37:39 +02:00
parent 294c8b03be
commit 6d9f0ea3ec
4 changed files with 64 additions and 45 deletions

View File

@ -546,6 +546,7 @@ main(int argc, char** argv)
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.);
int ads_index = param.getDefault("ads_index", 2);
std::vector<double> c_vals_visc(2, -1e100);
c_vals_visc[0] = 0.0;
c_vals_visc[1] = 7.0;
@ -562,8 +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, res_factor, c_max_ads,
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,
ads_index, c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
}
// Initialize polymer inflow function.

View File

@ -43,16 +43,19 @@ namespace Opm
double dead_pore_vol,
double res_factor,
double c_max_ads,
int ads_index,
const std::vector<double>& c_vals_visc,
const std::vector<double>& visc_mult_vals,
const std::vector<double>& c_vals_ads,
const std::vector<double>& ads_vals)
const std::vector<double>& ads_vals
)
: c_max_(c_max),
mix_param_(mix_param),
rock_density_(rock_density),
dead_pore_vol_(dead_pore_vol),
res_factor_(res_factor),
c_max_ads_(c_max_ads),
ads_index_(ads_index),
c_vals_visc_(c_vals_visc),
visc_mult_vals_(visc_mult_vals),
c_vals_ads_(c_vals_ads),
@ -71,10 +74,12 @@ namespace Opm
double dead_pore_vol,
double res_factor,
double c_max_ads,
int ads_index,
const std::vector<double>& c_vals_visc,
const std::vector<double>& visc_mult_vals,
const std::vector<double>& c_vals_ads,
const std::vector<double>& ads_vals)
const std::vector<double>& ads_vals
)
{
c_max_ = c_max;
mix_param_ = mix_param;
@ -86,6 +91,7 @@ namespace Opm
visc_mult_vals_ = visc_mult_vals;
c_vals_ads_ = c_vals_ads;
ads_vals_ = ads_vals;
ads_index_ = ads_index;
}
void readFromDeck(const EclipseGridParser& gridparser)
@ -102,6 +108,7 @@ namespace Opm
dead_pore_vol_ = plyrock[0];
res_factor_ = plyrock[2];
rock_density_ = plyrock[3];
ads_index_ = plyrock[4];
c_max_ads_ = plyrock[5];
// We assume NTPVT=1
@ -146,6 +153,11 @@ namespace Opm
return c_max_ads_;
}
int adsIndex() const
{
return ads_index_;
}
double viscMult(double c) const
{
return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c);
@ -157,17 +169,44 @@ namespace Opm
return Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c);
}
double adsorbtion(double c) const
double simpleAdsorbtion(double c) const
{
return Opm::linearInterpolation(c_vals_ads_, ads_vals_, c);
}
double adsorbtionWithDer(double c, double* der) const
double simpleAdsorbtionWithDer(double c, double* der) const
{
*der = Opm::linearInterpolationDerivative(c_vals_ads_, ads_vals_, c);
return Opm::linearInterpolation(c_vals_ads_, ads_vals_, c);
}
double adsorbtion(double c, double cmax) const
{
if (ads_index_ == 1) {
return simpleAdsorbtion(c);
} else if (ads_index_ == 2) {
return simpleAdsorbtion(std::max(c, cmax));
} else {
THROW("Invalid Adsoption index");
}
}
double adsorbtionWithDer(double c, double cmax, double* der) const
{
if (ads_index_ == 1) {
return simpleAdsorbtionWithDer(c, der);
} else if (ads_index_ == 2) {
if (c < cmax) {
*der = 0;
return simpleAdsorbtion(cmax);
} else {
return simpleAdsorbtionWithDer(c, der);
}
} else {
THROW("Invalid Adsoption index");
}
}
void effectiveInvVisc(const double c, const double* visc, double* inv_visc_eff) const
{
@ -238,6 +277,7 @@ namespace Opm
double dead_pore_vol_;
double res_factor_;
double c_max_ads_;
int ads_index_;
std::vector<double> c_vals_visc_;
std::vector<double> visc_mult_vals_;
std::vector<double> c_vals_ads_;

View File

@ -294,8 +294,8 @@ namespace Opm
ff = tm.fracFlow(s_arg, c_arg, cell);
mc = tm.computeMc(c_arg);
double rhor = tm.polyprops_.rockDensity();
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
double ads = tm.polyprops_.adsorbtion(std::max(c_arg, cmax0));
double ads0 = tm.polyprops_.adsorbtion(c0, cmax0);
double ads = tm.polyprops_.adsorbtion(c_arg, cmax0);
res_s = s_arg - s0 + dtpv*(outflux*ff + influx);
res_c = s_arg*(1 - dps)*c_arg - (s0 - dps)*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
@ -314,8 +314,8 @@ namespace Opm
double ff = tm.fracFlow(s, c, cell);
double mc = tm.computeMc(c);
double rhor = tm.polyprops_.rockDensity();
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
double ads0 = tm.polyprops_.adsorbtion(c0, cmax0);
double ads = tm.polyprops_.adsorbtion(c, cmax0);
double res = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+ dtpv*(outflux*ff*mc + influx_polymer);
@ -345,7 +345,7 @@ namespace Opm
cmax0 = tm.cmax_[cell];
dps = tm.polyprops_.deadPoreVol();
rhor = tm.polyprops_.rockDensity();
ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
ads0 = tm.polyprops_.adsorbtion(c0, cmax0);
res_factor = tm.polyprops_.resFactor();
c_max_ads = tm.polyprops_.cMaxAds();
double dflux = -tm.source_[cell];
@ -387,8 +387,8 @@ namespace Opm
double mc = tm.computeMc(c);
double dps = tm.polyprops_.deadPoreVol();
double rhor = tm.polyprops_.rockDensity();
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
double ads0 = tm.polyprops_.adsorbtion(c0, cmax0);
double ads = tm.polyprops_.adsorbtion(c, cmax0);
res[0] = s - s0 + dtpv*(outflux*ff + influx);
res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
@ -403,8 +403,8 @@ namespace Opm
mc = tm.computeMc(c);
double dps = tm.polyprops_.deadPoreVol();
double rhor = tm.polyprops_.rockDensity();
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
double ads0 = tm.polyprops_.adsorbtion(c0, cmax0);
double ads = tm.polyprops_.adsorbtion(c, cmax0);
res[0] = s - s0 + dtpv*(outflux*ff + influx);
res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
@ -426,7 +426,7 @@ namespace Opm
double c = x[1];
double ff = tm.fracFlow(s, c, cell);
double mc = tm.computeMc(c);
double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0));
double ads = tm.polyprops_.adsorbtion(c, cmax0);
return (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+ dtpv*(outflux*ff*mc + influx_polymer);
@ -456,14 +456,8 @@ namespace Opm
double ff = tm.fracFlowWithDer(s, c, cell, ff_ds_dc);
double mc_dc;
double mc = tm.computeMcWithDer(c, &mc_dc);
double ads;
double ads_dc;
if (c < cmax0) {
ads = tm.polyprops_.adsorbtion(cmax0);
ads_dc = 0;
} else {
ads = tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
}
double ads = tm.polyprops_.adsorbtionWithDer(c, cmax0, &ads_dc);
res[0] = s - s0 + dtpv*(outflux*ff + influx);
res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
@ -499,15 +493,9 @@ namespace Opm
double mc = tm.computeMcWithDer(c, &mc_dc);
double dps = tm.polyprops_.deadPoreVol();
double rhor = tm.polyprops_.rockDensity();
double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0));
double ads;
double ads0 = tm.polyprops_.adsorbtion(c0, cmax0);
double ads_dc;
if (c < cmax0) {
ads = tm.polyprops_.adsorbtion(cmax0);
ads_dc = 0;
} else {
ads = tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
}
double ads = tm.polyprops_.adsorbtionWithDer(c, cmax0, &ads_dc);
res[0] = s - s0 + dtpv*(outflux*ff + influx);
res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
@ -551,11 +539,7 @@ namespace Opm
double mc_dc;
double mc = tm.computeMcWithDer(c, &mc_dc);
double ads_dc;
if (c < cmax0) {
ads_dc = 0;
} else {
tm.polyprops_.adsorbtionWithDer(c, &ads_dc);
}
tm.polyprops_.adsorbtionWithDer(c, cmax0, &ads_dc);
res_s_ds_dc[0] = 1 + dtpv*outflux*ff_ds_dc[0];
res_s_ds_dc[1] = dtpv*outflux*ff_ds_dc[1];
res_c_ds_dc[0] = (1 - dps)*c + dtpv*outflux*(ff_ds_dc[0])*mc;
@ -849,7 +833,7 @@ 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_ads = polyprops_.adsorbtion(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;
@ -874,15 +858,9 @@ 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_ads = polyprops_.adsorbtionWithDer(c, c_max, &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;

View File

@ -204,7 +204,7 @@ namespace Opm
const double* poro = props.porosity();
double abs_mass = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
abs_mass += polyprops.adsorbtion(cmax[cell])*pv[cell]*((1.0 - poro[cell])/poro[cell])*rhor;
abs_mass += polyprops.simpleAdsorbtion(cmax[cell])*pv[cell]*((1.0 - poro[cell])/poro[cell])*rhor;
}
return abs_mass;
}