Merge pull request #3869 from plgbrts/vapoilwat

adding three-phase simulator including water evaporation and salt precipitation
This commit is contained in:
Bård Skaflestad
2022-05-25 15:33:53 +02:00
committed by GitHub
15 changed files with 316 additions and 21 deletions

View File

@@ -388,7 +388,7 @@ add_dependencies(moduleVersion opmsimulators)
set(FLOW_MODELS blackoil brine energy extbo foam gasoil gaswater
oilwater oilwater_brine gaswater_brine oilwater_polymer
oilwater_polymer_injectivity micp polymer solvent
gasoil_energy brine_saltprecipitation gaswater_saltprec_vapwat)
gasoil_energy brine_saltprecipitation gaswater_saltprec_vapwat brine_precsalt_vapwat)
set(FLOW_VARIANT_MODELS brine_energy onephase onephase_energy)
set(FLOW_TGTS)

View File

@@ -279,7 +279,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

@@ -0,0 +1,24 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <flow/flow_ebos_brine_precsalt_vapwat.hpp>
int main(int argc, char** argv)
{
return Opm::flowEbosBrinePrecsaltVapwatMainStandalone(argc, argv);
}

View File

@@ -0,0 +1,85 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <flow/flow_ebos_brine_precsalt_vapwat.hpp>
#include <opm/material/common/ResetLocale.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/simulators/flow/Main.hpp>
namespace Opm {
namespace Properties {
namespace TTag {
struct EclFlowBrinePrecsaltVapwatProblem {
using InheritsFrom = std::tuple<EclFlowProblem>;
};
}
template<class TypeTag>
struct EnableBrine<TypeTag, TTag::EclFlowBrinePrecsaltVapwatProblem> {
static constexpr bool value = true;
};
template<class TypeTag>
struct EnableSaltPrecipitation<TypeTag, TTag::EclFlowBrinePrecsaltVapwatProblem> {
static constexpr bool value = true;
};
template<class TypeTag>
struct EnableEvaporation<TypeTag, TTag::EclFlowBrinePrecsaltVapwatProblem> {
static constexpr bool value = true;
};
}}
namespace Opm {
void flowEbosBrinePrecsaltVapwatSetDeck(double setupTime, std::shared_ptr<Deck> deck,
std::shared_ptr<EclipseState> eclState,
std::shared_ptr<Schedule> schedule,
std::shared_ptr<SummaryConfig> summaryConfig)
{
using TypeTag = Properties::TTag::EclFlowBrinePrecsaltVapwatProblem;
using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
Vanguard::setExternalSetupTime(setupTime);
Vanguard::setExternalDeck(std::move(deck));
Vanguard::setExternalEclState(std::move(eclState));
Vanguard::setExternalSchedule(std::move(schedule));
Vanguard::setExternalSummaryConfig(std::move(summaryConfig));
}
// ----------------- Main program -----------------
int flowEbosBrinePrecsaltVapwatMain(int argc, char** argv, bool outputCout, bool outputFiles)
{
// we always want to use the default locale, and thus spare us the trouble
// with incorrect locale settings.
resetLocale();
FlowMainEbos<Properties::TTag::EclFlowBrinePrecsaltVapwatProblem>
mainfunc {argc, argv, outputCout, outputFiles};
return mainfunc.execute();
}
int flowEbosBrinePrecsaltVapwatMainStandalone(int argc, char** argv)
{
using TypeTag = Properties::TTag::EclFlowBrinePrecsaltVapwatProblem;
auto mainObject = Opm::Main(argc, argv);
return mainObject.runStatic<TypeTag>();
}
}

View File

