diff --git a/opm/polymer/PolymerProperties.cpp b/opm/polymer/PolymerProperties.cpp index 3afba933a..572c5f538 100644 --- a/opm/polymer/PolymerProperties.cpp +++ b/opm/polymer/PolymerProperties.cpp @@ -243,6 +243,52 @@ namespace Opm } } + void PolymerProperties::effectiveInvPolyVisc(const double c, + const double* visc, + double& inv_mu_p_eff) const + { + double dummy; + effectiveInvPolyViscBoth(c, visc, inv_mu_p_eff, dummy, false); + + } + + void PolymerProperties::effectiveInvPolyViscWithDer(const double c, + const double* visc, + double& inv_mu_p_eff, + double& d_inv_mu_p_eff_dc) const + { + effectiveInvPolyViscBoth(c, visc, inv_mu_p_eff, d_inv_mu_p_eff_dc, true); + + } + + void PolymerProperties::effectiveInvPolyViscBoth(const double c, + const double* visc, + double& inv_mu_p_eff, + double& dinv_mu_p_eff_dc, + const bool if_with_der) const + { + const double omega = mix_param_; + const double mu_w = visc[0]; + + double mu_m = 0.0; + double dmu_m_dc = 0.0; + + if (if_with_der) { + mu_m = viscMultWithDer(c, &dmu_m_dc)*mu_w; + dmu_m_dc *= mu_w; + } else { + mu_m = viscMult(c)*mu_w; + } + + const double inv_mu_m_omega = std::pow(mu_m, -omega); + const double mu_p = viscMult(c_max_) * mu_w; + inv_mu_p_eff = inv_mu_m_omega * std::pow(mu_p, omega - 1.); + + if (if_with_der) { + dinv_mu_p_eff_dc = -omega * dmu_m_dc * std::pow(mu_m, -omega - 1) * std::pow(mu_p, omega - 1); + } + } + void PolymerProperties::effectiveRelperm(const double c, const double cmax, const double* relperm, diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index cbdc9ed22..141073e86 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -292,6 +292,15 @@ namespace Opm const double* visc, double& inv_mu_w_eff, double& dinv_mu_w_eff_dc) const; + void effectiveInvPolyVisc(const double c, + const double* visc, + double& inv_mu_p_eff) const; + + void effectiveInvPolyViscWithDer(const double c, + const double* visc, + double& inv_mu_p_eff, + double& d_inv_mu_p_eff_dc) const; + void effectiveRelperm(const double c, const double cmax, const double* relperm, @@ -401,6 +410,13 @@ namespace Opm void effectiveInvViscBoth(const double c, const double* visc, double& inv_mu_w_eff, double& dinv_mu_w_eff_dc, bool if_with_der) const; + + void effectiveInvPolyViscBoth(const double c, + const double* visc, + double& inv_mu_p_eff, + double& dinv_mu_p_eff_dc, + const bool if_with_der) const; + void effectiveRelpermBoth(const double c, const double cmax, const double* relperm, diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 87b67e8fc..c8f8d2a52 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -464,7 +464,8 @@ namespace Opm { // Reduce mobility of water phase by relperm reduction and effective viscosity increase. rq_[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc; // Compute polymer mobility. - rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_wat_eff_visc; + const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value().data()); + rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_poly_eff_visc; rq_[poly_pos_].b = rq_[actph].b; rq_[poly_pos_].dh = rq_[actph].dh; UpwindSelector upwind(grid_, ops_, rq_[poly_pos_].dh.value()); diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp index 5db5abb2e..2deac485f 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -192,6 +192,34 @@ namespace Opm { + ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const double* visc) const + { + const int nc = c.size(); + V inv_mu_p_eff(nc); + V dinv_mu_p_eff(nc); + for (int i = 0; i < nc; ++i) { + double im = 0; + double dim = 0; + // TODO: the usage of visc can be likely wrong, while more investigation will be requried. + polymer_props_.effectiveInvPolyViscWithDer(c.value()(i), visc, im, dim); + inv_mu_p_eff(i) = im; + dinv_mu_p_eff(i) = dim; + } + + ADB::M dim_diag(dinv_mu_p_eff.matrix().asDiagonal()); + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dim_diag * c.derivative()[block]; + } + return ADB::function(std::move(inv_mu_p_eff), std::move(jacs)); + } + + + + + V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const { const int nc = c.size(); diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index cc54ce282..c970946dc 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -91,6 +91,12 @@ namespace Opm { /// \return value of inverse effective water viscosity. ADB effectiveInvWaterVisc(const ADB& c,const double* visc) const; + + /// \param[in] c ADB of polymer concentraion values. + /// \param[in] visc Array of water viscosity values + /// \return ADB of inverse effective polymer viscosity. + ADB + effectiveInvPolymerVisc(const ADB& c, const double* visc) const; /// \param[in] c Array of n polymer concentraion values. /// \return Array of n mc values, here mc means m(c) * c.