mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #577 from totto82/solventMiscible
Add miscible effects to the solvent model
This commit is contained in:
Atgeirr Flø Rasmussen
commit
eb228835c2
@ -454,6 +454,8 @@ namespace Opm {
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const ADB& kr ,
|
||||
const ADB& mu ,
|
||||
const ADB& rho ,
|
||||
const ADB& p ,
|
||||
const SolutionState& state );
|
||||
|
||||
|
||||
@ -752,7 +752,7 @@ namespace detail {
|
||||
for (int phase = 0; phase < maxnp; ++phase) {
|
||||
if (active_[ phase ]) {
|
||||
const int pos = pu.phase_pos[ phase ];
|
||||
rq_[pos].b = fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
|
||||
rq_[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
|
||||
rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
|
||||
// OPM_AD_DUMP(rq_[pos].b);
|
||||
// OPM_AD_DUMP(rq_[pos].accum[aix]);
|
||||
@ -989,7 +989,10 @@ namespace detail {
|
||||
const std::vector<ADB> kr = asImpl().computeRelPerm(state);
|
||||
#pragma omp parallel for schedule(static)
|
||||
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
|
||||
asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
const ADB mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
|
||||
const ADB rho = asImpl().fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv);
|
||||
asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
|
||||
|
||||
residual_.material_balance_eq[ phaseIdx ] =
|
||||
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
|
||||
@ -2404,25 +2407,21 @@ namespace detail {
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template <class Grid, class Implementation>
|
||||
void
|
||||
BlackoilModelBase<Grid, Implementation>::computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const ADB& kr ,
|
||||
const ADB& mu ,
|
||||
const ADB& rho ,
|
||||
const ADB& phasePressure,
|
||||
const SolutionState& state)
|
||||
{
|
||||
// Compute and store mobilities.
|
||||
const int canonicalPhaseIdx = canph_[ actph ];
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
const ADB tr_mult = transMult(state.pressure);
|
||||
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv, cond);
|
||||
rq_[ actph ].mob = tr_mult * kr / mu;
|
||||
|
||||
// Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST.
|
||||
const ADB rho = fluidDensity(canonicalPhaseIdx, rq_[actph].b, state.rs, state.rv);
|
||||
const ADB rhoavg = ops_.caver * rho;
|
||||
rq_[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
|
||||
if (use_threshold_pressure_) {
|
||||
|
||||
@ -63,6 +63,7 @@ namespace Opm {
|
||||
/// \param[in] has_vapoil turn on vaporized oil feature
|
||||
/// \param[in] terminal_output request output to cout/cerr
|
||||
/// \param[in] has_solvent turn on solvent feature
|
||||
/// \param[in] is_miscible turn on miscible feature
|
||||
BlackoilSolventModel(const typename Base::ModelParameters& param,
|
||||
const Grid& grid,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
@ -75,7 +76,8 @@ namespace Opm {
|
||||
const bool has_disgas,
|
||||
const bool has_vapoil,
|
||||
const bool terminal_output,
|
||||
const bool has_solvent);
|
||||
const bool has_solvent,
|
||||
const bool is_miscible);
|
||||
|
||||
/// Apply an update to the primary variables, chopped if appropriate.
|
||||
/// \param[in] dx updates to apply to primary variables
|
||||
@ -106,6 +108,10 @@ namespace Opm {
|
||||
const bool has_solvent_;
|
||||
const int solvent_pos_;
|
||||
const SolventPropsAdFromDeck& solvent_props_;
|
||||
const bool is_miscible_;
|
||||
std::vector<ADB> mu_eff_;
|
||||
std::vector<ADB> b_eff_;
|
||||
|
||||
|
||||
// Need to declare Base members we want to use here.
|
||||
using Base::grid_;
|
||||
@ -141,9 +147,6 @@ namespace Opm {
|
||||
using Base::computePressures;
|
||||
using Base::computeGasPressure;
|
||||
using Base::applyThresholdPressures;
|
||||
using Base::fluidViscosity;
|
||||
using Base::fluidReciprocFVF;
|
||||
using Base::fluidDensity;
|
||||
using Base::fluidRsSat;
|
||||
using Base::fluidRvSat;
|
||||
using Base::poroMult;
|
||||
@ -162,6 +165,28 @@ namespace Opm {
|
||||
std::vector<ADB>
|
||||
computeRelPerm(const SolutionState& state) const;
|
||||
|
||||
ADB
|
||||
fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
const ADB& temp ,
|
||||
const ADB& rs ,
|
||||
const ADB& rv ,
|
||||
const std::vector<PhasePresence>& cond) const;
|
||||
|
||||
ADB
|
||||
fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
const ADB& temp ,
|
||||
const ADB& rs ,
|
||||
const ADB& rv ,
|
||||
const std::vector<PhasePresence>& cond) const;
|
||||
|
||||
ADB
|
||||
fluidDensity(const int phase,
|
||||
const ADB& b,
|
||||
const ADB& rs,
|
||||
const ADB& rv) const;
|
||||
|
||||
void
|
||||
makeConstantState(SolutionState& state) const;
|
||||
|
||||
@ -198,12 +223,19 @@ namespace Opm {
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const ADB& kr ,
|
||||
const ADB& mu ,
|
||||
const ADB& rho ,
|
||||
const ADB& p ,
|
||||
const SolutionState& state );
|
||||
|
||||
const std::vector<PhasePresence>
|
||||
phaseCondition() const {return this->phaseCondition_;}
|
||||
|
||||
// compute effective viscosities (mu_eff_) and effective b factors (b_eff_) using the ToddLongstaff model
|
||||
void computeEffectiveProperties(const SolutionState& state);
|
||||
|
||||
// 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);
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -82,12 +82,15 @@ namespace Opm {
|
||||
const bool has_disgas,
|
||||
const bool has_vapoil,
|
||||
const bool terminal_output,
|
||||
const bool has_solvent)
|
||||
const bool has_solvent,
|
||||
const bool is_miscible)
|
||||
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
|
||||
eclState, has_disgas, has_vapoil, terminal_output),
|
||||
has_solvent_(has_solvent),
|
||||
solvent_pos_(detail::solventPos(fluid.phaseUsage())),
|
||||
solvent_props_(solvent_props)
|
||||
solvent_props_(solvent_props),
|
||||
is_miscible_(is_miscible)
|
||||
|
||||
{
|
||||
if (has_solvent_) {
|
||||
|
||||
@ -101,7 +104,10 @@ namespace Opm {
|
||||
}
|
||||
|
||||
residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
|
||||
|
||||
}
|
||||
if (is_miscible_) {
|
||||
mu_eff_.resize(fluid_.numPhases() + 1, ADB::null());
|
||||
b_eff_.resize(fluid_.numPhases() + 1, ADB::null());
|
||||
}
|
||||
}
|
||||
|
||||
@ -203,7 +209,8 @@ namespace Opm {
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
|
||||
const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
|
||||
rq_[solvent_pos_].b = solvent_props_.bSolvent(pg,cells_);
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
rq_[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond);
|
||||
rq_[solvent_pos_].accum[aix] = pv_mult * rq_[solvent_pos_].b * ss;
|
||||
}
|
||||
}
|
||||
@ -354,33 +361,42 @@ namespace Opm {
|
||||
DataBlock b(nperf, pu.num_phases);
|
||||
std::vector<double> rsmax_perf(nperf, 0.0);
|
||||
std::vector<double> rvmax_perf(nperf, 0.0);
|
||||
const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
|
||||
if (pu.phase_used[BlackoilPhases::Aqua]) {
|
||||
const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
|
||||
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
|
||||
}
|
||||
|
||||
assert(active_[Oil]);
|
||||
assert(active_[Gas]);
|
||||
const ADB perf_rv = subset(state.rv, well_cells);
|
||||
const ADB perf_rs = subset(state.rs, well_cells);
|
||||
const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
|
||||
if (pu.phase_used[BlackoilPhases::Liquid]) {
|
||||
const ADB perf_rs = subset(state.rs, well_cells);
|
||||
const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
|
||||
//const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value();
|
||||
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
|
||||
const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
|
||||
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
|
||||
}
|
||||
V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
|
||||
for (int phase = 1; phase < pu.num_phases; ++phase) {
|
||||
if ( phase != pu.phase_pos[BlackoilPhases::Vapour]) {
|
||||
if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
|
||||
continue; // the gas surface density is added after the solvent is accounted for.
|
||||
}
|
||||
surf_dens_copy += superset(fluid_.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
|
||||
}
|
||||
|
||||
if (pu.phase_used[BlackoilPhases::Vapour]) {
|
||||
const ADB perf_rv = subset(state.rv, well_cells);
|
||||
// Unclear wether the effective or the pure values should be used for the wells
|
||||
// the current usage of unmodified properties values gives best match.
|
||||
//V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value();
|
||||
V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
|
||||
V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
|
||||
if (has_solvent_) {
|
||||
|
||||
const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
|
||||
//const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
|
||||
|
||||
// A weighted sum of the b-factors of gas and solvent are used.
|
||||
const int nc = Opm::AutoDiffGrid::numCells(grid_);
|
||||
|
||||
@ -404,6 +420,7 @@ namespace Opm {
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
|
||||
|
||||
bg = bg * (ones - F_solvent);
|
||||
@ -510,53 +527,114 @@ namespace Opm {
|
||||
BlackoilSolventModel<Grid>::computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const ADB& kr ,
|
||||
const ADB& mu ,
|
||||
const ADB& rho ,
|
||||
const ADB& phasePressure,
|
||||
const SolutionState& state)
|
||||
{
|
||||
Base::computeMassFlux(actph, transi, kr, phasePressure, state);
|
||||
|
||||
const int canonicalPhaseIdx = canph_[ actph ];
|
||||
// make a copy to make it possible to modify it
|
||||
ADB kr_mod = kr;
|
||||
if (canonicalPhaseIdx == Gas) {
|
||||
if (has_solvent_) {
|
||||
const int nc = Opm::UgGridHelpers::numCells(grid_);
|
||||
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
const ADB zero = ADB::constant(V::Zero(nc));
|
||||
const V ones = V::Constant(nc, 1.0);
|
||||
|
||||
const ADB& ss = state.solvent_saturation;
|
||||
const ADB& sg = (active_[ Gas ]
|
||||
? state.saturation[ pu.phase_pos[ Gas ] ]
|
||||
: zero);
|
||||
|
||||
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
|
||||
ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
|
||||
V ones = V::Constant(nc, 1.0);
|
||||
const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg));
|
||||
|
||||
const ADB tr_mult = transMult(state.pressure);
|
||||
const ADB mu = solvent_props_.muSolvent(phasePressure,cells_);
|
||||
// Compute solvent properties
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond);
|
||||
ADB rho_s = fluidDensity(Solvent,rq_[solvent_pos_].b, state.rs, state.rv);
|
||||
|
||||
rq_[solvent_pos_].mob = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * tr_mult * kr / mu;
|
||||
// Compute solvent relperm and mass flux
|
||||
ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod;
|
||||
Base::computeMassFlux(solvent_pos_, transi, krs, mu_s, rho_s, phasePressure, state);
|
||||
|
||||
const ADB rho_solvent = solvent_props_.solventSurfaceDensity(cells_) * rq_[solvent_pos_].b;
|
||||
const ADB rhoavg_solvent = ops_.caver * rho_solvent;
|
||||
rq_[ solvent_pos_ ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg_solvent * (ops_.ngrad * geo_.z().matrix()));
|
||||
// Modify gas relperm
|
||||
kr_mod = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * kr_mod;
|
||||
|
||||
UpwindSelector<double> upwind_solvent(grid_, ops_, rq_[solvent_pos_].dh.value());
|
||||
// Compute solvent flux.
|
||||
rq_[solvent_pos_].mflux = upwind_solvent.select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh);
|
||||
|
||||
// Update gas mobility and flux
|
||||
rq_[actph].mob = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * rq_[actph].mob;
|
||||
|
||||
const ADB& b = rq_[ actph ].b;
|
||||
const ADB& mob = rq_[ actph ].mob;
|
||||
const ADB& dh = rq_[ actph ].dh;
|
||||
UpwindSelector<double> upwind_gas(grid_, ops_, dh.value());
|
||||
rq_[ actph ].mflux = upwind_gas.select(b * mob) * (transi * dh);
|
||||
}
|
||||
}
|
||||
// Compute mobility and flux
|
||||
Base::computeMassFlux(actph, transi, kr_mod, mu, rho, phasePressure, state);
|
||||
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
ADB
|
||||
BlackoilSolventModel<Grid>::fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
const ADB& temp ,
|
||||
const ADB& rs ,
|
||||
const ADB& rv ,
|
||||
const std::vector<PhasePresence>& cond) const
|
||||
{
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
if (phase == Solvent) {
|
||||
if (!is_miscible_) {
|
||||
return solvent_props_.muSolvent(p, cells_);
|
||||
} else {
|
||||
return mu_eff_[solvent_pos_];
|
||||
}
|
||||
|
||||
} else {
|
||||
if (!is_miscible_) {
|
||||
return Base::fluidViscosity(phase, p, temp, rs, rv, cond);
|
||||
} else {
|
||||
return mu_eff_[pu.phase_pos[ phase ]];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
ADB
|
||||
BlackoilSolventModel<Grid>::fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
const ADB& temp ,
|
||||
const ADB& rs ,
|
||||
const ADB& rv ,
|
||||
const std::vector<PhasePresence>& cond) const
|
||||
{
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
if (phase == Solvent) {
|
||||
if (!is_miscible_) {
|
||||
return solvent_props_.bSolvent(p, cells_);
|
||||
} else {
|
||||
return b_eff_[solvent_pos_];
|
||||
}
|
||||
|
||||
} else {
|
||||
if (!is_miscible_) {
|
||||
return Base::fluidReciprocFVF(phase, p, temp, rs, rv, cond);
|
||||
} else {
|
||||
return b_eff_[pu.phase_pos[ phase ]];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
ADB
|
||||
BlackoilSolventModel<Grid>::fluidDensity(const int phase,
|
||||
const ADB& b,
|
||||
const ADB& rs,
|
||||
const ADB& rv) const
|
||||
{
|
||||
if (phase == Solvent && has_solvent_) {
|
||||
return solvent_props_.solventSurfaceDensity(cells_) * b;
|
||||
} else {
|
||||
return Base::fluidDensity(phase, b, rs, rv);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
std::vector<ADB>
|
||||
BlackoilSolventModel<Grid>::computeRelPerm(const SolutionState& state) const
|
||||
@ -581,13 +659,227 @@ namespace Opm {
|
||||
|
||||
if (has_solvent_) {
|
||||
const ADB& ss = state.solvent_saturation;
|
||||
return fluid_.relperm(sw, so, sg+ss, cells_);
|
||||
if (is_miscible_) {
|
||||
|
||||
std::vector<ADB> relperm = fluid_.relperm(sw, so, sg+ss, cells_);
|
||||
|
||||
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_);
|
||||
|
||||
assert(active_[ Oil ]);
|
||||
assert(active_[ Gas ]);
|
||||
|
||||
const ADB sn = ss + so + sg;
|
||||
|
||||
// adjust endpoints
|
||||
const V sgcr = fluid_.scaledCriticalGasSaturations(cells_);
|
||||
const V sogcr = fluid_.scaledCriticalOilinGasSaturations(cells_);
|
||||
const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
|
||||
const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
|
||||
|
||||
const V ones = V::Constant(nc, 1.0);
|
||||
ADB sor = misc * sorwmis + (ones - misc) * sogcr;
|
||||
ADB sgc = misc * sgcwmis + (ones - misc) * sgcr;
|
||||
|
||||
// avoid negative values
|
||||
Selector<double> negative_selector(sgc.value(), Selector<double>::LessEqualZero);
|
||||
sgc = negative_selector.select(zero, sgc);
|
||||
|
||||
const ADB ssg = ss + sg - sgc;
|
||||
const ADB sn_eff = sn - sor - sgc;
|
||||
|
||||
// avoid negative values
|
||||
Selector<double> zeroSn_selector(sn_eff.value(), Selector<double>::Zero);
|
||||
const ADB F_totalGas = zeroSn_selector.select( zero, ssg / sn_eff);
|
||||
|
||||
const ADB mkrgt = solvent_props_.miscibleSolventGasRelPermMultiplier(F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
|
||||
const ADB mkro = solvent_props_.miscibleOilRelPermMultiplier(ones - F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
|
||||
|
||||
const V eps = V::Constant(nc, 1e-7);
|
||||
Selector<double> noOil_selector(so.value()-eps, Selector<double>::LessEqualZero);
|
||||
relperm[Gas] = noOil_selector.select(relperm[Gas], (ones - misc) * relperm[Gas] + misc * mkrgt);
|
||||
relperm[Oil] = noOil_selector.select(relperm[Oil], (ones - misc) * relperm[Oil] + misc * mkro);
|
||||
return relperm;
|
||||
} else {
|
||||
return fluid_.relperm(sw, so, sg+ss, cells_);
|
||||
}
|
||||
} else {
|
||||
return fluid_.relperm(sw, so, sg, cells_);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
BlackoilSolventModel<Grid>::computeEffectiveProperties(const SolutionState& state)
|
||||
{
|
||||
// Viscosity
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
const int np = fluid_.numPhases();
|
||||
const int nc = Opm::UgGridHelpers::numCells(grid_);
|
||||
const ADB zero = ADB::constant(V::Zero(nc));
|
||||
|
||||
const ADB& pw = state.canonical_phase_pressures[pu.phase_pos[Water]];
|
||||
const ADB& po = state.canonical_phase_pressures[pu.phase_pos[Oil]];
|
||||
const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
|
||||
const ADB mu_w = fluid_.muWat(pw, state.temperature, cells_);
|
||||
const ADB mu_o = fluid_.muOil(po, state.temperature, state.rs, cond, cells_);
|
||||
const ADB mu_g = fluid_.muGas(pg, state.temperature, state.rv, cond, cells_);
|
||||
const ADB mu_s = solvent_props_.muSolvent(pg,cells_);
|
||||
std::vector<ADB> viscosity(np + 1, ADB::null());
|
||||
viscosity[pu.phase_pos[Oil]] = mu_o;
|
||||
viscosity[pu.phase_pos[Gas]] = mu_g;
|
||||
viscosity[pu.phase_pos[Water]] = mu_w;
|
||||
viscosity[solvent_pos_] = mu_s;
|
||||
|
||||
// Density
|
||||
const ADB bw = fluid_.bWat(pw, state.temperature, cells_);
|
||||
const ADB bo = fluid_.bOil(po, state.temperature, state.rs, cond, cells_);
|
||||
const ADB bg = fluid_.bGas(pg, state.temperature, state.rv, cond, cells_);
|
||||
const ADB bs = solvent_props_.bSolvent(pg, cells_);
|
||||
|
||||
std::vector<ADB> density(np + 1, ADB::null());
|
||||
density[pu.phase_pos[Oil]] = fluidDensity(Oil, bo, state.rs, state.rv);
|
||||
density[pu.phase_pos[Gas]] = fluidDensity(Gas, bg, state.rs, state.rv);
|
||||
density[pu.phase_pos[Water]] = fluidDensity(Water, bw, state.rs, state.rv);
|
||||
density[solvent_pos_] = fluidDensity(Solvent, bs, state.rs, state.rv);
|
||||
|
||||
// Effective saturations
|
||||
const ADB& ss = state.solvent_saturation;
|
||||
const ADB& so = state.saturation[ pu.phase_pos[ Oil ] ];
|
||||
const ADB& sg = (active_[ Gas ]
|
||||
? state.saturation[ pu.phase_pos[ Gas ] ]
|
||||
: zero);
|
||||
const ADB& sw = (active_[ Water ]
|
||||
? state.saturation[ pu.phase_pos[ Water ] ]
|
||||
: zero);
|
||||
|
||||
const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
|
||||
const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
|
||||
|
||||
std::vector<ADB> effective_saturations (np + 1, ADB::null());
|
||||
effective_saturations[pu.phase_pos[Oil]] = so - sorwmis;
|
||||
effective_saturations[pu.phase_pos[Gas]] = sg - sgcwmis;
|
||||
effective_saturations[pu.phase_pos[Water]] = sw;
|
||||
effective_saturations[solvent_pos_] = ss - sgcwmis;
|
||||
|
||||
// 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_);
|
||||
|
||||
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_];
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
BlackoilSolventModel<Grid>::computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const Opm::PhaseUsage pu)
|
||||
{
|
||||
const int nc = cells_.size();
|
||||
const V ones = V::Constant(nc, 1.0);
|
||||
const ADB zero = ADB::constant(V::Zero(nc));
|
||||
|
||||
// Calculation of effective saturations
|
||||
ADB so_eff = saturations[pu.phase_pos[ Oil ]];
|
||||
ADB sg_eff = saturations[pu.phase_pos[ Gas ]];
|
||||
ADB ss_eff = saturations[solvent_pos_];
|
||||
|
||||
// Avoid negative values
|
||||
Selector<double> negative_selectorSo(so_eff.value(), Selector<double>::LessZero);
|
||||
Selector<double> negative_selectorSg(sg_eff.value(), Selector<double>::LessZero);
|
||||
Selector<double> negative_selectorSs(ss_eff.value(), Selector<double>::LessZero);
|
||||
so_eff = negative_selectorSo.select(zero, so_eff);
|
||||
sg_eff = negative_selectorSg.select(zero, sg_eff);
|
||||
ss_eff = negative_selectorSs.select(zero, ss_eff);
|
||||
|
||||
// Make copy of the pure viscosities
|
||||
const ADB mu_o = viscosity[pu.phase_pos[ Oil ]];
|
||||
const ADB mu_g = viscosity[pu.phase_pos[ Gas ]];
|
||||
const ADB mu_s = viscosity[solvent_pos_];
|
||||
|
||||
const ADB sn_eff = so_eff + sg_eff + ss_eff;
|
||||
const ADB sos_eff = so_eff + ss_eff;
|
||||
const ADB ssg_eff = ss_eff + sg_eff;
|
||||
|
||||
// Avoid division by zero
|
||||
Selector<double> zero_selectorSos(sos_eff.value(), Selector<double>::Zero);
|
||||
Selector<double> zero_selectorSsg(ssg_eff.value(), Selector<double>::Zero);
|
||||
Selector<double> zero_selectorSn(sn_eff.value(), Selector<double>::Zero);
|
||||
|
||||
const ADB mu_s_pow = pow(mu_s,0.25);
|
||||
const ADB mu_o_pow = pow(mu_o,0.25);
|
||||
const ADB mu_g_pow = pow(mu_g,0.25);
|
||||
|
||||
const ADB mu_mos = zero_selectorSos.select(mu_o + mu_s, mu_o * mu_s / pow( ( (so_eff / sos_eff) * mu_s_pow) + ( (ss_eff / sos_eff) * mu_o_pow) , 4.0));
|
||||
const ADB mu_msg = zero_selectorSsg.select(mu_g + mu_s , mu_g * mu_s / pow( ( (sg_eff / ssg_eff) * mu_s_pow) + ( (ss_eff / ssg_eff) * mu_g_pow) , 4.0));
|
||||
const ADB mu_m = zero_selectorSn.select(mu_s + mu_o + mu_g, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow * mu_g_pow)
|
||||
+ ( (ss_eff / sn_eff) * mu_o_pow * mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
|
||||
// Mixing parameter for viscosity
|
||||
const double mix_param_mu = solvent_props_.mixingParameterViscosity();
|
||||
|
||||
// Update viscosities
|
||||
viscosity[pu.phase_pos[ Oil ]] = pow(mu_o,1.0 - mix_param_mu) * pow(mu_mos,mix_param_mu);
|
||||
viscosity[pu.phase_pos[ Gas ]] = pow(mu_g,1.0 - mix_param_mu) * pow(mu_msg,mix_param_mu);
|
||||
viscosity[solvent_pos_] = pow(mu_s,1.0 - mix_param_mu) * pow(mu_m,mix_param_mu);
|
||||
|
||||
// Density
|
||||
ADB& rho_o = density[pu.phase_pos[ Oil ]];
|
||||
ADB& rho_g = density[pu.phase_pos[ Gas ]];
|
||||
ADB& rho_s = density[solvent_pos_];
|
||||
|
||||
// mixing parameter for density
|
||||
const double mix_param_rho = solvent_props_.mixingParameterDensity();
|
||||
|
||||
// compute effective viscosities for density calculations. These have to
|
||||
// be recomputed as a different mixing parameter may be used.
|
||||
const ADB mu_o_eff = pow(mu_o,1.0 - mix_param_rho) * pow(mu_mos,mix_param_rho);
|
||||
const ADB mu_g_eff = pow(mu_g,1.0 - mix_param_rho) * pow(mu_msg,mix_param_rho);
|
||||
const ADB mu_s_eff = pow(mu_s,1.0 - mix_param_rho) * pow(mu_m,mix_param_rho);
|
||||
|
||||
const ADB sog_eff = so_eff + sg_eff;
|
||||
const ADB sof = so_eff / sog_eff;
|
||||
const ADB sgf = sg_eff / sog_eff;
|
||||
|
||||
// Effective densities
|
||||
const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) );
|
||||
const ADB mu_o_eff_pow = pow(mu_o_eff,0.25);
|
||||
const ADB mu_g_eff_pow = pow(mu_g_eff,0.25);
|
||||
const ADB mu_s_eff_pow = pow(mu_s_eff,0.25);
|
||||
const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow));
|
||||
const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow));
|
||||
const ADB sfraction_se = (mu_sog_pow - ( mu_o_pow * mu_g_pow * mu_s_pow / mu_s_eff_pow) ) / ( mu_sog_pow - (mu_o_pow * mu_g_pow));
|
||||
const ADB rho_o_eff = (rho_o * sfraction_oe) + (rho_s * (ones - sfraction_oe));
|
||||
const ADB rho_g_eff = (rho_g * sfraction_ge) + (rho_s * (ones - sfraction_ge));
|
||||
const ADB rho_s_eff = (rho_s * sfraction_se) + (rho_g * sgf * (ones - sfraction_se)) + (rho_o * sof * (ones - sfraction_se));
|
||||
|
||||
// Avoid division by zero for equal mobilities. For equal mobilities the effecitive density is calculated
|
||||
// based on the saturation fraction directly.
|
||||
Selector<double> unitGasSolventMobilityRatio_selector(mu_s.value() - mu_g.value(), Selector<double>::Zero);
|
||||
Selector<double> unitOilSolventMobilityRatio_selector(mu_s.value() - mu_o.value(), Selector<double>::Zero);
|
||||
|
||||
// Effective densities when the mobilities are equal
|
||||
const ADB rho_m = zero_selectorSn.select(zero, (rho_o * so_eff / sn_eff) + (rho_g * sg_eff / sn_eff) + (rho_s * ss_eff / sn_eff));
|
||||
const ADB rho_o_eff_simple = ((ones - mix_param_rho) * rho_o) + (mix_param_rho * rho_m);
|
||||
const ADB rho_g_eff_simple = ((ones - mix_param_rho) * rho_g) + (mix_param_rho * rho_m);
|
||||
//const ADB rho_s_eff_simple = ((ones - mix_param_rho) * rho_s) + (mix_param_rho * rho_m);
|
||||
|
||||
// Update densities
|
||||
rho_o = unitOilSolventMobilityRatio_selector.select(rho_o_eff_simple, rho_o_eff);
|
||||
rho_g = unitGasSolventMobilityRatio_selector.select(rho_g_eff_simple, rho_g_eff);
|
||||
rho_s = rho_s_eff;
|
||||
|
||||
}
|
||||
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
@ -611,9 +903,16 @@ namespace Opm {
|
||||
makeConstantState(state0);
|
||||
// Compute initial accumulation contributions
|
||||
// and well connection pressures.
|
||||
if (is_miscible_) {
|
||||
computeEffectiveProperties(state0);
|
||||
}
|
||||
|
||||
computeAccum(state0, 0);
|
||||
computeWellConnectionPressures(state0, well_state);
|
||||
}
|
||||
if (is_miscible_) {
|
||||
computeEffectiveProperties(state);
|
||||
}
|
||||
|
||||
// -------- Mass balance equations --------
|
||||
assembleMassBalanceEq(state);
|
||||
@ -661,14 +960,15 @@ namespace Opm {
|
||||
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
|
||||
ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
|
||||
V ones = V::Constant(nperf,1.0);
|
||||
|
||||
b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
|
||||
b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells));
|
||||
|
||||
}
|
||||
if (param_.solve_welleq_initially_ && initial_assembly) {
|
||||
// solve the well equations as a pre-processing step
|
||||
Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
|
||||
}
|
||||
|
||||
Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
|
||||
Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
|
||||
Base::addWellFluxEq(cq_s, state);
|
||||
|
||||
@ -131,6 +131,7 @@ namespace Opm
|
||||
bool has_solvent_;
|
||||
DeckConstPtr deck_;
|
||||
SolventPropsAdFromDeck solvent_props_;
|
||||
bool is_miscible_;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -54,9 +54,10 @@ namespace Opm
|
||||
, has_solvent_(has_solvent)
|
||||
, deck_(deck)
|
||||
, solvent_props_(solvent_props)
|
||||
, is_miscible_(false)
|
||||
{
|
||||
if(deck->hasKeyword("MISCIBLE")) {
|
||||
std::cerr << "MISICIBLE keyword is present. Mixing is not currently supported" << std::endl;
|
||||
is_miscible_ = true;
|
||||
}
|
||||
}
|
||||
|
||||
@ -80,7 +81,8 @@ namespace Opm
|
||||
BaseType::has_disgas_,
|
||||
BaseType::has_vapoil_,
|
||||
BaseType::terminal_output_,
|
||||
has_solvent_));
|
||||
has_solvent_,
|
||||
is_miscible_));
|
||||
|
||||
if (!BaseType::threshold_pressures_by_face_.empty()) {
|
||||
model->setThresholdPressures(BaseType::threshold_pressures_by_face_);
|
||||
|
||||
@ -23,6 +23,7 @@
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/Tables/PvdsTable.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
|
||||
|
||||
|
||||
namespace Opm
|
||||
@ -96,7 +97,7 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
|
||||
// relative permeabilty multiplier
|
||||
if (!ssfnTables.empty()) {
|
||||
|
||||
int numRegions = pvdsTables.size();
|
||||
int numRegions = ssfnTables.size();
|
||||
|
||||
if(numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only single table saturation function supported for SSFN");
|
||||
@ -119,7 +120,157 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "SSFN must be specified in SOLVENT runs\n");
|
||||
}
|
||||
|
||||
|
||||
if (deck->hasKeyword("MISCIBLE") ) {
|
||||
// misicible hydrocabon relative permeability wrt water
|
||||
const TableContainer& sof2Tables = tables->getSof2Tables();
|
||||
if (!sof2Tables.empty()) {
|
||||
|
||||
int numRegions = sof2Tables.size();
|
||||
|
||||
if(numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only single table saturation function supported for SOF2");
|
||||
}
|
||||
// resize the attributes of the object
|
||||
krn_.resize(numRegions);
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
const Opm::Sof2Table& sof2Table = sof2Tables.getTable<Sof2Table>(regionIdx);
|
||||
|
||||
// Copy data
|
||||
// Sn = So + Sg + Ss;
|
||||
const auto& sn = sof2Table.getSoColumn();
|
||||
const auto& krn = sof2Table.getKroColumn();
|
||||
|
||||
krn_[regionIdx] = NonuniformTableLinear<double>(sn, krn);
|
||||
}
|
||||
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "SOF2 must be specified in MISCIBLE (SOLVENT) runs\n");
|
||||
}
|
||||
|
||||
const TableContainer& miscTables = tables->getMiscTables();
|
||||
if (!miscTables.empty()) {
|
||||
|
||||
int numRegions = miscTables.size();
|
||||
|
||||
if(numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only single table miscibility function supported for MISC");
|
||||
}
|
||||
// resize the attributes of the object
|
||||
misc_.resize(numRegions);
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
const Opm::MiscTable& miscTable = miscTables.getTable<MiscTable>(regionIdx);
|
||||
|
||||
// Copy data
|
||||
// solventFraction = Ss / (Ss + Sg);
|
||||
const auto& solventFraction = miscTable.getSolventFractionColumn();
|
||||
const auto& misc = miscTable.getMiscibilityColumn();
|
||||
|
||||
misc_[regionIdx] = NonuniformTableLinear<double>(solventFraction, misc);
|
||||
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n");
|
||||
}
|
||||
|
||||
// miscible relative permeability multipleiers
|
||||
const TableContainer& msfnTables = tables->getMsfnTables();
|
||||
if (!msfnTables.empty()) {
|
||||
|
||||
int numRegions = msfnTables.size();
|
||||
|
||||
if(numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only single table saturation function supported for MSFN");
|
||||
}
|
||||
// resize the attributes of the object
|
||||
mkrsg_.resize(numRegions);
|
||||
mkro_.resize(numRegions);
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
const Opm::MsfnTable& msfnTable = msfnTables.getTable<MsfnTable>(regionIdx);
|
||||
|
||||
// Copy data
|
||||
// Ssg = Ss + Sg;
|
||||
const auto& Ssg = msfnTable.getGasPhaseFractionColumn();
|
||||
const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn();
|
||||
const auto& kro = msfnTable.getOilRelpermMultiplierColumn();
|
||||
|
||||
mkrsg_[regionIdx] = NonuniformTableLinear<double>(Ssg, krsg);
|
||||
mkro_[regionIdx] = NonuniformTableLinear<double>(Ssg, kro);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
const TableContainer& sorwmisTables = tables->getSorwmisTables();
|
||||
if (!sorwmisTables.empty()) {
|
||||
|
||||
int numRegions = sorwmisTables.size();
|
||||
|
||||
if(numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SORWMIS");
|
||||
}
|
||||
// resize the attributes of the object
|
||||
sorwmis_.resize(numRegions);
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
const Opm::SorwmisTable& sorwmisTable = sorwmisTables.getTable<SorwmisTable>(regionIdx);
|
||||
|
||||
// Copy data
|
||||
const auto& sw = sorwmisTable.getWaterSaturationColumn();
|
||||
const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn();
|
||||
|
||||
sorwmis_[regionIdx] = NonuniformTableLinear<double>(sw, sorwmis);
|
||||
}
|
||||
}
|
||||
|
||||
const TableContainer& sgcwmisTables = tables->getSgcwmisTables();
|
||||
if (!sgcwmisTables.empty()) {
|
||||
|
||||
int numRegions = sgcwmisTables.size();
|
||||
|
||||
if(numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SGCWMIS");
|
||||
}
|
||||
// resize the attributes of the object
|
||||
sgcwmis_.resize(numRegions);
|
||||
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
const Opm::SgcwmisTable& sgcwmisTable = sgcwmisTables.getTable<SgcwmisTable>(regionIdx);
|
||||
|
||||
// Copy data
|
||||
const auto& sw = sgcwmisTable.getWaterSaturationColumn();
|
||||
const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn();
|
||||
|
||||
sgcwmis_[regionIdx] = NonuniformTableLinear<double>(sw, sgcwmis);
|
||||
}
|
||||
}
|
||||
|
||||
if (deck->hasKeyword("TLMIXPAR")) {
|
||||
const auto tlmixparRecord = deck->getKeyword("TLMIXPAR")->getRecord(0);
|
||||
std::vector<double> mix_params_viscosity = tlmixparRecord->getItem("TL_VISCOSITY_PARAMETER")->getSIDoubleData();
|
||||
const int numRegions = mix_params_viscosity.size();
|
||||
if (numRegions > 1) {
|
||||
OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
|
||||
}
|
||||
mix_param_viscosity_ = mix_params_viscosity[0];
|
||||
|
||||
std::vector<double> mix_params_density = tlmixparRecord->getItem("TL_DENSITY_PARAMETER")->getSIDoubleData();
|
||||
const int numDensityItems = mix_params_density.size();
|
||||
if (numDensityItems == 0) {
|
||||
mix_param_density_ = mix_param_viscosity_;
|
||||
} else if (numDensityItems == 1) {
|
||||
mix_param_density_ = mix_params_density[0];
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Only singel miscibility region is supported for TLMIXPAR.");
|
||||
}
|
||||
} else {
|
||||
mix_param_viscosity_ = 0.0;
|
||||
mix_param_density_ = 0.0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
|
||||
@ -149,75 +300,104 @@ ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg,
|
||||
const Cells& cells) const
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int n = cells.size();
|
||||
assert(pg.size() == n);
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, b_);
|
||||
|
||||
V b(n);
|
||||
V dbdp(n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const double& pg_i = pg.value()[i];
|
||||
int regionIdx = cellPvtRegionIdx_[i];
|
||||
b[i] = b_[regionIdx](pg_i);
|
||||
dbdp[i] = b_[regionIdx].derivative(pg_i);
|
||||
}
|
||||
|
||||
ADB::M dbdp_diag(dbdp.matrix().asDiagonal());
|
||||
const int num_blocks = pg.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]);
|
||||
}
|
||||
return ADB::function(std::move(b), std::move(jacs));
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction,
|
||||
const Cells& cells) const
|
||||
const Cells& cells) const
|
||||
{
|
||||
const int n = cells.size();
|
||||
assert(solventFraction.value().size() == n);
|
||||
V krg(n);
|
||||
V dkrgdsf(n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const double& solventFraction_i = solventFraction.value()[i];
|
||||
int regionIdx = 0; // TODO add mapping from cells to sat function table
|
||||
krg[i] = krg_[regionIdx](solventFraction_i);
|
||||
dkrgdsf[i] = krg_[regionIdx].derivative(solventFraction_i);
|
||||
}
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krg_);
|
||||
|
||||
ADB::M dkrgdsf_diag(dkrgdsf.matrix().asDiagonal());
|
||||
const int num_blocks = solventFraction.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
fastSparseProduct(dkrgdsf_diag, solventFraction.derivative()[block], jacs[block]);
|
||||
}
|
||||
return ADB::function(std::move(krg), std::move(jacs));
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction,
|
||||
const Cells& cells) const
|
||||
const Cells& cells) const
|
||||
{
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, krs_);
|
||||
}
|
||||
|
||||
|
||||
ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn,
|
||||
const Cells& cells) const
|
||||
{
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, krn_);
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (mkrsg_.size() > 0) {
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, mkrsg_);
|
||||
}
|
||||
// trivial function if not specified
|
||||
return Ssg;
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So,
|
||||
const Cells& cells) const
|
||||
{
|
||||
if (mkro_.size() > 0) {
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(So, cells, mkro_);
|
||||
}
|
||||
// trivial function if not specified
|
||||
return So;
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction,
|
||||
const Cells& cells) const
|
||||
{
|
||||
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, misc_);
|
||||
}
|
||||
|
||||
|
||||
ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw,
|
||||
const Cells& cells) const {
|
||||
if (sgcwmis_.size()>0) {
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sgcwmis_);
|
||||
}
|
||||
// return zeros if not specified
|
||||
return ADB::constant(V::Zero(Sw.size()));
|
||||
}
|
||||
|
||||
|
||||
ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw,
|
||||
const Cells& cells) const {
|
||||
if (sorwmis_.size()>0) {
|
||||
return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, sorwmis_);
|
||||
}
|
||||
// return zeros if not specified
|
||||
return ADB::constant(V::Zero(Sw.size()));
|
||||
}
|
||||
|
||||
ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
|
||||
const Cells& cells,
|
||||
const std::vector<NonuniformTableLinear<double>>& tables) const {
|
||||
const int n = cells.size();
|
||||
assert(solventFraction.value().size() == n);
|
||||
V krs(n);
|
||||
V dkrsdsf(n);
|
||||
assert(X_AD.value().size() == n);
|
||||
V x(n);
|
||||
V dx(n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
const double& solventFraction_i = solventFraction.value()[i];
|
||||
const double& X_i = X_AD.value()[i];
|
||||
int regionIdx = 0; // TODO add mapping from cells to sat function table
|
||||
krs[i] = krs_[regionIdx](solventFraction_i);
|
||||
dkrsdsf[i] = krs_[regionIdx].derivative(solventFraction_i);
|
||||
x[i] = tables[regionIdx](X_i);
|
||||
dx[i] = tables[regionIdx].derivative(X_i);
|
||||
}
|
||||
|
||||
ADB::M dkrsdsf_diag(dkrsdsf.matrix().asDiagonal());
|
||||
const int num_blocks = solventFraction.numBlocks();
|
||||
ADB::M dx_diag(dx.matrix().asDiagonal());
|
||||
const int num_blocks = X_AD.numBlocks();
|
||||
std::vector<ADB::M> jacs(num_blocks);
|
||||
for (int block = 0; block < num_blocks; ++block) {
|
||||
fastSparseProduct(dkrsdsf_diag, solventFraction.derivative()[block], jacs[block]);
|
||||
fastSparseProduct(dx_diag, X_AD.derivative()[block], jacs[block]);
|
||||
}
|
||||
return ADB::function(std::move(krs), std::move(jacs));
|
||||
return ADB::function(std::move(x), std::move(jacs));
|
||||
}
|
||||
|
||||
|
||||
|
||||
V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
|
||||
const int n = cells.size();
|
||||
V density(n);
|
||||
@ -228,4 +408,13 @@ V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
|
||||
return density;
|
||||
}
|
||||
|
||||
double SolventPropsAdFromDeck::mixingParameterViscosity() const {
|
||||
return mix_param_viscosity_;
|
||||
}
|
||||
|
||||
double SolventPropsAdFromDeck::mixingParameterDensity() const {
|
||||
return mix_param_density_;
|
||||
}
|
||||
|
||||
|
||||
} //namespace OPM
|
||||
|
||||
@ -56,35 +56,95 @@ public:
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n formation volume factor values.
|
||||
ADB bSolvent(const ADB& pg,
|
||||
const Cells& cells) const;
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Solvent viscosity.
|
||||
/// \param[in] pg Array of n gas pressure values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
||||
/// \return Array of n viscosity values.
|
||||
ADB muSolvent(const ADB& pg,
|
||||
const Cells& cells) const;
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Gas relPerm multipliers
|
||||
/// \param[in] solventFraction Array of n solvent fraction values.
|
||||
/// \param[in] gasFraction Array of n gas fraction Sg / (sg + Ss) values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n gas relPerm multiplier values.
|
||||
ADB gasRelPermMultiplier(const ADB& solventFraction,
|
||||
const Cells& cells) const;
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Solvent relPerm multipliers
|
||||
/// \param[in] solventFraction Array of n solvent fraction values.
|
||||
/// \param[in] solventFraction Array of n solvent fraction Ss / (Sg + Ss) values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n solvent relPerm multiplier values.
|
||||
ADB solventRelPermMultiplier(const ADB& solventFraction,
|
||||
const Cells& cells) const;
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Miscible hydrocrabon relPerm wrt water
|
||||
/// \param[in] Sn Array of n total hyrdrocarbon saturation values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n miscible hyrdrocabon wrt water relPerm values.
|
||||
ADB misicibleHydrocarbonWaterRelPerm(const ADB& Sn,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Miscible Solvent + Gas relPerm multiplier
|
||||
/// \param[in] Ssg Array of n total gas fraction (Sgas + Ssolvent) / Sn values, where
|
||||
/// Sn = Sgas + Ssolvent + Soil.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n solvent gas relperm multiplier.
|
||||
ADB miscibleSolventGasRelPermMultiplier(const ADB& Ssg,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Miscible Oil relPerm multiplier
|
||||
/// \param[in] So Array of n oil fraction values. Soil / Sn values, where Sn = Sgas + Ssolvent + Soil.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n oil relperm multiplier.
|
||||
ADB miscibleOilRelPermMultiplier(const ADB& So,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Miscible function
|
||||
/// \param[in] solventFraction Array of n solvent fraction Ss / (Sg + Ss) values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n miscibility values
|
||||
ADB miscibilityFunction(const ADB& solventFraction,
|
||||
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.
|
||||
/// \return Array of n miscible critical gas saturation values
|
||||
ADB miscibleCriticalGasSaturationFunction(const ADB& Sw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Miscible residual oil 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.
|
||||
/// \return Array of n miscible residual oil saturation values
|
||||
ADB miscibleResidualOilSaturationFunction(const ADB& Sw,
|
||||
const Cells& cells) const;
|
||||
|
||||
/// Solvent surface density
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \return Array of n solvent density values.
|
||||
V solventSurfaceDensity(const Cells& cells) const;
|
||||
|
||||
/// Todd-Longstaff mixing parameter for viscosity calculation
|
||||
double mixingParameterViscosity() const;
|
||||
|
||||
/// Todd-Longstaff mixing parameter for density calculation
|
||||
double mixingParameterDensity() const;
|
||||
|
||||
|
||||
private:
|
||||
|
||||
/// Makes ADB from table values
|
||||
/// \param[in] X Array of n table lookup values.
|
||||
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
|
||||
/// \param[in] tables Vector of tables, one for each PVT region.
|
||||
/// \return Array of n solvent density values.
|
||||
ADB makeADBfromTables(const ADB& X,
|
||||
const Cells& cells,
|
||||
const std::vector<NonuniformTableLinear<double>>& tables) const;
|
||||
|
||||
// The PVT region which is to be used for each cell
|
||||
std::vector<int> cellPvtRegionIdx_;
|
||||
std::vector<NonuniformTableLinear<double> > b_;
|
||||
@ -93,6 +153,14 @@ private:
|
||||
std::vector<double> solvent_surface_densities_;
|
||||
std::vector<NonuniformTableLinear<double> > krg_;
|
||||
std::vector<NonuniformTableLinear<double> > krs_;
|
||||
std::vector<NonuniformTableLinear<double> > krn_;
|
||||
std::vector<NonuniformTableLinear<double> > mkro_;
|
||||
std::vector<NonuniformTableLinear<double> > mkrsg_;
|
||||
std::vector<NonuniformTableLinear<double> > misc_;
|
||||
std::vector<NonuniformTableLinear<double> > sorwmis_;
|
||||
std::vector<NonuniformTableLinear<double> > sgcwmis_;
|
||||
double mix_param_viscosity_;
|
||||
double mix_param_density_;
|
||||
};
|
||||
|
||||
} // namespace OPM
|
||||
|
||||
@ -238,6 +238,8 @@ namespace Opm {
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const ADB& kr ,
|
||||
const ADB& mu ,
|
||||
const ADB& rho ,
|
||||
const ADB& p ,
|
||||
const SolutionState& state );
|
||||
|
||||
|
||||
@ -300,7 +300,10 @@ namespace Opm {
|
||||
}
|
||||
|
||||
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
|
||||
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
const ADB mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
|
||||
const ADB rho = fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv);
|
||||
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
|
||||
|
||||
residual_.material_balance_eq[ phaseIdx ] =
|
||||
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
|
||||
@ -413,18 +416,18 @@ namespace Opm {
|
||||
BlackoilPolymerModel<Grid>::computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const ADB& kr ,
|
||||
const ADB& mu ,
|
||||
const ADB& rho ,
|
||||
const ADB& phasePressure,
|
||||
const SolutionState& state)
|
||||
{
|
||||
Base::computeMassFlux(actph, transi, kr, phasePressure, state);
|
||||
Base::computeMassFlux(actph, transi, kr, mu, rho, phasePressure, state);
|
||||
|
||||
// Polymer treatment.
|
||||
const int canonicalPhaseIdx = canph_[ actph ];
|
||||
if (canonicalPhaseIdx == Water) {
|
||||
if (has_polymer_) {
|
||||
const std::vector<PhasePresence>& cond = phaseCondition();
|
||||
const ADB tr_mult = transMult(state.pressure);
|
||||
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv, cond);
|
||||
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
|
||||
const ADB mc = computeMc(state);
|
||||
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user