OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #427 from totto82/co2pvt

Browse Source 
Add a co2brine pvt model
This commit is contained in:
Tor Harald Sandve 2020-10-02 12:26:19 +02:00 committed by GitHub
commit 28dd5a2b36
23 changed files with 1445 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
opm/material/binarycoefficients/Brine_CO2.hpp
View File

@ -28,9 +28,6 @@
#ifndef OPM_BINARY_COEFF_BRINE_CO2_HPP
#define OPM_BINARY_COEFF_BRINE_CO2_HPP
#include <opm/material/components/Brine.hpp>
#include <opm/material/components/H2O.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/IdealGas.hpp>
namespace Opm {
@ -40,10 +37,8 @@ namespace BinaryCoeff {
* \ingroup Binarycoefficients
* \brief Binary coefficients for brine and CO2.
*/
template<class Scalar, class CO2Tables, bool verbose = true>
template<class Scalar, class H2O, class CO2, bool verbose = true>
class Brine_CO2 {
typedef Opm::H2O<Scalar> H2O;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
typedef Opm::IdealGas<Scalar> IdealGas;
static const int liquidPhaseIdx = 0; // index of the liquid phase
static const int gasPhaseIdx = 1; // index of the gas phase

399
opm/material/components/SimpleHuDuanH2O.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <opm/material/IdealGas.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Unused.hpp>
#include <cmath>
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 Scalar>
class SimpleHuDuanH2O : public Component<Scalar, SimpleHuDuanH2O<Scalar> >
{
typedef Opm::IdealGas<Scalar> IdealGas;
typedef IAPWS::Common<Scalar> 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 <class Evaluation>
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 <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& /*pressure*/)
{ return 1.976e3*temperature + 40.65e3/molarMass(); }
/*!
* \copydoc Component::gasHeatCapacity
*/
template <class Evaluation>
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 <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature,
const Evaluation& /*pressure*/)
{ return 4180*temperature; }
/*!
* \copydoc Component::liquidHeatCapacity
*/
template <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Scalar>
const Scalar SimpleHuDuanH2O<Scalar>::R = Opm::Constants<Scalar>::R / 18e-3;
} // namespace Opm
#endif

39
opm/material/components/co2tables.inc
View File

@ -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 */

20
opm/material/fluidsystems/BlackOilFluidSystem.hpp
View File

@ -31,6 +31,7 @@
#include "blackoilpvt/OilPvtMultiplexer.hpp"
#include "blackoilpvt/GasPvtMultiplexer.hpp"
#include "blackoilpvt/WaterPvtMultiplexer.hpp"
#include "blackoilpvt/BrineCo2Pvt.hpp"
#include <opm/material/fluidsystems/BaseFluidSystem.hpp>
#include <opm/material/Constants.hpp>
@ -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>();
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

12
opm/material/fluidsystems/BrineCO2FluidSystem.hpp
View File

@ -72,8 +72,13 @@ public:
struct ParameterCache : public Opm::NullParameterCache<Evaluation>
{};
//! 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<Scalar, CO2Tables> CO2;
//! The binary coefficients for brine and CO2 used by this fluid system
typedef Opm::BinaryCoeff::Brine_CO2<Scalar, CO2Tables> BinaryCoeffBrineCO2;
typedef Opm::BinaryCoeff::Brine_CO2<Scalar, H2O, 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<Scalar, CO2Tables> CO2;
/*!
* \copydoc BaseFluidSystem::componentName
*/

427
opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <opm/material/Constants.hpp>
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/components/Brine.hpp>
#include <opm/material/components/SimpleHuDuanH2O.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>
#include <opm/material/components/TabulatedComponent.hpp>
#include <opm/material/binarycoefficients/H2O_CO2.hpp>
#include <opm/material/binarycoefficients/Brine_CO2.hpp>
#include <opm/material/components/co2tables.inc>
#if HAVE_ECL_INPUT
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <vector>
namespace Opm {
/*!
* \brief This class represents the Pressure-Volume-Temperature relations of the liquid phase
* for a CO2-Brine system
*/
template <class Scalar>
class BrineCo2Pvt
{
typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
typedef Opm::SimpleHuDuanH2O<Scalar> H2O;
typedef Opm::Brine<Scalar, H2O> Brine;
//typedef Opm::H2O<Scalar> H2O_IAPWS;
//typedef Opm::Brine<Scalar, H2O_IAPWS> Brine_IAPWS;
//typedef Opm::TabulatedComponent<Scalar, H2O_IAPWS> H2O_Tabulated;
//typedef Opm::TabulatedComponent<Scalar, Brine_IAPWS> Brine_Tabulated;
//typedef H2O_Tabulated H2O;
//typedef Brine_Tabulated Brine;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
public:
typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
//! The binary coefficients for brine and CO2 used by this fluid system
typedef Opm::BinaryCoeff::Brine_CO2<Scalar, H2O, CO2> BinaryCoeffBrineCO2;
explicit BrineCo2Pvt() = default;
BrineCo2Pvt(const std::vector<Scalar>& brineReferenceDensity,
const std::vector<Scalar>& co2ReferenceDensity,
const std::vector<Scalar>& 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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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<Scalar>& data) const
{
return co2ReferenceDensity_ == data.co2ReferenceDensity_ &&
brineReferenceDensity_ == data.brineReferenceDensity_;
}
private:
std::vector<Scalar> brineReferenceDensity_;
std::vector<Scalar> co2ReferenceDensity_;
std::vector<Scalar> salinity_;
template <class LhsEval>
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 <class LhsEval>
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 "<<T<<"K)";
throw NumericalIssue(oss.str());
}
if(pl >= 2.5e8) {
std::ostringstream oss;
oss << "Liquid density for Brine and CO2 is only "
"defined below 250MPa (is "<<pl<<"Pa)";
throw NumericalIssue(oss.str());
}
const LhsEval& rho_brine = Brine::liquidDensity(T, pl);
const LhsEval& rho_pure = H2O::liquidDensity(T, pl);
const LhsEval& rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlCO2);
const LhsEval& contribCO2 = rho_lCO2 - rho_pure;
return rho_brine + contribCO2;
}
template <class LhsEval>
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 <class LhsEval>
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 <class LhsEval>
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 <class LhsEval>
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 <class LhsEval>
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 <class LhsEval>
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

