diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh index 83cb541df..d2ed3b212 100644 --- a/ebos/equil/initstateequil.hh +++ b/ebos/equil/initstateequil.hh @@ -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()) { diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 72469cfbb..0e034b806 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -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, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 112f44e24..e7c7bfc06 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -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. diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 68e3dec66..2f73f4bde 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -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]; diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp index 3fd99aa77..167efcf09 100644 --- a/opm/simulators/wells/SingleWellState.hpp +++ b/opm/simulators/wells/SingleWellState.hpp @@ -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 well_potentials; std::vector productivity_index; std::vector surface_rates; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 9cefbfdca..561acc4b1 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -280,6 +280,7 @@ namespace Opm std::vector& 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& b_perfcells_dense, const double Tw, const int perf, @@ -306,6 +308,7 @@ namespace Opm std::vector& 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, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 362d13813..f317a0e63 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -94,12 +94,15 @@ namespace Opm std::vector& 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 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 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 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& b_perfcells_dense, const double Tw, const int perf, @@ -239,6 +249,7 @@ namespace Opm std::vector& 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(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 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(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(); diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 16d632e43..919cd5e8e 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -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;