mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Support for liveoil in combination with solvent
- a solvent specific updateState is used to assure that the correct oil saturation is used to detect phase transision - presence of gas is compensated for in the oil phase
This commit is contained in:
Tor Harald Sandve
parent
1f4b9e56f8
commit
8be3ca7557
@ -99,12 +99,7 @@ namespace Opm {
|
|||||||
residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
|
residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
|
||||||
Base::material_name_.push_back("Solvent");
|
Base::material_name_.push_back("Solvent");
|
||||||
assert(solvent_pos_ == fluid_.numPhases());
|
assert(solvent_pos_ == fluid_.numPhases());
|
||||||
if (has_vapoil_) {
|
|
||||||
OPM_THROW(std::runtime_error, "Solvent option only works with dead gas\n");
|
|
||||||
}
|
|
||||||
|
|
||||||
residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
|
residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
|
||||||
|
|
||||||
stdWells().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
|
stdWells().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
|
||||||
}
|
}
|
||||||
if (is_miscible_) {
|
if (is_miscible_) {
|
||||||
@ -365,12 +360,12 @@ namespace Opm {
|
|||||||
}
|
}
|
||||||
|
|
||||||
const ADB& rs_perfcells = subset(state.rs, well_cells);
|
const ADB& rs_perfcells = subset(state.rs, well_cells);
|
||||||
|
const ADB& rv_perfcells = subset(state.rv, well_cells);
|
||||||
int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
|
int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
|
||||||
int oil_pos = fluid_.phaseUsage().phase_pos[Oil];
|
int oil_pos = fluid_.phaseUsage().phase_pos[Oil];
|
||||||
|
|
||||||
// remove contribution from the dissolved gas.
|
// remove contribution from the dissolved gas.
|
||||||
// TODO compensate for gas in the oil phase
|
const ADB cq_s_solvent = (isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction) * (cq_s[gas_pos] - (rs_perfcells * cq_s[oil_pos] / (ones - rs_perfcells * rv_perfcells)));
|
||||||
assert(!has_vapoil_);
|
|
||||||
const ADB cq_s_solvent = (isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction) * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]);
|
|
||||||
|
|
||||||
// Solvent contribution to the mass balance equation is given as a fraction
|
// Solvent contribution to the mass balance equation is given as a fraction
|
||||||
// of the gas contribution.
|
// of the gas contribution.
|
||||||
@ -382,62 +377,246 @@ namespace Opm {
|
|||||||
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
template <class Grid>
|
template <class Grid>
|
||||||
void BlackoilSolventModel<Grid>::updateState(const V& dx,
|
void
|
||||||
|
BlackoilSolventModel<Grid>::
|
||||||
|
updateState(const V& dx,
|
||||||
ReservoirState& reservoir_state,
|
ReservoirState& reservoir_state,
|
||||||
WellState& well_state)
|
WellState& well_state)
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
if (has_solvent_) {
|
|
||||||
// Extract solvent change.
|
|
||||||
const int np = fluid_.numPhases();
|
const int np = fluid_.numPhases();
|
||||||
const int nc = Opm::AutoDiffGrid::numCells(grid_);
|
const int nc = numCells(grid_);
|
||||||
|
const V null;
|
||||||
|
assert(null.size() == 0);
|
||||||
const V zero = V::Zero(nc);
|
const V zero = V::Zero(nc);
|
||||||
const int solvent_start = nc * np;
|
|
||||||
const V dss = subset(dx, Span(nc, 1, solvent_start));
|
|
||||||
|
|
||||||
// Create new dx with the dss part deleted.
|
// Extract parts of dx corresponding to each part.
|
||||||
V modified_dx = V::Zero(dx.size() - nc);
|
const V dp = subset(dx, Span(nc));
|
||||||
modified_dx.head(solvent_start) = dx.head(solvent_start);
|
int varstart = nc;
|
||||||
const int tail_len = dx.size() - solvent_start - nc;
|
const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
|
||||||
modified_dx.tail(tail_len) = dx.tail(tail_len);
|
varstart += dsw.size();
|
||||||
|
|
||||||
// Call base version.
|
const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
|
||||||
Base::updateState(modified_dx, reservoir_state, well_state);
|
varstart += dxvar.size();
|
||||||
|
|
||||||
// Update solvent.
|
const V dss = has_solvent_ ? subset(dx, Span(nc, 1, varstart)) : null;
|
||||||
|
varstart += dss.size();
|
||||||
|
|
||||||
|
// Extract well parts np phase rates + bhp
|
||||||
|
const V dwells = subset(dx, Span(Base::numWellVars(), 1, varstart));
|
||||||
|
varstart += dwells.size();
|
||||||
|
|
||||||
|
assert(varstart == dx.size());
|
||||||
|
|
||||||
|
// Pressure update.
|
||||||
|
const double dpmaxrel = dpMaxRel();
|
||||||
|
const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
|
||||||
|
const V absdpmax = dpmaxrel*p_old.abs();
|
||||||
|
const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
|
||||||
|
const V p = (p_old - dp_limited).max(zero);
|
||||||
|
std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
|
||||||
|
|
||||||
|
|
||||||
|
// Saturation updates.
|
||||||
|
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||||
|
const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
|
||||||
|
const double dsmax = dsMax();
|
||||||
|
|
||||||
|
// initialize with zeros
|
||||||
|
// if the phase is active the saturation are overwritten.
|
||||||
|
V so = zero;
|
||||||
|
V sw = zero;
|
||||||
|
V sg = zero;
|
||||||
|
V ss = zero;
|
||||||
|
|
||||||
|
// Appleyard chop process.
|
||||||
|
// We chop to large updates in saturations
|
||||||
|
{
|
||||||
|
V maxVal = zero;
|
||||||
|
V dso = zero;
|
||||||
|
if (active_[Water]){
|
||||||
|
maxVal = dsw.abs().max(maxVal);
|
||||||
|
dso = dso - dsw;
|
||||||
|
}
|
||||||
|
|
||||||
|
V dsg;
|
||||||
|
if (active_[Gas]){
|
||||||
|
dsg = Base::isSg_ * dxvar - Base::isRv_ * dsw;
|
||||||
|
maxVal = dsg.abs().max(maxVal);
|
||||||
|
dso = dso - dsg;
|
||||||
|
}
|
||||||
|
if (has_solvent_){
|
||||||
|
maxVal = dss.abs().max(maxVal);
|
||||||
|
}
|
||||||
|
|
||||||
|
maxVal = dso.abs().max(maxVal);
|
||||||
|
|
||||||
|
V step = dsmax/maxVal;
|
||||||
|
step = step.min(1.);
|
||||||
|
|
||||||
|
if (active_[Water]) {
|
||||||
|
const int pos = pu.phase_pos[ Water ];
|
||||||
|
const V sw_old = s_old.col(pos);
|
||||||
|
sw = sw_old - step * dsw;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (active_[Gas]) {
|
||||||
|
const int pos = pu.phase_pos[ Gas ];
|
||||||
|
const V sg_old = s_old.col(pos);
|
||||||
|
sg = sg_old - step * dsg;
|
||||||
|
}
|
||||||
|
if (has_solvent_) {
|
||||||
auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
|
auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
|
||||||
const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
|
const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
|
||||||
const V ss = (ss_old - dss).max(zero);
|
ss = ss_old - step * dss;
|
||||||
std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
|
}
|
||||||
|
|
||||||
// adjust oil saturation
|
const int pos = pu.phase_pos[ Oil ];
|
||||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
const V so_old = s_old.col(pos);
|
||||||
const int oilpos = pu.phase_pos[ Oil ];
|
so = so_old - step * dso;
|
||||||
|
}
|
||||||
|
|
||||||
|
auto ixg = sg < 0;
|
||||||
for (int c = 0; c < nc; ++c) {
|
for (int c = 0; c < nc; ++c) {
|
||||||
reservoir_state.saturation()[c*np + oilpos] = 1 - ss[c];
|
if (ixg[c]) {
|
||||||
if (pu.phase_used[ Gas ]) {
|
sw[c] = sw[c] / (1-sg[c]);
|
||||||
const int gaspos = pu.phase_pos[ Gas ];
|
so[c] = so[c] / (1-sg[c]);
|
||||||
reservoir_state.saturation()[c*np + oilpos] -= reservoir_state.saturation()[c*np + gaspos];
|
sg[c] = 0;
|
||||||
}
|
|
||||||
if (pu.phase_used[ Water ]) {
|
|
||||||
const int waterpos = pu.phase_pos[ Water ];
|
|
||||||
reservoir_state.saturation()[c*np + oilpos] -= reservoir_state.saturation()[c*np + waterpos];
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
auto ixo = so < 0;
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
if (ixo[c]) {
|
||||||
|
sw[c] = sw[c] / (1-so[c]);
|
||||||
|
sg[c] = sg[c] / (1-so[c]);
|
||||||
|
so[c] = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
auto ixw = sw < 0;
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
if (ixw[c]) {
|
||||||
|
so[c] = so[c] / (1-sw[c]);
|
||||||
|
sg[c] = sg[c] / (1-sw[c]);
|
||||||
|
sw[c] = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// The oil saturation is defined to
|
||||||
|
// fill the rest of the pore space.
|
||||||
|
// For convergence reasons oil saturations
|
||||||
|
// must be included in the appelyard copping
|
||||||
|
so = V::Constant(nc,1.0) - sw - sg - ss;
|
||||||
|
|
||||||
|
// Update rs and rv
|
||||||
|
const double drmaxrel = drMaxRel();
|
||||||
|
V rs;
|
||||||
|
if (has_disgas_) {
|
||||||
|
const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
|
||||||
|
const V drs = Base::isRs_ * dxvar;
|
||||||
|
const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
|
||||||
|
rs = rs_old - drs_limited;
|
||||||
|
}
|
||||||
|
V rv;
|
||||||
|
if (has_vapoil_) {
|
||||||
|
const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
|
||||||
|
const V drv = Base::isRv_ * dxvar;
|
||||||
|
const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
|
||||||
|
rv = rv_old - drv_limited;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Sg is used as primal variable for water only cells.
|
||||||
|
const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
|
||||||
|
auto watOnly = sw > (1 - epsilon);
|
||||||
|
|
||||||
|
// phase translation sg <-> rs
|
||||||
|
std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
|
||||||
|
|
||||||
|
if (has_disgas_) {
|
||||||
|
const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
|
||||||
|
const V rsSat = fluidRsSat(p, so, cells_);
|
||||||
|
// The obvious case
|
||||||
|
auto hasGas = (sg > 0 && Base::isRs_ == 0);
|
||||||
|
|
||||||
|
// Set oil saturated if previous rs is sufficiently large
|
||||||
|
const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
|
||||||
|
auto gasVaporized = ( (rs > rsSat * (1+epsilon) && Base::isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
|
||||||
|
auto useSg = watOnly || hasGas || gasVaporized;
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
if (useSg[c]) {
|
||||||
|
rs[c] = rsSat[c];
|
||||||
} else {
|
} else {
|
||||||
// Just forward call to base version.
|
primalVariable_[c] = PrimalVariables::RS;
|
||||||
Base::updateState(dx, reservoir_state, well_state);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
// phase transitions so <-> rv
|
||||||
|
if (has_vapoil_) {
|
||||||
|
|
||||||
|
// The gas pressure is needed for the rvSat calculations
|
||||||
|
const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
|
||||||
|
const V gaspress = computeGasPressure(p, sw, so, sg);
|
||||||
|
const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
|
||||||
|
const V rvSat = fluidRvSat(gaspress, so, cells_);
|
||||||
|
|
||||||
|
// The obvious case
|
||||||
|
auto hasOil = (so > 0 && Base::isRv_ == 0);
|
||||||
|
|
||||||
|
// Set oil saturated if previous rv is sufficiently large
|
||||||
|
const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
|
||||||
|
auto oilCondensed = ( (rv > rvSat * (1+epsilon) && Base::isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
|
||||||
|
auto useSg = watOnly || hasOil || oilCondensed;
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
if (useSg[c]) {
|
||||||
|
rv[c] = rvSat[c];
|
||||||
|
} else {
|
||||||
|
primalVariable_[c] = PrimalVariables::RV;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Update the reservoir_state
|
||||||
|
if (has_solvent_) {
|
||||||
|
auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
|
||||||
|
std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
|
||||||
|
}
|
||||||
|
|
||||||
|
if (active_[ Oil ]) {
|
||||||
|
const int pos = pu.phase_pos[ Oil ];
|
||||||
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
reservoir_state.saturation()[c*np + pos] = so[c];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Update the reservoir_state
|
||||||
|
if (has_disgas_) {
|
||||||
|
std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
|
||||||
|
}
|
||||||
|
|
||||||
|
if (has_vapoil_) {
|
||||||
|
std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
|
||||||
|
}
|
||||||
|
|
||||||
|
// TODO: gravity should be stored as a member
|
||||||
|
// const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
|
||||||
|
// asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state);
|
||||||
|
Base::updateWellState(dwells,well_state);
|
||||||
|
|
||||||
|
// Update phase conditions used for property calculations.
|
||||||
|
updatePhaseCondFromPrimalVariable();
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user