222
opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <opm/material/Constants.hpp>
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>
#include <opm/material/components/co2tables.inc>
#if HAVE_ECL_INPUT
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
#include <vector>
namespace Opm {
/*!
* \brief This class represents the Pressure-Volume-Temperature relations of the gas phase
* for CO2
*/
template <class Scalar>
class Co2GasPvt
{
typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
public:
typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
explicit Co2GasPvt() = default;
Co2GasPvt(const std::vector<Scalar>& 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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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 <class Evaluation>
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<Scalar>& data) const
{
return gasReferenceDensity_ == data.gasReferenceDensity_;
}
private:
std::vector<Scalar> gasReferenceDensity_;
};
} // namespace Opm
#endif

4
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp
View File

@ -202,8 +202,8 @@ public:
}
const std::vector<Scalar>& waterReferenceDensity() const
{ return waterReferenceDensity_; }
const Scalar waterReferenceDensity(unsigned regionIdx) const
{ return waterReferenceDensity_[regionIdx]; }
const std::vector<Scalar>& referencePressure() const
{ return referencePressure_; }

7
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp
View File

@ -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<Scalar>& oilReferenceDensity() const
{ return oilReferenceDensity_; }
const std::vector<Scalar>& oilReferencePressure() const
{ return oilReferencePressure_; }
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return oilReferenceDensity_[regionIdx]; }
const std::vector<Scalar>& oilReferenceFormationVolumeFactor() const
{ return oilReferenceFormationVolumeFactor_; }

4
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp
View File

@ -210,8 +210,8 @@ public:
return (1.0 + X*(1.0 + X/2.0))/BwRef;
}
const std::vector<Scalar>& waterReferenceDensity() const
{ return waterReferenceDensity_; }
const Scalar waterReferenceDensity(unsigned regionIdx) const
{ return waterReferenceDensity_[regionIdx]; }
const std::vector<Scalar>& waterReferencePressure() const
{ return waterReferencePressure_; }

4
opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp
View File

@ -278,8 +278,8 @@ public:
const Evaluation& /*pressure*/) const
{ return 0.0; /* this is dead oil! */ }
const std::vector<Scalar>& oilReferenceDensity() const
{ return oilReferenceDensity_; }
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return oilReferenceDensity_[regionIdx]; }
const std::vector<TabulatedOneDFunction>& inverseOilB() const
{ return inverseOilB_; }

4
opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp
View File

@ -288,8 +288,8 @@ public:
const Evaluation& /*pressure*/) const
{ return 0.0; /* this is dry gas! */ }
const std::vector<Scalar>& gasReferenceDensity() const
{ return gasReferenceDensity_; }
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }
const std::vector<TabulatedOneDFunction>& inverseGasB() const
{ return inverseGasB_; }

50
opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp
View File

@ -30,6 +30,7 @@
#include "DryGasPvt.hpp"
#include "WetGasPvt.hpp"
#include "GasPvtThermal.hpp"
#include "Co2GasPvt.hpp"
#if HAVE_ECL_INPUT
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@ -53,9 +54,14 @@ namespace Opm {
codeToCall; \
break; \
} \
case Co2GasPvt: { \
auto& pvtImpl = getRealPvt<Co2GasPvt>(); \
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<ThermalGasPvt>();
break;
}
case Co2GasPvt: {
delete &getRealPvt<Co2GasPvt>();
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<Scalar>;
break;
case Co2GasPvt:
realGasPvt_ = new Opm::Co2GasPvt<Scalar>;
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<const Opm::GasPvtThermal<Scalar>* >(realGasPvt_);
}
template <GasPvtApproach approachV>
typename std::enable_if<approachV == Co2GasPvt, Opm::Co2GasPvt<Scalar> >::type& getRealPvt()
{
assert(gasPvtApproach() == approachV);
return *static_cast<Opm::Co2GasPvt<Scalar>* >(realGasPvt_);
}
template <GasPvtApproach approachV>
typename std::enable_if<approachV == Co2GasPvt, const Opm::Co2GasPvt<Scalar> >::type& getRealPvt() const
{
assert(gasPvtApproach() == approachV);
return *static_cast<const Opm::Co2GasPvt<Scalar>* >(realGasPvt_);
}
const void* realGasPvt() const { return realGasPvt_; }
bool operator==(const GasPvtMultiplexer<Scalar,enableThermal>& data) const
@ -320,6 +356,9 @@ public:
case ThermalGasPvt:
return *static_cast<const Opm::GasPvtThermal<Scalar>*>(realGasPvt_) ==
*static_cast<const Opm::GasPvtThermal<Scalar>*>(data.realGasPvt_);
case Co2GasPvt:
return *static_cast<const Opm::Co2GasPvt<Scalar>*>(realGasPvt_) ==
*static_cast<const Opm::Co2GasPvt<Scalar>*>(data.realGasPvt_);
default:
return true;
}
@ -338,6 +377,9 @@ public:
case ThermalGasPvt:
realGasPvt_ = new Opm::GasPvtThermal<Scalar>(*static_cast<const Opm::GasPvtThermal<Scalar>*>(data.realGasPvt_));
break;
case Co2GasPvt:
realGasPvt_ = new Opm::Co2GasPvt<Scalar>(*static_cast<const Opm::Co2GasPvt<Scalar>*>(data.realGasPvt_));
break;
default:
break;
}

