only use d factor for gas wells

This commit is contained in:
Tor Harald Sandve 2023-11-14 12:45:25 +01:00
parent e33d27ca90
commit 0f12977b17
6 changed files with 47 additions and 40 deletions

View File

@ -188,7 +188,7 @@ namespace Opm
template<class Value>
void computePerfRate(const IntensiveQuantities& int_quants,
const std::vector<Value>& mob_perfcells,
const double Tw,
const std::vector<Scalar>& Tw,
const int seg,
const int perf,
const Value& segment_pressure,
@ -204,7 +204,7 @@ namespace Opm
const Value& rv,
const std::vector<Value>& b_perfcells,
const std::vector<Value>& mob_perfcells,
const double Tw,
const std::vector<Scalar>& Tw,
const int perf,
const Value& segment_pressure,
const Value& segment_density,

View File

@ -349,7 +349,7 @@ namespace Opm
getMobility(ebosSimulator, perf, mob, deferred_logger);
const double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
const Scalar seg_pressure = segment_pressure[seg];
std::vector<Scalar> cq_s(this->num_components_, 0.);
Scalar perf_press = 0.0;
@ -768,7 +768,7 @@ namespace Opm
const Value& rv,
const std::vector<Value>& b_perfcells,
const std::vector<Value>& mob_perfcells,
const double Tw,
const std::vector<Scalar>& Tw,
const int perf,
const Value& segment_pressure,
const Value& segment_density,
@ -804,7 +804,7 @@ namespace Opm
// compute component volumetric rates at standard conditions
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown);
cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
}
@ -828,9 +828,6 @@ namespace Opm
total_mob += mob_perfcells[comp_idx];
}
// injection perforations total volume rates
const Value cqt_i = - Tw * (total_mob * drawdown);
// compute volume ratio between connection and at standard conditions
Value volume_ratio = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
@ -871,9 +868,9 @@ namespace Opm
}
}
// injecting connections total volumerates at standard conditions
Value cqt_is = cqt_i / volume_ratio;
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
const Value factor = total_mob / volume_ratio;
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
cq_s[componentIdx] = - Tw[componentIdx] * cmix_s[componentIdx] * factor * drawdown;
}
} // end for injection perforations
@ -907,7 +904,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
computePerfRate(const IntensiveQuantities& int_quants,
const std::vector<Value>& mob_perfcells,
const double Tw,
const std::vector<Scalar>& Tw,
const int seg,
const int perf,
const Value& segment_pressure,
@ -1182,11 +1179,11 @@ namespace Opm
// the well index associated with the connection
const double trans_mult = ebos_simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
const std::vector<Scalar> tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> ipr_b_perf(this->ipr_b_.size());
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
@ -1682,7 +1679,7 @@ namespace Opm
getMobility(ebosSimulator, perf, mob, deferred_logger);
const double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(int_quants, cell_idx);
const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
const std::vector<Scalar> Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
std::vector<EvalWell> cq_s(this->num_components_, 0.0);
EvalWell perf_press;
PerforationRates perfRates;
@ -1995,7 +1992,7 @@ namespace Opm
getMobility(ebosSimulator, perf, mob, deferred_logger);
const double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(int_quants, cell_idx);
const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
const std::vector<Scalar> Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
std::vector<Scalar> cq_s(this->num_components_, 0.0);
Scalar perf_press = 0.0;
PerforationRates perf_rates;

View File

@ -286,7 +286,7 @@ namespace Opm
void computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<Value>& mob,
const Value& bhp,
const double Tw,
const std::vector<Scalar>& Tw,
const int perf,
const bool allow_cf,
std::vector<Value>& cq_s,
@ -302,7 +302,7 @@ namespace Opm
const Value& rvw,
const Value& rsw,
std::vector<Value>& b_perfcells_dense,
const double Tw,
const std::vector<Scalar>& Tw,
const int perf,
const bool allow_cf,
const Value& skin_pressure,

View File

@ -119,7 +119,7 @@ namespace Opm
computePerfRate(const IntensiveQuantities& intQuants,
const std::vector<Value>& mob,
const Value& bhp,
const double Tw,
const std::vector<Scalar>& Tw,
const int perf,
const bool allow_cf,
std::vector<Value>& cq_s,
@ -224,7 +224,7 @@ namespace Opm
const Value& rvw,
const Value& rsw,
std::vector<Value>& b_perfcells_dense,
const double Tw,
const std::vector<Scalar>& Tw,
const int perf,
const bool allow_cf,
const Value& skin_pressure,
@ -249,7 +249,7 @@ namespace Opm
// compute component volumetric rates at standard conditions
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
}
@ -270,9 +270,6 @@ namespace Opm
total_mob_dense += mob[componentIdx];
}
// injection perforations total volume rates
const Value cqt_i = - Tw * (total_mob_dense * drawdown);
// compute volume ratio between connection at standard conditions
Value volumeRatio = bhp * 0.0; // initialize it with the correct type
;
@ -306,9 +303,9 @@ namespace Opm
}
// injecting connections total volumerates at standard conditions
Value cqt_is = cqt_i / volumeRatio;
const Value factor = total_mob_dense / volumeRatio;
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
cq_s[componentIdx] = - Tw[componentIdx] * cmix_s[componentIdx] * factor * drawdown;
}
// calculating the perforation solution gas rate and solution oil rates
@ -496,7 +493,7 @@ namespace Opm
PerforationRates perf_rates;
double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_rates, deferred_logger);
@ -790,12 +787,13 @@ namespace Opm
}
// the well index associated with the connection
const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
double trans_mult = ebos_simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const std::vector<Scalar> tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> ipr_b_perf(this->ipr_b_.size());
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
@ -1366,7 +1364,7 @@ namespace Opm
getMobility(ebosSimulator, perf, mob, deferred_logger);
double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
std::vector<Scalar> cq_s(this->num_components_, 0.);
PerforationRates perf_rates;
@ -1669,7 +1667,8 @@ namespace Opm
std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
PerforationRates perf_rates;
double trans_mult = ebos_simulator.problem().template wellTransMultiplier<double>(int_quant, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);
computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
perf_rates, deferred_logger);
// TODO: make area a member
@ -2274,7 +2273,7 @@ namespace Opm
std::vector<Scalar> cq_s(this->num_components_, 0.);
double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
PerforationRates perf_rates;
computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
cq_s, perf_rates, deferred_logger);

