From 6d9f0ea3ec295426dcdd8a338175483a4a9347f4 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 26 Mar 2012 16:37:39 +0200 Subject: [PATCH] Added support for adsorbtion index (desorbtion allowed or not). --- examples/polymer_reorder.cpp | 5 ++- opm/polymer/PolymerProperties.hpp | 48 ++++++++++++++++++++++-- opm/polymer/TransportModelPolymer.cpp | 54 ++++++++------------------- opm/polymer/polymerUtilities.cpp | 2 +- 4 files changed, 64 insertions(+), 45 deletions(-) diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 8a8d8bcb3..d12424624 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -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 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. diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index 81e4c0f1e..7c09f4bb5 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -43,16 +43,19 @@ namespace Opm double dead_pore_vol, double res_factor, double c_max_ads, + int ads_index, const std::vector& c_vals_visc, const std::vector& visc_mult_vals, const std::vector& c_vals_ads, - const std::vector& ads_vals) + const std::vector& 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& c_vals_visc, const std::vector& visc_mult_vals, const std::vector& c_vals_ads, - const std::vector& ads_vals) + const std::vector& 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 c_vals_visc_; std::vector visc_mult_vals_; std::vector c_vals_ads_; diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index b5b2c052a..71eb2669c 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -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; diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp index e9f0dd76f..f8b2c9f45 100644 --- a/opm/polymer/polymerUtilities.cpp +++ b/opm/polymer/polymerUtilities.cpp @@ -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; }