@@ -0,0 +1,42 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef FLOW_EBOS_BRINE_PRECSALT_VAPWAT_HPP
#define FLOW_EBOS_BRINE_PRECSALT_VAPWAT_HPP
#include <memory>
namespace Opm {
class Deck;
class EclipseState;
class Schedule;
class SummaryConfig;
void flowEbosBrinePrecsaltVapwatSetDeck(double setupTime, std::shared_ptr<Deck> deck,
std::shared_ptr<EclipseState> eclState,
std::shared_ptr<Schedule> schedule,
std::shared_ptr<SummaryConfig> summaryConfig);
//! \brief Main function used in flow binary.
int flowEbosBrinePrecsaltVapwatMain(int argc, char** argv, bool outputCout, bool outputFiles);
//! \brief Main function used in flow_brine binary.
int flowEbosBrinePrecsaltVapwatMainStandalone(int argc, char** argv);
}
#endif // FLOW_EBOS_BRINE_PRECSALT_VAPWAT_HPP

View File

@@ -35,6 +35,7 @@
#include <flow/flow_ebos_brine.hpp>
#include <flow/flow_ebos_brine_saltprecipitation.hpp>
#include <flow/flow_ebos_gaswater_saltprec_vapwat.hpp>
#include <flow/flow_ebos_brine_precsalt_vapwat.hpp>
#include <flow/flow_ebos_onephase.hpp>
#include <flow/flow_ebos_onephase_energy.hpp>
#include <flow/flow_ebos_oilwater_brine.hpp>
@@ -727,9 +728,17 @@ private:
}
}
else if (eclipseState_->getSimulationConfig().hasPRECSALT()) {
flowEbosBrineSaltPrecipitationSetDeck(
setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
return flowEbosBrineSaltPrecipitationMain(argc_, argv_, outputCout_, outputFiles_);
if (eclipseState_->getSimulationConfig().hasVAPWAT()) {
//case with water vaporization into gas phase and salt precipitation
flowEbosBrinePrecsaltVapwatSetDeck(
setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
return flowEbosBrinePrecsaltVapwatMain(argc_, argv_, outputCout_, outputFiles_);
}
else {
flowEbosBrineSaltPrecipitationSetDeck(
setupTime_, deck_, eclipseState_, schedule_, summaryConfig_);
return flowEbosBrineSaltPrecipitationMain(argc_, argv_, outputCout_, outputFiles_);
}
}
else {
flowEbosBrineSetDeck(

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

@@ -72,6 +72,10 @@ namespace Opm
if constexpr (Base::has_brine) {
OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet");
}
if constexpr (Base::has_watVapor) {
OPM_THROW(std::runtime_error, "water evaporation is not supported by multisegment well yet");
}
}
@@ -1529,6 +1533,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

@@ -257,6 +257,7 @@ namespace Opm
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& rvwmax_perf,
std::vector<double>& surf_dens_perf) const;
void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator,
@@ -264,6 +265,7 @@ namespace Opm
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& rvwmax_perf,
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger);
@@ -280,6 +282,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 +300,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 +310,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

@@ -831,6 +831,7 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& rvwmax_perf,
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger)
{
@@ -939,7 +940,7 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
// Compute volume ratio.
x = mix;
// Subtract dissolved gas from oil phase and vapporized oil from gas phase
// Subtract dissolved gas from oil phase and vapporized oil from gas phase and vaporized water from gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
@@ -971,7 +972,37 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d;
}
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
//matrix system: (mix[oilpos] = q_os, x[oilpos] = bo*q_or, etc...)
//┌ ┐ ┌ ┐ ┌ ┐
//│mix[oilpos] │ | 1 Rv 0 | |x[oilpos] |
//│mix[gaspos] │ = │ Rs 1 0 │ │x[gaspos] │
//│mix[waterpos]│ │ 0 Rvw 1 │ │x[waterpos │
//└ ┘ └ ┘ └ ┘
const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
double rvw = 0.0;
if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) {
rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]);
}
if (rvw > 0.0) {
// Subtract water in gas from water mixture
x[waterpos] = mix[waterpos] - x[gaspos] * rvw;
}
}
} else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
//no oil
const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
double rvw = 0.0;
if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) {
rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]);
}
if (rvw > 0.0) {
// Subtract water in gas from water mixture
x[waterpos] = mix[waterpos] - mix[gaspos] * rvw;
}
}
double volrat = 0.0;
for (int component = 0; component < num_comp; ++component) {
volrat += x[component] / b_perf[perf*num_comp+ component];

View File

@@ -144,6 +144,7 @@ protected:
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& rvwmax_perf,
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger);

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,13 @@ namespace Opm
perf_dis_gas_rate = getValue(dis_gas);
perf_vap_oil_rate = getValue(vap_oil);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
const Value vap_wat = rvw * cq_sGas;
cq_s[waterCompIdx] += vap_wat;
if (this->isProducer())
perf_vap_wat_rate = getValue(vap_wat);
}
}
} else {
@@ -353,7 +371,7 @@ 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
// 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)
@@ -377,6 +395,18 @@ namespace Opm
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
// q_ws = q_wr * b_w + rvw * q_gr * b_g
// q_wr = 1 / b_w * (q_ws - rvw * q_gr * b_g) = 1 / b_w * (q_ws - rvw * 1 / d (q_gs - rs * q_os))
// vaporized water in gas
// rvw * q_gr * b_g = q_ws -q_wr *b_w = rvw * (q_gs -rs *q_os) / d
perf_vap_wat_rate = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
}
}
else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
//no oil
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
perf_vap_wat_rate = getValue(rvw) * getValue(cq_s[gasCompIdx]);
}
}
}
@@ -426,6 +456,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 +521,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 +580,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 +604,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 +730,8 @@ 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];
// Correction salt rate; evaporated water does not contain salt
EvalWell cq_s_sm = cq_s[waterCompIdx] - perf_vap_wat_rate;
if (this->isInjector()) {
cq_s_sm *= this->wsalt();
} else {
@@ -1249,6 +1284,7 @@ namespace Opm
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& rvwmax_perf,
std::vector<double>& surf_dens_perf) const
{
const int nperf = this->number_of_perforations_;
@@ -1266,6 +1302,10 @@ namespace Opm
rsmax_perf.resize(nperf);
rvmax_perf.resize(nperf);
}
//rvw is only used if both water and gas is present
if (waterPresent && gasPresent) {
rvwmax_perf.resize(nperf);
}
// Compute the average pressure in each well block
const auto& perf_press = ws.perf_data.pressure;
@@ -1292,7 +1332,35 @@ namespace Opm
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const int gaspos = gasCompIdx + perf * this->num_components_;
if (oilPresent) {
if (oilPresent && waterPresent) {
const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers
rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
double rv = 0.0;
double rvw = 0.0;
if (oilrate > 0) {
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
if (gasrate > 0) {
rv = oilrate / gasrate;
}
rv = std::min(rv, rvmax_perf[perf]);
}
if (waterrate > 0) {
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
if (gasrate > 0) {
rvw = waterrate / gasrate;
}
rvw = std::min(rvw, rvwmax_perf[perf]);
}
if (rv > 0.0 || rvw > 0.0){
b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, rvw);
}
else {
b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
}
} else if (oilPresent) {
//no water
const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
if (oilrate > 0) {
@@ -1303,7 +1371,24 @@ 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);
}
} else if (waterPresent) {
//no oil
const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers
rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
if (waterrate > 0) {
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
double rvw = 0.0;
if (gasrate > 0) {
rvw = waterrate / gasrate;
}
rvw = std::min(rvw, rvwmax_perf[perf]);
b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, 0.0 /*Rv*/, rvw);
}
else {
b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
@@ -1472,6 +1557,7 @@ namespace Opm
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& rvwmax_perf,
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger)
{
@@ -1528,7 +1614,7 @@ namespace Opm
}
}
this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, deferred_logger);
this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
this->computeConnectionPressureDelta();
}
@@ -1550,9 +1636,10 @@ namespace Opm
std::vector<double> b_perf;
std::vector<double> rsmax_perf;
std::vector<double> rvmax_perf;
std::vector<double> rvwmax_perf;
std::vector<double> surf_dens_perf;
computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, deferred_logger);
computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf);
computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
}
@@ -2017,10 +2104,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;