diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 0f45bc9df..3d5b37847 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -694,7 +694,7 @@ main(int argc, char** argv) // Report volume balances. Opm::computeSaturatedVol(porevol, state.saturation(), satvol); - polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration()); + polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), polydata.deadPoreVol()); Opm::computeInjectedProduced(*props, polydata, state.saturation(), state.concentration(), src, simtimer.currentStepLength(), inflow_c, injected, produced, polyinj, polyprod); diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp index 8585d1a43..3470eada4 100644 --- a/opm/polymer/polymerUtilities.cpp +++ b/opm/polymer/polymerUtilities.cpp @@ -167,10 +167,12 @@ namespace Opm /// @param[in] pv the pore volume by cell. /// @param[in] s saturation values (for all P phases) /// @param[in] c polymer concentration + /// @param[in] dps dead pore space /// @return total polymer mass in grid. double computePolymerMass(const std::vector& pv, const std::vector& s, - const std::vector& c) + const std::vector& c, + const double dps) { const int num_cells = pv.size(); const int np = s.size()/pv.size(); @@ -179,10 +181,40 @@ namespace Opm } double polymass = 0.0; for (int cell = 0; cell < num_cells; ++cell) { - polymass += c[cell]*pv[cell]*s[np*cell + 0]; + polymass += c[cell]*pv[cell]*s[np*cell + 0]*(1.0 - dps); } return polymass; } + + + /// @brief Computes total absorbed polymer mass over all grid cells. + /// @param[in] polyprops polymer properties + /// @param[in] pv the pore volume by cell. + /// @param[in] s saturation values (for all P phases) + /// @param[in] cmax max polymer concentration for cell + /// @param[in] dps dead pore space + /// @return total absorbed polymer mass. + double computePolymerAbsorbed(const Opm::PolymerProperties& polyprops, + const std::vector& pv, + const std::vector& s, + const std::vector& cmax, + const double dps) + { + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + double abs_mass = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + const double max_polymass = cmax[cell]*pv[cell]*s[np*cell + 0]*(1.0 - dps); + abs_mass += polyprops.adsorbtion(cmax[cell])*max_polymass; + } + return abs_mass; + } + + + } // namespace Opm diff --git a/opm/polymer/polymerUtilities.hpp b/opm/polymer/polymerUtilities.hpp index 78184304f..5837c2e53 100644 --- a/opm/polymer/polymerUtilities.hpp +++ b/opm/polymer/polymerUtilities.hpp @@ -91,14 +91,29 @@ namespace Opm double& polyinj, double& polyprod); - /// @brief Computes total polymer mass over all grid cells. + /// @brief Computes total (free) polymer mass over all grid cells. /// @param[in] pv the pore volume by cell. /// @param[in] s saturation values (for all P phases) /// @param[in] c polymer concentration + /// @param[in] dps dead pore space /// @return total polymer mass in grid. double computePolymerMass(const std::vector& pv, const std::vector& s, - const std::vector& c); + const std::vector& c, + const double dps); + + /// @brief Computes total absorbed polymer mass over all grid cells. + /// @param[in] polyprops polymer properties + /// @param[in] pv the pore volume by cell. + /// @param[in] s saturation values (for all P phases) + /// @param[in] cmax max polymer concentration for cell + /// @param[in] dps dead pore space + /// @return total absorbed polymer mass. + double computePolymerAbsorbed(const Opm::PolymerProperties& polyprops, + const std::vector& pv, + const std::vector& s, + const std::vector& cmax, + const double dps); } // namespace Opm