OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Add support for PMISC

Browse Source 
Pressure effects are added to relative permeability, capillary pressure,
viscosity and density miscibility
This commit is contained in:
Tor Harald Sandve 2016-03-04 11:15:46 +01:00
parent a19607fabc
commit 948d985f56
4 changed files with 152 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -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<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const Opm::PhaseUsage pu);
// compute phase pressures.
std::vector<ADB>
computePressures(const ADB& po,
const ADB& sw,
const ADB& so,
const ADB& sg,
const ADB& ss) const;
};

122
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -177,17 +177,72 @@ namespace Opm {
BlackoilSolventModel<Grid>::variableStateExtractVars(const ReservoirState& x,
const std::vector<int>& indices,
std::vector<ADB>& 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<const V>(& 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<double> zero_selector(ss.value() + sg.value(), Selector<double>::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 <class Grid>
@ -881,6 +949,30 @@ namespace Opm {
}
template <class Grid>
std::vector<ADB>
BlackoilSolventModel<Grid>::
computePressures(const ADB& po,
const ADB& sw,
const ADB& so,
const ADB& sg,
const ADB& ss) const
{
std::vector<ADB> pressures = Base::computePressures(po, sw, so, sg);
// The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss)
std::vector<ADB> 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 <class Grid>
void
BlackoilSolventModel<Grid>::assemble(const ReservoirState& reservoir_state,

29
opm/autodiff/SolventPropsAdFromDeck.cpp
View File

@ -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<PmiscTable>(regionIdx);
// Copy data
const auto& po = pmiscTable.getOilPhasePressureColumn();
const auto& pmisc = pmiscTable.getMiscibilityColumn();
pmisc_[regionIdx] = NonuniformTableLinear<double>(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 {

8
opm/autodiff/SolventPropsAdFromDeck.hpp
View File

@ -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<NonuniformTableLinear<double> > mkro_;
std::vector<NonuniformTableLinear<double> > mkrsg_;
std::vector<NonuniformTableLinear<double> > misc_;
std::vector<NonuniformTableLinear<double> > pmisc_;
std::vector<NonuniformTableLinear<double> > sorwmis_;
std::vector<NonuniformTableLinear<double> > sgcwmis_;
std::vector<double> mix_param_viscosity_;