diff --git a/opm/material/binarycoefficients/Brine_CO2.hpp b/opm/material/binarycoefficients/Brine_CO2.hpp index f4691898b..808934a7f 100644 --- a/opm/material/binarycoefficients/Brine_CO2.hpp +++ b/opm/material/binarycoefficients/Brine_CO2.hpp @@ -28,9 +28,6 @@ #ifndef OPM_BINARY_COEFF_BRINE_CO2_HPP #define OPM_BINARY_COEFF_BRINE_CO2_HPP -#include -#include -#include #include namespace Opm { @@ -40,10 +37,8 @@ namespace BinaryCoeff { * \ingroup Binarycoefficients * \brief Binary coefficients for brine and CO2. */ -template +template class Brine_CO2 { - typedef Opm::H2O H2O; - typedef Opm::CO2 CO2; typedef Opm::IdealGas IdealGas; static const int liquidPhaseIdx = 0; // index of the liquid phase static const int gasPhaseIdx = 1; // index of the gas phase diff --git a/opm/material/components/SimpleHuDuanH2O.hpp b/opm/material/components/SimpleHuDuanH2O.hpp new file mode 100644 index 000000000..7388bd51a --- /dev/null +++ b/opm/material/components/SimpleHuDuanH2O.hpp @@ -0,0 +1,399 @@ +// -*- 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::SimpleHuDuanH2O + */ +#ifndef OPM_SIMPLE_HU_DUAN_H2O_HPP +#define OPM_SIMPLE_HU_DUAN_H2O_HPP + +#include "Component.hpp" +#include "iapws/Common.hpp" + + +#include + +#include + +#include + +#include + +namespace Opm { + + /*! + * \ingroup Components + * + * \brief A simple version of pure water with density from Hu et al. + * + * Compared to the water formulation of IAPWS'97, this class provides + * a much simpler component that represents the thermodynamic + * properties of pure water. This implies that the likelyhood for + * bugs in this class is reduced and the numerical performance is + * increased. (At the cost of accuracy for the representation of the + * physical quantities, of course.) + * + * Density from Hu, Duan, Zhu and Chou: PVTx properties of the CO2-H2O and CO2-H2O-NaCl + * systems below 647 K: Assessment of experimental data and + * thermodynamics models, Chemical Geology, 2007. + * + * \tparam Scalar The type used for representing scalar values + */ +template +class SimpleHuDuanH2O : public Component > +{ + typedef Opm::IdealGas IdealGas; + typedef IAPWS::Common Common; + + static const Scalar R; // specific gas constant of water + +public: + /*! + * \brief A human readable name for the water. + */ + static const char* name() + { return "H2O"; } + + /*! + * \brief Returns true iff the gas phase is assumed to be compressible + */ + static bool gasIsCompressible() + { return true; } + + /*! + * \brief Returns true iff the liquid phase is assumed to be compressible + */ + static bool liquidIsCompressible() + { return false; } + + /*! + * \brief Returns true iff the gas phase is assumed to be ideal + */ + static bool gasIsIdeal() + { return true; } + + /*! + * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of water. + */ + static Scalar molarMass() + { return 18e-3; } + + /*! + * \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of water. + */ + static Scalar criticalTemperature() + { return 647.096; /* [K] */ } + + /*! + * \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of water. + */ + static Scalar criticalPressure() + { return 22.064e6; /* [N/m^2] */ } + + /*! + * \brief Returns the temperature \f$\mathrm{[K]}\f$ at water's triple point. + */ + static Scalar tripleTemperature() + { return 273.16; /* [K] */ } + + /*! + * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at water's triple point. + */ + static Scalar triplePressure() + { return 611.657; /* [N/m^2] */ } + + /*! + * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure water + * at a given temperature. + * + *\param T temperature of component in \f$\mathrm{[K]}\f$ + * + * See: + * + * IAPWS: "Revised Release on the IAPWS Industrial Formulation + * 1997 for the Thermodynamic Properties of Water and Steam", + * http://www.iapws.org/relguide/IF97-Rev.pdf + */ + template + static Evaluation vaporPressure(const Evaluation& T) + { + if (T > criticalTemperature()) + return criticalPressure(); + if (T < tripleTemperature()) + return 0; // water is solid: We don't take sublimation into account + + static const Scalar n[10] = { + 0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2, + 0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2, + -0.48232657361591e4, 0.40511340542057e6, -0.23855557567849, + 0.65017534844798e3 + }; + + Evaluation sigma = T + n[8]/(T - n[9]); + + Evaluation A = (sigma + n[0])*sigma + n[1]; + Evaluation B = (n[2]*sigma + n[3])*sigma + n[4]; + Evaluation C = (n[5]*sigma + n[6])*sigma + n[7]; + + Evaluation tmp = 2.0*C/(Opm::sqrt(B*B - 4.0*A*C) - B); + tmp *= tmp; + tmp *= tmp; + + return 1e6*tmp; + } + + /*! + * \brief Specific enthalpy of water steam \f$\mathrm{[J/kg]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation gasEnthalpy(const Evaluation& temperature, + const Evaluation& /*pressure*/) + { return 1.976e3*temperature + 40.65e3/molarMass(); } + + + /*! + * \copydoc Component::gasHeatCapacity + */ + template + static Evaluation gasHeatCapacity(const Evaluation& temperature OPM_UNUSED, + const Evaluation& pressure OPM_UNUSED) + { return 1.976e3; } + + /*! + * \brief Specific enthalpy of liquid water \f$\mathrm{[J/kg]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation liquidEnthalpy(const Evaluation& temperature, + const Evaluation& /*pressure*/) + { return 4180*temperature; } + + /*! + * \copydoc Component::liquidHeatCapacity + */ + template + static Evaluation liquidHeatCapacity(const Evaluation& temperature OPM_UNUSED, + const Evaluation& pressure OPM_UNUSED) + { return 4.184e3; } + + /*! + * \brief Specific internal energy of steam \f$\mathrm{[J/kg]}\f$. + * + * Definition of enthalpy: \f$h= u + pv = u + p / \rho\f$. + * + * Rearranging for internal energy yields: \f$u = h - pv\f$. + * + * Exploiting the Ideal Gas assumption (\f$pv = R_{\textnormal{specific}} T\f$)gives: \f$u = h - R / M T \f$. + * + * The universal gas constant can only be used in the case of molar formulations. + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation gasInternalEnergy(const Evaluation& temperature, + const Evaluation& pressure) + { + return + gasEnthalpy(temperature, pressure) - + 1/molarMass()* // conversion from [J/(mol K)] to [J/(kg K)] + IdealGas::R*temperature; // = pressure *spec. volume for an ideal gas + } + + /*! + * \brief Specific internal energy of liquid water \f$\mathrm{[J/kg]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation liquidInternalEnergy(const Evaluation& temperature, + const Evaluation& pressure) + { + return + liquidEnthalpy(temperature, pressure) - + pressure/liquidDensity(temperature, pressure); + } + + /*! + * \brief Specific heat conductivity of liquid water \f$\mathrm{[W/(m K)]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation liquidThermalConductivity(const Evaluation& /*temperature*/, + const Evaluation& /*pressure*/) + { + return 0.578078; // conductivity of liquid water [W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C + } + + /*! + * \brief Specific heat conductivity of steam \f$\mathrm{[W/(m K)]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation gasThermalConductivity(const Evaluation& /*temperature*/, + const Evaluation& /*pressure*/) + { + return 0.028224; // conductivity of steam [W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8°C + } + + /*! + * \brief The density \f$\mathrm{[kg/m^3]}\f$ of steam at a given pressure and temperature. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure) + { + // Assume an ideal gas + return molarMass()*IdealGas::molarDensity(temperature, pressure); + } + + /*! + * \brief The pressure of steam in \f$\mathrm{[Pa]}\f$ at a given density and temperature. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param density density of component in \f$\mathrm{[kg/m^3]}\f$ + */ + template + static Evaluation gasPressure(const Evaluation& temperature, const Evaluation& density) + { + // Assume an ideal gas + return IdealGas::pressure(temperature, density/molarMass()); + } + + /*! + * \brief The density of pure water at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure) + { + return liquidDensity_(temperature, pressure); + } + + /*! + * \brief The pressure of water in \f$\mathrm{[Pa]}\f$ at a given density and temperature. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param density density of component in \f$\mathrm{[kg/m^3]}\f$ + */ + template + static Evaluation liquidPressure(const Evaluation& /*temperature*/, const Evaluation& /*density*/) + { + throw std::logic_error("The liquid pressure is undefined for incompressible fluids"); + } + + /*! + * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of steam. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + * \param regularize defines, if the functions is regularized or not, set to true by default + */ + template + static Evaluation gasViscosity(const Evaluation& /*temperature*/, + const Evaluation& /*pressure*/) + { + return 1e-05; + } + + /*! + * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure water. + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure) + { + if (temperature > 570) { + std::ostringstream oss; + oss << "Viscosity of water based on Hu et al is too different from IAPWS for T above 570K and " + << "(T = " << temperature << ")"; + throw NumericalIssue(oss.str()); + } + + const Evaluation& rho = liquidDensity(temperature, pressure); + return Common::viscosity(temperature, rho); + } + +private: + + /*! + * \brief The density of pure water in \f$\mathrm{[kg/m3]}\f$ + * + * \param temperature temperature of component in \f$\mathrm{[K]}\f$ + * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ + */ + template + static Evaluation liquidDensity_(const Evaluation& T, const Evaluation& pressure) { + // Hu, Duan, Zhu and Chou: PVTx properties of the CO2-H2O and CO2-H2O-NaCl + // systems below 647 K: Assessment of experimental data and + // thermodynamics models, Chemical Geology, 2007. + if (T > 647 || pressure > 100e6) { + std::ostringstream oss; + oss << "Density of water is only implemented for temperatures below 647K and " + << "pressures below 100MPa. (T = " << T << ", p=" << pressure; + throw NumericalIssue(oss.str()); + } + + Evaluation p = pressure / 1e6; // to MPa + Scalar Mw = molarMass() * 1e3; //kg/kmol + + static const Scalar k0[5] = { 3.27225e-07, -4.20950e-04, 2.32594e-01, -4.16920e+01, 5.71292e+03 }; + static const Scalar k1[5] = { -2.32306e-10, 2.91138e-07, -1.49662e-04, 3.59860e-02, -3.55071 }; + static const Scalar k2[3] = { 2.57241e-14, -1.24336e-11, 5.42707e-07 }; + static const Scalar k3[3] = { -4.42028e-18, 2.10007e-15, -8.11491e-11 }; + Evaluation k0_eval = 1e-3 * (((k0[0]*T + k0[1])*T + k0[2])*T + k0[3] + k0[4]/T); + Evaluation k1_eval = 1e-2 * (((k1[0]*T + k1[1])*T + k1[2])*T + k1[3] + k1[4]/T); + Evaluation k2_eval = 1e-1 * ((k2[0]*T + k2[1])*T*T + k2[2]); + Evaluation k3_eval = (k3[0]*T + k3[1])*T*T + k3[2]; + + // molar volum (m³/kmol): + Evaluation vw = ((k3_eval*p + k2_eval)*p + k1_eval)*p + k0_eval; + + // density kg/m3 + return Mw / vw; + + } + +}; + +template +const Scalar SimpleHuDuanH2O::R = Opm::Constants::R / 18e-3; + +} // namespace Opm + +#endif diff --git a/opm/material/components/co2tables.inc b/opm/material/components/co2tables.inc index 1c327aae9..e94562aed 100644 --- a/opm/material/components/co2tables.inc +++ b/opm/material/components/co2tables.inc @@ -1,3 +1,6 @@ +#ifndef CO2TABLES_INC +#define CO2TABLES_INC + /* Tables for CO2 fluid properties calculated according to Span and * Wagner (1996). * @@ -23,13 +26,13 @@ struct TabulatedDensityTraits { static const Scalar vals[numX][numY]; }; -const double TabulatedDensityTraits::xMin = 2.800000000000000e+02; -const double TabulatedDensityTraits::xMax = 4.000000000000000e+02; -const double TabulatedDensityTraits::yMin = 1.000000000000000e+05; -const double TabulatedDensityTraits::yMax = 1.000000000000000e+08; -const char *TabulatedDensityTraits::name = "density"; +inline const double TabulatedDensityTraits::xMin = 2.800000000000000e+02; +inline const double TabulatedDensityTraits::xMax = 4.000000000000000e+02; +inline const double TabulatedDensityTraits::yMin = 1.000000000000000e+05; +inline const double TabulatedDensityTraits::yMax = 1.000000000000000e+08; +inline const char *TabulatedDensityTraits::name = "density"; -const double TabulatedDensityTraits::vals[200][500] = +inline const double TabulatedDensityTraits::vals[200][500] = { { 1.902062465274410e+00, 5.782408947482105e+00, 9.764909779046345e+00, 1.385694970929571e+01, 1.806682704790397e+01, @@ -20445,13 +20448,13 @@ struct TabulatedEnthalpyTraits { static const Scalar vals[numX][numY]; }; -const double TabulatedEnthalpyTraits::xMin = 2.800000000000000e+02; -const double TabulatedEnthalpyTraits::xMax = 4.000000000000000e+02; -const double TabulatedEnthalpyTraits::yMin = 1.000000000000000e+05; -const double TabulatedEnthalpyTraits::yMax = 1.000000000000000e+08; -const char *TabulatedEnthalpyTraits::name = "enthalpy"; +inline const double TabulatedEnthalpyTraits::xMin = 2.800000000000000e+02; +inline const double TabulatedEnthalpyTraits::xMax = 4.000000000000000e+02; +inline const double TabulatedEnthalpyTraits::yMin = 1.000000000000000e+05; +inline const double TabulatedEnthalpyTraits::yMax = 1.000000000000000e+08; +inline const char *TabulatedEnthalpyTraits::name = "enthalpy"; -const double TabulatedEnthalpyTraits::vals[200][500] = +inline const double TabulatedEnthalpyTraits::vals[200][500] = { { 5.700580815972213e+03, 3.531328941638982e+03, 1.311308233079838e+03, -9.631993272332195e+02, -3.296365473729317e+03, @@ -40864,13 +40867,13 @@ struct CO2Tables { static const double brineSalinity; }; -TabulatedFunction CO2Tables::tabulatedEnthalpy; -TabulatedFunction CO2Tables::tabulatedDensity; -const double CO2Tables::brineSalinity = 1.000000000000000e-01; +inline TabulatedFunction CO2Tables::tabulatedEnthalpy; +inline TabulatedFunction CO2Tables::tabulatedDensity; +inline const double CO2Tables::brineSalinity = 1.000000000000000e-01; // initialize the static tables once. this is a bit hacky in so far as it uses some // advanced C++ features (static initializer functions) -int initCO2Tables_(); +inline int initCO2Tables_(); int initCO2Tables_() { CO2Tables::tabulatedEnthalpy.resize(TabulatedEnthalpyTraits::xMin, @@ -40899,4 +40902,6 @@ int initCO2Tables_() } extern int co2TablesInitDummy__; -int co2TablesInitDummy__ = initCO2Tables_(); +inline int co2TablesInitDummy__ = initCO2Tables_(); + +#endif /* CO2TABLES_INC */ diff --git a/opm/material/fluidsystems/BlackOilFluidSystem.hpp b/opm/material/fluidsystems/BlackOilFluidSystem.hpp index 16b15afca..ffbd12361 100644 --- a/opm/material/fluidsystems/BlackOilFluidSystem.hpp +++ b/opm/material/fluidsystems/BlackOilFluidSystem.hpp @@ -31,6 +31,7 @@ #include "blackoilpvt/OilPvtMultiplexer.hpp" #include "blackoilpvt/GasPvtMultiplexer.hpp" #include "blackoilpvt/WaterPvtMultiplexer.hpp" +#include "blackoilpvt/BrineCo2Pvt.hpp" #include #include @@ -183,8 +184,7 @@ public: */ static void initFromState(const EclipseState& eclState, const Schedule& schedule) { - const auto& densityTable = eclState.getTableManager().getDensityTable(); - size_t numRegions = densityTable.size(); + size_t numRegions = eclState.runspec().tabdims().getNumPVTTables(); initBegin(numRegions); numActivePhases_ = 0; @@ -220,14 +220,6 @@ public: setEnableDissolvedGas(eclState.getSimulationConfig().hasDISGAS()); setEnableVaporizedOil(eclState.getSimulationConfig().hasVAPOIL()); - // set the reference densities of all PVT regions - for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - setReferenceDensities(densityTable[regionIdx].oil, - densityTable[regionIdx].water, - densityTable[regionIdx].gas, - regionIdx); - } - if (phaseIsActive(gasPhaseIdx)) { gasPvt_ = std::make_shared(); gasPvt_->initFromState(eclState, schedule); @@ -243,6 +235,14 @@ public: waterPvt_->initFromState(eclState, schedule); } + // set the reference densities of all PVT regions + for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + setReferenceDensities(phaseIsActive(oilPhaseIdx)? oilPvt_->oilReferenceDensity(regionIdx):700., + phaseIsActive(waterPhaseIdx)? waterPvt_->waterReferenceDensity(regionIdx):1000., + phaseIsActive(gasPhaseIdx)? gasPvt_->gasReferenceDensity(regionIdx):2., + regionIdx); + } + initEnd(); } #endif // HAVE_ECL_INPUT diff --git a/opm/material/fluidsystems/BrineCO2FluidSystem.hpp b/opm/material/fluidsystems/BrineCO2FluidSystem.hpp index 24553a9f1..739381125 100644 --- a/opm/material/fluidsystems/BrineCO2FluidSystem.hpp +++ b/opm/material/fluidsystems/BrineCO2FluidSystem.hpp @@ -72,8 +72,13 @@ public: struct ParameterCache : public Opm::NullParameterCache {}; + //! The type of the component for brine used by the fluid system + typedef Brine_Tabulated Brine; + //! The type of the component for pure CO2 used by the fluid system + typedef Opm::CO2 CO2; + //! The binary coefficients for brine and CO2 used by this fluid system - typedef Opm::BinaryCoeff::Brine_CO2 BinaryCoeffBrineCO2; + typedef Opm::BinaryCoeff::Brine_CO2 BinaryCoeffBrineCO2; /**************************************** * Fluid phase related static parameters @@ -154,11 +159,6 @@ public: //! The index of the CO2 component static const int CO2Idx = 1; - //! The type of the component for brine used by the fluid system - typedef Brine_Tabulated Brine; - //! The type of the component for pure CO2 used by the fluid system - typedef Opm::CO2 CO2; - /*! * \copydoc BaseFluidSystem::componentName */ diff --git a/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp b/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp new file mode 100644 index 000000000..4e1af295c --- /dev/null +++ b/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp @@ -0,0 +1,427 @@ +// -*- 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::BrineCo2Pvt + */ +#ifndef OPM_BRINE_CO2_PVT_HPP +#define OPM_BRINE_CO2_PVT_HPP + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#if HAVE_ECL_INPUT +#include +#include +#include +#endif + +#include + +namespace Opm { +/*! + * \brief This class represents the Pressure-Volume-Temperature relations of the liquid phase + * for a CO2-Brine system + */ +template +class BrineCo2Pvt +{ + typedef std::vector > SamplingPoints; + typedef Opm::SimpleHuDuanH2O H2O; + typedef Opm::Brine Brine; + + //typedef Opm::H2O H2O_IAPWS; + //typedef Opm::Brine Brine_IAPWS; + //typedef Opm::TabulatedComponent H2O_Tabulated; + //typedef Opm::TabulatedComponent Brine_Tabulated; + + //typedef H2O_Tabulated H2O; + //typedef Brine_Tabulated Brine; + + typedef Opm::CO2 CO2; + +public: + typedef Opm::Tabulated1DFunction TabulatedOneDFunction; + + //! The binary coefficients for brine and CO2 used by this fluid system + typedef Opm::BinaryCoeff::Brine_CO2 BinaryCoeffBrineCO2; + + explicit BrineCo2Pvt() = default; + BrineCo2Pvt(const std::vector& brineReferenceDensity, + const std::vector& co2ReferenceDensity, + const std::vector& salinity) + : brineReferenceDensity_(brineReferenceDensity), + co2ReferenceDensity_(co2ReferenceDensity), + salinity_(salinity) + { + } +#if HAVE_ECL_INPUT + /*! + * \brief Initialize the parameters for Brine-CO2 system using an ECL deck. + * + */ + void initFromState(const EclipseState& eclState, const Schedule&) + { + if( !eclState.getTableManager().getDensityTable().empty()) { + std::cerr << "WARNING: CO2STOR is enabled but DENSITY is in the deck. \n" << + "The surface density is computed based on CO2-BRINE PVT at standard conditions (STCOND) and DENSITY is ignored " << std::endl; + } + + if( eclState.getTableManager().hasTables("PVDO") || !eclState.getTableManager().getPvtoTables().empty()) { + std::cerr << "WARNING: CO2STOR is enabled but PVDO or PVTO is in the deck. \n" << + "BRINE PVT properties are computed based on the Hu et al. pvt model and PVDO/PVTO input is ignored. " << std::endl; + } + + // We only supported single pvt region for the co2-brine module + size_t numRegions = 1; + setNumRegions(numRegions); + size_t regionIdx = 0; + // Currently we only support constant salinity + const Scalar molality = eclState.getTableManager().salinity(); // mol/kg + const Scalar MmNaCl = 58e-3; // molar mass of NaCl [kg/mol] + // convert to mass fraction + Brine::salinity = 1 / ( 1 + 1 / (molality*MmNaCl)); // + salinity_[numRegions] = Brine::salinity; + // set the surface conditions using the STCOND keyword + Scalar T_ref = eclState.getTableManager().stCond().temperature; + Scalar P_ref = eclState.getTableManager().stCond().pressure; + + brineReferenceDensity_[regionIdx] = Brine::liquidDensity(T_ref, P_ref); + co2ReferenceDensity_[regionIdx] = CO2::gasDensity(T_ref, P_ref); + } +#endif + + void setNumRegions(size_t numRegions) + { + brineReferenceDensity_.resize(numRegions); + co2ReferenceDensity_.resize(numRegions); + salinity_.resize(numRegions); + } + + + /*! + * \brief Initialize the reference densities of all fluids for a given PVT region + */ + void setReferenceDensities(unsigned regionIdx, + Scalar rhoRefBrine, + Scalar rhoRefCO2, + Scalar /*rhoRefWater*/) + { + brineReferenceDensity_[regionIdx] = rhoRefBrine; + co2ReferenceDensity_[regionIdx] = rhoRefCO2; + } + + + /*! + * \brief Finish initializing the oil phase PVT properties. + */ + void initEnd() + { + + } + + /*! + * \brief Return the number of PVT regions which are considered by this PVT-object. + */ + unsigned numRegions() const + { return brineReferenceDensity_.size(); } + + /*! + * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters. + */ + template + Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED, + const Evaluation& temperature OPM_UNUSED, + const Evaluation& pressure OPM_UNUSED, + const Evaluation& Rs OPM_UNUSED) const + { + throw std::runtime_error("Requested the enthalpy of brine 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& /*Rs*/) const + { + //TODO: The viscosity does not yet depend on the composition + return saturatedViscosity(regionIdx, temperature, pressure); + } + + /*! + * \brief Returns the dynamic viscosity [Pa s] of oil saturated gas at given pressure. + */ + template + Evaluation saturatedViscosity(unsigned /*regionIdx*/, + const Evaluation& temperature, + const Evaluation& pressure) const + { + return Brine::liquidViscosity(temperature, pressure); + } + + /*! + * \brief Returns the formation volume factor [-] of the fluid phase. + */ + template + Evaluation inverseFormationVolumeFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& Rs) const + { + return density_(regionIdx, temperature, pressure, Rs)/brineReferenceDensity_[regionIdx]; + } + + /*! + * \brief Returns the formation volume factor [-] of brine saturated with CO2 at a given pressure. + */ + template + Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const + { + Evaluation rsSat = rsSat_(regionIdx, temperature, pressure); + return density_(regionIdx, temperature, pressure, rsSat)/brineReferenceDensity_[regionIdx]; + } + + /*! + * \brief Returns the saturation pressure of the brine phase [Pa] + * depending on its mass fraction of the gas component + * + * \param Rs + */ + template + Evaluation saturationPressure(unsigned /*regionIdx*/, + const Evaluation& /*temperature*/, + const Evaluation& /*Rs*/) const + { + throw std::runtime_error("Requested the saturation pressure for the brine-co2 pvt module. Not yet implemented."); + } + + /*! + * \brief Returns the gas dissoluiton factor \f$R_s\f$ [m^3/m^3] of the liquid phase. + */ + template + Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& /*oilSaturation*/, + const Evaluation& /*maxOilSaturation*/) const + { + //TODO support VAPPARS + return rsSat_(regionIdx, temperature, pressure); + } + + /*! + * \brief Returns thegas dissoluiton factor \f$R_s\f$ [m^3/m^3] of the liquid phase. + */ + template + Evaluation saturatedGasDissolutionFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const + { + return rsSat_(regionIdx, temperature, pressure); + } + + const Scalar oilReferenceDensity(unsigned regionIdx) const + { return brineReferenceDensity_[regionIdx]; } + + const Scalar gasReferenceDensity(unsigned regionIdx) const + { return co2ReferenceDensity_[regionIdx]; } + + const Scalar salinity(unsigned regionIdx) const + { return salinity_[regionIdx]; } + + bool operator==(const BrineCo2Pvt& data) const + { + return co2ReferenceDensity_ == data.co2ReferenceDensity_ && + brineReferenceDensity_ == data.brineReferenceDensity_; + } + +private: + std::vector brineReferenceDensity_; + std::vector co2ReferenceDensity_; + std::vector salinity_; + + template + LhsEval density_(unsigned regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& Rs) const + { + LhsEval xlCO2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx)); + LhsEval result = liquidDensity_(temperature, + pressure, + xlCO2); + + Valgrind::CheckDefined(result); + return result; + } + + + template + LhsEval liquidDensity_(const LhsEval& T, + const LhsEval& pl, + const LhsEval& xlCO2) const + { + Valgrind::CheckDefined(T); + Valgrind::CheckDefined(pl); + Valgrind::CheckDefined(xlCO2); + + if(T < 273.15) { + std::ostringstream oss; + oss << "Liquid density for Brine and CO2 is only " + "defined above 273.15K (is "<= 2.5e8) { + std::ostringstream oss; + oss << "Liquid density for Brine and CO2 is only " + "defined below 250MPa (is "< + LhsEval liquidDensityWaterCO2_(const LhsEval& temperature, + const LhsEval& pl, + const LhsEval& xlCO2) const + { + Scalar M_CO2 = CO2::molarMass(); + Scalar M_H2O = H2O::molarMass(); + + const LhsEval& tempC = temperature - 273.15; /* tempC : temperature in °C */ + const LhsEval& rho_pure = H2O::liquidDensity(temperature, pl); + // calculate the mole fraction of CO2 in the liquid. note that xlH2O is available + // as a function parameter, but in the case of a pure gas phase the value of M_T + // for the virtual liquid phase can become very large + const LhsEval xlH2O = 1.0 - xlCO2; + const LhsEval& M_T = M_H2O * xlH2O + M_CO2 * xlCO2; + const LhsEval& V_phi = + (37.51 + + tempC*(-9.585e-2 + + tempC*(8.74e-4 - + tempC*5.044e-7))) / 1.0e6; + return 1/ (xlCO2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T)); + } + + /*! + * \brief Convert a gas dissolution factor to the the corresponding mass fraction + * of the gas component in the oil phase. + */ + template + LhsEval convertRsToXoG_(const LhsEval& Rs, unsigned regionIdx) const + { + Scalar rho_oRef = brineReferenceDensity_[regionIdx]; + Scalar rho_gRef = co2ReferenceDensity_[regionIdx]; + + const LhsEval& rho_oG = Rs*rho_gRef; + return rho_oG/(rho_oRef + rho_oG); + } + + + /*! + * \brief Convert a gas mass fraction in the oil phase the corresponding mole fraction. + */ + template + LhsEval convertXoGToxoG_(const LhsEval& XoG) const + { + Scalar M_CO2 = CO2::molarMass(); + Scalar M_H2O = H2O::molarMass(); + return XoG*M_H2O / (M_CO2*(1 - XoG) + XoG*M_H2O); + } + + + /*! + * \brief Convert a gas mole fraction in the oil phase the corresponding mass fraction. + */ + template + LhsEval convertxoGToXoG(const LhsEval& xoG) const + { + Scalar M_CO2 = CO2::molarMass(); + Scalar M_H2O = H2O::molarMass(); + + return xoG*M_CO2 / (xoG*(M_CO2 - M_H2O) + M_H2O); + } + + + /*! + * \brief Convert the mass fraction of the gas component in the oil phase to the + * corresponding gas dissolution factor. + */ + template + LhsEval convertXoGToRs(const LhsEval& XoG, unsigned regionIdx) const + { + Scalar rho_oRef = brineReferenceDensity_[regionIdx]; + Scalar rho_gRef = co2ReferenceDensity_[regionIdx]; + + return XoG/(1.0 - XoG)*(rho_oRef/rho_gRef); + } + + + template + LhsEval rsSat_(unsigned regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + // calulate the equilibrium composition for the given + // temperature and pressure. + LhsEval xgH2O; + LhsEval xlCO2; + BinaryCoeffBrineCO2::calculateMoleFractions(temperature, + pressure, + salinity_[regionIdx], + /*knownPhaseIdx=*/-1, + xlCO2, + xgH2O); + + // normalize the phase compositions + xlCO2 = Opm::max(0.0, Opm::min(1.0, xlCO2)); + + return convertXoGToRs(convertxoGToXoG(xlCO2), regionIdx); + } + +}; + +} // namespace Opm + +#endif diff --git a/opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp b/opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp new file mode 100644 index 000000000..eac2b6826 --- /dev/null +++ b/opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp @@ -0,0 +1,222 @@ +// -*- 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::Co2GasPvt + */ +#ifndef OPM_CO2_GAS_PVT_HPP +#define OPM_Co2_GAS_PVT_HPP + +#include + +#include +#include +#include +#include + +#if HAVE_ECL_INPUT +#include +#include +#include +#endif + +#include + +namespace Opm { +/*! + * \brief This class represents the Pressure-Volume-Temperature relations of the gas phase + * for CO2 + */ +template +class Co2GasPvt +{ + typedef std::vector > SamplingPoints; + typedef Opm::CO2 CO2; + +public: + typedef Opm::Tabulated1DFunction TabulatedOneDFunction; + + explicit Co2GasPvt() = default; + Co2GasPvt(const std::vector& gasReferenceDensity) + : gasReferenceDensity_(gasReferenceDensity) + { + } +#if HAVE_ECL_INPUT + /*! + * \brief Initialize the parameters for co2 gas using an ECL deck. + */ + void initFromState(const EclipseState& eclState, const Schedule&) + { + if( !eclState.getTableManager().getDensityTable().empty()) { + std::cerr << "WARNING: CO2STOR is enabled but DENSITY is in the deck. \n" << + "The surface density is computed based on CO2-BRINE PVT at standard conditions (STCOND) and DENSITY is ignored " << std::endl; + } + + if( eclState.getTableManager().hasTables("PVDG") || !eclState.getTableManager().getPvtgTables().empty()) { + std::cerr << "WARNING: CO2STOR is enabled but PVDG or PVTG is in the deck. \n" << + "CO2 PVT properties are computed based on the Span-Wagner pvt model and PVDG/PVTG input is ignored. " << std::endl; + } + + // We only supported single pvt region for the co2-brine module + size_t numRegions = 1; + setNumRegions(numRegions); + size_t regionIdx = 0; + Scalar T_ref = eclState.getTableManager().stCond().temperature; + Scalar P_ref = eclState.getTableManager().stCond().pressure; + gasReferenceDensity_[regionIdx] = CO2::gasDensity(T_ref, P_ref); + initEnd(); + } +#endif + + void setNumRegions(size_t numRegions) + { + gasReferenceDensity_.resize(numRegions); + } + + + /*! + * \brief Initialize the reference densities of all fluids for a given PVT region + */ + void setReferenceDensities(unsigned regionIdx, + Scalar /*rhoRefOil*/, + Scalar rhoRefGas, + Scalar /*rhoRefWater*/) + { + gasReferenceDensity_[regionIdx] = rhoRefGas; + } + + /*! + * \brief Finish initializing the oil phase PVT properties. + */ + void initEnd() + { + + } + + /*! + * \brief Return the number of PVT regions which are considered by this PVT-object. + */ + unsigned numRegions() const + { return gasReferenceDensity_.size(); } + + /*! + * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters. + */ + template + Evaluation internalEnergy(unsigned regionIdx OPM_UNUSED, + const Evaluation& temperature OPM_UNUSED, + const Evaluation& pressure OPM_UNUSED, + const Evaluation& Rv OPM_UNUSED) const + { + throw std::runtime_error("Requested the enthalpy of gas 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& /*Rv*/) const + { return saturatedViscosity(regionIdx, temperature, pressure); } + + /*! + * \brief Returns the dynamic viscosity [Pa s] of oil saturated gas at given pressure. + */ + template + Evaluation saturatedViscosity(unsigned /*regionIdx*/, + const Evaluation& temperature, + const Evaluation& pressure) const + { + return CO2::gasViscosity(temperature, pressure); + } + + /*! + * \brief Returns the formation volume factor [-] of the fluid phase. + */ + template + Evaluation inverseFormationVolumeFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& /*Rv*/) const + { return saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure); } + + /*! + * \brief Returns the formation volume factor [-] of oil saturated gas at given pressure. + */ + template + Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const + { + return CO2::gasDensity(temperature, pressure)/gasReferenceDensity_[regionIdx]; + } + + /*! + * \brief Returns the saturation pressure of the gas phase [Pa] + * depending on its mass fraction of the oil component + * + * \param Rv The surface volume of oil component dissolved in what will yield one cubic meter of gas at the surface [-] + */ + template + Evaluation saturationPressure(unsigned /*regionIdx*/, + const Evaluation& /*temperature*/, + const Evaluation& /*Rv*/) const + { return 0.0; /* this is dry gas! */ } + + /*! + * \brief Returns the oil vaporization factor \f$R_v\f$ [m^3/m^3] of the oil phase. + */ + template + Evaluation saturatedOilVaporizationFactor(unsigned /*regionIdx*/, + const Evaluation& /*temperature*/, + const Evaluation& /*pressure*/, + const Evaluation& /*oilSaturation*/, + const Evaluation& /*maxOilSaturation*/) const + { return 0.0; /* this is dry gas! */ } + + /*! + * \brief Returns the oil vaporization factor \f$R_v\f$ [m^3/m^3] of the oil phase. + */ + template + Evaluation saturatedOilVaporizationFactor(unsigned /*regionIdx*/, + const Evaluation& /*temperature*/, + const Evaluation& /*pressure*/) const + { return 0.0; /* this is dry gas! */ } + + const Scalar gasReferenceDensity(unsigned regionIdx) const + { return gasReferenceDensity_[regionIdx]; } + + bool operator==(const Co2GasPvt& data) const + { + return gasReferenceDensity_ == data.gasReferenceDensity_; + } + +private: + std::vector gasReferenceDensity_; +}; + +} // namespace Opm + +#endif diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp index dfc153df9..49fa6e9cd 100644 --- a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp @@ -202,8 +202,8 @@ public: } - const std::vector& waterReferenceDensity() const - { return waterReferenceDensity_; } + const Scalar waterReferenceDensity(unsigned regionIdx) const + { return waterReferenceDensity_[regionIdx]; } const std::vector& referencePressure() const { return referencePressure_; } diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp index 4b256cb81..470c51160 100644 --- a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp @@ -272,11 +272,8 @@ public: const Evaluation& /*Rs*/) const { return 0.0; /* this is dead oil, so there isn't any meaningful saturation pressure! */ } - const std::vector& oilReferenceDensity() const - { return oilReferenceDensity_; } - - const std::vector& oilReferencePressure() const - { return oilReferencePressure_; } + const Scalar oilReferenceDensity(unsigned regionIdx) const + { return oilReferenceDensity_[regionIdx]; } const std::vector& oilReferenceFormationVolumeFactor() const { return oilReferenceFormationVolumeFactor_; } diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp index 5ce3e9668..f14175ef0 100644 --- a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp @@ -210,8 +210,8 @@ public: return (1.0 + X*(1.0 + X/2.0))/BwRef; } - const std::vector& waterReferenceDensity() const - { return waterReferenceDensity_; } + const Scalar waterReferenceDensity(unsigned regionIdx) const + { return waterReferenceDensity_[regionIdx]; } const std::vector& waterReferencePressure() const { return waterReferencePressure_; } diff --git a/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp b/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp index 5d7057ee0..513b7f23d 100644 --- a/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp @@ -278,8 +278,8 @@ public: const Evaluation& /*pressure*/) const { return 0.0; /* this is dead oil! */ } - const std::vector& oilReferenceDensity() const - { return oilReferenceDensity_; } + const Scalar oilReferenceDensity(unsigned regionIdx) const + { return oilReferenceDensity_[regionIdx]; } const std::vector& inverseOilB() const { return inverseOilB_; } diff --git a/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp b/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp index 4160a1b22..acfba2683 100644 --- a/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp @@ -288,8 +288,8 @@ public: const Evaluation& /*pressure*/) const { return 0.0; /* this is dry gas! */ } - const std::vector& gasReferenceDensity() const - { return gasReferenceDensity_; } + const Scalar gasReferenceDensity(unsigned regionIdx) const + { return gasReferenceDensity_[regionIdx]; } const std::vector& inverseGasB() const { return inverseGasB_; } diff --git a/opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp b/opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp index 16079ead5..845c18674 100644 --- a/opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp +++ b/opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp @@ -30,6 +30,7 @@ #include "DryGasPvt.hpp" #include "WetGasPvt.hpp" #include "GasPvtThermal.hpp" +#include "Co2GasPvt.hpp" #if HAVE_ECL_INPUT #include @@ -53,9 +54,14 @@ namespace Opm { codeToCall; \ break; \ } \ + case Co2GasPvt: { \ + auto& pvtImpl = getRealPvt(); \ + codeToCall; \ + break; \ + } \ case NoGasPvt: \ throw std::logic_error("Not implemented: Gas PVT of this deck!"); \ - } + } \ /*! @@ -78,7 +84,8 @@ public: NoGasPvt, DryGasPvt, WetGasPvt, - ThermalGasPvt + ThermalGasPvt, + Co2GasPvt }; GasPvtMultiplexer() @@ -112,6 +119,10 @@ public: delete &getRealPvt(); break; } + case Co2GasPvt: { + delete &getRealPvt(); + break; + } case NoGasPvt: break; } @@ -127,8 +138,9 @@ public: { if (!eclState.runspec().phases().active(Phase::GAS)) return; - - if (enableThermal && eclState.getSimulationConfig().isThermal()) + if (eclState.runspec().co2Storage()) + setApproach(Co2GasPvt); + else if (enableThermal && eclState.getSimulationConfig().isThermal()) setApproach(ThermalGasPvt); else if (!eclState.getTableManager().getPvtgTables().empty()) setApproach(WetGasPvt); @@ -154,6 +166,10 @@ public: realGasPvt_ = new Opm::GasPvtThermal; break; + case Co2GasPvt: + realGasPvt_ = new Opm::Co2GasPvt; + break; + case NoGasPvt: throw std::logic_error("Not implemented: Gas PVT of this deck!"); } @@ -170,6 +186,12 @@ public: unsigned numRegions() const { OPM_GAS_PVT_MULTIPLEXER_CALL(return pvtImpl.numRegions()); return 1; } + /*! + * \brief Return the reference density which are considered by this PVT-object. + */ + const Scalar gasReferenceDensity(unsigned regionIdx) + { OPM_GAS_PVT_MULTIPLEXER_CALL(return pvtImpl.gasReferenceDensity(regionIdx)); return 2.; } + /*! * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters. */ @@ -303,6 +325,20 @@ public: return *static_cast* >(realGasPvt_); } + template + typename std::enable_if >::type& getRealPvt() + { + assert(gasPvtApproach() == approachV); + return *static_cast* >(realGasPvt_); + } + + template + typename std::enable_if >::type& getRealPvt() const + { + assert(gasPvtApproach() == approachV); + return *static_cast* >(realGasPvt_); + } + const void* realGasPvt() const { return realGasPvt_; } bool operator==(const GasPvtMultiplexer& data) const @@ -320,6 +356,9 @@ public: case ThermalGasPvt: return *static_cast*>(realGasPvt_) == *static_cast*>(data.realGasPvt_); + case Co2GasPvt: + return *static_cast*>(realGasPvt_) == + *static_cast*>(data.realGasPvt_); default: return true; } @@ -338,6 +377,9 @@ public: case ThermalGasPvt: realGasPvt_ = new Opm::GasPvtThermal(*static_cast*>(data.realGasPvt_)); break; + case Co2GasPvt: + realGasPvt_ = new Opm::Co2GasPvt(*static_cast*>(data.realGasPvt_)); + break; default: break; } diff --git a/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp index f88b88f23..01c98eb9c 100644 --- a/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp +++ b/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp @@ -370,6 +370,9 @@ public: const IsothermalPvt* isoThermalPvt() const { return isothermalPvt_; } + const Scalar gasReferenceDensity(unsigned regionIdx) const + { return isothermalPvt_->gasReferenceDensity(regionIdx); } + const std::vector& gasvisctCurves() const { return gasvisctCurves_; } diff --git a/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp b/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp index a987b93e0..fd51a89eb 100644 --- a/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp @@ -102,7 +102,7 @@ public: for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { Scalar rhoRefO = densityTable[regionIdx].oil; - Scalar rhoRefG = densityTable[regionIdx].gas;; + Scalar rhoRefG = densityTable[regionIdx].gas; Scalar rhoRefW = densityTable[regionIdx].water; setReferenceDensities(regionIdx, rhoRefO, rhoRefG, rhoRefW); @@ -602,11 +602,11 @@ public: throw NumericalIssue(errlog.str()); } - const std::vector& gasReferenceDensity() const - { return gasReferenceDensity_; } + const Scalar gasReferenceDensity(unsigned regionIdx) const + { return gasReferenceDensity_[regionIdx]; } - const std::vector& oilReferenceDensity() const - { return oilReferenceDensity_; } + const Scalar oilReferenceDensity(unsigned regionIdx) const + { return oilReferenceDensity_[regionIdx]; } const std::vector& inverseOilBTable() const { return inverseOilBTable_; } diff --git a/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp b/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp index 30bbf0d70..c68f64352 100644 --- a/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp +++ b/opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp @@ -31,6 +31,7 @@ #include "DeadOilPvt.hpp" #include "LiveOilPvt.hpp" #include "OilPvtThermal.hpp" +#include "BrineCo2Pvt.hpp" #if HAVE_ECL_INPUT #include @@ -60,9 +61,14 @@ namespace Opm { codeToCall; \ break; \ } \ + case BrineCo2Pvt: { \ + auto& pvtImpl = getRealPvt(); \ + codeToCall; \ + break; \ + } \ case NoOilPvt: \ throw std::logic_error("Not implemented: Oil PVT of this deck!"); \ - } + } \ /*! * \brief This class represents the Pressure-Volume-Temperature relations of the oil @@ -87,7 +93,8 @@ public: LiveOilPvt, DeadOilPvt, ConstantCompressibilityOilPvt, - ThermalOilPvt + ThermalOilPvt, + BrineCo2Pvt }; OilPvtMultiplexer() @@ -125,6 +132,10 @@ public: delete &getRealPvt(); break; } + case BrineCo2Pvt: { + delete &getRealPvt(); + break; + } case NoOilPvt: break; @@ -141,8 +152,11 @@ public: { if (!eclState.runspec().phases().active(Phase::OIL)) return; - - if (enableThermal && eclState.getSimulationConfig().isThermal()) + // TODO move the BrineCo2 approach to the waterPvtMultiplexer + // when a proper gas-water simulator is supported + if (eclState.runspec().co2Storage()) + setApproach(BrineCo2Pvt); + else if (enableThermal && eclState.getSimulationConfig().isThermal()) setApproach(ThermalOilPvt); else if (!eclState.getTableManager().getPvcdoTable().empty()) setApproach(ConstantCompressibilityOilPvt); @@ -165,6 +179,12 @@ public: unsigned numRegions() const { OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.numRegions()); return 1; } + /*! + * \brief Return the reference density which are considered by this PVT-object. + */ + const Scalar oilReferenceDensity(unsigned regionIdx) + { OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.oilReferenceDensity(regionIdx)); return 700.; } + /*! * \brief Returns the specific enthalpy [J/kg] oil given a set of parameters. */ @@ -265,6 +285,10 @@ public: realOilPvt_ = new Opm::OilPvtThermal; break; + case BrineCo2Pvt: + realOilPvt_ = new Opm::BrineCo2Pvt; + break; + case NoOilPvt: throw std::logic_error("Not implemented: Oil PVT of this deck!"); } @@ -337,6 +361,20 @@ public: return *static_cast* >(realOilPvt_); } + template + typename std::enable_if >::type& getRealPvt() + { + assert(approach() == approachV); + return *static_cast* >(realOilPvt_); + } + + template + typename std::enable_if >::type& getRealPvt() const + { + assert(approach() == approachV); + return *static_cast* >(realOilPvt_); + } + const void* realOilPvt() const { return realOilPvt_; } bool operator==(const OilPvtMultiplexer& data) const @@ -357,6 +395,9 @@ public: case ThermalOilPvt: return *static_cast*>(realOilPvt_) == *static_cast*>(data.realOilPvt_); + case BrineCo2Pvt: + return *static_cast*>(realOilPvt_) == + *static_cast*>(data.realOilPvt_); default: return true; } @@ -378,6 +419,9 @@ public: case ThermalOilPvt: realOilPvt_ = new Opm::OilPvtThermal(*static_cast*>(data.realOilPvt_)); break; + case BrineCo2Pvt: + realOilPvt_ = new Opm::BrineCo2Pvt(*static_cast*>(data.realOilPvt_)); + break; default: break; } diff --git a/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp index 7f6e04d8a..8abf2de11 100644 --- a/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp +++ b/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp @@ -380,6 +380,9 @@ public: const IsothermalPvt* isoThermalPvt() const { return isothermalPvt_; } + const Scalar oilReferenceDensity(unsigned regionIdx) const + { return isothermalPvt_->oilReferenceDensity(regionIdx); } + const std::vector& oilvisctCurves() const { return oilvisctCurves_; } diff --git a/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp b/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp index d1b4b248e..addf68bb5 100644 --- a/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp @@ -142,6 +142,12 @@ public: unsigned numRegions() const { OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.numRegions()); return 1; } + /*! + * \brief Return the reference density which are considered by this PVT-object. + */ + const Scalar waterReferenceDensity(unsigned regionIdx) + { OPM_WATER_PVT_MULTIPLEXER_CALL(return pvtImpl.waterReferenceDensity(regionIdx)); return 1000.; } + /*! * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters. */ diff --git a/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp b/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp index 3c87ff6a7..68ca91d57 100644 --- a/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp @@ -305,6 +305,9 @@ public: const IsothermalPvt* isoThermalPvt() const { return isothermalPvt_; } + const Scalar waterReferenceDensity(unsigned regionIdx) const + { return isothermalPvt_->waterReferenceDensity(regionIdx); } + const std::vector& viscrefPress() const { return viscrefPress_; } diff --git a/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp b/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp index 5bda3868c..51d6e6bb4 100644 --- a/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp @@ -628,13 +628,11 @@ public: throw NumericalIssue(errlog.str()); } - const std::vector& gasReferenceDensity() const { - return gasReferenceDensity_; - } + const Scalar gasReferenceDensity(unsigned regionIdx) const + { return gasReferenceDensity_[regionIdx]; } - const std::vector& oilReferenceDensity() const { - return oilReferenceDensity_; - } + const Scalar oilReferenceDensity(unsigned regionIdx) const + { return oilReferenceDensity_[regionIdx]; } const std::vector& inverseGasB() const { return inverseGasB_; diff --git a/tests/test_co2brinepvt.cpp b/tests/test_co2brinepvt.cpp new file mode 100644 index 000000000..d0eef52d1 --- /dev/null +++ b/tests/test_co2brinepvt.cpp @@ -0,0 +1,196 @@ +// -*- 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 + * + * \brief This is the unit test for the co2 brine PVT model + * + * This test requires the presence of opm-common. + */ +#include "config.h" + +#if !HAVE_ECL_INPUT +#error "The test for the co2 brine PVT classes requires eclipse input support in opm-common" +#endif + +//#include +//#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +// values of strings based on the first SPE1 test case of opm-data. note that in the +// real world it does not make much sense to specify a fluid phase using more than a +// single keyword, but for a unit test, this saves a lot of boiler-plate code. +static const char* deckString1 = + "RUNSPEC\n" + "\n" + "DIMENS\n" + " 10 10 3 /\n" + "\n" + "TABDIMS\n" + " * 1 /\n" + "\n" + "OIL\n" + "GAS\n" + "CO2STOR\n" + "\n" + "DISGAS\n" + "\n" + "METRIC\n" + "\n" + "GRID\n" + "\n" + "DX\n" + " 300*1000 /\n" + "DY\n" + " 300*1000 /\n" + "DZ\n" + " 100*20 100*30 100*50 /\n" + "\n" + "TOPS\n" + " 100*1234 /\n" + "\n" + "PORO\n" + " 300*0.15 /\n" + "PROPS\n" + "\n"; + +template +void ensurePvtApi(const BrinePvt& brinePvt, const Co2Pvt& co2Pvt) +{ + // we don't want to run this, we just want to make sure that it compiles + while (0) { + Evaluation temperature = 273.15 + 20.0; + Evaluation pressure = 1e5; + Evaluation Rs = 0.0; + Evaluation Rv = 0.0; + Evaluation So = 0.5; + Evaluation maxSo = 1.0; + Evaluation tmp; + + ///// + // brine PVT API + ///// + tmp = brinePvt.viscosity(/*regionIdx=*/0, + temperature, + pressure, + Rs); + tmp = brinePvt.inverseFormationVolumeFactor(/*regionIdx=*/0, + temperature, + pressure, + Rs); + tmp = brinePvt.saturatedViscosity(/*regionIdx=*/0, + temperature, + pressure); + tmp = brinePvt.saturatedInverseFormationVolumeFactor(/*regionIdx=*/0, + temperature, + pressure); + tmp = brinePvt.saturationPressure(/*regionIdx=*/0, + temperature, + Rs); + tmp = brinePvt.saturatedGasDissolutionFactor(/*regionIdx=*/0, + temperature, + pressure); + tmp = brinePvt.saturatedGasDissolutionFactor(/*regionIdx=*/0, + temperature, + pressure, + So, + maxSo); + + ///// + // co2 PVT API + ///// + tmp = co2Pvt.viscosity(/*regionIdx=*/0, + temperature, + pressure, + Rv); + tmp = co2Pvt.inverseFormationVolumeFactor(/*regionIdx=*/0, + temperature, + pressure, + Rv); + tmp = co2Pvt.saturatedViscosity(/*regionIdx=*/0, + temperature, + pressure); + tmp = co2Pvt.saturatedInverseFormationVolumeFactor(/*regionIdx=*/0, + temperature, + pressure); + tmp = co2Pvt.saturationPressure(/*regionIdx=*/0, + temperature, + Rv); + tmp = co2Pvt.saturatedOilVaporizationFactor(/*regionIdx=*/0, + temperature, + pressure); + tmp = co2Pvt.saturatedOilVaporizationFactor(/*regionIdx=*/0, + temperature, + pressure, + So, + maxSo); + + // prevent GCC from producing a "variable assigned but unused" warning + tmp = 2.0*tmp; + } +} + +template +inline void testAll() +{ + Opm::Parser parser; + auto python = std::make_shared(); + + auto deck = parser.parseString(deckString1); + Opm::EclipseState eclState(deck); + Opm::Schedule schedule(deck, eclState, python); + + Opm::GasPvtMultiplexer co2Pvt; + Opm::OilPvtMultiplexer brinePvt; + + co2Pvt.initFromState(eclState, schedule); + brinePvt.initFromState(eclState, schedule); + + typedef Opm::DenseAd::Evaluation FooEval; + ensurePvtApi(brinePvt, co2Pvt); + ensurePvtApi(brinePvt, co2Pvt); +} + + +int main(int argc, char **argv) +{ + Dune::MPIHelper::instance(argc, argv); + + testAll(); + testAll(); + + return 0; +} diff --git a/tests/test_components.cpp b/tests/test_components.cpp index 2e198eaa4..7efa8263d 100644 --- a/tests/test_components.cpp +++ b/tests/test_components.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -62,6 +63,36 @@ namespace ComponentsTest { #include +template +void testSimpleH2O() +{ + typedef Opm::H2O H2O; + typedef Opm::SimpleHuDuanH2O SimpleHuDuanH2O; + typedef Opm::MathToolbox EvalToolbox; + + int numT = 67; + int numP = 45; + Evaluation T = 280; + Evaluation p = 1e6; + + for (int iT = 0; iT < numT; ++iT) { + p = 1e6; + T += 5; + for (int iP = 0; iP < numP; ++iP) { + p *= 1.1; + if (!EvalToolbox::isSame(H2O::liquidDensity(T,p), SimpleHuDuanH2O::liquidDensity(T,p), /*tolerance=*/1e-3*H2O::liquidDensity(T,p).value())) + throw std::logic_error("oops: the water density based on Hu-Duan has more then 1e-3 deviation from IAPWS'97"); + + if (T >= 570) // for temperature larger then 570 the viscosity based on HuDuan is too far from IAPWS. + continue; + + if (!EvalToolbox::isSame(H2O::liquidViscosity(T,p), SimpleHuDuanH2O::liquidViscosity(T,p), /*tolerance=*/5.e-2*H2O::liquidViscosity(T,p).value())){ + throw std::logic_error("oops: the water viscosity based on Hu-Duan has more then 5e-2 deviation from IAPWS'97"); + } + } + } +} + template void testAllComponents() { @@ -91,6 +122,8 @@ inline void testAll() // ensure that all components are API-compliant testAllComponents(); testAllComponents(); + testSimpleH2O(); + } diff --git a/tests/test_fluidsystems.cpp b/tests/test_fluidsystems.cpp index ca0d1e7cb..5cb0a3c9b 100644 --- a/tests/test_fluidsystems.cpp +++ b/tests/test_fluidsystems.cpp @@ -60,12 +60,6 @@ #include - -namespace Opm { -namespace FluidSystemsTest { -#include -} } - #include // check that the blackoil fluid system implements all non-standard functions @@ -248,7 +242,7 @@ void testAllFluidSystems() } // Brine -- CO2 - { typedef Opm::BrineCO2FluidSystem FluidSystem; + { typedef Opm::BrineCO2FluidSystem FluidSystem; checkFluidSystem(); } // H2O -- N2