From 948d985f56a87f57e811e81e4280d9a4652340dc Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 4 Mar 2016 11:15:46 +0100 Subject: [PATCH] Add support for PMISC Pressure effects are added to relative permeability, capillary pressure, viscosity and density miscibility --- opm/autodiff/BlackoilSolventModel.hpp | 10 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 122 ++++++++++++++++++--- opm/autodiff/SolventPropsAdFromDeck.cpp | 29 +++++ opm/autodiff/SolventPropsAdFromDeck.hpp | 8 ++ 4 files changed, 152 insertions(+), 17 deletions(-) diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index e8cc9a063..bc2a72ebf 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -144,7 +144,6 @@ namespace Opm { using Base::wellsActive; using Base::wells; using Base::variableState; - using Base::computePressures; using Base::computeGasPressure; using Base::applyThresholdPressures; using Base::fluidRsSat; @@ -157,7 +156,6 @@ namespace Opm { using Base::dsMax; using Base::drMaxRel; using Base::maxResidualAllowed; - using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; @@ -237,6 +235,14 @@ namespace Opm { // compute density and viscosity using the ToddLongstaff mixing model void computeToddLongstaffMixing(std::vector& viscosity, std::vector& density, const std::vector& saturations, const Opm::PhaseUsage pu); + // compute phase pressures. + std::vector + computePressures(const ADB& po, + const ADB& sw, + const ADB& so, + const ADB& sg, + const ADB& ss) const; + }; diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index f7fda17b0..3f703a045 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -177,17 +177,72 @@ namespace Opm { BlackoilSolventModel::variableStateExtractVars(const ReservoirState& x, const std::vector& indices, std::vector& vars) const - { - SolutionState state = Base::variableStateExtractVars(x, indices, vars); - if (has_solvent_) { - state.solvent_saturation = std::move(vars[indices[Solvent]]); + { + // This is more or less a copy of the base class. Refactoring is needed in the base class + // to avoid unnecessary copying. + + //using namespace Opm::AutoDiffGrid; + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const Opm::PhaseUsage pu = fluid_.phaseUsage(); + + SolutionState state(fluid_.numPhases()); + + // Pressure. + state.pressure = std::move(vars[indices[Pressure]]); + + // Temperature cannot be a variable at this time (only constant). + const V temp = Eigen::Map(& x.temperature()[0], x.temperature().size()); + state.temperature = ADB::constant(temp); + + // Saturations + { + ADB so = ADB::constant(V::Ones(nc, 1)); + + if (active_[ Water ]) { + state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]); + const ADB& sw = state.saturation[pu.phase_pos[ Water ]]; + so -= sw; + } + if (has_solvent_) { + state.solvent_saturation = std::move(vars[indices[Solvent]]); + so -= state.solvent_saturation; + } + + if (active_[ Gas ]) { + // Define Sg Rs and Rv in terms of xvar. + // Xvar is only defined if gas phase is active + const ADB& xvar = vars[indices[Xvar]]; + ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ]; + sg = Base::isSg_*xvar + Base::isRv_*so; + so -= sg; + + if (active_[ Oil ]) { + // RS and RV is only defined if both oil and gas phase are active. + state.canonical_phase_pressures = computePressures(state.pressure, state.saturation[pu.phase_pos[ Water ]], so, sg, state.solvent_saturation); + const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_); + if (has_disgas_) { + state.rs = (1-Base::isRs_)*rsSat + Base::isRs_*xvar; + } else { + state.rs = rsSat; + } + const ADB rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_); + if (has_vapoil_) { + state.rv = (1-Base::isRv_)*rvSat + Base::isRv_*xvar; + } else { + state.rv = rvSat; + } + } + } + if (active_[ Oil ]) { // Note that so is never a primary variable. - const Opm::PhaseUsage pu = fluid_.phaseUsage(); - state.saturation[pu.phase_pos[ Oil ]] -= state.solvent_saturation; + state.saturation[pu.phase_pos[ Oil ]] = std::move(so); } } + // wells + Base::variableStateExtractWellsVars(indices, vars, state); return state; + } @@ -665,7 +720,9 @@ namespace Opm { Selector zero_selector(ss.value() + sg.value(), Selector::Zero); ADB F_solvent = zero_selector.select(ss, ss / (ss + sg)); - const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_); + const ADB& po = state.canonical_phase_pressures[ Oil ]; + const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_) + * solvent_props_.pressureMiscibilityFunction(po, cells_); assert(active_[ Oil ]); assert(active_[ Gas ]); @@ -769,16 +826,27 @@ namespace Opm { // Compute effective viscosities and densities computeToddLongstaffMixing(viscosity, density, effective_saturations, pu); - // Store the computed volume factors and viscosities - b_eff_[pu.phase_pos[ Water ]] = bw; - b_eff_[pu.phase_pos[ Oil ]] = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs); - b_eff_[pu.phase_pos[ Gas ]] = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv); - b_eff_[solvent_pos_] = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_); + // compute the volume factors from the densities + const ADB b_eff_o = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs); + const ADB b_eff_g = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv); + const ADB b_eff_s = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_); + // account for pressure effects and store the computed volume factors and viscosities + const V ones = V::Constant(nc, 1.0); + const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_); + + b_eff_[pu.phase_pos[ Oil ]] = pmisc * b_eff_o + (ones - pmisc) * bo; + b_eff_[pu.phase_pos[ Gas ]] = pmisc * b_eff_g + (ones - pmisc) * bg; + b_eff_[solvent_pos_] = pmisc * b_eff_s + (ones - pmisc) * bs; + + // keep the mu*b interpolation + mu_eff_[pu.phase_pos[ Oil ]] = b_eff_[pu.phase_pos[ Oil ]] / (pmisc * b_eff_o / viscosity[pu.phase_pos[ Oil ]] + (ones - pmisc) * bo / mu_o); + mu_eff_[pu.phase_pos[ Gas ]] = b_eff_[pu.phase_pos[ Gas ]] / (pmisc * b_eff_g / viscosity[pu.phase_pos[ Gas ]] + (ones - pmisc) * bg / mu_g); + mu_eff_[solvent_pos_] = b_eff_[solvent_pos_] / (pmisc * b_eff_s / viscosity[solvent_pos_] + (ones - pmisc) * bs / mu_s); + + // for water the pure values are used mu_eff_[pu.phase_pos[ Water ]] = mu_w; - mu_eff_[pu.phase_pos[ Oil ]] = viscosity[pu.phase_pos[ Oil ]]; - mu_eff_[pu.phase_pos[ Gas ]] = viscosity[pu.phase_pos[ Gas ]]; - mu_eff_[solvent_pos_] = viscosity[solvent_pos_]; + b_eff_[pu.phase_pos[ Water ]] = bw; } template @@ -881,6 +949,30 @@ namespace Opm { } + template + std::vector + BlackoilSolventModel:: + computePressures(const ADB& po, + const ADB& sw, + const ADB& so, + const ADB& sg, + const ADB& ss) const + { + std::vector pressures = Base::computePressures(po, sw, so, sg); + // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss) + std::vector pressures_imisc = Base::computePressures(po, sw, so, sg + ss); + + // Pressure effects on capillary pressure miscibility + const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_); + // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent + // to changing the gas pressure directly. + const int nc = cells_.size(); + const V ones = V::Constant(nc, 1.0); + pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas])); + return pressures; + } + + template void BlackoilSolventModel::assemble(const ReservoirState& reservoir_state, diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index 2fd967021..b94d1a64c 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -172,6 +172,25 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n"); } + const TableContainer& pmiscTables = tables->getPmiscTables(); + if (!pmiscTables.empty()) { + + int numRegions = pmiscTables.size(); + + // resize the attributes of the object + pmisc_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::PmiscTable& pmiscTable = pmiscTables.getTable(regionIdx); + + // Copy data + const auto& po = pmiscTable.getOilPhasePressureColumn(); + const auto& pmisc = pmiscTable.getMiscibilityColumn(); + + pmisc_[regionIdx] = NonuniformTableLinear(po, pmisc); + + } + } + // miscible relative permeability multipleiers const TableContainer& msfnTables = tables->getMsfnTables(); if (!msfnTables.empty()) { @@ -339,6 +358,16 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction, return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_); } +ADB SolventPropsAdFromDeck::pressureMiscibilityFunction(const ADB& po, + const Cells& cells) const +{ + if (pmisc_.size() > 0) { + return SolventPropsAdFromDeck::makeADBfromTables(po, cells, cellMiscRegionIdx_, pmisc_); + } + // return ones if not specified i.e. no effect. + return ADB::constant(V::Constant(po.size(), 1.0)); +} + ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, const Cells& cells) const { diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index 4aa69b4c4..cdc94022e 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/opm/autodiff/SolventPropsAdFromDeck.hpp @@ -108,6 +108,13 @@ public: ADB miscibilityFunction(const ADB& solventFraction, const Cells& cells) const; + /// Pressure dependent miscibility function + /// \param[in] solventFraction Array of n oil phase pressure . + /// \param[in] cells Array of n cell indices to be associated with the pressure values. + /// \return Array of n miscibility values + ADB pressureMiscibilityFunction(const ADB& po, + const Cells& cells) const; + /// Miscible critical gas saturation function /// \param[in] Sw Array of n water saturation values. /// \param[in] cells Array of n cell indices to be associated with the saturation values. @@ -177,6 +184,7 @@ private: std::vector > mkro_; std::vector > mkrsg_; std::vector > misc_; + std::vector > pmisc_; std::vector > sorwmis_; std::vector > sgcwmis_; std::vector mix_param_viscosity_;