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

Merge pull request #3704 from GitPaean/ptflash_simpletest

Browse Source 
various correction to prepare for the PTFlash simulation
This commit is contained in:
Bård Skaflestad
2023-10-13 17:56:03 +02:00
committed by GitHub
parent 49be8698cc f807281c15
commit 8f5ec261ec
12 changed files with 72 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
opm/material/components/C1.hpp
View File

@@ -75,9 +75,9 @@ public:
{ return 4.60e6; /* [N/m^2] */ }
/*!
* \brief Critical volume of \f$C_1\f$ [m2/kmol].
* \brief Critical volume of \f$C_1\f$ [m3/kmol].
*/
static Scalar criticalVolume() {return 9.863e-5; }
static Scalar criticalVolume() {return 9.863e-2; }
/*!
* \brief Acentric factor of \f$C_1\f$.

4
opm/material/components/C10.hpp
View File

@@ -75,9 +75,9 @@ public:
{ return 2.10e6; /* [N/m^2] */ }
/*!
* \brief Critical volume of \f$C_10\f$ [m2/kmol].
* \brief Critical volume of \f$C_10\f$ [m3/kmol].
*/
static Scalar criticalVolume() {return 6.098e-4; }
static Scalar criticalVolume() {return 6.098e-1; }
/*!
* \brief Acentric factor of \f$C_10\f$.

2
opm/material/components/SimpleCO2.hpp
View File

@@ -89,7 +89,7 @@ public:
* \brief Critical volume of \f$CO_2\f$ [m2/kmol].
*/
// Critical volume [m3/kmol]
static Scalar criticalVolume() {return 9.412e-5; }
static Scalar criticalVolume() {return 9.412e-2; }
/*!
* \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at the triple point of \f$CO_2\f$.

8
opm/material/constraintsolvers/PTFlash.hpp
View File

@@ -79,7 +79,6 @@ public:
template <class FluidState>
static void solve(FluidState& fluid_state,
const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
int spatialIdx,
std::string twoPhaseMethod,
Scalar tolerance = -1.,
int verbosity = 0)
@@ -104,7 +103,6 @@ public:
// Print header
if (verbosity >= 1) {
std::cout << "********" << std::endl;
std::cout << "Flash calculations on Cell " << spatialIdx << std::endl;
std::cout << "Inputs are K = [" << K << "], L = [" << L << "], z = [" << z << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl;
}
@@ -1127,7 +1125,7 @@ protected:
const bool newton_afterwards, const int verbosity)
{
// Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method.
const int maxIterations = newton_afterwards ? 3 : 10;
const int maxIterations = newton_afterwards ? 3 : 100;
// Store cout format before manipulation
std::ios_base::fmtflags f(std::cout.flags());
@@ -1169,7 +1167,7 @@ protected:
}
// Print iteration info
if (verbosity == 2 || verbosity == 4) {
if (verbosity >= 2) {
int prec = 5;
int fugWidth = (prec + 3);
int convWidth = prec + 9;
@@ -1184,7 +1182,7 @@ protected:
}
// Check convergence
if (convFugRatio.two_norm() < 1e-6){
if (convFugRatio.two_norm() < 1e-6) {
// Restore cout format
std::cout.flags(f);

38
opm/material/eos/PengRobinson.hpp
View File

@@ -54,7 +54,7 @@ namespace Opm {
* R. Reid, et al.: The Properties of Gases and Liquids,
* 4th edition, McGraw-Hill, 1987, pp. 42-44, 82
*/
template <class Scalar>
template <class Scalar, bool UseLegacy=true>
class PengRobinson
{
//! The ideal gas constant [Pa * m^3/mol/K]
@@ -205,27 +205,25 @@ public:
Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
Vm = VmCubic;
// find the extrema (if they are present)
Evaluation Vmin, Vmax, pmin, pmax;
if (findExtrema_(Vmin, Vmax,
pmin, pmax,
a, b, T))
{
if (isGasPhase)
Vm = std::max(Vmax, VmCubic);
else {
if (Vmin > 0)
Vm = std::min(Vmin, VmCubic);
else
Vm = VmCubic;
if (UseLegacy) {
// find the extrema (if they are present)
Evaluation Vmin, Vmax, pmin, pmax;
if (findExtrema_(Vmin, Vmax, pmin, pmax, a, b, T)) {
if (isGasPhase)
Vm = std::max(Vmax, VmCubic);
else {
if (Vmin > 0)
Vm = std::min(Vmin, VmCubic);
else
Vm = VmCubic;
}
} else {
// the EOS does not exhibit any physically meaningful
// extrema, and the fluid is critical...
Vm = VmCubic;
handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
}
}
else {
// the EOS does not exhibit any physically meaningful
// extrema, and the fluid is critical...
Vm = VmCubic;
handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase);
}
}
Valgrind::CheckDefined(Vm);

15
opm/material/eos/PengRobinsonMixture.hpp
View File

@@ -40,7 +40,6 @@ template <class Scalar, class StaticParameters>
class PengRobinsonMixture
{
enum { numComponents = StaticParameters::numComponents };
typedef ::Opm::PengRobinson<Scalar> PengRobinson;
// this class cannot be instantiated!
PengRobinsonMixture() = default;
@@ -53,20 +52,6 @@ class PengRobinsonMixture
static const Scalar w;
public:
/*!
* \brief Computes molar volumes where the Peng-Robinson EOS is
* true.
*
* \return Number of solutions.
*/
template <class MutableParams, class FluidState>
static int computeMolarVolumes(Scalar* Vm,
const MutableParams& params,
unsigned phaseIdx,
const FluidState& fs)
{
return PengRobinson::computeMolarVolumes(Vm, params, phaseIdx, fs);
}
/*!
* \brief Returns the fugacity coefficient of an individual

4
opm/material/fluidsystems/Co2BrineFluidSystem.hh
View File

@@ -34,13 +34,12 @@
#include <opm/material/fluidsystems/PTFlashParameterCache.hpp>
#include <opm/material/viscositymodels/LBC.hpp>
//#include <opm/material/viscositymodels/LBC.hpp>
namespace Opm {
/*!
* \ingroup FluidSystem
*
* \brief A two phase two component system, co2 brine
* \brief A two phase two component system with components co2 brine
*/
template<class Scalar>
@@ -64,7 +63,6 @@ namespace Opm {
template <class ValueType>
using ParameterCache = Opm::PTFlashParameterCache<ValueType, Co2BrineFluidSystem<Scalar>>;
using ViscosityModel = typename Opm::ViscosityModels<Scalar, Co2BrineFluidSystem<Scalar>>;
//using ViscosityModel = typename Opm::ViscosityModels<Scalar, Co2BrineFluidSystem<Scalar>>;
using PengRobinsonMixture = typename Opm::PengRobinsonMixture<Scalar, Co2BrineFluidSystem<Scalar>>;

2
opm/material/fluidsystems/PTFlashParameterCache.hpp
View File

@@ -48,7 +48,7 @@ class PTFlashParameterCache
{
using ThisType = PTFlashParameterCache<Scalar, FluidSystem>;
using ParentType = Opm::ParameterCacheBase<ThisType>;
using PengRobinson = Opm::PengRobinson<Scalar>;
using PengRobinson = Opm::PengRobinson<Scalar, false>; // false refer to discard some old code
enum { numPhases = FluidSystem::numPhases };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };

5
opm/material/fluidsystems/ThreeComponentFluidSystem.hh
View File

@@ -153,7 +153,7 @@ namespace Opm {
static const char* name[] = {"o", // oleic phase
"g"}; // gas phase
assert(0 <= phaseIdx && phaseIdx < 2);
assert(phaseIdx < 2);
return name[phaseIdx];
}
@@ -166,7 +166,7 @@ namespace Opm {
Comp2::name(),
};
assert(0 <= compIdx && compIdx < 3);
assert(compIdx < 3);
return name[compIdx];
}
@@ -196,6 +196,7 @@ namespace Opm {
// Use LBC method to calculate viscosity
LhsEval mu;
mu = ViscosityModel::LBC(fluidState, paramCache, phaseIdx);
return mu;
}

41
opm/material/viscositymodels/LBC.hpp
View File

@@ -49,54 +49,51 @@ public:
unsigned phaseIdx)
{
const Scalar MPa_atm = 0.101325;
const Scalar R = 8.3144598e-3;//Mj/kmol*K
const Scalar R = Opm::Constants<Scalar>::R;
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& rho = Opm::decay<LhsEval>(fluidState.density(phaseIdx));
const auto& P = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& Z = Opm::decay<LhsEval>(fluidState.compressFactor(phaseIdx));
LhsEval sumMm = 0.0;
LhsEval sumVolume = 0.0;
for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
const Scalar v_c = FluidSystem::criticalVolume(compIdx); // in m3/kmol
sumMm += x*Mm;
const Scalar v_c = FluidSystem::criticalVolume(compIdx) / 1000; // converting to m3/mol from m3/kmol
sumVolume += x*v_c;
}
LhsEval rho_pc = sumMm/sumVolume; //mixture pseudocritical density
LhsEval rho_r = rho/rho_pc;
LhsEval rho_pc = 1.0 / sumVolume;
LhsEval V = (R * T * Z)/P;
LhsEval rho = 1.0 / V;
LhsEval rho_r = rho / rho_pc;
LhsEval xsum_T_c = 0.0; //mixture pseudocritical temperature
LhsEval xsum_Mm = 0.0; //mixture molar mass
LhsEval xsum_p_ca = 0.0; //mixture pseudocritical pressure
LhsEval xsum_T_c = 0.0; // mixture pseudocritical temperature
LhsEval xsum_Mm = 0.0; // mixture molar mass
LhsEval xsum_p_ca = 0.0; // mixture pseudocritical pressure
for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6; // converting to Mpa from pascal
const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; // converting to kg/kmol from kg/mol;
const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
Scalar p_ca = p_c / MPa_atm;
xsum_T_c += x*T_c;
xsum_Mm += x*Mm;
xsum_p_ca += x*p_ca;
xsum_T_c += x * T_c;
xsum_Mm += x * Mm;
xsum_p_ca += x * p_ca;
}
LhsEval zeta_tot = Opm::pow(xsum_T_c / (Opm::pow(xsum_Mm,3.0) * Opm::pow(xsum_p_ca,4.0)),1./6);
LhsEval my0 = 0.0;
LhsEval sumxrM = 0.0;
for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6; // converting to Mpa from pa;
const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; // converting to kg/kmol from kg/mol;
const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
Scalar p_ca = p_c / MPa_atm;
Scalar zeta = std::pow(T_c / (std::pow(Mm,3.0) * std::pow(p_ca,4.0)),1./6);
LhsEval T_r = T/T_c;
LhsEval xrM = x * std::pow(Mm,0.5);
LhsEval mys = 0.0;
if (T_r <=1.5) {
if (T_r <= 1.5) {
mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta;
} else {
mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta;