adjustments to account for vaporized water

This commit is contained in:
Paul Egberts 2022-04-08 22:31:21 +02:00
parent 6f29bf715c
commit ab3be6dce4
8 changed files with 58 additions and 12 deletions

View File

@ -273,7 +273,7 @@ private:
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
}
else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv);
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=Rvw*/);
}
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
if (FluidSystem::enableVaporizedOil()) {

View File

@ -677,6 +677,7 @@ computeSegmentFluidProperties(const EvalWell& temperature,
}
EvalWell rv(0.0);
EvalWell rvw(0.0);
// gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
@ -692,9 +693,9 @@ computeSegmentFluidProperties(const EvalWell& temperature,
rv = rvmax;
}
b[gasCompIdx] =
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv, rvw);
visc[gasCompIdx] =
FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv);
FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv, rvw);
phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx]
+ rv * b[gasCompIdx] * surf_dens[oilCompIdx];
} else { // no oil exists
@ -1114,6 +1115,7 @@ getSegmentSurfaceVolume(const EvalWell& temperature,
}
EvalWell rv(0.0);
EvalWell rvw(0.0);
// gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
@ -1137,7 +1139,8 @@ getSegmentSurfaceVolume(const EvalWell& temperature,
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
rv);
rv,
rvw);
} else { // no oil exists
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,

View File

@ -1509,6 +1509,7 @@ namespace Opm
auto& ws = well_state.well(this->index_of_well_);
ws.dissolved_gas_rate = 0;
ws.vaporized_oil_rate = 0;
ws.vaporized_wat_rate = 0;
// for the black oil cases, there will be four equations,
// the first three of them are the mass balance equations, the last one is the pressure equations.

View File

@ -318,7 +318,7 @@ namespace Opm {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, 0.0 /*=Rvw*/);
const double den = bg * detR;
coeff[ig] += 1.0 / den;
@ -359,7 +359,7 @@ namespace Opm {
}
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, 0.0);
coeff[ig] += 1.0 / bg;
}
}
@ -447,7 +447,7 @@ namespace Opm {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, 0.0 /*=Rvw*/);
const double den = bg * detR;
voidage_rates[ig] = surface_rates[ig];

View File

@ -56,6 +56,7 @@ public:
double temperature{0};
double dissolved_gas_rate{0};
double vaporized_oil_rate{0};
double vaporized_wat_rate{0};
std::vector<double> well_potentials;
std::vector<double> productivity_index;
std::vector<double> surface_rates;

View File

@ -280,6 +280,7 @@ namespace Opm
std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
double& perf_vap_wat_rate,
DeferredLogger& deferred_logger) const;
void computePerfRateScalar(const IntensiveQuantities& intQuants,
@ -297,6 +298,7 @@ namespace Opm
const Value& bhp,
const Value& rs,
const Value& rv,
const Value& rvw,
std::vector<Value>& b_perfcells_dense,
const double Tw,
const int perf,
@ -306,6 +308,7 @@ namespace Opm
std::vector<Value>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
double& perf_vap_wat_rate,
DeferredLogger& deferred_logger) const;
void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,

View File

