From 50cf80910b57d56119e99c984ee5cfb7f6b42992 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 17 Sep 2021 15:26:46 +0200 Subject: [PATCH] support aquct and co2store --- .../aquifers/AquiferCarterTracy.hpp | 53 ++++++++++++++++--- opm/simulators/aquifers/AquiferFetkovich.hpp | 3 -- opm/simulators/aquifers/AquiferInterface.hpp | 37 +++++++++---- 3 files changed, 73 insertions(+), 20 deletions(-) diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index 1af0db128..2a654507a 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -50,8 +50,6 @@ public: using typename Base::Simulator; using typename Base::ElementMapper; - using Base::waterCompIdx; - using Base::waterPhaseIdx; AquiferCarterTracy(const std::vector& connections, const Simulator& ebosSimulator, const AquiferCT::AQUCT_data& aquct_data) @@ -83,12 +81,21 @@ public: data.initPressure = this->pa0_; auto* aquCT = data.typeData.template create(); - aquCT->timeConstant = this->aquct_data_.timeConstant(); - aquCT->influxConstant = this->aquct_data_.influxConstant(); - aquCT->waterDensity = this->aquct_data_.waterDensity(); - aquCT->waterViscosity = this->aquct_data_.waterViscosity(); + aquCT->dimensionless_time = this->dimensionless_time_; aquCT->dimensionless_pressure = this->dimensionless_pressure_; + aquCT->influxConstant = this->aquct_data_.influxConstant(); + + if (!this->co2store_()) { + aquCT->timeConstant = this->aquct_data_.timeConstant(); + aquCT->waterDensity = this->aquct_data_.waterDensity(); + aquCT->waterViscosity = this->aquct_data_.waterViscosity(); + } else { + aquCT->waterDensity = this->rhow_; + aquCT->timeConstant = this->Tc_; + const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius; + aquCT->waterViscosity = this->Tc_ * this->aquct_data_.permeability / x; + } return data; } @@ -159,6 +166,11 @@ protected: return std::make_pair(a, b); } + std::size_t pvtRegionIdx() const + { + return this->aquct_data_.pvttableID - 1; + } + // This function implements Eq 5.7 of the EclipseTechnicalDescription inline void calculateInflowRate(int idx, const Simulator& simulator) override { @@ -170,7 +182,19 @@ protected: inline void calculateAquiferConstants() override { - this->Tc_ = this->aquct_data_.timeConstant(); + if(this->co2store_()) { + const auto press = this->aquct_data_.initial_pressure.value(); + Scalar temp = FluidSystem::reservoirTemperature(); + if (this->aquct_data_.initial_temperature.has_value()) + temp = this->aquct_data_.initial_temperature.value(); + + Scalar rs = 0.0; // no dissolved CO2 + Scalar waterViscosity = FluidSystem::oilPvt().viscosity(pvtRegionIdx(), temp, press, rs); + const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius; + this->Tc_ = waterViscosity * x / this->aquct_data_.permeability; + } else { + this->Tc_ = this->aquct_data_.timeConstant(); + } this->beta_ = this->aquct_data_.influxConstant(); } @@ -191,7 +215,20 @@ protected: } this->pa0_ = this->aquct_data_.initial_pressure.value(); - this->rhow_ = this->aquct_data_.waterDensity(); + if(this->co2store_()) { + const auto press = this->aquct_data_.initial_pressure.value(); + + Scalar temp = FluidSystem::reservoirTemperature(); + if (this->aquct_data_.initial_temperature.has_value()) + temp = this->aquct_data_.initial_temperature.value(); + + Scalar rs = 0.0; // no dissolved CO2 + Scalar waterDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rs) + * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIdx()); + this->rhow_ = waterDensity; + } else { + this->rhow_ = this->aquct_data_.waterDensity(); + } } virtual Scalar aquiferDepth() const override diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index 68bb7e7be..0fad5501c 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -50,9 +50,6 @@ public: using typename Base::Simulator; using typename Base::ElementMapper; - using Base::waterCompIdx; - using Base::waterPhaseIdx; - AquiferFetkovich(const std::vector& connections, const Simulator& ebosSimulator, const Aquifetp::AQUFETP_data& aqufetp_data) diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 9f0a0104f..abea05af4 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -72,9 +72,6 @@ public: BlackoilIndices::numPhases> FluidState; - static const auto waterCompIdx = FluidSystem::waterCompIdx; - static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; - // Constructor AquiferInterface(int aqID, const std::vector& connections, @@ -125,7 +122,7 @@ public: elemCtx.updateIntensiveQuantities(0); const auto& iq = elemCtx.intensiveQuantities(0, 0); - pressure_previous_[idx] = getValue(iq.fluidState().pressure(waterPhaseIdx)); + pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_())); } } @@ -150,7 +147,7 @@ public: this->updateCellPressure(this->pressure_current_, idx, intQuants); this->calculateInflowRate(idx, context.simulator()); - rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] + rates[BlackoilIndices::conti0EqIdx + compIdx_()] += this->Qai_[idx] / context.dofVolume(spaceIdx, timeIdx); } @@ -166,6 +163,28 @@ protected: return ebos_simulator_.problem().gravity()[2]; } + inline bool co2store_() const + { + return ebos_simulator_.vanguard().eclState().runspec().co2Storage(); + } + + inline bool phaseIdx_() const + { + if(co2store_()) + return FluidSystem::oilPhaseIdx; + + return FluidSystem::waterPhaseIdx; + } + + inline bool compIdx_() const + { + if(co2store_()) + return FluidSystem::oilCompIdx; + + return FluidSystem::waterCompIdx; + } + + inline void initQuantities() { // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 @@ -188,14 +207,14 @@ protected: updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - pressure_water.at(idx) = fs.pressure(waterPhaseIdx); + pressure_water.at(idx) = fs.pressure(phaseIdx_()); } inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); + pressure_water.at(idx) = fs.pressure(phaseIdx_()).value(); } virtual void endTimeStep() = 0; @@ -380,8 +399,8 @@ protected: const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = iq0.fluidState(); - water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); - const auto water_density = fs.density(waterPhaseIdx); + water_pressure_reservoir = fs.pressure(phaseIdx_()).value(); + const auto water_density = fs.density(phaseIdx_()); const auto gdz = this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());