3
opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
View File

@ -370,6 +370,9 @@ public:
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return isothermalPvt_->gasReferenceDensity(regionIdx); }
const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
{ return gasvisctCurves_; }

10
opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp
View File

@ -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<Scalar>& gasReferenceDensity() const
{ return gasReferenceDensity_; }
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }
const std::vector<Scalar>& oilReferenceDensity() const
{ return oilReferenceDensity_; }
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return oilReferenceDensity_[regionIdx]; }
const std::vector<TabulatedTwoDFunction>& inverseOilBTable() const
{ return inverseOilBTable_; }

52
opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp
View File

@ -31,6 +31,7 @@
#include "DeadOilPvt.hpp"
#include "LiveOilPvt.hpp"
#include "OilPvtThermal.hpp"
#include "BrineCo2Pvt.hpp"
#if HAVE_ECL_INPUT
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@ -60,9 +61,14 @@ namespace Opm {
codeToCall; \
break; \
} \
case BrineCo2Pvt: { \
auto& pvtImpl = getRealPvt<BrineCo2Pvt>(); \
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<ThermalOilPvt>();
break;
}
case BrineCo2Pvt: {
delete &getRealPvt<BrineCo2Pvt>();
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<Scalar>;
break;
case BrineCo2Pvt:
realOilPvt_ = new Opm::BrineCo2Pvt<Scalar>;
break;
case NoOilPvt:
throw std::logic_error("Not implemented: Oil PVT of this deck!");
}
@ -337,6 +361,20 @@ public:
return *static_cast<const Opm::OilPvtThermal<Scalar>* >(realOilPvt_);
}
template <OilPvtApproach approachV>
typename std::enable_if<approachV == BrineCo2Pvt, Opm::BrineCo2Pvt<Scalar> >::type& getRealPvt()
{
assert(approach() == approachV);
return *static_cast<Opm::BrineCo2Pvt<Scalar>* >(realOilPvt_);
}
template <OilPvtApproach approachV>
typename std::enable_if<approachV == BrineCo2Pvt, const Opm::BrineCo2Pvt<Scalar> >::type& getRealPvt() const
{
assert(approach() == approachV);
return *static_cast<const Opm::BrineCo2Pvt<Scalar>* >(realOilPvt_);
}
const void* realOilPvt() const { return realOilPvt_; }
bool operator==(const OilPvtMultiplexer<Scalar,enableThermal>& data) const
@ -357,6 +395,9 @@ public:
case ThermalOilPvt:
return *static_cast<const Opm::OilPvtThermal<Scalar>*>(realOilPvt_) ==
*static_cast<const Opm::OilPvtThermal<Scalar>*>(data.realOilPvt_);
case BrineCo2Pvt:
return *static_cast<const Opm::BrineCo2Pvt<Scalar>*>(realOilPvt_) ==
*static_cast<const Opm::BrineCo2Pvt<Scalar>*>(data.realOilPvt_);
default:
return true;
}
@ -378,6 +419,9 @@ public:
case ThermalOilPvt:
realOilPvt_ = new Opm::OilPvtThermal<Scalar>(*static_cast<const Opm::OilPvtThermal<Scalar>*>(data.realOilPvt_));
break;
case BrineCo2Pvt:
realOilPvt_ = new Opm::BrineCo2Pvt<Scalar>(*static_cast<const Opm::BrineCo2Pvt<Scalar>*>(data.realOilPvt_));
break;
default:
break;
}

3
opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp
View File

@ -380,6 +380,9 @@ public:
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return isothermalPvt_->oilReferenceDensity(regionIdx); }
const std::vector<TabulatedOneDFunction>& oilvisctCurves() const
{ return oilvisctCurves_; }

6
opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp
View File

@ -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.
*/

3
opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp
View File

