diff --git a/CMakeLists.txt b/CMakeLists.txt index 2de078e61..bef3bc997 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -372,7 +372,8 @@ set(FLOW_MODELS blackoil brine energy extbo foam gasoil gaswater oilwater_polymer_injectivity micp polymer solvent gasoil_energy brine_saltprecipitation gaswater_saltprec_vapwat brine_precsalt_vapwat - blackoil_legacyassembly gasoildiffuse) + blackoil_legacyassembly gasoildiffuse gaswater_dissolution + gaswater_dissolution_diffuse) set(FLOW_VARIANT_MODELS brine_energy onephase onephase_energy) set(FLOW_TGTS) diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index f7fb08ab3..79256c39b 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -77,6 +77,8 @@ class EclEquilInitializer enum { enableBrine = getPropValue() }; enum { enableEvaporation = getPropValue() }; enum { enableSaltPrecipitation = getPropValue() }; + enum { has_disgas_in_water = getPropValue() }; + public: // NB: setting the enableEnergy argument to true enables storage of enthalpy and @@ -89,6 +91,7 @@ public: enableEvaporation, enableBrine, enableSaltPrecipitation, + has_disgas_in_water, Indices::numPhases>; diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc index efbbf2ef7..dca10d1f9 100644 --- a/ebos/eclgenericoutputblackoilmodule.cc +++ b/ebos/eclgenericoutputblackoilmodule.cc @@ -685,6 +685,7 @@ assignToSolution(data::Solution& sol) {"PCOG", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, pcog_}, {"PRES_OVB", UnitSystem::measure::pressure, data::TargetType::RESTART_SOLUTION, overburdenPressure_}, {"RS", UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_SOLUTION, rs_}, + {"RSW", UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_SOLUTION, rsw_}, {"RSSAT", UnitSystem::measure::gas_oil_ratio, data::TargetType::RESTART_AUXILIARY, gasDissolutionFactor_}, {"RV", UnitSystem::measure::oil_gas_ratio, data::TargetType::RESTART_SOLUTION, rv_}, {"RVSAT", UnitSystem::measure::oil_gas_ratio, data::TargetType::RESTART_AUXILIARY, oilVaporizationFactor_}, @@ -791,6 +792,8 @@ setRestart(const data::Solution& sol, temperature_[elemIdx] = sol.data("TEMP")[globalDofIndex]; if (!rs_.empty() && sol.has("RS")) rs_[elemIdx] = sol.data("RS")[globalDofIndex]; + if (!rsw_.empty() && sol.has("RSW")) + rsw_[elemIdx] = sol.data("RSW")[globalDofIndex]; if (!rv_.empty() && sol.has("RV")) rv_[elemIdx] = sol.data("RV")[globalDofIndex]; if (!rvw_.empty() && sol.has("RVW")) @@ -985,6 +988,10 @@ doAllocBuffers(unsigned bufferSize, rs_.resize(bufferSize, 0.0); rstKeywords["RS"] = 0; } + if (FluidSystem::enableDissolvedGasInWater()) { + rsw_.resize(bufferSize, 0.0); + rstKeywords["RSW"] = 0; + } if (FluidSystem::enableVaporizedOil()) { rv_.resize(bufferSize, 0.0); rstKeywords["RV"] = 0; diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh index a4c0b9bcd..6016a9ae6 100644 --- a/ebos/eclgenericoutputblackoilmodule.hh +++ b/ebos/eclgenericoutputblackoilmodule.hh @@ -414,6 +414,7 @@ protected: ScalarBuffer oilPressure_; ScalarBuffer temperature_; ScalarBuffer rs_; + ScalarBuffer rsw_; ScalarBuffer rv_; ScalarBuffer rvw_; ScalarBuffer overburdenPressure_; diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index ec28c91f1..7a3430f04 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -285,6 +285,10 @@ public: this->rs_[globalDofIdx] = getValue(fs.Rs()); Valgrind::CheckDefined(this->rs_[globalDofIdx]); } + if (!this->rsw_.empty()) { + this->rsw_[globalDofIdx] = getValue(fs.Rsw()); + Valgrind::CheckDefined(this->rsw_[globalDofIdx]); + } if (!this->rv_.empty()) { this->rv_[globalDofIdx] = getValue(fs.Rv()); @@ -498,6 +502,9 @@ public: if (!this->rs_.empty()) this->rs_[globalDofIdx] = fsInitial.Rs(); + if (!this->rsw_.empty()) + this->rsw_[globalDofIdx] = fsInitial.Rsw(); + if (!this->rvw_.empty()) this->rvw_[globalDofIdx] = fsInitial.Rvw(); @@ -823,6 +830,8 @@ public: fs.setTemperature(this->temperature_[elemIdx]); if (!this->rs_.empty()) fs.setRs(this->rs_[elemIdx]); + if (!this->rsw_.empty()) + fs.setRsw(this->rsw_[elemIdx]); if (!this->rv_.empty()) fs.setRv(this->rv_[elemIdx]); if (!this->rvw_.empty()) diff --git a/ebos/equil/initstateequil.cc b/ebos/equil/initstateequil.cc index 3e23a3f24..da93a3e44 100644 --- a/ebos/equil/initstateequil.cc +++ b/ebos/equil/initstateequil.cc @@ -280,7 +280,7 @@ density(const double depth, { // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); - double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, saltConcentration); + double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0 /*=Rsw*/, saltConcentration); rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); return rho; } @@ -1583,8 +1583,9 @@ applyNumericalAquifers_(const GridView& gridView, const auto search = num_aqu_cells.find(cartIx); if (search != num_aqu_cells.end()) { // numerical aquifer cells are filled with water initially - // for co2store the oilphase is used for brine - const auto watPos = co2store? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx; + // for co2store the oilphase may be used for the brine + bool co2store_oil_as_brine = co2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const auto watPos = co2store_oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx; if (FluidSystem::phaseIsActive(watPos)) { this->sat_[watPos][elemIdx] = 1.; } else { diff --git a/flow/flow_ebos_gaswater_dissolution.cpp b/flow/flow_ebos_gaswater_dissolution.cpp new file mode 100644 index 000000000..e8e83ae63 --- /dev/null +++ b/flow/flow_ebos_gaswater_dissolution.cpp @@ -0,0 +1,102 @@ +/* + 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 . +*/ +#include "config.h" + +// Define making clear that the simulator supports AMG +#define FLOW_SUPPORT_AMG 1 + +#include + +#include +#include + +#include +#include +#include + +#include +#include + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlowGasWaterDissolutionProblem { + using InheritsFrom = std::tuple; +}; +} + +template +struct Linearizer { using type = TpfaLinearizer; }; + +template +struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; + +template +struct EnableDiffusion { static constexpr bool value = false; }; + +template +struct EnableDisgasInWater { + static constexpr bool value = true; +}; + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef BlackOilTwoPhaseIndices(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*disabledCompIdx=*/FluidSystem::oilCompIdx, + getPropValue()> type; +}; +}} + +namespace Opm { + + +// ----------------- Main program ----------------- +int flowEbosGasWaterDissolutionMain(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 + mainfunc {argc, argv, outputCout, outputFiles} ; + return mainfunc.execute(); +} + +int flowEbosGasWaterDissolutionMainStandalone(int argc, char** argv) +{ + using TypeTag = Properties::TTag::EclFlowGasWaterDissolutionProblem; + auto mainObject = Opm::Main(argc, argv); + return mainObject.runStatic(); +} + +} diff --git a/flow/flow_ebos_gaswater_dissolution.hpp b/flow/flow_ebos_gaswater_dissolution.hpp new file mode 100644 index 000000000..3bdd8f21f --- /dev/null +++ b/flow/flow_ebos_gaswater_dissolution.hpp @@ -0,0 +1,30 @@ +/* + 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 . +*/ +#ifndef FLOW_EBOS_GASWATER_DISSOLUTION_HPP +#define FLOW_EBOS_GASWATER_DISSOLUTION_HPP + +namespace Opm { + +//! \brief Main function used in flow binary. +int flowEbosGasWaterDissolutionMain(int argc, char** argv, bool outputCout, bool outputFiles); + +//! \brief Main function used in flow_gaswater_dissolution binary. +int flowEbosGasWaterDissolutionMainStandalone(int argc, char** argv); + +} + +#endif // FLOW_EBOS_GASWATER_DISSOLUTION_HPP diff --git a/flow/flow_ebos_gaswater_dissolution_diffuse.cpp b/flow/flow_ebos_gaswater_dissolution_diffuse.cpp new file mode 100644 index 000000000..95c8e4a9c --- /dev/null +++ b/flow/flow_ebos_gaswater_dissolution_diffuse.cpp @@ -0,0 +1,102 @@ +/* + 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 . +*/ +#include "config.h" + +// Define making clear that the simulator supports AMG +#define FLOW_SUPPORT_AMG 1 + +#include + +#include +#include + +#include +#include +#include + +#include +#include + +namespace Opm { +namespace Properties { +namespace TTag { +struct EclFlowGasWaterDissolutionDiffuseProblem { + using InheritsFrom = std::tuple; +}; +} + +//template +//struct Linearizer { using type = TpfaLinearizer; }; + +//template +//struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; + +template +struct EnableDiffusion { static constexpr bool value = true; }; + +template +struct EnableDisgasInWater { + static constexpr bool value = true; +}; + +//! The indices required by the model +template +struct Indices +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + using BaseTypeTag = TTag::EclFlowProblem; + using FluidSystem = GetPropType; + +public: + typedef BlackOilTwoPhaseIndices(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + getPropValue(), + /*PVOffset=*/0, + /*disabledCompIdx=*/FluidSystem::oilCompIdx, + getPropValue()> type; +}; +}} + +namespace Opm { + + +// ----------------- Main program ----------------- +int flowEbosGasWaterDissolutionDiffuseMain(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 + mainfunc {argc, argv, outputCout, outputFiles} ; + return mainfunc.execute(); +} + +int flowEbosGasWaterDissolutionDiffuseMainStandalone(int argc, char** argv) +{ + using TypeTag = Properties::TTag::EclFlowGasWaterDissolutionDiffuseProblem; + auto mainObject = Opm::Main(argc, argv); + return mainObject.runStatic(); +} + +} diff --git a/flow/flow_ebos_gaswater_dissolution_diffuse.hpp b/flow/flow_ebos_gaswater_dissolution_diffuse.hpp new file mode 100644 index 000000000..120530c05 --- /dev/null +++ b/flow/flow_ebos_gaswater_dissolution_diffuse.hpp @@ -0,0 +1,30 @@ +/* + 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 . +*/ +#ifndef FLOW_EBOS_GASWATER_DISSOLUTION_DIFFUSE_HPP +#define FLOW_EBOS_GASWATER_DISSOLUTION_DIFFUSE_HPP + +namespace Opm { + +//! \brief Main function used in flow binary. +int flowEbosGasWaterDissolutionDiffuseMain(int argc, char** argv, bool outputCout, bool outputFiles); + +//! \brief Main function used in flow_gaswater_dissolution binary. +int flowEbosGasWaterDissolutionDiffuseMainStandalone(int argc, char** argv); + +} + +#endif // FLOW_EBOS_GASWATER_DISSOLUTION_DIFFUSE_HPP diff --git a/flow/flow_gaswater_dissolution.cpp b/flow/flow_gaswater_dissolution.cpp new file mode 100644 index 000000000..43f11a1d3 --- /dev/null +++ b/flow/flow_gaswater_dissolution.cpp @@ -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 . +*/ +#include "config.h" +#include + + +int main(int argc, char** argv) +{ + return Opm::flowEbosGasWaterDissolutionMainStandalone(argc, argv); +} diff --git a/flow/flow_gaswater_dissolution_diffuse.cpp b/flow/flow_gaswater_dissolution_diffuse.cpp new file mode 100644 index 000000000..1ba4a5967 --- /dev/null +++ b/flow/flow_gaswater_dissolution_diffuse.cpp @@ -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 . +*/ +#include "config.h" +#include + + +int main(int argc, char** argv) +{ + return Opm::flowEbosGasWaterDissolutionDiffuseMainStandalone(argc, argv); +} diff --git a/opm/simulators/aquifers/AquiferAnalytical.hpp b/opm/simulators/aquifers/AquiferAnalytical.hpp index fc048b91b..3cc9894fa 100644 --- a/opm/simulators/aquifers/AquiferAnalytical.hpp +++ b/opm/simulators/aquifers/AquiferAnalytical.hpp @@ -62,6 +62,8 @@ public: enum { enableEnergy = getPropValue() }; enum { enableBrine = getPropValue() }; enum { enableEvaporation = getPropValue() }; + enum { has_disgas_in_water = getPropValue() }; + enum { enableSaltPrecipitation = getPropValue() }; static constexpr int numEq = BlackoilIndices::numEq; @@ -77,6 +79,7 @@ public: enableEvaporation, enableBrine, enableSaltPrecipitation, + has_disgas_in_water, BlackoilIndices::numPhases>; // Constructor diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 132fbeccd..f2a5b6b8e 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -79,7 +79,8 @@ protected: int phaseIdx_() const { - if (co2store_()) + // If OIL is used to model brine the aquifer should do the same + if (co2store_() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) return FluidSystem::oilPhaseIdx; return FluidSystem::waterPhaseIdx; diff --git a/opm/simulators/flow/Main.hpp b/opm/simulators/flow/Main.hpp index 3d6544183..b351342c2 100644 --- a/opm/simulators/flow/Main.hpp +++ b/opm/simulators/flow/Main.hpp @@ -42,6 +42,8 @@ #include #include #include +#include +#include #include #include #include @@ -632,6 +634,7 @@ private: int runTwoPhase(const Phases& phases) { const bool diffusive = eclipseState_->getSimulationConfig().isDiffusive(); + const bool disgasw = eclipseState_->getSimulationConfig().hasDISGASW(); // oil-gas if (phases.active( Phase::OIL ) && phases.active( Phase::GAS )) { @@ -655,12 +658,19 @@ private: // gas-water else if ( phases.active( Phase::GAS ) && phases.active( Phase::WATER ) ) { + if (disgasw) { + if (diffusive) { + return flowEbosGasWaterDissolutionDiffuseMain(argc_, argv_, outputCout_, outputFiles_); + } + return flowEbosGasWaterDissolutionMain(argc_, argv_, outputCout_, outputFiles_); + } if (diffusive) { if (outputCout_) { - std::cerr << "The DIFFUSE option is not available for the two-phase gas/water model." << std::endl; + std::cerr << "The DIFFUSE option is not available for the two-phase gas/water model without disgasw." << std::endl; } return EXIT_FAILURE; } + return flowEbosGasWaterMain(argc_, argv_, outputCout_, outputFiles_); } else { diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 98f8ed7e5..981a19115 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -305,6 +305,8 @@ copyToWellState(const MultisegmentWellGeneric& mswell, (pvtReg, segment_pressure[seg], std::max(0.0, Rs), std::max(0.0, Rv), + 0.0, // Rsw + 0.0, // Rvw temperature, saltConc, surf_rates, resv_rates); } @@ -340,7 +342,7 @@ copyToWellState(const MultisegmentWellGeneric& mswell, if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Water]] = - FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], saltConc); + FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], 0.0 /*Rsw*/, saltConc); } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 6c8e15b0a..1bf43d9d0 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -143,11 +143,12 @@ computeFluidProperties(const EvalWell& temperature, const EvalWell seg_pressure = primary_variables.getSegmentPressure(seg); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + EvalWell rsw(0.0); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); b[waterCompIdx] = - FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration); + FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rsw, saltConcentration); visc[waterCompIdx] = - FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration); + FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, rsw, saltConcentration); // TODO: double check here // TODO: should not we use phaseIndex here? phase_densities[waterCompIdx] = b[waterCompIdx] * surf_dens[waterCompIdx]; @@ -348,10 +349,12 @@ getSurfaceVolume(const EvalWell& temperature, std::vector b(well_.numComponents(), 0.); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + EvalWell rsw(0.0); b[waterCompIdx] = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, + rsw, saltConcentration); } diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 19b30d8dd..47c3cde2b 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -84,6 +84,10 @@ namespace Opm OPM_THROW(std::runtime_error, "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet." << " \n See (WCONINJE item 10 / WCONHIST item 8)"); } + if constexpr (!Indices::oilEnabled && Indices::numPhases > 1) { + OPM_THROW(std::runtime_error, "water + gas case not supported by multisegment well yet"); + } + } @@ -1546,6 +1550,7 @@ namespace Opm auto& ws = well_state.well(this->index_of_well_); ws.dissolved_gas_rate = 0; + ws.dissolved_gas_rate_in_water = 0; ws.vaporized_oil_rate = 0; ws.vaporized_wat_rate = 0; diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index a0ca88f27..28cbf2225 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -104,6 +104,8 @@ namespace Opm { ra.rv = 0.0; ra.pv = 0.0; ra.saltConcentration = 0.0; + ra.rsw = 0.0; + ra.rvw = 0.0; } @@ -162,6 +164,12 @@ namespace Opm { attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV; } attr.saltConcentration += fs.saltConcentration().value() * hydrocarbonPV; + if (FluidSystem::enableDissolvedGasInWater()) { + attr.rsw += fs.Rsw().value() * hydrocarbonPV; // scale with total volume? + } + if (FluidSystem::enableVaporizedWater()) { + attr.rvw += fs.Rvw().value() * hydrocarbonPV; // scale with total volume? + } } if (pv_cell > 0.) { @@ -183,6 +191,12 @@ namespace Opm { attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell; } attr.saltConcentration += fs.saltConcentration().value() * pv_cell; + if (FluidSystem::enableDissolvedGasInWater()) { + attr.rsw += fs.Rsw().value() * pv_cell; + } + if (FluidSystem::enableVaporizedWater()) { + attr.rvw += fs.Rvw().value() * pv_cell; + } } } @@ -199,11 +213,15 @@ namespace Opm { const double rs_hpv_sum = comm.sum(attri_hpv.rs); const double rv_hpv_sum = comm.sum(attri_hpv.rv); const double sc_hpv_sum = comm.sum(attri_hpv.saltConcentration); + const double rsw_hpv_sum = comm.sum(attri_hpv.rsw); + const double rvw_hpv_sum = comm.sum(attri_hpv.rvw); ra.pressure = p_hpv_sum / hpv_sum; ra.temperature = T_hpv_sum / hpv_sum; ra.rs = rs_hpv_sum / hpv_sum; ra.rv = rv_hpv_sum / hpv_sum; + ra.rsw = rsw_hpv_sum / hpv_sum; + ra.rvw = rvw_hpv_sum / hpv_sum; ra.pv = hpv_sum; ra.saltConcentration = sc_hpv_sum / hpv_sum; } else { @@ -215,12 +233,16 @@ namespace Opm { const double T_pv_sum = comm.sum(attri_pv.temperature); const double rs_pv_sum = comm.sum(attri_pv.rs); const double rv_pv_sum = comm.sum(attri_pv.rv); + const double rsw_pv_sum = comm.sum(attri_pv.rsw); + const double rvw_pv_sum = comm.sum(attri_pv.rvw); const double sc_pv_sum = comm.sum(attri_pv.saltConcentration); ra.pressure = p_pv_sum / pv_sum; ra.temperature = T_pv_sum / pv_sum; ra.rs = rs_pv_sum / pv_sum; ra.rv = rv_pv_sum / pv_sum; + ra.rsw = rsw_pv_sum / pv_sum; + ra.rvw = rvw_pv_sum / pv_sum; ra.pv = pv_sum; ra.saltConcentration = sc_pv_sum / pv_sum; } @@ -277,12 +299,24 @@ namespace Opm { std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); + // Actual Rsw and Rvw: + double Rsw = ra.rsw; + double Rvw = ra.rvw; + // Determinant of 'R' matrix + const double detRw = 1.0 - (Rsw * Rvw); + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - // q[w]_r = q[w]_s / bw + // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s) - const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration); + const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rsw, saltConcentration); - coeff[iw] = 1.0 / bw; + const double den = bw * detRw; + + coeff[iw] += 1.0 / den; + + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + coeff[ig] -= ra.rvw / den; + } } // Actual Rs and Rv: @@ -292,6 +326,13 @@ namespace Opm { // Determinant of 'R' matrix const double detR = 1.0 - (Rs * Rv); + // Currently we only support either gas in water or gas in oil + // not both + if(detR != 1 && detRw != 1 ) { + std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR); + throw(msg); + } + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) @@ -307,14 +348,19 @@ 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, 0.0 /*=Rvw*/); - const double den = bg * detR; - - coeff[ig] += 1.0 / den; - - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { - coeff[io] -= ra.rs / den; + const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, Rvw); + if (FluidSystem::enableDissolvedGasInWater()) { + const double denw = bg * detRw; + coeff[ig] += 1.0 / denw; + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + coeff[iw] -= ra.rsw / denw; + } + } else { + const double den = bg * detR; + coeff[ig] += 1.0 / den; + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + coeff[io] -= ra.rs / den; + } } } } @@ -338,7 +384,7 @@ namespace Opm { if (RegionAttributeHelpers::PhaseUsed::water(pu)) { // q[w]_r = q[w]_s / bw - const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration); + const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, saltConcentration); coeff[iw] = 1.0 / bw; } @@ -385,7 +431,8 @@ namespace Opm { const auto& ra = this->attr_.attributes(r); this->calcReservoirVoidageRates(pvtRegionIdx, - ra.pressure, ra.rs, ra.rv, + ra.pressure, ra.rs, ra.rv, + ra.rsw, ra.rvw, ra.temperature, ra.saltConcentration, surface_rates, @@ -409,6 +456,10 @@ namespace Opm { * \param[in] rs Dissolved gas/oil ratio. * * \param[in] rv Vaporised oil/gas ratio. + * + * \param[in] rsw Dissolved gas/water ratio. + * + * \param[in] rwv Vaporised water/gas ratio. * * \param[in] T Temperature. Unused in non-thermal simulation * runs. @@ -427,6 +478,8 @@ namespace Opm { const double p, const double rs, const double rv, + const double rsw, + const double rvw, const double T, const double saltConcentration, const SurfaceRates& surface_rates, @@ -440,15 +493,29 @@ namespace Opm { const auto [Rs, Rv] = this-> dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates); + const auto [Rsw, Rvw] = this-> + dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates); + + std::fill_n(&voidage_rates[0], pu.num_phases, 0.0); + + // Determinant of 'R' matrix + const auto detRw = 1.0 - (Rsw * Rvw); + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - // q[w]_r = q[w]_s / bw + // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s) + voidage_rates[iw] = surface_rates[iw]; + const auto bw = FluidSystem::waterPvt() .inverseFormationVolumeFactor(pvtRegionIdx, T, p, + Rsw, saltConcentration); - voidage_rates[iw] = surface_rates[iw] / bw; + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + voidage_rates[iw] -= Rvw * surface_rates[ig]; + } + voidage_rates[iw] /= bw * detRw; } // Determinant of 'R' matrix @@ -467,18 +534,32 @@ namespace Opm { voidage_rates[io] /= bo * detR; } + // we only support either gas in water + // or gas in oil + if(detR != 1 && detRw != 1 ) { + std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR); + throw(msg); + } if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) voidage_rates[ig] = surface_rates[ig]; if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { voidage_rates[ig] -= Rs * surface_rates[io]; } + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + voidage_rates[ig] -= Rsw * surface_rates[iw]; + } const auto bg = FluidSystem::gasPvt() .inverseFormationVolumeFactor(pvtRegionIdx, T, p, - Rv, 0.0 /*=Rvw*/); + Rv, Rvw); - voidage_rates[ig] /= bg * detR; + // we only support either gas in water or gas in oil + if (detRw == 1) { + voidage_rates[ig] /= bg * detR; + } else { + voidage_rates[ig] /= bg * detRw; + } } } @@ -537,6 +618,8 @@ namespace Opm { , temperature(0.0) , rs(0.0) , rv(0.0) + , rsw(0.0) + , rvw(0.0) , pv(0.0) , saltConcentration(0.0) {} @@ -545,6 +628,8 @@ namespace Opm { double temperature; double rs; double rv; + double rsw; + double rvw; double pv; double saltConcentration; }; diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp index a6ac2c02f..9caccb3d8 100644 --- a/opm/simulators/wells/SingleWellState.hpp +++ b/opm/simulators/wells/SingleWellState.hpp @@ -55,6 +55,7 @@ public: double thp{0}; double temperature{0}; double dissolved_gas_rate{0}; + double dissolved_gas_rate_in_water{0}; double vaporized_oil_rate{0}; double vaporized_wat_rate{0}; std::vector well_potentials; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 4e59ebdcb..97b9b1557 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -266,6 +266,7 @@ namespace Opm std::vector& rsmax_perf, std::vector& rvmax_perf, std::vector& rvwmax_perf, + std::vector& rswmax_perf, std::vector& surf_dens_perf) const; void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator, @@ -274,6 +275,7 @@ namespace Opm const std::vector& rsmax_perf, const std::vector& rvmax_perf, const std::vector& rvwmax_perf, + const std::vector& rswmax_perf, const std::vector& surf_dens_perf, DeferredLogger& deferred_logger); @@ -289,6 +291,7 @@ namespace Opm const bool allow_cf, std::vector& cq_s, double& perf_dis_gas_rate, + double& perf_dis_gas_rate_in_water, double& perf_vap_oil_rate, double& perf_vap_wat_rate, DeferredLogger& deferred_logger) const; @@ -309,6 +312,7 @@ namespace Opm const Value& rs, const Value& rv, const Value& rvw, + const Value& rsw, std::vector& b_perfcells_dense, const double Tw, const int perf, @@ -317,6 +321,7 @@ namespace Opm const std::vector& cmix_s, std::vector& cq_s, double& perf_dis_gas_rate, + double& perf_dis_gas_rate_in_water, double& perf_vap_oil_rate, double& perf_vap_wat_rate, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index 50486e687..1bee8a869 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -92,6 +92,7 @@ computeDensities(const std::vector& perfComponentRates, const std::vector& rsmax_perf, const std::vector& rvmax_perf, const std::vector& rvwmax_perf, + const std::vector& rswmax_perf, const std::vector& surf_dens_perf, DeferredLogger& deferred_logger) { @@ -216,7 +217,7 @@ computeDensities(const std::vector& perfComponentRates, if (d <= 0.0) { std::ostringstream sstr; sstr << "Problematic d value " << d << " obtained for well " << well_.name() - << " during ccomputeConnectionDensities with rs " << rs + << " during computeConnectionDensities with rs " << rs << ", rv " << rv << " obtaining d " << d << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) " @@ -232,34 +233,44 @@ computeDensities(const std::vector& 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); - Scalar 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); + } + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + //matrix system: (mix[oilpos] = q_os, x[oilpos] = bo*q_or, etc...) + //┌ ┐ ┌ ┐ ┌ ┐ + //│mix[oilpos] │ | 1 Rv 0 | |x[oilpos] | + //│mix[gaspos] │ = │ Rs 1 Rsw│ │x[gaspos] │ + //│mix[waterpos]│ │ 0 Rvw 1 │ │x[waterpos │ + //└ ┘ └ ┘ └ ┘ const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); Scalar 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; + Scalar rsw = 0.0; + if (!rswmax_perf.empty() && mix[waterpos] > 1e-12) { + rsw = std::min(mix[gaspos]/mix[waterpos], rswmax_perf[perf]); + } + const Scalar d = 1.0 - rsw*rvw; + if (d <= 0.0) { + std::ostringstream sstr; + sstr << "Problematic d value " << d << " obtained for well " << well_.name() + << " during computeConnectionDensities with rsw " << rsw + << ", rvw " << rvw + << " obtaining d " << d + << " Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) " + << " for this connection."; + deferred_logger.debug(sstr.str()); + } else { + if (rsw > 0.0) { + // Subtract gas in water from gas mixture + x[gaspos] = (mix[gaspos] - mix[waterpos]*rsw)/d; + } + if (rvw > 0.0) { + // Subtract water in gas from water mixture + x[waterpos] = (mix[waterpos] - mix[gaspos]*rvw)/d; + } } } @@ -288,6 +299,7 @@ computePropertiesForPressures(const WellState& well_state, std::vector& rsmax_perf, std::vector& rvmax_perf, std::vector& rvwmax_perf, + std::vector& rswmax_perf, std::vector& surf_dens_perf) const { const int nperf = well_.numPerfs(); @@ -312,6 +324,7 @@ computePropertiesForPressures(const WellState& well_state, //rvw is only used if both water and gas is present if (waterPresent && gasPresent) { rvwmax_perf.resize(nperf); + rswmax_perf.resize(nperf); } // Compute the average pressure in each well block @@ -330,8 +343,21 @@ computePropertiesForPressures(const WellState& well_state, if (waterPresent) { const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - b_perf[ waterCompIdx + perf * well_.numComponents()] = - FluidSystem::waterPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, saltConcentration); + double rsw = 0.0; + if (FluidSystem::enableDissolvedGasInWater()) { + // TODO support mutual solubility in water and oil + assert(!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)); + const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); + rswmax_perf[perf] = FluidSystem::waterPvt().saturatedGasDissolutionFactor(region_idx, temperature, p_avg, saltConcentration); + if (waterrate > 0) { + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0); + if (gasrate > 0) { + rsw = waterrate / gasrate; + } + rsw = std::min(rsw, rswmax_perf[perf]); + } + } + b_perf[ waterCompIdx + perf * well_.numComponents()] = FluidSystem::waterPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rsw, saltConcentration); } if (gasPresent) { @@ -456,6 +482,7 @@ computeProperties(const WellState& well_state, const std::vector& rsmax_perf, const std::vector& rvmax_perf, const std::vector& rvwmax_perf, + const std::vector& rswmax_perf, const std::vector& surf_dens_perf, DeferredLogger& deferred_logger) { @@ -521,7 +548,7 @@ computeProperties(const WellState& well_state, } } - this->computeDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger); + this->computeDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf, deferred_logger); this->computePressureDelta(); } diff --git a/opm/simulators/wells/StandardWellConnections.hpp b/opm/simulators/wells/StandardWellConnections.hpp index b29638dbd..127637f01 100644 --- a/opm/simulators/wells/StandardWellConnections.hpp +++ b/opm/simulators/wells/StandardWellConnections.hpp @@ -49,6 +49,7 @@ public: std::vector& rsmax_perf, std::vector& rvmax_perf, std::vector& rvwmax_perf, + std::vector& rswmax_perf, std::vector& surf_dens_perf) const; //! \brief Compute connection properties (densities, pressure drop, ...) @@ -61,6 +62,7 @@ public: const std::vector& rsmax_perf, const std::vector& rvmax_perf, const std::vector& rvwmax_perf, + const std::vector& rswmax_perf, const std::vector& surf_dens_perf, DeferredLogger& deferred_logger); @@ -84,6 +86,7 @@ private: const std::vector& rsmax_perf, const std::vector& rvmax_perf, const std::vector& rvwmax_perf, + const std::vector& rswmax_perf, const std::vector& surf_dens_perf, DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 4afb4374e..11bfb4eda 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -97,6 +97,7 @@ namespace Opm const bool allow_cf, std::vector& cq_s, double& perf_dis_gas_rate, + double& perf_dis_gas_rate_in_water, double& perf_vap_oil_rate, double& perf_vap_wat_rate, DeferredLogger& deferred_logger) const @@ -106,6 +107,8 @@ namespace Opm const EvalWell rs = this->extendEval(fs.Rs()); const EvalWell rv = this->extendEval(fs.Rv()); const EvalWell rvw = this->extendEval(fs.Rvw()); + const EvalWell rsw = this->extendEval(fs.Rsw()); + std::vector b_perfcells_dense(this->num_components_, EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0}); for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { @@ -147,6 +150,7 @@ namespace Opm rs, rv, rvw, + rsw, b_perfcells_dense, Tw, perf, @@ -155,6 +159,7 @@ namespace Opm cmix_s, cq_s, perf_dis_gas_rate, + perf_dis_gas_rate_in_water, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger); @@ -177,6 +182,7 @@ namespace Opm const Scalar rs = fs.Rs().value(); const Scalar rv = fs.Rv().value(); const Scalar rvw = fs.Rvw().value(); + const Scalar rsw = fs.Rsw().value(); std::vector b_perfcells_dense(this->num_components_, 0.0); for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { @@ -208,6 +214,7 @@ namespace Opm Scalar perf_dis_gas_rate = 0.0; Scalar perf_vap_oil_rate = 0.0; Scalar perf_vap_wat_rate = 0.0; + Scalar perf_dis_gas_rate_in_water = 0.0; // surface volume fraction of fluids within wellbore std::vector cmix_s(this->numComponents(), 0.0); @@ -221,6 +228,7 @@ namespace Opm rs, rv, rvw, + rsw, b_perfcells_dense, Tw, perf, @@ -229,6 +237,7 @@ namespace Opm cmix_s, cq_s, perf_dis_gas_rate, + perf_dis_gas_rate_in_water, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger); @@ -244,6 +253,7 @@ namespace Opm const Value& rs, const Value& rv, const Value& rvw, + const Value& rsw, std::vector& b_perfcells_dense, const double Tw, const int perf, @@ -252,6 +262,7 @@ namespace Opm const std::vector& cmix_s, std::vector& cq_s, double& perf_dis_gas_rate, + double& perf_dis_gas_rate_in_water, double& perf_vap_oil_rate, double& perf_vap_wat_rate, DeferredLogger& deferred_logger) const @@ -292,12 +303,20 @@ 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); + } + + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const Value cq_sWat = cq_s[waterCompIdx]; + const Value cq_sGas = cq_s[gasCompIdx]; + const Value vap_wat = rvw * cq_sGas; + const Value dis_gas_wat = rsw * cq_sWat; + cq_s[waterCompIdx] += vap_wat; + cq_s[gasCompIdx] += dis_gas_wat; + if (this->isProducer()) { + perf_vap_wat_rate = getValue(vap_wat); + perf_dis_gas_rate_in_water = getValue(dis_gas_wat); } } @@ -319,9 +338,32 @@ namespace Opm // compute volume ratio between connection at standard conditions Value volumeRatio = bhp * 0.0; // initialize it with the correct type ; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) { const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx]; + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + // Incorporate RSW/RVW factors if both water and gas active + const Value d = 1.0 - rvw * rsw; + + if (d <= 0.0) { + std::ostringstream sstr; + sstr << "Problematic d value " << d << " obtained for well " << this->name() + << " during computePerfRate calculations with rsw " << rsw + << ", rvw " << rvw << " and pressure " << pressure + << " obtaining d " << d + << " Continue as if no dissolution (rsw = 0) and vaporization (rvw = 0) " + << " for this connection."; + deferred_logger.debug(sstr.str()); + } + const Value tmp_wat = d > 0.0? (cmix_s[waterCompIdx] - rvw * cmix_s[gasCompIdx]) / d : cmix_s[waterCompIdx]; + volumeRatio += tmp_wat / b_perfcells_dense[waterCompIdx]; + + const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rsw * cmix_s[waterCompIdx]) / d : cmix_s[waterCompIdx]; + volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx]; + } else { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx]; + } } if constexpr (Indices::enableSolvent) { @@ -407,10 +449,12 @@ namespace Opm 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)) { + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { //no oil const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); perf_vap_wat_rate = getValue(rvw) * getValue(cq_s[gasCompIdx]); + perf_dis_gas_rate_in_water = getValue(rsw) * getValue(cq_s[waterCompIdx]); } } } @@ -458,6 +502,7 @@ namespace Opm ws.vaporized_oil_rate = 0; ws.dissolved_gas_rate = 0; + ws.dissolved_gas_rate_in_water = 0; ws.vaporized_wat_rate = 0; const int np = this->number_of_phases_; @@ -517,6 +562,7 @@ namespace Opm { const auto& comm = this->parallel_well_info_.communication(); ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate); + ws.dissolved_gas_rate_in_water = comm.sum(ws.dissolved_gas_rate_in_water); ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate); ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate); } @@ -584,12 +630,13 @@ namespace Opm getMobilityEval(ebosSimulator, perf, mob, deferred_logger); double perf_dis_gas_rate = 0.; + double perf_dis_gas_rate_in_water = 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, perf_vap_wat_rate, deferred_logger); + cq_s, perf_dis_gas_rate, perf_dis_gas_rate_in_water, 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; @@ -608,6 +655,7 @@ namespace Opm // updating the solution gas rate and solution oil rate if (this->isProducer()) { ws.dissolved_gas_rate += perf_dis_gas_rate; + ws.dissolved_gas_rate_in_water += perf_dis_gas_rate_in_water; ws.vaporized_oil_rate += perf_vap_oil_rate; ws.vaporized_wat_rate += perf_vap_wat_rate; } @@ -1288,6 +1336,7 @@ namespace Opm std::vector& rsmax_perf, std::vector& rvmax_perf, std::vector& rvwmax_perf, + std::vector& rswmax_perf, std::vector& surf_dens_perf) const { std::function getTemperature = @@ -1326,6 +1375,7 @@ namespace Opm rsmax_perf, rvmax_perf, rvwmax_perf, + rswmax_perf, surf_dens_perf); } @@ -1448,6 +1498,7 @@ namespace Opm const std::vector& rsmax_perf, const std::vector& rvmax_perf, const std::vector& rvwmax_perf, + const std::vector& rswmax_perf, const std::vector& surf_dens_perf, DeferredLogger& deferred_logger) { @@ -1481,6 +1532,7 @@ namespace Opm rsmax_perf, rvmax_perf, rvwmax_perf, + rswmax_perf, surf_dens_perf, deferred_logger); } @@ -1503,9 +1555,10 @@ namespace Opm std::vector rsmax_perf; std::vector rvmax_perf; std::vector rvwmax_perf; + std::vector rswmax_perf; std::vector surf_dens_perf; - 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); + computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf); + computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, rswmax_perf, surf_dens_perf, deferred_logger); } @@ -1927,12 +1980,13 @@ namespace Opm std::vector cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.}); double perf_dis_gas_rate = 0.; + double perf_dis_gas_rate_in_water = 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, perf_vap_wat_rate, deferred_logger); + cq_s, perf_dis_gas_rate, perf_dis_gas_rate_in_water, 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/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 34967810f..9a0151158 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -113,6 +113,7 @@ public: static constexpr bool has_foam = getPropValue(); static constexpr bool has_brine = getPropValue(); static constexpr bool has_watVapor = getPropValue(); + static constexpr bool has_disgas_in_water = getPropValue(); static constexpr bool has_saltPrecip = getPropValue(); static constexpr bool has_micp = getPropValue(); @@ -125,6 +126,7 @@ public: has_watVapor, has_brine, has_saltPrecip, + has_disgas_in_water, Indices::numPhases >; /// Constructor WellInterface(const Well& well, diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index a779ccba3..ce1a322f0 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -495,7 +495,7 @@ WellState::report(const int* globalCellIdxMap, well.rates.set(rt::alq, 0.0); } - well.rates.set(rt::dissolved_gas, ws.dissolved_gas_rate); + well.rates.set(rt::dissolved_gas, ws.dissolved_gas_rate + ws.dissolved_gas_rate_in_water); well.rates.set(rt::vaporized_oil, ws.vaporized_oil_rate); well.rates.set(rt::vaporized_water, ws.vaporized_wat_rate);