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

Merge ebd70f4d22 into 9a0ccfd853

Browse Source
This commit is contained in:
Tor Harald Sandve 2023-11-26 11:16:32 +00:00 committed by GitHub
commit 778f4f68a6
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
15 changed files with 42001 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -444,6 +444,7 @@ list (APPEND TEST_SOURCE_FILES
tests/material/test_2dtables.cpp
tests/material/test_blackoilfluidstate.cpp
tests/material/test_components.cpp
tests/material/test_binarycoefficients.cpp
tests/material/test_fluidmatrixinteractions.cpp
tests/material/test_fluidsystems.cpp
tests/material/test_spline.cpp

4
opm/input/eclipse/EclipseState/Tables/TableManager.hpp
View File

@ -201,6 +201,8 @@ namespace Opm {
double salinity() const;
int actco2s() const;
bool operator==(const TableManager& data) const;
template<class Serializer>
@ -261,6 +263,7 @@ namespace Opm {
serializer(m_gas_comp_index);
serializer(m_rtemp);
serializer(m_salinity);
serializer(m_actco2s);
serializer(m_tlmixpar);
serializer(m_ppcwmax);
if (!serializer.isSerializing()) {
@ -412,6 +415,7 @@ namespace Opm {
std::size_t m_gas_comp_index = 77;
double m_rtemp {288.7056}; // 60 Fahrenheit in Kelvin
double m_salinity {0.0};
int m_actco2s {1};
struct SplitSimpleTables {
size_t plyshMax = 0;

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

@ -30,6 +30,9 @@
#include <opm/material/IdealGas.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/common/TimingMacros.hpp>
#include <array>
namespace Opm {
namespace BinaryCoeff {
@ -100,37 +103,70 @@ public:
const int knownPhaseIdx,
Evaluation& xlCO2,
Evaluation& ygH2O,
const int& activityModel,
bool extrapolate = false)
{
OPM_TIMEFUNCTION_LOCAL();
Evaluation A = computeA_(temperature, pg, extrapolate);
/* salinity: conversion from mass fraction to mol fraction */
Evaluation x_NaCl = salinityToMolFrac_(salinity);
// if both phases are present the mole fractions in each phase can be calculate
// with the mutual solubility function
// High- or low-temperature case?
bool highTemp;
if (temperature > 372.15) {
highTemp = true;
}
else {
highTemp = false;
}
// If both phases are present the mole fractions in each phase can be calculate with the mutual solubility
// function
if (knownPhaseIdx < 0) {
Evaluation molalityNaCl = moleFracToMolality_(x_NaCl); // molality of NaCl //CHANGED
Evaluation m0_CO2 = molalityCO2inPureWater_(temperature, pg, extrapolate); // molality of CO2 in pure water
Evaluation gammaStar = activityCoefficient_(temperature, pg, molalityNaCl);// activity coefficient of CO2 in brine
Evaluation m_CO2 = m0_CO2 / gammaStar; // molality of CO2 in brine
xlCO2 = m_CO2 / (molalityNaCl + 55.508 + m_CO2); // mole fraction of CO2 in brine
ygH2O = A * (1 - xlCO2 - x_NaCl); // mole fraction of water in the gas phase
// Duan-Sun model as given in Spycher & Pruess (2005) have a different fugacity coefficient formula and
// activity coefficient definition (not a true activity coefficient but a ratio).
// Technically only valid below T = 100 C, but we use low-temp. parameters and formulas even above 100 C as
// an approximation.
if (activityModel == 3) {
auto [xCO2, yH2O] = mutualSolubilitySpycherPruess2005_(temperature, pg, molalityNaCl, extrapolate);
xlCO2 = xCO2;
ygH2O = yH2O;
}
else {
// High-temp. cases need iterations
if (highTemp) {
auto [xCO2, yH2O] = highTempSolubility_(temperature, pg, molalityNaCl, activityModel, extrapolate);
xlCO2 = xCO2;
ygH2O = yH2O;
}
// While low-temp. cases don't need it
else {
auto [xCO2, yH2O] = lowTempSolubility_(temperature, pg, molalityNaCl, activityModel, extrapolate);
xlCO2 = xCO2;
ygH2O = yH2O;
}
}
}
// if only liquid phase is present the mole fraction of CO2 in brine is given and
// and the virtual equilibrium mole fraction of water in the non-existing gas phase can be estimated
// with the mutual solubility function
if (knownPhaseIdx == liquidPhaseIdx)
if (knownPhaseIdx == liquidPhaseIdx && activityModel == 3) {
const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
ygH2O = A * (1 - xlCO2 - x_NaCl);
}
// if only gas phase is present the mole fraction of water in the gas phase is given and
// and the virtual equilibrium mole fraction of CO2 in the non-existing liquid phase can be estimated
// with the mutual solubility function
if (knownPhaseIdx == gasPhaseIdx)
if (knownPhaseIdx == gasPhaseIdx && activityModel == 3) {
//y_H2o = fluidstate.
const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
xlCO2 = 1 - x_NaCl - ygH2O / A;
}
}
/*!
@ -149,7 +185,12 @@ public:
* \param pg the gas phase pressure [Pa]
*/
template <class Evaluation>
static Evaluation fugacityCoefficientCO2(const Evaluation& temperature, const Evaluation& pg, bool extrapolate = false)
static Evaluation fugacityCoefficientCO2(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& yH2O,
const bool highTemp,
bool extrapolate = false,
bool spycherPruess2005 = false)
{
OPM_TIMEFUNCTION_LOCAL();
Valgrind::CheckDefined(temperature);
@ -157,24 +198,36 @@ public:
Evaluation V = 1 / (CO2::gasDensity(temperature, pg, extrapolate) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
Evaluation pg_bar = pg / 1.e5; // gas phase pressure in bar
Evaluation a_CO2 = (7.54e7 - 4.13e4 * temperature); // mixture parameter of Redlich-Kwong equation
Scalar b_CO2 = 27.8; // mixture parameter of Redlich-Kwong equation
Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
// Parameters in Redlich-Kwong equation
Evaluation a_CO2 = aCO2_(temperature, highTemp);
Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
Scalar b_CO2 = bCO2_(highTemp);
Evaluation b_mix = bMix_(yH2O, highTemp);
Evaluation lnPhiCO2;
lnPhiCO2 = log(V / (V - b_CO2));
lnPhiCO2 += b_CO2 / (V - b_CO2);
lnPhiCO2 -= 2 * a_CO2 / (R * pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
lnPhiCO2 +=
a_CO2 * b_CO2
/ (R
* pow(temperature, 1.5)
* b_CO2
* b_CO2)
* (log((V + b_CO2) / V)
- b_CO2 / (V + b_CO2));
lnPhiCO2 -= log(pg_bar * V / (R * temperature));
if (spycherPruess2005) {
lnPhiCO2 = log(V / (V - b_CO2));
lnPhiCO2 += b_CO2 / (V - b_CO2);
lnPhiCO2 -= 2 * a_CO2 / (R * pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
lnPhiCO2 +=
a_CO2 * b_CO2
/ (R
* pow(temperature, 1.5)
* b_CO2
* b_CO2)
* (log((V + b_CO2) / V)
- b_CO2 / (V + b_CO2));
lnPhiCO2 -= log(pg_bar * V / (R * temperature));
}
else {
lnPhiCO2 = (b_CO2 / b_mix) * (pg_bar * V / (R * temperature) - 1);
lnPhiCO2 -= log(pg_bar * (V - b_mix) / (R * temperature));
lnPhiCO2 += (2 * (yH2O * a_CO2_H2O + (1 - yH2O) * a_CO2) / a_mix - (b_CO2 / b_mix)) *
a_mix / (b_mix * R * pow(temperature, 1.5)) * log(V / (V + b_mix));
}
return exp(lnPhiCO2); // fugacity coefficient of CO2
}
@ -187,29 +240,261 @@ public:
* \param pg the gas phase pressure [Pa]
*/
template <class Evaluation>
static Evaluation fugacityCoefficientH2O(const Evaluation& temperature, const Evaluation& pg, bool extrapolate = false)
static Evaluation fugacityCoefficientH2O(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& yH2O,
const bool highTemp,
bool extrapolate = false,
bool spycherPruess2005 = false)
{
OPM_TIMEFUNCTION_LOCAL();
Valgrind::CheckDefined(temperature);
Valgrind::CheckDefined(pg);
const Evaluation& V = 1 / (CO2::gasDensity(temperature, pg, extrapolate) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
const Evaluation& pg_bar = pg / 1.e5; // gas phase pressure in bar
const Evaluation& a_CO2 = (7.54e7 - 4.13e4 * temperature);// mixture parameter of Redlich-Kwong equation
Scalar a_CO2_H2O = 7.89e7;// mixture parameter of Redlich-Kwong equation
Scalar b_CO2 = 27.8;// mixture parameter of Redlich-Kwong equation
Scalar b_H2O = 18.18;// mixture parameter of Redlich-Kwong equation
Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
Evaluation lnPhiH2O;
lnPhiH2O =
log(V/(V - b_CO2))
+ b_H2O/(V - b_CO2) - 2*a_CO2_H2O
/ (R*pow(temperature, 1.5)*b_CO2)*log((V + b_CO2)/V)
+ a_CO2*b_H2O/(R*pow(temperature, 1.5)*b_CO2*b_CO2)
*(log((V + b_CO2)/V) - b_CO2/(V + b_CO2))
- log(pg_bar*V/(R*temperature));
// Mixture parameter of Redlich-Kwong equation
Evaluation a_H2O = aH2O_(temperature, highTemp);
Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
Scalar b_H2O = bH2O_(highTemp);
Evaluation b_mix = bMix_(yH2O, highTemp);
Evaluation lnPhiH2O;
if (spycherPruess2005) {
lnPhiH2O =
log(V/(V - b_mix))
+ b_H2O/(V - b_mix) - 2*a_CO2_H2O
/ (R*pow(temperature, 1.5)*b_mix)*log((V + b_mix)/V)
+ a_mix*b_H2O/(R*pow(temperature, 1.5)*b_mix*b_mix)
*(log((V + b_mix)/V) - b_mix/(V + b_mix))
- log(pg_bar*V/(R*temperature));
}
else {
lnPhiH2O = (b_H2O / b_mix) * (pg_bar * V / (R * temperature) - 1);
lnPhiH2O -= log(pg_bar * (V - b_mix) / (R * temperature));
lnPhiH2O += (2 * (yH2O * a_H2O + (1 - yH2O) * a_CO2_H2O) / a_mix - (b_H2O / b_mix)) *
a_mix / (b_mix * R * pow(temperature, 1.5)) * log(V / (V + b_mix));
}
return exp(lnPhiH2O); // fugacity coefficient of H2O
}
private:
/*!
* \brief
*/
template <class Evaluation>
static Evaluation aCO2_(const Evaluation& temperature, const bool& highTemp)
{
if (highTemp) {
return 8.008e7 - 4.984e4 * temperature;
}
else {
return 7.54e7 - 4.13e4 * temperature;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation aH2O_(const Evaluation& temperature, const bool& highTemp)
{
if (highTemp) {
return 1.337e8 - 1.4e4 * temperature;
}
else {
return 0.0;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation aCO2_H2O_(const Evaluation& temperature, const Evaluation& yH2O, const bool& highTemp)
{
if (highTemp) {
// Pure parameters
Evaluation aCO2 = aCO2_(temperature, highTemp);
Evaluation aH2O = aH2O_(temperature, highTemp);
// Mixture Eq. (A-6)
Evaluation K_CO2_H2O = 0.4228 - 7.422e-4 * temperature;
Evaluation K_H2O_CO2 = 1.427e-2 - 4.037e-4 * temperature;
Evaluation k_CO2_H2O = yH2O * K_H2O_CO2 + (1 - yH2O) * K_CO2_H2O;
// Eq. (A-5)
return sqrt(aCO2 * aH2O) * (1 - k_CO2_H2O);
}
else {
return 7.89e7;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation aMix_(const Evaluation& temperature, const Evaluation& yH2O, const bool& highTemp)
{
if (highTemp) {
// Parameters
Evaluation aCO2 = aCO2_(temperature, highTemp);
Evaluation aH2O = aH2O_(temperature, highTemp);
Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
return yH2O * yH2O * aH2O + 2 * yH2O * (1 - yH2O) * a_CO2_H2O + (1 - yH2O) * (1 - yH2O) * aCO2;
}
else {
return aCO2_(temperature, highTemp);
}
}
/*!
* \brief
*/
static Scalar bCO2_(const bool& highTemp)
{
if (highTemp) {
return 28.25;
}
else {
return 27.8;
}
}
/*!
* \brief
*/
static Scalar bH2O_(const bool& highTemp)
{
if (highTemp) {
return 15.7;
}
else {
return 18.18;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation bMix_(const Evaluation& yH2O, const bool& highTemp)
{
if (highTemp) {
// Parameters
Scalar bCO2 = bCO2_(highTemp);
Scalar bH2O = bH2O_(highTemp);
return yH2O * bH2O + (1 - yH2O) * bCO2;
}
else {
return bCO2_(highTemp);
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation V_avg_CO2_(const Evaluation& temperature, const bool& highTemp)
{
if (highTemp && (temperature > 373.15)) {
return 32.6 + 3.413e-2 * (temperature - 373.15);
}
else {
return 32.6;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation V_avg_H2O_(const Evaluation& temperature, const bool& highTemp)
{
if (highTemp && (temperature > 373.15)) {
return 18.1 + 3.137e-2 * (temperature - 373.15);
}
else {
return 18.1;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation AM_(const Evaluation& temperature, const bool& highTemp)
{
if (highTemp && temperature > 373.15) {
Evaluation deltaTk = temperature - 373.15;
return deltaTk * (-3.084e-2 + 1.927e-5 * deltaTk);
}
else {
return 0.0;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation Pref_(const Evaluation& temperature, const bool& highTemp)
{
if (highTemp && temperature > 373.15) {
const Evaluation& temperatureCelcius = temperature - 273.15;
static const Scalar c[5] = { -1.9906e-1, 2.0471e-3, 1.0152e-4, -1.4234e-6, 1.4168e-8 };
return c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
}
else {
return 1.0;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation activityCoefficientCO2_(const Evaluation& temperature,
const Evaluation& xCO2,
const bool& highTemp)
{
if (highTemp) {
// Eq. (13)
Evaluation AM = AM_(temperature, highTemp);
Evaluation lnGammaCO2 = 2 * AM * xCO2 * (1 - xCO2) * (1 - xCO2);
return exp(lnGammaCO2);
}
else {
return 1.0;
}
}
/*!
* \brief
*/
template <class Evaluation>
static Evaluation activityCoefficientH2O_(const Evaluation& temperature,
const Evaluation& xCO2,
const bool& highTemp)
{
if (highTemp) {
// Eq. (12)
Evaluation AM = AM_(temperature, highTemp);
Evaluation lnGammaH2O = (1 - 2 * (1 - xCO2)) * AM * xCO2 * xCO2;
return exp(lnGammaH2O);
}
else {
return 1.0;
}
}
/*!
* \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction
*
@ -228,10 +513,10 @@ private:
}
/*!
* \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction (mol NaCl / mol solution)
*
* \param x_NaCl mole fraction of NaCL in brine [mol/mol]
*/
* \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction (mol NaCl / mol solution)
*
* \param x_NaCl mole fraction of NaCL in brine [mol/mol]
*/
template <class Evaluation>
static Evaluation moleFracToMolality_(const Evaluation& x_NaCl)
{
@ -240,45 +525,164 @@ private:
}
/*!
* \brief Returns the equilibrium molality of CO2 (mol CO2 / kg water) for a
* CO2-water mixture at a given pressure and temperature
*
* \param temperature The temperature [K]
* \param pg The gas phase pressure [Pa]
*/
* \brief Returns the mole fraction NaCl; inverse of moleFracToMolality
*
* \param x_NaCl mole fraction of NaCL in brine [mol/mol]
*/
template <class Evaluation>
static Evaluation molalityCO2inPureWater_(const Evaluation& temperature, const Evaluation& pg, bool extrapolate = false)
static Evaluation molalityToMoleFrac_(const Evaluation& m_NaCl)
{
OPM_TIMEFUNCTION_LOCAL();
const Evaluation& A = computeA_(temperature, pg, extrapolate); // according to Spycher, Pruess and Ennis-King (2003)
const Evaluation& B = computeB_(temperature, pg, extrapolate); // according to Spycher, Pruess and Ennis-King (2003)
const Evaluation& yH2OinGas = (1 - B) / (1. / A - B); // equilibrium mol fraction of H2O in the gas phase
const Evaluation& xCO2inWater = B * (1 - yH2OinGas); // equilibrium mol fraction of CO2 in the water phase
return (xCO2inWater * 55.508) / (1 - xCO2inWater); // CO2 molality
// conversion from molality to mole fractio (dissolved CO2 neglected)
return m_NaCl / (55.508 + m_NaCl);
}
/*!
* \brief Returns the activity coefficient of CO2 in brine for a
* molal description. According to "Duan and Sun 2003"
* given in "Spycher and Pruess 2005"
*
* \param temperature the temperature [K]
* \param pg the gas phase pressure [Pa]
* \param molalityNaCl molality of NaCl (mol NaCl / kg water)
*/
* \brief Fixed-point iterations for high-temperature cases
*/
template <class Evaluation>
static Evaluation activityCoefficient_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& molalityNaCl)
static std::pair<Evaluation, Evaluation> highTempSolubility_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& m_NaCl,
const int& activityModel,
bool extrapolate = false)
{
OPM_TIMEFUNCTION_LOCAL();
const Evaluation& lambda = computeLambda_(temperature, pg); // lambda_{CO2-Na+}
const Evaluation& xi = computeXi_(temperature, pg); // Xi_{CO2-Na+-Cl-}
const Evaluation& lnGammaStar =
2*molalityNaCl*lambda + xi*molalityNaCl*molalityNaCl;
return exp(lnGammaStar);
// Start point for fixed-point iterations as recommended below in section 2.2
Evaluation yH2O = H2O::vaporPressure(temperature) / pg; // ideal mixing
Evaluation xCO2 = 0.009; // same as ~0.5 mol/kg
Evaluation gammaNaCl = 1.0; // default salt activity coeff = 1.0
// We can pre-calculate Duan-Sun, Spycher & Pruess (2009) salt activity coeff.
if (m_NaCl > 0.0 && activityModel == 2) {
gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
}
// Options
int max_iter = 100;
Scalar tol = 1e-8;
const bool highTemp = true;
// Fixed-point loop x_i+1 = F(x_i)
for (int i = 0; i < max_iter; ++i) {
// Calculate activity coefficient for Rumpf et al (1994) model
if (m_NaCl > 0.0 && activityModel == 1) {
gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, xCO2, activityModel);
}
// F(x_i) is the mutual solubilities
auto [xCO2_new, yH2O_new] = mutualSolubility_(temperature, pg, xCO2, yH2O, m_NaCl, gammaNaCl, highTemp,
extrapolate);
// Check for convergence
if (abs(xCO2_new - xCO2) < tol && abs(yH2O_new - yH2O) < tol) {
xCO2 = xCO2_new;
yH2O = yH2O_new;
break;
}
// Else update mole fractions for next iteration
else {
xCO2 = xCO2_new;
yH2O = yH2O_new;
}
}
return {xCO2, yH2O};
}
/*!
* \brief Fixed-point iterations for high-temperature cases
*/
template <class Evaluation>
static std::pair<Evaluation, Evaluation> lowTempSolubility_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& m_NaCl,
const int& activityModel,
bool extrapolate = false)
{
// Calculate activity coefficient for salt
Evaluation gammaNaCl = 1.0;
if (m_NaCl > 0.0 && activityModel > 0 && activityModel < 3) {
gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
}
// Calculate mutual solubility.
// Note that we don't use xCO2 and yH2O input in low-temperature case, so we set them to 0.0
const bool highTemp = false;
auto [xCO2, yH2O] = mutualSolubility_(temperature, pg, Evaluation(0.0), Evaluation(0.0), m_NaCl, gammaNaCl,
highTemp, extrapolate);
return {xCO2, yH2O};
}
/*!
* \brief Mutual solubility according to Spycher & Pruess (2009)
*/
template <class Evaluation>
static std::pair<Evaluation, Evaluation> mutualSolubility_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& xCO2,
const Evaluation& yH2O,
const Evaluation& m_NaCl,
const Evaluation& gammaNaCl,
const bool& highTemp,
bool extrapolate = false)
{
// Calculate A and B (without salt effect); Eqs. (8) and (9)
const Evaluation& A = computeA_(temperature, pg, yH2O, xCO2, highTemp, extrapolate);
Evaluation B = computeB_(temperature, pg, yH2O, xCO2, highTemp, extrapolate);
// Add salt effect to B, Eq. (17)
B /= gammaNaCl;
// Compute yH2O and xCO2, Eqs. (B-7) and (B-2)
Evaluation yH2O_new = (1. - B) * 55.508 / ((1. / A - B) * (m_NaCl + 55.508) + m_NaCl * B);
Evaluation xCO2_new;
if (highTemp) {
xCO2_new = B * (1 - yH2O);
}
else {
xCO2_new = B * (1 - yH2O_new);
}
return {xCO2_new, yH2O_new};
}
/*!
* \brief Mutual solubility according to Spycher & Pruess (2009)
*/
template <class Evaluation>
static std::pair<Evaluation, Evaluation>
mutualSolubilitySpycherPruess2005_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& m_NaCl,
bool extrapolate = false)
{
// Calculate A and B (without salt effect); Eqs. (8) and (9)
const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
const Evaluation& B = computeB_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
// Mole fractions and molality in pure water
Evaluation yH2O = (1 - B) / (1. / A - B);
Evaluation xCO2 = B * (1 - yH2O);
// Modifiy mole fractions with Duan-Sun "activity coefficient" if salt is involved
if (m_NaCl > 0.0) {
const Evaluation& gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), 3);
// Molality with salt
Evaluation mCO2 = (xCO2 * 55.508) / (1 - xCO2); // pure water
mCO2 /= gammaNaCl;
xCO2 = mCO2 / (m_NaCl + 55.508 + mCO2);
// new yH2O with salt
const Evaluation& xNaCl = molalityToMoleFrac_(m_NaCl);
yH2O = A * (1 - xCO2 - xNaCl);
}
return {xCO2, yH2O};
}
/*!
* \brief Returns the paramater A for the calculation of
* them mutual solubility in the water-CO2 system.
@ -288,16 +692,36 @@ private:
* \param pg the gas phase pressure [Pa]
*/
template <class Evaluation>
static Evaluation computeA_(const Evaluation& temperature, const Evaluation& pg, bool extrapolate = false)
static Evaluation computeA_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& yH2O,
const Evaluation& xCO2,
const bool& highTemp,
bool extrapolate = false,
bool spycherPruess2005 = false)
{
OPM_TIMEFUNCTION_LOCAL();
const Evaluation& deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
Scalar v_av_H2O = 18.1; // average partial molar volume of H2O [cm^3/mol]
Scalar R = IdealGas::R * 10;
const Evaluation& k0_H2O = equilibriumConstantH2O_(temperature); // equilibrium constant for H2O at 1 bar
const Evaluation& phi_H2O = fugacityCoefficientH2O(temperature, pg, extrapolate); // fugacity coefficient of H2O for the water-CO2 system
// Intermediate calculations
const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp); // pressure range [bar] from pref to pg[bar]
Evaluation v_av_H2O = V_avg_H2O_(temperature, highTemp); // average partial molar volume of H2O [cm^3/mol]
Evaluation k0_H2O = equilibriumConstantH2O_(temperature, highTemp); // equilibrium constant for H2O at 1 bar
Evaluation phi_H2O = fugacityCoefficientH2O(temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005); // fugacity coefficient of H2O for the water-CO2 system
Evaluation gammaH2O = activityCoefficientH2O_(temperature, xCO2, highTemp);
// In the intermediate temperature range 99-109 C, equilibrium constants and fugacity coeff. are linearly
// weighted
if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
const Evaluation weight = (382.15 - temperature) / 10.;
const Evaluation& k0_H2O_low = equilibriumConstantH2O_(temperature, false);
const Evaluation& phi_H2O_low = fugacityCoefficientH2O(temperature, pg, Evaluation(0.0), false, extrapolate);
k0_H2O = k0_H2O * (1 - weight) + k0_H2O_low * weight;
phi_H2O = phi_H2O * (1 - weight) + phi_H2O_low * weight;
}
// Eq. (10)
const Evaluation& pg_bar = pg / 1.e5;
return k0_H2O/(phi_H2O*pg_bar)*exp(deltaP*v_av_H2O/(R*temperature));
Scalar R = IdealGas::R * 10;
return k0_H2O * gammaH2O / (phi_H2O * pg_bar) * exp(deltaP * v_av_H2O / (R * temperature));
}
/*!
@ -309,16 +733,119 @@ private:
* \param pg the gas phase pressure [Pa]
*/
template <class Evaluation>
static Evaluation computeB_(const Evaluation& temperature, const Evaluation& pg, bool extrapolate = false)
static Evaluation computeB_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& yH2O,
const Evaluation& xCO2,
const bool& highTemp,
bool extrapolate = false,
bool spycherPruess2005 = false)
{
OPM_TIMEFUNCTION_LOCAL();
const Evaluation& deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
const Scalar v_av_CO2 = 32.6; // average partial molar volume of CO2 [cm^3/mol]
const Scalar R = IdealGas::R * 10;
const Evaluation& k0_CO2 = equilibriumConstantCO2_(temperature); // equilibrium constant for CO2 at 1 bar
const Evaluation& phi_CO2 = fugacityCoefficientCO2(temperature, pg, extrapolate); // fugacity coefficient of CO2 for the water-CO2 system
// Intermediate calculations
const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp); // pressure range [bar] from pref to pg[bar]
Evaluation v_av_CO2 = V_avg_CO2_(temperature, highTemp); // average partial molar volume of CO2 [cm^3/mol]
Evaluation k0_CO2 = equilibriumConstantCO2_(temperature, pg, highTemp, spycherPruess2005); // equilibrium constant for CO2 at 1 bar
Evaluation phi_CO2 = fugacityCoefficientCO2(temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005); // fugacity coefficient of CO2 for the water-CO2 system
Evaluation gammaCO2 = activityCoefficientCO2_(temperature, xCO2, highTemp);
// In the intermediate temperature range 99-109 C, equilibrium constants and fugacity coeff. are linearly
// weighted
if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
const Evaluation weight = (382.15 - temperature) / 10.;
const Evaluation& k0_CO2_low = equilibriumConstantCO2_(temperature, pg, false, spycherPruess2005);
const Evaluation& phi_CO2_low = fugacityCoefficientCO2(temperature, pg, Evaluation(0.0), false, extrapolate);
k0_CO2 = k0_CO2 * (1 - weight) + k0_CO2_low * weight;
phi_CO2 = phi_CO2 * (1 - weight) + phi_CO2_low * weight;
}
// Eq. (11)
const Evaluation& pg_bar = pg / 1.e5;
return phi_CO2*pg_bar/(55.508*k0_CO2)*exp(-(deltaP*v_av_CO2)/(R*temperature));
const Scalar R = IdealGas::R * 10;
return phi_CO2 * pg_bar / (55.508 * k0_CO2 * gammaCO2) * exp(-deltaP * v_av_CO2 / (R * temperature));
}
/*!
* \brief Activity model of salt in Spycher & Pruess (2009)
*/
template <class Evaluation>
static Evaluation activityCoefficientSalt_(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& m_NaCl,
const Evaluation& xCO2,
const int& activityModel)
{
OPM_TIMEFUNCTION_LOCAL();
// Lambda and xi parameter for either Rumpf et al (1994) (activityModel = 1) or Duan-Sun as modified by Spycher
// & Pruess (2009) (activityModel = 2) or Duan & Sun (2003) as given in Spycher & Pruess (2005) (activityModel =
// 3)
Evaluation lambda;
Evaluation xi;
Evaluation convTerm;
if (activityModel == 1) {
lambda = computeLambdaRumpfetal_(temperature);
xi = -0.0028 * 3.0;
Evaluation m_CO2 = xCO2 * (m_NaCl + 55.508) / (1 - xCO2);
convTerm = (1 + (m_CO2 + m_NaCl) / 55.508) / (1 + m_CO2 / 55.508);
}
else if (activityModel == 2) {
lambda = computeLambdaSpycherPruess2009_(temperature);
xi = computeXiSpycherPruess2009_(temperature);
convTerm = 1 + m_NaCl / 55.508;
}
else if (activityModel == 3) {
lambda = computeLambdaDuanSun_(temperature, pg);
xi = computeXiDuanSun_(temperature, pg);
convTerm = 1.0;
}
else {
throw std::runtime_error("Activity model for salt-out effect has not been implemented!");
}
// Eq. (18)
const Evaluation& lnGamma = 2 * lambda * m_NaCl + xi * m_NaCl * m_NaCl;
// Eq. (18), return activity coeff. on mole-fraction scale
return convTerm * exp(lnGamma);
}
/*!
* \brief Lambda parameter in Duan & Sun model, as modified and detailed in Spycher & Pruess (2009)
*/
template <class Evaluation>
static Evaluation computeLambdaSpycherPruess2009_(const Evaluation& temperature)
{
// Table 1
static const Scalar c[3] = { 2.217e-4, 1.074, 2648. };
// Eq. (19)
return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
}
/*!
* \brief Xi parameter in Duan & Sun model, as modified and detailed in Spycher & Pruess (2009)
*/
template <class Evaluation>
static Evaluation computeXiSpycherPruess2009_(const Evaluation& temperature)
{
// Table 1
static const Scalar c[3] = { 1.3e-5, -20.12, 5259. };
// Eq. (19)
return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
}
/*!
* \brief Lambda parameter in Rumpf et al. (1994), as detailed in Spycher & Pruess (2005)
*/
template <class Evaluation>
static Evaluation computeLambdaRumpfetal_(const Evaluation& temperature)
{
// B^(0) below Eq. (A-6)
static const Scalar c[4] = { 0.254, -76.82, -10656, 6312e3 };
return c[0] + c[1] / temperature + c[2] / (temperature * temperature) +
c[3] / (temperature * temperature * temperature);
}
/*!
@ -329,19 +856,13 @@ private:
* \param pg the gas phase pressure [Pa]
*/
template <class Evaluation>
static Evaluation computeLambda_(const Evaluation& temperature, const Evaluation& pg)
static Evaluation computeLambdaDuanSun_(const Evaluation& temperature, const Evaluation& pg)
{
OPM_TIMEFUNCTION_LOCAL();
static const Scalar c[6] =
{ -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
return
c[0]
+ c[1]*temperature
+ c[2]/temperature
+ c[3]*pg_bar/temperature
+ c[4]*pg_bar/(630.0 - temperature)
return c[0] + c[1]*temperature + c[2]/temperature + c[3]*pg_bar/temperature + c[4]*pg_bar/(630.0 - temperature)
+ c[5]*temperature*log(pg_bar);
}
@ -353,9 +874,8 @@ private:
* \param pg the gas phase pressure [Pa]
*/
template <class Evaluation>
static Evaluation computeXi_(const Evaluation& temperature, const Evaluation& pg)
static Evaluation computeXiDuanSun_(const Evaluation& temperature, const Evaluation& pg)
{
OPM_TIMEFUNCTION_LOCAL();
static const Scalar c[4] =
{ 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
@ -370,12 +890,29 @@ private:
* \param temperature the temperature [K]
*/
template <class Evaluation>
static Evaluation equilibriumConstantCO2_(const Evaluation& temperature)
static Evaluation equilibriumConstantCO2_(const Evaluation& temperature,
const Evaluation& pg,
const bool& highTemp,
bool spycherPruess2005 = false)
{
OPM_TIMEFUNCTION_LOCAL();
Evaluation temperatureCelcius = temperature - 273.15;
static const Scalar c[3] = { 1.189, 1.304e-2, -5.446e-5 };
Evaluation logk0_CO2 = c[0] + temperatureCelcius*(c[1] + temperatureCelcius*c[2]);
std::array<Scalar, 4> c;
if (highTemp) {
c = { 1.668, 3.992e-3, -1.156e-5, 1.593e-9 };
}
else {
// For temperature below 31 C and pressures above saturation pressure, separate parameters are needed
Evaluation psat = CO2::vaporPressure(temperature);
if (temperatureCelcius < 31 && pg > psat && !spycherPruess2005) {
c = { 1.169, 1.368e-2, -5.38e-5, 0.0 };
}
else {
c = { 1.189, 1.304e-2, -5.446e-5, 0.0 };
}
}
Evaluation logk0_CO2 = c[0] + temperatureCelcius * (c[1] + temperatureCelcius *
(c[2] + temperatureCelcius * c[3]));
Evaluation k0_CO2 = pow(10.0, logk0_CO2);
return k0_CO2;
}
@ -387,13 +924,18 @@ private:
* \param temperature the temperature [K]
*/
template <class Evaluation>
static Evaluation equilibriumConstantH2O_(const Evaluation& temperature)
static Evaluation equilibriumConstantH2O_(const Evaluation& temperature, const bool& highTemp)
{
OPM_TIMEFUNCTION_LOCAL();
Evaluation temperatureCelcius = temperature - 273.15;
static const Scalar c[4] = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7 };
Evaluation logk0_H2O =
c[0] + temperatureCelcius*(c[1] + temperatureCelcius*(c[2] + temperatureCelcius*c[3]));
std::array<Scalar, 5> c;
if (highTemp){
c = { -2.1077, 2.8127e-2, -8.4298e-5, 1.4969e-7, -1.1812e-10 };
}
else {
c = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7, 0.0 };
}
Evaluation logk0_H2O = c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
return pow(10.0, logk0_H2O);
}

63
opm/material/components/BrineDynamic.hpp
View File

@ -260,24 +260,31 @@ public:
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure, const Evaluation& salinity, bool extrapolate = false)
{
Evaluation tempC = temperature - 273.15;
Evaluation pMPa = pressure/1.0E6;
// Evaluation tempC = temperature - 273.15;
// Evaluation pMPa = pressure/1.0E6;
const Evaluation rhow = H2O::liquidDensity(temperature, pressure, extrapolate);
return
rhow +
1000*salinity*(
0.668 +
0.44*salinity +
1.0E-6*(
300*pMPa -
2400*pMPa*salinity +
tempC*(
80.0 +
3*tempC -
3300*salinity -
13*pMPa +
47*pMPa*salinity)));
static constexpr Scalar a[3] = {0.30960, -0.000069, 0.0};
const Evaluation& T = temperature - 273.15;
const Evaluation coeff = a[0] + T * (a[1] + a[2] * T);
const Evaluation exponent = coeff * salinity;
return rhow * pow(10.0, exponent);
// return
// rhow +
// 1000*salinity*(
// 0.668 +
// 0.44*salinity +
// 1.0E-6*(
// 300*pMPa -
// 2400*pMPa*salinity +
// tempC*(
// 80.0 +
// 3*tempC -
// 3300*salinity -
// 13*pMPa +
// 47*pMPa*salinity)));
}
/*!
@ -335,16 +342,24 @@ public:
* "Equations of State for basin geofluids"
*/
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& /*pressure*/, const Evaluation& salinity)
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure, const Evaluation& salinity, bool extrapolate = false)
{
Evaluation T_C = temperature - 273.15;
if(temperature <= 275.) // regularization
T_C = 275.0;
const Evaluation muw = H2O::liquidViscosity(temperature, pressure, extrapolate);
static constexpr Scalar a[3] = {0.718000, 0.003590, 0.0};
Evaluation A = (0.42*Opm::pow((Opm::pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
Evaluation mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*exp(-A);
const Evaluation& T = temperature - 273.15;
const Evaluation coeff = a[0] + T * (a[1] + a[2] * T);
const Evaluation exponent = coeff * salinity;
return mu_brine/1000.0; // convert to [Pa s] (todo: check if correct cP->Pa s is times 10...)
return muw * pow(10.0, exponent);
// Evaluation T_C = temperature - 273.15;
// if(temperature <= 275.) // regularization
// T_C = 275.0;
// Evaluation A = (0.42*Opm::pow((Opm::pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
// Evaluation mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*exp(-A);
// return mu_brine/1000.0; // convert to [Pa s] (todo: check if correct cP->Pa s is times 10...)
}
//Molar mass salt (assumes pure NaCl) [kg/mol]

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

@ -331,6 +331,10 @@ public:
assert(temperature > 0);
assert(pressure > 0);
// Activity model for salt-out effect in brine-CO2 mutual solubility
// 2 = Duan & Sun as modified in Spycher & Pruess, Trans. Porous Media, (2009)
const int activityModel = 3;
// calulate the equilibrium composition for the given
// temperature and pressure. TODO: calculateMoleFractions()
// could use some cleanup.
@ -341,7 +345,8 @@ public:
LhsEval(Brine_IAPWS::salinity),
/*knownPhaseIdx=*/-1,
xlCO2,
xgH2O);
xgH2O,
activityModel);
// normalize the phase compositions
xlCO2 = max(0.0, min(1.0, xlCO2));

150
opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp
View File

@ -78,15 +78,12 @@ public:
explicit BrineCo2Pvt() = default;
BrineCo2Pvt(const std::vector<Scalar>& salinity,
int activityModel = 1,
Scalar T_ref = 288.71, //(273.15 + 15.56)
Scalar P_ref = 101325)
: salinity_(salinity)
{
// Throw an error if reference state is not (T, p) = (15.56 C, 1 atm) = (288.71 K, 1.01325e5 Pa)
if (T_ref != Scalar(288.71) || P_ref != Scalar(1.01325e5)) {
OPM_THROW(std::runtime_error,
"BrineCo2Pvt class can only be used with default reference state (T, P) = (288.71 K, 1.01325e5 Pa)!");
}
setActivityModelSalt(activityModel);
int num_regions = salinity_.size();
co2ReferenceDensity_.resize(num_regions);
brineReferenceDensity_.resize(num_regions);
@ -154,6 +151,18 @@ public:
void setEnableSaltConcentration(bool yesno)
{ enableSaltConcentration_ = yesno; }
/*!
* \brief Set activity coefficient model for salt in solubility model
*/
void setActivityModelSalt(int activityModel)
{
// Only 0, 1, 2 and 3 are allowed
if (activityModel > 3 || activityModel < 0) {
throw std::runtime_error("ACTCO2S options are 0, 1, 2 or 3; see manual!");
}
activityModel_ = activityModel;
}
/*!
* \brief Return the number of PVT regions which are considered by this PVT-object.
*/
@ -177,7 +186,7 @@ public:
pressure,
salinity,
xlCO2)
- pressure / density_(regionIdx, temperature, pressure, Rs, salinity ));
- pressure / density(regionIdx, temperature, pressure, Rs, salinity ));
}
/*!
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
@ -194,7 +203,7 @@ public:
pressure,
Evaluation(salinity_[regionIdx]),
xlCO2)
- pressure / density_(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx])));
- pressure / density(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx])));
}
/*!
@ -263,8 +272,8 @@ public:
{
OPM_TIMEFUNCTION_LOCAL();
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltconcentration);
Evaluation rsSat = rsSat_(regionIdx, temperature, pressure, salinity);
return (1.0 - convertRsToXoG_(rsSat,regionIdx)) * density_(regionIdx, temperature, pressure, rsSat, salinity)/brineReferenceDensity_[regionIdx];
Evaluation rs_sat = rsSat(regionIdx, temperature, pressure, salinity);
return (1.0 - convertRsToXoG_(rs_sat,regionIdx)) * density(regionIdx, temperature, pressure, rs_sat, salinity)/brineReferenceDensity_[regionIdx];
}
/*!
* \brief Returns the formation volume factor [-] of the fluid phase.
@ -278,7 +287,7 @@ public:
{
OPM_TIMEFUNCTION_LOCAL();
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
return (1.0 - convertRsToXoG_(Rs,regionIdx)) * density_(regionIdx, temperature, pressure, Rs, salinity)/brineReferenceDensity_[regionIdx];
return (1.0 - convertRsToXoG_(Rs,regionIdx)) * density(regionIdx, temperature, pressure, Rs, salinity)/brineReferenceDensity_[regionIdx];
}
/*!
* \brief Returns the formation volume factor [-] of the fluid phase.
@ -289,7 +298,7 @@ public:
const Evaluation& pressure,
const Evaluation& Rs) const
{
return (1.0 - convertRsToXoG_(Rs,regionIdx)) * density_(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx]))/brineReferenceDensity_[regionIdx];
return (1.0 - convertRsToXoG_(Rs,regionIdx)) * density(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx]))/brineReferenceDensity_[regionIdx];
}
/*!
@ -300,9 +309,9 @@ public:
const Evaluation& temperature,
const Evaluation& pressure) const
{
OPM_TIMEFUNCTION_LOCAL();
Evaluation rsSat = rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
return (1.0 - convertRsToXoG_(rsSat,regionIdx)) * density_(regionIdx, temperature, pressure, rsSat, Evaluation(salinity_[regionIdx]))/brineReferenceDensity_[regionIdx];
OPM_TIMEFUNCTION_LOCAL();
Evaluation rs_sat = rsSat(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
return (1.0 - convertRsToXoG_(rs_sat,regionIdx)) * density(regionIdx, temperature, pressure, rs_sat, Evaluation(salinity_[regionIdx]))/brineReferenceDensity_[regionIdx];
}
/*!
@ -345,7 +354,7 @@ public:
const Evaluation& /*maxOilSaturation*/) const
{
//TODO support VAPPARS
return rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
return rsSat(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
/*!
@ -358,7 +367,7 @@ public:
const Evaluation& saltConcentration) const
{
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
return rsSat_(regionIdx, temperature, pressure, salinity);
return rsSat(regionIdx, temperature, pressure, salinity);
}
/*!
@ -369,7 +378,7 @@ public:
const Evaluation& temperature,
const Evaluation& pressure) const
{
return rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
return rsSat(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
const Scalar oilReferenceDensity(unsigned regionIdx) const
@ -402,23 +411,16 @@ public:
return pow(Evaluation(10), log_D_Brine) * 1e-4; // convert from cm2/s to m2/s
}
private:
std::vector<Scalar> brineReferenceDensity_;
std::vector<Scalar> co2ReferenceDensity_;
std::vector<Scalar> salinity_;
bool enableDissolution_ = true;
bool enableSaltConcentration_ = false;
template <class LhsEval>
LhsEval density_(unsigned regionIdx,
const LhsEval& temperature,
const LhsEval& pressure,
const LhsEval& Rs,
const LhsEval& salinity) const
template <class Evaluation>
Evaluation density(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rs,
const Evaluation& salinity) const
{
OPM_TIMEFUNCTION_LOCAL();
LhsEval xlCO2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx), salinity);
LhsEval result = liquidDensity_(temperature,
OPM_TIMEFUNCTION_LOCAL();
Evaluation xlCO2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx), salinity);
Evaluation result = liquidDensity_(temperature,
pressure,
xlCO2,
salinity);
@ -427,6 +429,42 @@ private:
return result;
}
template <class Evaluation>
Evaluation rsSat(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& salinity) const
{
OPM_TIMEFUNCTION_LOCAL();
if (!enableDissolution_)
return 0.0;
// calulate the equilibrium composition for the given
// temperature and pressure.
Evaluation xgH2O;
Evaluation xlCO2;
BinaryCoeffBrineCO2::calculateMoleFractions(temperature,
pressure,
salinity,
/*knownPhaseIdx=*/-1,
xlCO2,
xgH2O,
activityModel_,
extrapolate);
// normalize the phase compositions
xlCO2 = max(0.0, min(1.0, xlCO2));
return convertXoGToRs(convertxoGToXoG(xlCO2, salinity), regionIdx);
}
private:
std::vector<Scalar> brineReferenceDensity_;
std::vector<Scalar> co2ReferenceDensity_;
std::vector<Scalar> salinity_;
bool enableDissolution_ = true;
bool enableSaltConcentration_ = false;
int activityModel_;
template <class LhsEval>
LhsEval liquidDensity_(const LhsEval& T,
@ -455,11 +493,19 @@ private:
}
const LhsEval& rho_brine = Brine::liquidDensity(T, pl, salinity, extrapolate);
const LhsEval& rho_pure = H2O::liquidDensity(T, pl, extrapolate);
const LhsEval& rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlCO2);
const LhsEval& contribCO2 = rho_lCO2 - rho_pure;
static constexpr Scalar a[3] = {0.103300, -0.000023, -0.000002};
const LhsEval& T_c = T - 273.15;
const LhsEval coeff = a[0] + T_c * (a[1] + a[2] * T_c);
const LhsEval& XCO2 = convertxoGToXoG(xlCO2, salinity);
const LhsEval exponent = coeff * XCO2;
return rho_brine + contribCO2;
return rho_brine * pow(10.0, exponent);
// const LhsEval& rho_pure = H2O::liquidDensity(T, pl, extrapolate);
// const LhsEval& rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlCO2);
// const LhsEval& contribCO2 = rho_lCO2 - rho_pure;
// return rho_brine + contribCO2;
}
template <class LhsEval>
@ -542,42 +588,12 @@ private:
return XoG/(1.0 - XoG)*(rho_oRef/rho_gRef);
}
template <class LhsEval>
LhsEval rsSat_(unsigned regionIdx,
const LhsEval& temperature,
const LhsEval& pressure,
const LhsEval& salinity) const
{
OPM_TIMEFUNCTION_LOCAL();
if (!enableDissolution_)
return 0.0;
// calulate the equilibrium composition for the given
// temperature and pressure.
LhsEval xgH2O;
LhsEval xlCO2;
BinaryCoeffBrineCO2::calculateMoleFractions(temperature,
pressure,
salinity,
/*knownPhaseIdx=*/-1,
xlCO2,
xgH2O,
extrapolate);
// normalize the phase compositions
xlCO2 = max(0.0, min(1.0, xlCO2));
return convertXoGToRs(convertxoGToXoG(xlCO2, salinity), regionIdx);
}
template <class LhsEval>
static LhsEval liquidEnthalpyBrineCO2_(const LhsEval& T,
const LhsEval& p,
const LhsEval& salinity,
const LhsEval& X_CO2_w)
{
OPM_TIMEFUNCTION_LOCAL();
/* X_CO2_w : mass fraction of CO2 in brine */
/* same function as enthalpy_brine, only extended by CO2 content */

20
opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp
View File

@ -64,15 +64,12 @@ public:
explicit Co2GasPvt() = default;
Co2GasPvt(const std::vector<Scalar>& salinity,
int activityModel = 1,
Scalar T_ref = 288.71, //(273.15 + 15.56)
Scalar P_ref = 101325)
: salinity_(salinity)
{
// Throw an error if reference state is not (T, p) = (15.56 C, 1 atm) = (288.71 K, 1.01325e5 Pa)
if (T_ref != Scalar(288.71) || P_ref != Scalar(1.01325e5)) {
OPM_THROW(std::runtime_error,
"BrineCo2Pvt class can only be used with default reference state (T, P) = (288.71 K, 1.01325e5 Pa)!");
}
setActivityModelSalt(activityModel);
int num_regions = salinity_.size();
setNumRegions(num_regions);
for (int i = 0; i < num_regions; ++i) {
@ -119,6 +116,17 @@ public:
void setEnableVaporizationWater(bool yesno)
{ enableVaporization_ = yesno; }
/*!
* \brief Set activity coefficient model for salt in solubility model
*/
void setActivityModelSalt(int activityModel)
{
// Only 0, 1, 2, 3 are allowed
if (activityModel > 3 || activityModel < 0) {
throw std::runtime_error("ACTCO2S options are 0, 1, 2 or 3; see manual!");
}
activityModel_ = activityModel;}
/*!
* \brief Finish initializing the co2 phase PVT properties.
*/
@ -314,6 +322,7 @@ private:
/*knownPhaseIdx=*/-1,
xlCO2,
xgH2O,
activityModel_,
extrapolate);
// normalize the phase compositions
@ -371,6 +380,7 @@ private:
std::vector<Scalar> gasReferenceDensity_;
std::vector<Scalar> salinity_;
bool enableVaporization_ = true;
int activityModel_;
};
} // namespace Opm

11
src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp
View File

@ -147,6 +147,8 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
m_salinity = ParserKeywords::SALINITY::MOLALITY::defaultValue;
m_actco2s = ParserKeywords::ACTCO2S::ACTIVITY_MODEL::defaultValue;
initDims( deck );
initSimpleTables( deck );
initFullTables(deck, "PVTG", m_pvtgTables);
@ -205,6 +207,9 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
if( deck.hasKeyword( "SALINITY" ) )
m_salinity = deck["SALINITY"].back().getRecord(0).getItem("MOLALITY").get<double>( 0 ); //unit independent of unit systems
if ( deck.hasKeyword( "ACTCO2S" ) )
m_actco2s = deck["ACTCO2S"].back().getRecord(0).getItem("ACTIVITY_MODEL").get<int>( 0 );
if ( deck.hasKeyword( "ROCK2D") )
initRockTables(deck, "ROCK2D", m_rock2dTables );
@ -322,6 +327,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
result.m_gas_comp_index = 77;
result.m_rtemp = 1.0;
result.m_salinity = 1.0;
result.m_actco2s = 1;
result.m_tlmixpar = TLMixpar::serializationTestObject();
result.m_ppcwmax = Ppcwmax::serializationTestObject();
return result;
@ -1247,6 +1253,10 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
return this->m_salinity;
}
int TableManager::actco2s() const {
return this->m_actco2s;
}
std::size_t TableManager::gas_comp_index() const {
return this->m_gas_comp_index;
}
@ -1300,6 +1310,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
jfunc == data.jfunc &&
m_rtemp == data.m_rtemp &&
m_salinity == data.m_salinity &&
m_actco2s == data.m_actco2s &&
m_gas_comp_index == data.m_gas_comp_index;
}

14
src/opm/input/eclipse/share/keywords/001_Eclipse300/A/ACTCO2S Normal file
View File

@ -0,0 +1,14 @@
{
"name": "ACTCO2S",
"sections": [
"PROPS"
],
"size": 1,
"items": [
{
"name": "ACTIVITY_MODEL",
"value_type": "INT",
"default": 3
}
]
}

1
src/opm/input/eclipse/share/keywords/keyword_list.cmake
View File

@ -1045,6 +1045,7 @@ set( keywords
000_Eclipse100/Z/ZIPPY2
000_Eclipse100/Z/ZIPP2OFF
001_Eclipse300/A/ACTCO2S
001_Eclipse300/B/BLOCK_PROBE300
001_Eclipse300/C/CIRCLE
001_Eclipse300/C/CNAMES

2
src/opm/material/components/CO2.cpp
View File

@ -27,7 +27,7 @@
#include <opm/material/common/quad.hpp>
#endif
#include "co2tables.inc"
#include "co2tables_rc.inc"
namespace Opm {

40896
src/opm/material/components/co2tables_rc.inc Normal file
View File

File diff suppressed because it is too large Load Diff

1
src/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.cpp
View File

@ -51,6 +51,7 @@ initFromState(const EclipseState& eclState, const Schedule&)
setEnableDissolvedGas(eclState.getSimulationConfig().hasDISGASW() || eclState.getSimulationConfig().hasDISGAS());
setEnableSaltConcentration(eclState.runspec().phases().active(Phase::BRINE));
setActivityModelSalt(eclState.getTableManager().actco2s());
// We only supported single pvt region for the co2-brine module
size_t numRegions = 1;
setNumRegions(numRegions);

1
src/opm/material/fluidsystems/blackoilpvt/Co2GasPvt.cpp
View File

@ -38,6 +38,7 @@ initFromState(const EclipseState& eclState, const Schedule&)
{
setEnableVaporizationWater(eclState.getSimulationConfig().hasVAPOIL() || eclState.getSimulationConfig().hasVAPWAT());
setActivityModelSalt(eclState.getTableManager().actco2s());
if (eclState.getTableManager().hasTables("PVDG") ||
!eclState.getTableManager().getPvtgTables().empty()) {

280
tests/material/test_binarycoefficients.cpp Normal file
View File

@ -0,0 +1,280 @@
// -*- 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 test makes sure that mandated API is adhered to by all component classes
*/
#include "config.h"
#define BOOST_TEST_MODULE Binarycoefficients
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/binarycoefficients/Brine_CO2.hpp>
#include <opm/material/components/SimpleHuDuanH2O.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp>
template<class Scalar>
bool close_at_tolerance(Scalar n1, Scalar n2, Scalar tolerance)
{
auto comp = boost::math::fpc::close_at_tolerance<Scalar>(tolerance);
return comp(n1, n2);
}
using Types = std::tuple<float,double>;
BOOST_AUTO_TEST_CASE_TEMPLATE(Brine_CO2, Scalar, Types)
{
using Evaluation = Opm::DenseAd::Evaluation<Scalar, 3>;
using H2O = Opm::SimpleHuDuanH2O<Scalar>;
using CO2 = Opm::CO2<Scalar>;
using BinaryCoeffBrineCO2 = Opm::BinaryCoeff::Brine_CO2<Scalar, H2O, CO2>;
// Molar mass of NaCl [kg/mol]
const Scalar MmNaCl = 58.44e-3;
// Extrapolate density?
const bool extrapolate = true;
// Activity model for salt
// 1 = Rumpf et al. (1994) as given in Spycher & Pruess (2005)
// 2 = Duan-Sun model as modified in Spycher & Pruess (2009)
// 3 = Duan-Sun model as given in Spycher & Pruess (2005)
const int activityModel = 2;
// Init. water gas mole fraction (which we don't care about here) and dissolved CO2 mole fraction
Evaluation xgH2O;
Evaluation xlCO2;
// Init pressure, temperature, and salinity variables
Evaluation T;
Evaluation p;
Evaluation s;
// Tolerance for Zhao et al (2015) data
const Scalar tol_zhao = 0.5e-2;
// Reference data Zhao et al., Geochimica et Cosmochimica Acta 149, 2015
// Experimental data in Table 3
// static constexpr Scalar xSol_Zhao[3][7] = {
// {1.245, 1.016, 0.859, 0.734, 0.623, 0.551, 0.498},
// {1.020, 0.836, 0.706, 0.618, 0.527, 0.469, 0.427},
// {1.001, 0.800, 0.647, 0.566, 0.498, 0.440, 0.390}
// };
// Spycher & Pruess (2009) (i.e., activityModel = 2 in our case) data given in Table 3
static constexpr Scalar xSol_Zhao[3][7] = {
{1.233, 1.012, 0.845, 0.716, 0.618, 0.541, 0.481},
{1.001, 0.819, 0.686, 0.588, 0.515, 0.463, 0.425},
{1.028, 0.824, 0.679, 0.577, 0.503, 0.450, 0.414}
};
// Temperature, pressure and salinity for Zhao et al (2015) data; Table 3
T = 323.15; // K
const int numT = 3;
p = 150e5; // Pa
const int numS = 7;
// Test against Zhao et al (2015) data
Evaluation salinity;
Evaluation mCO2;
for (int iT = 0; iT < numT; ++iT) {
// Salinity in mol/kg
s = 0.0;
for (int iS = 0; iS < numS; ++iS) {
// convert to mass fraction
salinity = 1 / ( 1 + 1 / (s * MmNaCl));
// Calculate solubility as mole fraction
BinaryCoeffBrineCO2::calculateMoleFractions(T, p, salinity, -1, xlCO2, xgH2O, activityModel, extrapolate);
// Convert to molality CO2
mCO2 = xlCO2 * (s + 55.508) / (1 - xlCO2);
// Compare to experimental data
// BOOST_CHECK_CLOSE(mCO2.value(), xSol_Zhao[iT][iS], tol);
BOOST_CHECK_MESSAGE(close_at_tolerance(mCO2.value(), xSol_Zhao[iT][iS], tol_zhao),
"relative difference between solubility {"<<mCO2.value()<<"} and Zhao et al. (2015) {"<<
xSol_Zhao[iT][iS]<<"} exceeds tolerance {"<<tol_zhao<<"} at (T, p, S) = ("<<T.value()<<", "<<
p.value()<<", "<<s.value()<<")");
// Increment salinity (mol/kg)
s += 1.0;
}
// Increment T
T += 50.0;
}
// Reference (experimental?) data Koschel et al., Fluid Phase Equilibria 247, 2006; Table 6
// Data for T = 323.1 K
static constexpr Scalar xSol_Koschel_1[2][4] = {
{0.0111, 0.0164, 0.0177, 0.0181},
{0.0083, 0.0114, 0.0124, 0.0137}
};
// Data for T = 373.1 K
static constexpr Scalar xSol_Koschel_2[2][3] = {
{0.006, 0.0112, 0.0161},
{0.0043, 0.008, 0.0112}
};
// Pressure data in Koschel et al. (2006) data
std::vector< std::vector<Evaluation> > p_1 = {
{5.1e6, 10.03e6, 14.38e6, 20.24e6},
{5.0e6, 10.04e6, 14.41e6, 20.24e6}
};
std::vector< std::vector<Evaluation> > p_2 = {
{5.07e6, 10.4e6, 19.14e6},
{5.04e6, 10.29e6, 19.02e6}
};
// Tolerance for Koschel et al. (2006) data
const Scalar tol_koschel = 14e-2;
// Test against reference data in Koschel et al. (2006)
s = 1.0;
for (int i = 0; i < 2; ++i) {
// Convert to mass fraction
salinity = 1 / ( 1 + 1 / (s * MmNaCl));
// First table
T = 323.1;
for (int j = 0; j < 4; ++j) {
// Get pressure
p = p_1[i][j];
// Calculate solubility as mole fraction
BinaryCoeffBrineCO2::calculateMoleFractions(T, p, salinity, -1, xlCO2, xgH2O, activityModel, extrapolate);
// Compare to experimental data
// BOOST_CHECK_CLOSE(xlCO2.value(), xSol_Koschel_1[i][j], tol);
BOOST_CHECK_MESSAGE(close_at_tolerance(xlCO2.value(), xSol_Koschel_1[i][j], tol_koschel),
"relative difference between solubility {"<<xlCO2.value()<<"} and Koschel et al. (2006) {"<<
xSol_Koschel_1[i][j]<<"} exceeds tolerance {"<<tol_koschel<<"} at (T, p, S) = ("<<T.value()<<", "<<
p.value()<<", "<<s.value()<<")");
}
// Second table
T = 373.1;
for (int j = 0; j < 3; ++j) {
// Get pressure
p = p_2[i][j];
// Calculate solubility as mole fraction
BinaryCoeffBrineCO2::calculateMoleFractions(T, p, salinity, -1, xlCO2, xgH2O, activityModel, extrapolate);
// Compare to experimental data
// BOOST_CHECK_CLOSE(xlCO2.value(), xSol_Koschel_2[i][j], tol);
BOOST_CHECK_MESSAGE(close_at_tolerance(xlCO2.value(), xSol_Koschel_2[i][j], tol_koschel),
"relative difference between solubility {"<<xlCO2.value()<<"} and Koschel et al. (2006) {"<<
xSol_Koschel_2[i][j]<<"} exceeds tolerance {"<<tol_koschel<<"} at (T, p, S) = ("<<T.value()<<", "<<
p.value()<<", "<<s.value()<<")");
}
// Increment salinity (to 3 mol/kg)
s += 2.0;
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(BrineDensityWithCO2, Scalar, Types)
{
using Evaluation = Opm::DenseAd::Evaluation<Scalar, 3>;
// Molar mass of NaCl [kg/mol]
const Scalar MmNaCl = 58.44e-3;
// Activity model for salt
// 1 = Rumpf et al. (1994) as given in Spycher & Pruess (2005)
// 2 = Duan-Sun model as modified in Spycher & Pruess (2009)
// 3 = Duan-Sun model as given in Spycher & Pruess (2005)
const int activityModel = 3;
// Tolerance for Yan et al. (2011) data
const Scalar tol_yan = 5e-3;
// Yan et al, Int. J. Greenh. Gas Control, 5, 2011; Table 4
static constexpr Scalar rho_Yan_1[3][6][3] = {{
{0.99722e3, 0.96370e3, 0.92928e3},
{1.00268e3, 0.96741e3, 0.93367e3},
{1.00528e3, 0.97062e3, 0.93760e3},
{1.00688e3, 0.97425e3, 0.94108e3},
{1.01293e3, 0.97962e3, 0.94700e3},
{1.01744e3, 0.98506e3, 0.95282e3}
},
{
{1.03116e3, 1.00026e3, 0.96883e3},
{1.03491e3, 1.00321e3, 0.97169e3},
{1.03968e3, 1.00667e3, 0.97483e3},
{1.04173e3, 1.00961e3, 0.97778e3},
{1.04602e3, 1.01448e3, 0.98301e3},
{1.05024e3, 1.01980e3, 0.98817e3}
},
{
{1.15824e3, 1.12727e3, 1.09559e3},
{1.16090e3, 1.12902e3, 1.10183e3},
{1.16290e3, 1.13066e3, 1.10349e3},
{1.16468e3, 1.13214e3, 1.10499e3},
{1.16810e3, 1.13566e3, 1.10882e3},
{1.17118e3, 1.13893e3, 1.11254e3}
}};
// Temperature, pressure and salinity for Yan et al (2011) data; Table 4
std::vector<Evaluation> T = {323.2, 373.2, 413.2}; // K
std::vector<Evaluation> p = {5e6, 10e6, 15e6, 20e6, 30e6, 40e6}; // Pa
std::vector<Scalar> s = {0.0, 1.0, 5.0};
// Test against Yan et al. (2011) data
Evaluation rs_sat;
Evaluation rho;
std::vector<Scalar> salinity = {0.0};
for (std::size_t iS = 0; iS < s.size(); ++iS) {
for (std::size_t iP = 0; iP < p.size(); ++iP) {
for (std::size_t iT = 0; iT < T.size(); ++iT) {
// Calculate salinity in mass fraction
salinity[0] = 1 / ( 1 + 1 / (s[iS] * MmNaCl));
// Instantiate BrineCo2Pvt class
Opm::BrineCo2Pvt<Scalar> brineCo2Pvt(salinity, activityModel);
// Calculate saturated Rs (dissolved CO2 in brine)
rs_sat = brineCo2Pvt.rsSat(/*regionIdx=*/0, T[iT], p[iP], Evaluation(salinity[0]));
// Calculate density of brine with dissolved CO2
rho = brineCo2Pvt.density(/*regionIdx=*/0, T[iT], p[iP], rs_sat, Evaluation(salinity[0]));
// Compare with data
BOOST_CHECK_MESSAGE(close_at_tolerance(rho.value(), rho_Yan_1[iS][iP][iT], tol_yan),
"relative difference between density {"<<rho.value()<<"} and Yan et al. (2011) {"<<
rho_Yan_1[iS][iP][iT]<<"} exceeds tolerance {"<<tol_yan<<"} at (T, p, S) = ("<<
T[iT].value()<<", "<<p[iP].value()<<", "<<s[iS]<<")");
}
}
}
}