@ -305,6 +305,9 @@ public:
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }
const Scalar waterReferenceDensity(unsigned regionIdx) const
{ return isothermalPvt_->waterReferenceDensity(regionIdx); }
const std::vector<Scalar>& viscrefPress() const
{ return viscrefPress_; }

10
opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp
View File

@ -628,13 +628,11 @@ public:
throw NumericalIssue(errlog.str());
}
const std::vector<Scalar>& gasReferenceDensity() const {
return gasReferenceDensity_;
}
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }
const std::vector<Scalar>& oilReferenceDensity() const {
return oilReferenceDensity_;
}
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return oilReferenceDensity_[regionIdx]; }
const std::vector<TabulatedTwoDFunction>& inverseGasB() const {
return inverseGasB_;

196
tests/test_co2brinepvt.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp>
//#include <opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp>
#include <opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp>
#include <opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp>
#include <opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/Python/Python.hpp>
#include <dune/common/parallel/mpihelper.hh>
// 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 <class Evaluation, class BrinePvt, class Co2Pvt>
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 <class Scalar>
inline void testAll()
{
Opm::Parser parser;
auto python = std::make_shared<Opm::Python>();
auto deck = parser.parseString(deckString1);
Opm::EclipseState eclState(deck);
Opm::Schedule schedule(deck, eclState, python);
Opm::GasPvtMultiplexer<Scalar> co2Pvt;
Opm::OilPvtMultiplexer<Scalar> brinePvt;
co2Pvt.initFromState(eclState, schedule);
brinePvt.initFromState(eclState, schedule);
typedef Opm::DenseAd::Evaluation<Scalar, 1> FooEval;
ensurePvtApi<Scalar>(brinePvt, co2Pvt);
ensurePvtApi<FooEval>(brinePvt, co2Pvt);
}
int main(int argc, char **argv)
{
Dune::MPIHelper::instance(argc, argv);
testAll<double>();
testAll<float>();
return 0;
}

33
tests/test_components.cpp
View File

@ -44,6 +44,7 @@
#include <opm/material/components/iapws/Common.hpp>
#include <opm/material/components/iapws/Region4.hpp>
#include <opm/material/components/H2O.hpp>
#include <opm/material/components/SimpleHuDuanH2O.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/components/Mesitylene.hpp>
#include <opm/material/components/TabulatedComponent.hpp>
@ -62,6 +63,36 @@ namespace ComponentsTest {
#include <dune/common/parallel/mpihelper.hh>
template <class Scalar, class Evaluation>
void testSimpleH2O()
{
typedef Opm::H2O<Scalar> H2O;
typedef Opm::SimpleHuDuanH2O<Scalar> SimpleHuDuanH2O;
typedef Opm::MathToolbox<Evaluation> 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 <class Scalar, class Evaluation>
void testAllComponents()
{
@ -91,6 +122,8 @@ inline void testAll()
// ensure that all components are API-compliant
testAllComponents<Scalar, Scalar>();
testAllComponents<Scalar, Evaluation>();
testSimpleH2O<Scalar, Evaluation>();
}

8
tests/test_fluidsystems.cpp
View File

@ -60,12 +60,6 @@
#include <opm/parser/eclipse/Python/Python.hpp>
namespace Opm {
namespace FluidSystemsTest {
#include <opm/material/components/co2tables.inc>
} }
#include <dune/common/parallel/mpihelper.hh>
// check that the blackoil fluid system implements all non-standard functions
@ -248,7 +242,7 @@ void testAllFluidSystems()
}
// Brine -- CO2
{ typedef Opm::BrineCO2FluidSystem<Scalar, Opm::FluidSystemsTest::CO2Tables> FluidSystem;
{ typedef Opm::BrineCO2FluidSystem<Scalar, CO2Tables> FluidSystem;
checkFluidSystem<Scalar, FluidSystem, FluidStateEval, LhsEval>(); }
// H2O -- N2