View File

@ -349,7 +349,7 @@ public:
return 0;
}
double wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const;
std::vector<double> wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const;
void updateConnectionDFactor(const Simulator& simulator, SingleWellState& ws) const;

View File

@ -1324,28 +1324,35 @@ namespace Opm
}
template <typename TypeTag>
double
std::vector<double>
WellInterface<TypeTag>::
wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const {
std::vector<Scalar> wi(this->num_components_, this->well_index_[perf] * trans_mult);
const auto& wdfac = this->well_ecl_.getWDFAC();
if (!wdfac.useDFactor()) {
return this->well_index_[perf] * trans_mult;
return wi;
}
if constexpr (! Indices::gasEnabled) {
return this->well_index_[perf] * trans_mult;
return wi;
}
// closed connection are still closed
if (this->well_index_[perf] == 0)
return 0.0;
return std::vector<Scalar>(this->num_components_, 0.0);
// for gas wells we may want to add a Forchheimer term if the WDFAC or WDFACCOR keyword is used
const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
// viscosity is evaluated at connection pressure
const double connection_pressure = ws.perf_data.pressure[perf];
const double mu = FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, getValue(intQuants.fluidState().Rv()), getValue(intQuants.fluidState().Rvw()));
const auto& rv = getValue(intQuants.fluidState().Rv());
const double psat = FluidSystem::gasPvt().saturationPressure(this->pvtRegionIdx(), ws.temperature, rv);
const double mu = connection_pressure < psat ?
FluidSystem::gasPvt().saturatedViscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure) :
FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, rv, getValue(intQuants.fluidState().Rvw()));
const double phi = getValue(intQuants.porosity());
//double k = connection.Kh()/h * trans_mult;
double Kh = connection.Kh()* trans_mult;
@ -1357,7 +1364,11 @@ namespace Opm
double d = wdfac.useConnectionDFactor()? connection.dFactor() : wdfac.getDFactor(rho, mu, Ke, phi, rw, h);
const PhaseUsage& pu = this->phaseUsage();
double Q = std::abs(ws.perf_data.phase_rates[perf*pu.num_phases + pu.phase_pos[Gas]]);
return 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (Q/2 * d / scaling));
const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
wi[gas_comp_idx] = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (Q/2 * d / scaling));
// solvent???
return wi;
}
template <typename TypeTag>