@ -94,12 +94,15 @@ namespace Opm
std::vector<EvalWell>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
double& perf_vap_wat_rate,
DeferredLogger& deferred_logger) const
{
const auto& fs = intQuants.fluidState();
const EvalWell pressure = this->extendEval(this->getPerfCellPressure(fs));
const EvalWell rs = this->extendEval(fs.Rs());
const EvalWell rv = this->extendEval(fs.Rv());
const EvalWell rvw = this->extendEval(fs.Rvw());
std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@ -139,6 +142,7 @@ namespace Opm
bhp,
rs,
rv,
rvw,
b_perfcells_dense,
Tw,
perf,
@ -148,6 +152,7 @@ namespace Opm
cq_s,
perf_dis_gas_rate,
perf_vap_oil_rate,
perf_vap_wat_rate,
deferred_logger);
}
@ -167,6 +172,7 @@ namespace Opm
const Scalar pressure = this->getPerfCellPressure(fs).value();
const Scalar rs = fs.Rs().value();
const Scalar rv = fs.Rv().value();
const Scalar rvw = fs.Rvw().value();
std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@ -197,6 +203,7 @@ namespace Opm
Scalar perf_dis_gas_rate = 0.0;
Scalar perf_vap_oil_rate = 0.0;
Scalar perf_vap_wat_rate = 0.0;
// surface volume fraction of fluids within wellbore
std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
@ -209,6 +216,7 @@ namespace Opm
bhp,
rs,
rv,
rvw,
b_perfcells_dense,
Tw,
perf,
@ -218,6 +226,7 @@ namespace Opm
cq_s,
perf_dis_gas_rate,
perf_vap_oil_rate,
perf_vap_wat_rate,
deferred_logger);
}
@ -230,6 +239,7 @@ namespace Opm
const Value& bhp,
const Value& rs,
const Value& rv,
const Value& rvw,
std::vector<Value>& b_perfcells_dense,
const double Tw,
const int perf,
@ -239,6 +249,7 @@ namespace Opm
std::vector<Value>& cq_s,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
double& perf_vap_wat_rate,
DeferredLogger& deferred_logger) const
{
// Pressure drawdown (also used to determine direction of flow)
@ -277,6 +288,12 @@ namespace Opm
perf_dis_gas_rate = getValue(dis_gas);
perf_vap_oil_rate = getValue(vap_oil);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const Value vap_wat = rvw * cq_sGas;
cq_s[waterCompIdx] += vap_wat;
if (this->isProducer())
perf_vap_wat_rate = getValue(vap_wat);
}
}
} else {
@ -353,10 +370,13 @@ namespace Opm
// TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
// s means standard condition, r means reservoir condition
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o
// q_gs = q_gr * b_g + rs * q_or * b_o
// q_ws = q_wr * b_w + rvw * q_gr * b_g
// d = 1.0 - rs * rv
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
// q_wr = 1/b_w * (rvw * q_gr * b_g -q_ws)= 1/b_w * (rvw * 1/d * (q_gs - rs * q_os) - q_ws) =
// 1 / (b_w * d) * (rvw * (q_gs - rs * q_os) - d * q_ws)
const double d = 1.0 - getValue(rv) * getValue(rs);
@ -378,6 +398,18 @@ namespace Opm
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
}
if (/*has_watVapor && */ FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// q_ws = q_wr * b_w + rvw * q_gr * b_g
// d = 1.0 - rs * rv
// q_wr = 1 / (b_w * d) * (rvw * (q_gs - rs * q_os) - d * q_ws)
const double d = 1.0 - getValue(rv) * getValue(rs);
// vaporized water in gas
// rvw * q_gr * b_g = q_ws -q_wr *b_w = q_ws - (rvw * (q_gs - rs * q_os) - d * q_ws)/d = rvw * (q_gs -rs *q_os) / d
perf_vap_wat_rate = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
}
}
}
}
@ -426,6 +458,7 @@ namespace Opm
ws.vaporized_oil_rate = 0;
ws.dissolved_gas_rate = 0;
ws.vaporized_wat_rate = 0;
const int np = this->number_of_phases_;
@ -490,6 +523,7 @@ namespace Opm
const auto& comm = this->parallel_well_info_.communication();
ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate);
ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate);
ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate);
}
// accumulate resWell_ and invDuneD_ in parallel to get effects of all perforations (might be distributed)
@ -548,10 +582,11 @@ namespace Opm
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
double perf_vap_wat_rate = 0.;
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
auto& ws = well_state.well(this->index_of_well_);
auto& perf_data = ws.perf_data;
@ -571,6 +606,7 @@ namespace Opm
if (this->isProducer()) {
ws.dissolved_gas_rate += perf_dis_gas_rate;
ws.vaporized_oil_rate += perf_vap_oil_rate;
ws.vaporized_wat_rate += perf_vap_wat_rate;
}
if constexpr (has_energy) {
@ -696,7 +732,7 @@ namespace Opm
if constexpr (has_brine) {
// TODO: the application of well efficiency factor has not been tested with an example yet
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
EvalWell cq_s_sm = cq_s[waterCompIdx];
EvalWell cq_s_sm = cq_s[waterCompIdx] - perf_vap_wat_rate;
if (this->isInjector()) {
cq_s_sm *= this->wsalt();
} else {
@ -1303,7 +1339,7 @@ namespace Opm
}
rv = std::min(rv, rvmax_perf[perf]);
b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, 0.0 /*Rvw*/);
}
else {
b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
@ -2017,10 +2053,11 @@ namespace Opm
std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
double perf_vap_wat_rate = 0.;
double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
// TODO: make area a member
const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
const auto& material_law_manager = ebos_simulator.problem().materialLawManager();

View File

@ -407,6 +407,7 @@ WellState::report(const int* globalCellIdxMap,
well.rates.set(rt::dissolved_gas, ws.dissolved_gas_rate);
well.rates.set(rt::vaporized_oil, ws.vaporized_oil_rate);
well.rates.set(rt::vaporized_water, ws.vaporized_wat_rate);
{
auto& curr = well.current_control;