opm-common/examples/co2brinepvt.cpp
2023-05-04 09:48:03 +02:00

195 lines
8.1 KiB
C++

// -*- 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
*
*/
#include "config.h"
#include <opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp>
#include <opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp>
#include <iostream>
template <class Co2Pvt>
double densityGas(const Co2Pvt& co2Pvt, const double p, const double T, const double Rv)
{
return co2Pvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
T,
p,
Rv,
/*Rvw=*/0.0) * co2Pvt.gasReferenceDensity(0);
}
template <class BrinePvt>
double densityBrine(const BrinePvt& brinePvt, const double p, const double T, const double Rs)
{
double bo = brinePvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
T,
p,
Rs);
return bo * (brinePvt.oilReferenceDensity(0) + Rs * brinePvt.gasReferenceDensity(0));
}
int main(int argc, char **argv)
{
bool help = false;
for (int i = 1; i < argc; ++i) {
std::string tmp = argv[i];
help = help || (tmp == "--h") || (tmp == "--help");
}
if (argc < 5 || help) {
std::cout << "USAGE:" << std::endl;
std::cout << "co2brinepvt <prop> <phase> <p> <T> <salinity> <rs> "<< std::endl;
std::cout << "prop = {density, invB, B, viscosity, rsSat, internalEnergy, enthalpy, diffusionCoefficient}" << std::endl;
std::cout << "phase = {CO2, brine}" << std::endl;
std::cout << "p: pressure in pascal" << std::endl;
std::cout << "T: temperature in kelvin" << std::endl;
std::cout << "salinity(optional): salt molality in mol/kg" << std::endl;
std::cout << "rs(optional): amount of dissolved CO2 in Brine in SM3/SM3" << std::endl;
std::cout << "OPTIONS:" << std::endl;
std::cout << "--h/--help Print help and exit." << std::endl;
std::cout << "DESCRIPTION:" << std::endl;
std::cout << "co2brinepvt computes PVT properties of a brine/co2 system " << std::endl;
std::cout << "for a given phase (oil or brine), pressure, temperature, salinity and rs." << std::endl;
std::cout << "The properties support are: density, the inverse phase formation volume factor (invB), viscosity, " << std::endl;
std::cout << "saturated dissolution factor (rsSat) " << std::endl;
std::cout << "See CO2STORE in the OPM manual for more details." << std::endl;
return EXIT_FAILURE;
}
std::string prop = argv[1];
std::string phase = argv[2];
double p = atof(argv[3]);
double T = atof(argv[4]);
double molality = 0.0;
double rs = 0.0;
double rv = 0.0; // only support 0.0 for now
if (argc > 5)
molality = atof(argv[5]);
if (argc > 6)
rs = atof(argv[6]);
const double MmNaCl = 58.44e-3; // molar mass of NaCl [kg/mol]
// convert to mass fraction
std::vector<double> salinity = {0.0};
if (molality > 0.0)
salinity[0] = 1 / ( 1 + 1 / (molality*MmNaCl));
Opm::BrineCo2Pvt<double> brineCo2Pvt(salinity);
Opm::Co2GasPvt<double> co2Pvt(salinity);
double value;
if (prop == "density") {
if (phase == "CO2") {
value = densityGas(co2Pvt, p, T, rv);
} else if (phase == "brine") {
value = densityBrine(brineCo2Pvt, p, T, rs);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
} else if (prop == "invB" || prop == "B") {
if (phase == "CO2") {
value = co2Pvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
T,
p,
rv,
/*Rvw=*/0.0);
} else if (phase == "brine") {
value = brineCo2Pvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
T,
p,
rs);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
if (prop == "B")
value = 1 / value;
} else if (prop == "viscosity") {
if (phase == "CO2") {
value = co2Pvt.viscosity(/*regionIdx=*/0,
T,
p,
rv,
/*Rvw=*/0.0);
} else if (phase == "brine") {
value = brineCo2Pvt.viscosity(/*regionIdx=*/0,
T,
p,
rs);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
} else if (prop == "rsSat") {
if (phase == "CO2") {
value = co2Pvt.saturatedWaterVaporizationFactor(/*regionIdx=*/0,
T,
p);
} else if (phase == "brine") {
value = brineCo2Pvt.saturatedGasDissolutionFactor(/*regionIdx=*/0,
T,
p);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
} else if (prop == "diffusionCoefficient") {
size_t comp_idx = 0; // not used
if (phase == "CO2") {
value = co2Pvt.diffusionCoefficient(T,p, comp_idx);
} else if (phase == "brine") {
value = brineCo2Pvt.diffusionCoefficient(T,p, comp_idx);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
} else if (prop == "internalEnergy") {
if (phase == "CO2") {
value = co2Pvt.internalEnergy(/*regionIdx=*/0 ,T,p, rv, 0.0);
} else if (phase == "brine") {
value = brineCo2Pvt.internalEnergy(/*regionIdx=*/0 ,T,p, rs);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
} else if (prop == "enthalpy") {
if (phase == "CO2") {
value = p / densityGas(co2Pvt, p, T, rv) + co2Pvt.internalEnergy(/*regionIdx=*/0 ,T,p, rv, 0.0);
} else if (phase == "brine") {
value = p / densityBrine(brineCo2Pvt, p, T, rs) + brineCo2Pvt.internalEnergy(/*regionIdx=*/0 ,T,p, rs);
} else {
throw std::runtime_error("phase " + phase + " not recognized. Use either CO2 or brine");
}
} else {
throw std::runtime_error("prop " + prop + " not recognized. "
+ "Use either density, visosity, invB, B, internalEnergy, enthalpy or diffusionCoefficient");
}
std::cout << value << std::endl;
return 0;
}