diff --git a/opm/material/fluidstates/BlackOilFluidState.hpp b/opm/material/fluidstates/BlackOilFluidState.hpp index e7102c1f9..a40c04d1f 100644 --- a/opm/material/fluidstates/BlackOilFluidState.hpp +++ b/opm/material/fluidstates/BlackOilFluidState.hpp @@ -88,6 +88,7 @@ template class BlackOilFluidState { @@ -132,6 +133,10 @@ public: Opm::Valgrind::CheckDefined(*Rv_); } + if (enableSaltWater) { + Opm::Valgrind::CheckDefined(*saltconcentration_); + } + if (enableTemperature || enableEnergy) Opm::Valgrind::CheckDefined(*temperature_); #endif // NDEBUG @@ -155,6 +160,9 @@ public: setRv(Opm::BlackOil::getRv_(fs, pvtRegionIdx)); } + if (enableSaltWater){ + setSaltconcentration(Opm::BlackOil::getSaltconcentration_(fs, pvtRegionIdx)); + } for (unsigned storagePhaseIdx = 0; storagePhaseIdx < numStoragePhases; ++storagePhaseIdx) { unsigned phaseIdx = storageToCanonicalPhaseIndex_(storagePhaseIdx); setSaturation(phaseIdx, fs.saturation(phaseIdx)); @@ -243,6 +251,12 @@ public: void setRv(const Scalar& newRv) { *Rv_ = newRv; } + /*! + * \brief Set the salt concentration. + */ + void setSaltconcentration(const Scalar& newSaltconcentration) + { *saltconcentration_ = newSaltconcentration; } + /*! * \brief Return the pressure of a fluid phase [Pa] */ @@ -311,6 +325,19 @@ public: return *Rv_; } + /*! + * \brief Return the concentration of salt in water + */ + const Scalar& saltconcentration() const + { + if (!enableSaltWater) { + static Scalar null = 0.0; + return null; + } + + return *saltconcentration_; + } + /*! * \brief Return the PVT region where the current fluid state is assumed to be part of. * @@ -517,6 +544,7 @@ private: std::array density_; Opm::ConditionalStorage Rs_; Opm::ConditionalStorage Rv_; + Opm::ConditionalStorage saltconcentration_; unsigned short pvtRegionIdx_; }; diff --git a/opm/material/fluidsystems/BlackOilFluidSystem.hpp b/opm/material/fluidsystems/BlackOilFluidSystem.hpp index d087f0088..6e91466d5 100644 --- a/opm/material/fluidsystems/BlackOilFluidSystem.hpp +++ b/opm/material/fluidsystems/BlackOilFluidSystem.hpp @@ -79,6 +79,21 @@ auto getRv_(typename std::enable_if::value, const Fluid -> decltype(Opm::decay(fluidState.Rv())) { return Opm::decay(fluidState.Rv()); } +template +LhsEval getSaltconcentration_(typename std::enable_if::value, const FluidState&>::type fluidState, + unsigned regionIdx) +{ + //const auto& XoG = + // Opm::decay(fluidState.massFraction(FluidSystem::oilPhaseIdx, FluidSystem::gasCompIdx)); + //return FluidSystem::convertXoGToRs(XoG, regionIdx); +} + +template +auto getSaltconcentration_(typename std::enable_if::value, const FluidState&>::type fluidState, + unsigned regionIdx OPM_UNUSED) + -> decltype(Opm::decay(fluidState.saltconcentration())) +{ return Opm::decay(fluidState.saltconcentration()); } + } /*! diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilitySaltWaterPvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilitySaltWaterPvt.hpp new file mode 100644 index 000000000..949e4592b --- /dev/null +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilitySaltWaterPvt.hpp @@ -0,0 +1,210 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::ConstantCompressibilitySaltWaterPvt + */ +#ifndef OPM_CONSTANT_COMPRESSIBILITY_SALTWATER_PVT_HPP +#define OPM_CONSTANT_COMPRESSIBILITY_SALTWATER_PVT_HPP + +#include + +#if HAVE_ECL_INPUT +#include +#include +#include +#include +#include +#include +#endif + +#include + +namespace Opm { +template +class WaterPvtMultiplexer; +/*! + * \brief This class represents the Pressure-Volume-Temperature relations of the gas phase + * without vaporized oil. + */ +template +class ConstantCompressibilitySaltWaterPvt +{ + typedef Opm::Tabulated1DFunction TabulatedOneDFunction; + typedef typename Opm::Tabulated1DFunction TabulatedFunction; + typedef std::vector > SamplingPoints; + +public: +#if HAVE_ECL_INPUT + /*! + * \brief Sets the pressure-dependent water viscosity and density + * using a table stemming from the Eclipse PVTWSALT keyword. + */ + void initFromDeck(const Deck& deck, const EclipseState& eclState) + { + const auto& tableManager = eclState.getTableManager(); + size_t numRegions = tableManager.getTabdims().getNumPVTTables(); + const auto& densityKeyword = deck.getKeyword("DENSITY"); + + formationVolumeTables_.resize(numRegions); + compressibilityTables_.resize(numRegions); + viscosityTables_.resize(numRegions); + viscosibilityTables_.resize(numRegions); + referencePressure_.resize(numRegions); + + const auto& pvtwsaltTables = tableManager.getPvtwSaltTables(); + if(!pvtwsaltTables.empty()){ + assert(numRegions == pvtwsaltTables.size()); + for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { + const auto& pvtwsaltTable = pvtwsaltTables[regionIdx]; + const auto& c = pvtwsaltTable.getSaltConcentrationColumn(); + + const auto& B = pvtwsaltTable.getFormationVolumeFactorColumn(); + formationVolumeTables_[regionIdx].setXYContainers(c, B); + + const auto& compressibility = pvtwsaltTable.getCompressibilityColumn(); + compressibilityTables_[regionIdx].setXYContainers(c, compressibility); + + const auto& viscositytable = pvtwsaltTable.getViscosityColumn(); + viscosityTables_[regionIdx].setXYContainers(c, viscositytable); + + const auto& viscosibility = pvtwsaltTable.getViscosibilityColumn(); + viscosibilityTables_[regionIdx].setXYContainers(c, viscosibility); + referencePressure_[regionIdx] = pvtwsaltTable.getReferencePressureValue(); + } + } + else { + throw std::runtime_error("PVTWSALT must be specified in SALTWATER runs\n"); + } + + + size_t numPvtwRegions = numRegions; + setNumRegions(numPvtwRegions); + + for (unsigned regionIdx = 0; regionIdx < numPvtwRegions; ++ regionIdx) { + auto densityRecord = densityKeyword.getRecord(regionIdx); + + waterReferenceDensity_[regionIdx] = + densityRecord.getItem("WATER").getSIDouble(0); + } + + initEnd(); + } +#endif + + void setNumRegions(size_t numRegions) + { + waterReferenceDensity_.resize(numRegions); + + for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + setReferenceDensities(regionIdx, 650.0, 1.0, 1000.0); + } + } + + /*! + * \brief Set the water reference density [kg / m^3] + */ + void setReferenceDensities(unsigned regionIdx, + Scalar /*rhoRefOil*/, + Scalar /*rhoRefGas*/, + Scalar rhoRefWater) + { waterReferenceDensity_[regionIdx] = rhoRefWater; } + + /*! + * \brief Finish initializing the water phase PVT properties. + */ + void initEnd() + { } + + + /*! + * \brief Return the number of PVT regions which are considered by this PVT-object. + */ + unsigned numRegions() const + { return waterReferenceDensity_.size(); } + + /*! + * \brief Returns the specific enthalpy [J/kg] of water given a set of parameters. + */ + template + Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED, + const Evaluation& temperature OPM_UNUSED, + const Evaluation& pressure OPM_UNUSED) const + { + throw std::runtime_error("Requested the enthalpy of water but the thermal option is not enabled"); + } + + /*! + * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. + */ + template + Evaluation viscosity(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& saltwaterconcentration) const + { + // cf. ECLiPSE 2013.2 technical description, p. 114 + Scalar pRef = referencePressure_[regionIdx]; + const Evaluation C = compressibilityTables_[regionIdx].eval(saltwaterconcentration, /*extrapolate=*/true); + const Evaluation Cv = viscosibilityTables_[regionIdx].eval(saltwaterconcentration, /*extrapolate=*/true); + const Evaluation BwRef = formationVolumeTables_[regionIdx].eval(saltwaterconcentration, /*extrapolate=*/true); + const Evaluation Y = (C-Cv)* (pressure - pRef); + Evaluation MuwRef = viscosityTables_[regionIdx].eval(saltwaterconcentration, /*extrapolate=*/true); + + const Evaluation& bw = inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltwaterconcentration); + + return MuwRef*BwRef*bw/(1 + Y*(1 + Y/2)); + } + + /*! + * \brief Returns the formation volume factor [-] of the fluid phase. + */ + template + Evaluation inverseFormationVolumeFactor(unsigned regionIdx, + const Evaluation& /*temperature*/, + const Evaluation& pressure, + const Evaluation& saltwaterconcentration) const + { + Scalar pRef = referencePressure_[regionIdx]; + + const Evaluation BwRef = formationVolumeTables_[regionIdx].eval(saltwaterconcentration, /*extrapolate=*/true); + const Evaluation C = compressibilityTables_[regionIdx].eval(saltwaterconcentration, /*extrapolate=*/true); + const Evaluation X = C * (pressure - pRef); + + return (1.0 + X*(1.0 + X/2.0))/BwRef; + + } + +private: + std::vector formationVolumeTables_; + std::vector compressibilityTables_; + std::vector viscosityTables_; + std::vector viscosibilityTables_; + std::vector referencePressure_; + std::vector waterReferenceDensity_; + +}; + +} // namespace Opm + +#endif diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp index 698334bc3..69c55168f 100644 --- a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp @@ -184,16 +184,18 @@ public: throw std::runtime_error("Requested the enthalpy of water but the thermal option is not enabled"); } + /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ template Evaluation viscosity(unsigned regionIdx, const Evaluation& temperature, - const Evaluation& pressure) const + const Evaluation& pressure, + const Evaluation& saltconcentration) const { Scalar BwMuwRef = waterViscosity_[regionIdx]*waterReferenceFormationVolumeFactor_[regionIdx]; - const Evaluation& bw = inverseFormationVolumeFactor(regionIdx, temperature, pressure); + const Evaluation& bw = inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration); Scalar pRef = waterReferencePressure_[regionIdx]; const Evaluation& Y = @@ -208,7 +210,8 @@ public: template Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation& /*temperature*/, - const Evaluation& pressure) const + const Evaluation& pressure, + const Evaluation& /*saltconcentration*/) const { // cf. ECLiPSE 2011 technical description, p. 116 Scalar pRef = waterReferencePressure_[regionIdx]; diff --git a/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp b/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp index 5154af280..e3435f1b8 100644 --- a/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp @@ -28,6 +28,7 @@ #define OPM_WATER_PVT_MULTIPLEXER_HPP #include "ConstantCompressibilityWaterPvt.hpp" +#include "ConstantCompressibilitySaltWaterPvt.hpp" #include "WaterPvtThermal.hpp" #define OPM_WATER_PVT_MULTIPLEXER_CALL(codeToCall) \ @@ -37,6 +38,11 @@ codeToCall; \ break; \ } \ + case ConstantCompressibilitySaltWaterPvt: { \ + auto& pvtImpl = getRealPvt(); \ + codeToCall; \ + break; \ + } \ case ThermalWaterPvt: { \ auto& pvtImpl = getRealPvt(); \ codeToCall; \ @@ -51,7 +57,7 @@ namespace Opm { * \brief This class represents the Pressure-Volume-Temperature relations of the water * phase in the black-oil model. */ -template +template class WaterPvtMultiplexer { public: @@ -59,6 +65,7 @@ public: enum WaterPvtApproach { NoWaterPvt, + ConstantCompressibilitySaltWaterPvt, ConstantCompressibilityWaterPvt, ThermalWaterPvt }; @@ -86,6 +93,10 @@ public: delete &getRealPvt(); break; } + case ConstantCompressibilitySaltWaterPvt: { + delete &getRealPvt(); + break; + } case ThermalWaterPvt: { delete &getRealPvt(); break; @@ -111,6 +122,8 @@ public: setApproach(ThermalWaterPvt); else if (deck.hasKeyword("PVTW")) setApproach(ConstantCompressibilityWaterPvt); + else if (enableSaltWater && deck.hasKeyword("PVTWSALT")) + setApproach(ConstantCompressibilitySaltWaterPvt); OPM_WATER_PVT_MULTIPLEXER_CALL(pvtImpl.initFromDeck(deck, eclState)); } @@ -141,7 +154,37 @@ public: Evaluation viscosity(unsigned regionIdx, const Evaluation& temperature, const Evaluation& pressure) const - { OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.viscosity(regionIdx, temperature, pressure)); return 0; } + { + // assert(realWaterPvt_ != ConstantCompressibilitySaltWaterPvt ); + const Evaluation saltconcentration = 0.0; + OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.viscosity(regionIdx, temperature, pressure, saltconcentration)); + return 0; + } + + /*! + * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. + */ + template + Evaluation viscosity(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& saltconcentration) const + { + OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.viscosity(regionIdx, temperature, pressure, saltconcentration)); + return 0; + } + + /*! + * \brief Returns the formation volume factor [-] of the fluid phase. + */ + template + Evaluation inverseFormationVolumeFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& saltconcentration) const + { OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration)); + return 0; + } /*! * \brief Returns the formation volume factor [-] of the fluid phase. @@ -150,7 +193,11 @@ public: Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation& temperature, const Evaluation& pressure) const - { OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.inverseFormationVolumeFactor(regionIdx, temperature, pressure)); return 0; } + { + const Evaluation saltconcentration = 0.0; + OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltconcentration)); + return 0; + } void setApproach(WaterPvtApproach appr) { @@ -159,6 +206,10 @@ public: realWaterPvt_ = new Opm::ConstantCompressibilityWaterPvt; break; + case ConstantCompressibilitySaltWaterPvt: + realWaterPvt_ = new Opm::ConstantCompressibilitySaltWaterPvt; + break; + case ThermalWaterPvt: realWaterPvt_ = new Opm::WaterPvtThermal; break; @@ -193,6 +244,20 @@ public: return *static_cast* >(realWaterPvt_); } + template + typename std::enable_if >::type& getRealPvt() + { + assert(approach() == approachV); + return *static_cast* >(realWaterPvt_); + } + + template + typename std::enable_if >::type& getRealPvt() const + { + assert(approach() == approachV); + return *static_cast* >(realWaterPvt_); + } + template typename std::enable_if >::type& getRealPvt() { diff --git a/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp index 73ded5c55..8cdcd4b61 100644 --- a/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp @@ -42,7 +42,7 @@ #endif namespace Opm { -template +template class WaterPvtMultiplexer; /*! @@ -56,7 +56,7 @@ class WaterPvtThermal { public: typedef Opm::Tabulated1DFunction TabulatedOneDFunction; - typedef WaterPvtMultiplexer IsothermalPvt; + typedef WaterPvtMultiplexer IsothermalPvt; WaterPvtThermal() { @@ -269,9 +269,10 @@ public: template Evaluation viscosity(unsigned regionIdx, const Evaluation& temperature, - const Evaluation& pressure) const + const Evaluation& pressure, + const Evaluation& saltwaterconcentration) const { - const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure); + const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, saltwaterconcentration); if (!enableThermalViscosity()) return isothermalMu; @@ -289,10 +290,11 @@ public: template Evaluation inverseFormationVolumeFactor(unsigned regionIdx, const Evaluation& temperature, - const Evaluation& pressure) const + const Evaluation& pressure, + const Evaluation& saltwaterconcentration) const { if (!enableThermalDensity()) - return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure); + return isothermalPvt_->inverseFormationVolumeFactor(regionIdx, temperature, pressure, saltwaterconcentration); Scalar BwRef = pvtwRefB_[regionIdx]; Scalar TRef = watdentRefTemp_[regionIdx]; diff --git a/tests/test_eclblackoilpvt.cpp b/tests/test_eclblackoilpvt.cpp index 8a9b8534f..e61aa9076 100644 --- a/tests/test_eclblackoilpvt.cpp +++ b/tests/test_eclblackoilpvt.cpp @@ -247,7 +247,8 @@ inline void testAll() refTmp = 1.1e-3; // the deck value is given in cP, while the SI units use Pa s... tmp = constCompWaterPvt.viscosity(/*regionIdx=*/0, /*temperature=*/273.15 + 20.0, - /*pressure=*/1e5); + /*pressure=*/1e5, + /*saltconcentration=*/0.0); if (std::abs(tmp - refTmp) > tolerance) throw std::logic_error("The reference water viscosity at region 0 is supposed to be "+std::to_string(refTmp) +". (is "+std::to_string(tmp)+")"); @@ -255,7 +256,8 @@ inline void testAll() refTmp = 1.2e-3; tmp = constCompWaterPvt.viscosity(/*regionIdx=*/1, /*temperature=*/273.15 + 20.0, - /*pressure=*/2e5); + /*pressure=*/2e5, + /*saltconcentration=*/0.0); if (std::abs(tmp - refTmp) > tolerance) throw std::logic_error("The reference water viscosity at region 1 is supposed to be "+std::to_string(refTmp) +". (is "+std::to_string(tmp)+")");