diff --git a/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp b/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp new file mode 100644 index 000000000..6a299597d --- /dev/null +++ b/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp @@ -0,0 +1,326 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::ChiParameterCache + */ +#ifndef OPM_CHI_PARAMETER_CACHE_HPP +#define OPM_CHI_PARAMETER_CACHE_HPP + +#include + +#include +#include + +#include +#include + +namespace Opm { + +/*! + * \ingroup Fluidsystems + * \brief Specifies the parameter cache used by the SPE-5 fluid system. + */ +template +class ChiParameterCache + : public Opm::ParameterCacheBase > +{ + using ThisType = ChiParameterCache; + using ParentType = Opm::ParameterCacheBase; + using PengRobinson = Opm::PengRobinson; + + enum { numPhases = FluidSystem::numPhases }; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx}; + +public: + //! The cached parameters for the oil phase + using OilPhaseParams = Opm::PengRobinsonParamsMixture; + //! The cached parameters for the gas phase + using GasPhaseParams = Opm::PengRobinsonParamsMixture; + + ChiParameterCache() + { + VmUpToDate_[oilPhaseIdx] = false; + Valgrind::SetUndefined(Vm_[oilPhaseIdx]); + VmUpToDate_[gasPhaseIdx] = false; + Valgrind::SetUndefined(Vm_[gasPhaseIdx]); + } + + //! \copydoc ParameterCacheBase::updatePhase + template + void updatePhase(const FluidState& fluidState, + unsigned phaseIdx, + int exceptQuantities = ParentType::None) + { + // if (phaseIdx != oilPhaseIdx) + // return; + + updateEosParams(fluidState, phaseIdx, exceptQuantities); + + // if we don't need to recalculate the molar volume, we exit + // here + // if (VmUpToDate_[phaseIdx]) + // return; + + // update the phase's molar volume + updateMolarVolume_(fluidState, phaseIdx); + } + + //! \copydoc ParameterCacheBase::updateSingleMoleFraction + template + void updateSingleMoleFraction(const FluidState& fluidState, + unsigned phaseIdx, + unsigned compIdx) + { + if (phaseIdx == oilPhaseIdx) + oilPhaseParams_.updateSingleMoleFraction(fluidState, compIdx); + if (phaseIdx == gasPhaseIdx) + gasPhaseParams_.updateSingleMoleFraction(fluidState, compIdx); + else + return; + + // update the phase's molar volume + updateMolarVolume_(fluidState, phaseIdx); + } + + /*! + * \brief The Peng-Robinson attractive parameter for a phase. + * + * \param phaseIdx The fluid phase of interest + */ + Scalar a(unsigned phaseIdx) const + { + switch (phaseIdx) + { + case oilPhaseIdx: return oilPhaseParams_.a(); + case gasPhaseIdx: return gasPhaseParams_.a(); + default: + throw std::logic_error("The a() parameter is only defined for " + "oil phases"); + }; + } + + /*! + * \brief The Peng-Robinson covolume for a phase. + * + * \param phaseIdx The fluid phase of interest + */ + Scalar b(unsigned phaseIdx) const + { + switch (phaseIdx) + { + case oilPhaseIdx: return oilPhaseParams_.b(); + case gasPhaseIdx: return gasPhaseParams_.b(); + default: + throw std::logic_error("The b() parameter is only defined for " + "oil phase"); + }; + } + + /*! + * \brief The Peng-Robinson attractive parameter for a pure + * component given the same temperature and pressure of the + * phase. + * + * \param phaseIdx The fluid phase of interest + * \param compIdx The component phase of interest + */ + Scalar aPure(unsigned phaseIdx, unsigned compIdx) const + { + switch (phaseIdx) + { + case oilPhaseIdx: return oilPhaseParams_.pureParams(compIdx).a(); + case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a(); + default: + throw std::logic_error("The a() parameter is only defined for " + "oil phase"); + }; + } + + /*! + * \brief The Peng-Robinson covolume for a pure component given + * the same temperature and pressure of the phase. + * + * \param phaseIdx The fluid phase of interest + * \param compIdx The component phase of interest + */ + Scalar bPure(unsigned phaseIdx, unsigned compIdx) const + { + switch (phaseIdx) + { + case oilPhaseIdx: return oilPhaseParams_.pureParams(compIdx).b(); + case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b(); + default: + throw std::logic_error("The b() parameter is only defined for " + "oil phase"); + }; + } + + /*! + * \brief Returns the molar volume of a phase [m^3/mol] + * + * \param phaseIdx The fluid phase of interest + */ + Scalar molarVolume(unsigned phaseIdx) const + { assert(VmUpToDate_[phaseIdx]); return Vm_[phaseIdx]; } + + + /*! + * \brief Returns the Peng-Robinson mixture parameters for the oil + * phase. + */ + const OilPhaseParams& oilPhaseParams() const + { return oilPhaseParams_; } + + /*! + * \brief Returns the Peng-Robinson mixture parameters for the gas + * phase. + */ + const GasPhaseParams& gasPhaseParams() const + // { throw std::invalid_argument("gas phase does not exist");} + { return gasPhaseParams_; } + + /*! + * \brief Update all parameters required by the equation of state to + * calculate some quantities for the phase. + * + * \param fluidState The representation of the thermodynamic system of interest. + * \param phaseIdx The index of the fluid phase of interest. + * \param exceptQuantities The quantities of the fluid state that have not changed since the last update. + */ + template + void updateEosParams(const FluidState& fluidState, + unsigned phaseIdx, + int exceptQuantities = ParentType::None) + { + // if (phaseIdx != oilPhaseIdx) + // return; + + if (!(exceptQuantities & ParentType::Temperature)) + { + updatePure_(fluidState, phaseIdx); + updateMix_(fluidState, phaseIdx); + VmUpToDate_[phaseIdx] = false; + } + else if (!(exceptQuantities & ParentType::Composition)) + { + updateMix_(fluidState, phaseIdx); + VmUpToDate_[phaseIdx] = false; + } + else if (!(exceptQuantities & ParentType::Pressure)) { + VmUpToDate_[phaseIdx] = false; + } + } + +protected: + /*! + * \brief Update all parameters of a phase which only depend on + * temperature and/or pressure. + * + * This usually means the parameters for the pure components. + */ + template + void updatePure_(const FluidState& fluidState, unsigned phaseIdx) + { + Scalar T = fluidState.temperature(phaseIdx); + Scalar p = fluidState.pressure(phaseIdx); + + switch (phaseIdx) + { + case oilPhaseIdx: oilPhaseParams_.updatePure(T, p); break; + case gasPhaseIdx: gasPhaseParams_.updatePure(T, p); break; + } + } + + /*! + * \brief Update all parameters of a phase which depend on the + * fluid composition. It is assumed that updatePure() has + * been called before this method. + * + * Here, the mixing rule kicks in. + */ + template + void updateMix_(const FluidState& fluidState, unsigned phaseIdx) + { + Valgrind::CheckDefined(fluidState.averageMolarMass(phaseIdx)); + switch (phaseIdx) + { + case oilPhaseIdx: + oilPhaseParams_.updateMix(fluidState); + break; + case gasPhaseIdx: + gasPhaseParams_.updateMix(fluidState); + break; + } + } + + template + void updateMolarVolume_(const FluidState& fluidState, + unsigned phaseIdx) + { + VmUpToDate_[phaseIdx] = true; + + // calculate molar volume of the phase (we will need this for the + // fugacity coefficients and the density anyway) + switch (phaseIdx) { + case gasPhaseIdx: { + // calculate molar volumes for the given composition. although + // this isn't a Peng-Robinson parameter strictly speaking, the + // molar volume appears in basically every quantity the fluid + // system can get queried, so it is okay to calculate it + // here... + Vm_[gasPhaseIdx] = + PengRobinson::computeMolarVolume(fluidState, + *this, + phaseIdx, + /*isGasPhase=*/true); + break; + } + case oilPhaseIdx: { + // calculate molar volumes for the given composition. although + // this isn't a Peng-Robinson parameter strictly speaking, the + // molar volume appears in basically every quantity the fluid + // system can get queried, so it is okay to calculate it + // here... + Vm_[oilPhaseIdx] = + PengRobinson::computeMolarVolume(fluidState, + *this, + phaseIdx, + /*isGasPhase=*/false); + + break; + } + }; + } + + bool VmUpToDate_[numPhases]; + Scalar Vm_[numPhases]; + + OilPhaseParams oilPhaseParams_; + GasPhaseParams gasPhaseParams_; +}; + +} // namespace Opm + +#endif diff --git a/opm/material/fluidsystems/chifluid/LBCviscosity.hpp b/opm/material/fluidsystems/chifluid/LBCviscosity.hpp new file mode 100644 index 000000000..8fdc75c00 --- /dev/null +++ b/opm/material/fluidsystems/chifluid/LBCviscosity.hpp @@ -0,0 +1,212 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::LBCviscosity + */ + +#ifndef LBC_VISCOSITY_HPP +#define LBC_VISCOSITY_HPP + +#include +#include + +namespace Opm +{ +template +class LBCviscosity +{ + +public: + + // Standard LBC model. (Lohrenz, Bray & Clark: "Calculating Viscosities of Reservoir + // fluids from Their Compositions", JPT 16.10 (1964). + template + static LhsEval LBC(const FluidState& fluidState, + const Params& /*paramCache*/, + unsigned phaseIdx) + { + const Scalar MPa_atm = 0.101325; + const Scalar R = 8.3144598e-3;//Mj/kmol*K + const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); + const auto& rho = Opm::decay(fluidState.density(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(fluidState.moleFraction(phaseIdx, compIdx)); + const Scalar v_c = FluidSystem::criticalVolume(compIdx); // in m3/kmol + sumMm += x*Mm; + sumVolume += x*v_c; + } + + LhsEval rho_pc = sumMm/sumVolume; //mixture pseudocritical density + 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 + 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(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; + } + 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& T_c = FluidSystem::criticalTemperature(compIdx); + const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol; + const auto& x = Opm::decay(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) { + 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; + } + my0 += xrM*mys; + sumxrM += xrM; + } + my0 /= sumxrM; + + std::vector LBC = {0.10230, + 0.023364, + 0.058533, + -0.040758, // trykkfeil i 1964-artikkel: -0.40758 + 0.0093324}; + + LhsEval sumLBC = 0.0; + for (int i = 0; i < 5; ++i) { + sumLBC += Opm::pow(rho_r,i)*LBC[i]; + } + + return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot)/1e3; // mPas-> Pas + } + + + // Improved LBC model for CO2 rich mixtures. (Lansangan, Taylor, Smith & Kovarik - 1993) + template + static LhsEval LBCmod(const FluidState& fluidState, + const Params& /*paramCache*/, + unsigned phaseIdx) + { + const Scalar MPa_atm = 0.101325; + const Scalar R = 8.3144598e-3;//Mj/kmol*K + const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); + const auto& rho = Opm::decay(fluidState.density(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(fluidState.moleFraction(phaseIdx, compIdx)); + const Scalar v_c = FluidSystem::criticalVolume(compIdx); // in m3/kmol + sumMm += x*Mm; + sumVolume += x*v_c; + } + + LhsEval rho_pc = sumMm/sumVolume; //mixture pseudocritical density + LhsEval rho_r = rho/rho_pc; + + LhsEval xxT_p = 0.0; // x*x*T_c/p_c + LhsEval xxT2_p = 0.0; // x*x*T^2_c/p_c + for (unsigned i_compIdx = 0; i_compIdx < FluidSystem::numComponents; ++i_compIdx) { + const Scalar& T_c_i = FluidSystem::criticalTemperature(i_compIdx); + const Scalar& p_c_i = FluidSystem::criticalPressure(i_compIdx)/1e6; // in Mpa; + const auto& x_i = Opm::decay(fluidState.moleFraction(phaseIdx, i_compIdx)); + for (unsigned j_compIdx = 0; j_compIdx < FluidSystem::numComponents; ++j_compIdx) { + const Scalar& T_c_j = FluidSystem::criticalTemperature(j_compIdx); + const Scalar& p_c_j = FluidSystem::criticalPressure(j_compIdx)/1e6; // in Mpa; + const auto& x_j = Opm::decay(fluidState.moleFraction(phaseIdx, j_compIdx)); + + const Scalar T_c_ij = std::sqrt(T_c_i*T_c_j); + const Scalar p_c_ij = 8.0*T_c_ij / Opm::pow(Opm::pow(T_c_i/p_c_i,1.0/3)+Opm::pow(T_c_j/p_c_j,1.0/3),3); + + xxT_p += x_i*x_j*T_c_ij/p_c_ij; + xxT2_p += x_i*x_j*T_c_ij*T_c_ij/p_c_ij; + } + } + + const LhsEval T_pc = xxT2_p/xxT_p; //mixture pseudocritical temperature + const LhsEval p_pc = T_pc/xxT_p; //mixture pseudocritical pressure + + LhsEval p_pca = p_pc / MPa_atm; + LhsEval zeta_tot = Opm::pow(T_pc / (Opm::pow(sumMm,3.0) * Opm::pow(p_pca,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& T_c = FluidSystem::criticalTemperature(compIdx); + const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol; + const auto& x = Opm::decay(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) { + 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; + } + my0 += xrM*mys; + sumxrM += xrM; + } + my0 /= sumxrM; + + std::vector LBC = {0.10230, + 0.023364, + 0.058533, + -0.040758, // trykkfeil i 1964-artikkel: -0.40758 + 0.0093324}; + + LhsEval sumLBC = 0.0; + for (int i = 0; i < 5; ++i) { + sumLBC += Opm::pow(rho_r,i)*LBC[i]; + } + + return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot -1.8366e-8*Opm::pow(rho_r,13.992))/1e3; // mPas-> Pas + } +}; + +} // namespace Opm + +#endif // LBC_VISCOSITY_HPP diff --git a/opm/material/fluidsystems/chifluid/chiwoms.h b/opm/material/fluidsystems/chifluid/chiwoms.h new file mode 100644 index 000000000..3c1db10be --- /dev/null +++ b/opm/material/fluidsystems/chifluid/chiwoms.h @@ -0,0 +1,31 @@ +#ifndef __CHIWOMS_H__ +#define __CHIWOMS_H__ + +// values are from the paper "Field-scale implications of density-driven +// convection in CO2-EOR reservoirs", to be presented at the Fifth CO2 +// Geological Storage Workshop, at 21–23 November 2018, in Utrecht, +// The Netherlands. + +constexpr double TEMPERATURE = 80; /* degree Celsius */ +constexpr double GRAVITYFACTOR = 1; /* fraction og gravity */ +constexpr double MIN_PRES = 75; /* bars */ +const double MAX_PRES = 220; /* bars */ +constexpr double SIM_TIME = 1; /* days */ +constexpr double Y_SIZE = 1.0; /* meter */ +constexpr double X_SIZE = 1.0; /* meter */ +constexpr double Z_SIZE = 1.0; /* meter */ +const unsigned NX = 5; /* number of cells x-dir */ +const unsigned NY = 5; /* number of cells y-dir */ +const unsigned NZ = 5; /* number of cells z-dir */ +const double POROSITY = 0.2; /* non-dimensional */ +const double PERMEABILITY = 100; /* milli-Darcy */ +const double DIFFUSIVITY = 1e-9; /* square meter per second */ +const double MFCOMP0 = 0.9999999; +const double MFCOMP1 = 0.0000001; +const double MFCOMP2 = 0.0; +constexpr double INFLOW_RATE = -1e-4; /* unit kg/s ? */ + +/* "random" fields will be equal as long as this is set the same */ +const double SEED = 5163166242092481088; + +#endif /* __CHIWOMS_H__ */ diff --git a/opm/material/fluidsystems/chifluid/components.hh b/opm/material/fluidsystems/chifluid/components.hh new file mode 100644 index 000000000..2a2b5fb0e --- /dev/null +++ b/opm/material/fluidsystems/chifluid/components.hh @@ -0,0 +1,221 @@ +#ifndef COMPONENTS_HH +#define COMPONENTS_HH + +#include "chiwoms.h" + +#include +#include +#include +#include +#include + +namespace Opm { +/*! + * \ingroup Components + * + * \brief A simple representation of linear octane + * + * \tparam Scalar The type used for scalar values + */ +template +class Octane : public Opm::Component > +{ +public: + /// Chemical name + static const char* name() { return "C8"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.11423; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 568.7; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 2.49e6; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.398; } + + // Critical volume [m3/kmol] (same as [L/mol]) + static Scalar criticalVolume() {return 4.92e-1; } +}; + +template +class NDekane : public Opm::Component > +{ +public: + /// Chemical name + static const char* name() { return "C10"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.1423; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 617.7; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 2.103e6; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.4884; } + + // Critical volume [m3/kmol] (same as [L/mol]) + static Scalar criticalVolume() {return 6.0976e-1; } +}; + +template +class Methane : public Opm::Component > +{ +public: + /// Chemical name + static const char* name() { return "CH4"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.0160; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 190.5640; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 4.599e6; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.0114; } + + // Critical volume [m3/kmol] + static Scalar criticalVolume() {return 9.8628e-2; } +}; + + +template +class Hydrogen : public Opm::Component > +{ +public: + /// Chemical name + static const char* name() { return "H2"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.0020156; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 33.2; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 1.297e6; } + + /// Acentric factor + static Scalar acentricFactor() { return -0.22; } + + // Critical volume [m3/kmol] + static Scalar criticalVolume() {return 6.45e-2; } + +}; + +template +class Nitrogen : public Opm::Component > +{ +public: + /// Chemical name + static const char* name() { return "N2"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.0280134; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 126.192; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 3.3958e6; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.039; } + + // Critical volume [m3/kmol] + static Scalar criticalVolume() {return 8.94e-2; } +}; + +template +class Water : public Opm::Component > +{ +public: + /// Chemical name + static const char* name() { return "H20"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.01801528; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 647; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 22.064e6; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.344; } + + // Critical volume [m3/kmol] + static Scalar criticalVolume() {return 5.595e-2; } +}; + +template +class ChiwomsCO2 : public Opm::SimpleCO2 +{ +public: + /// Chemical name + static const char* name() { return "CO2"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.0440095; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 304.1; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 7.38e6; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.225; } + + // Critical volume [m3/kmol] + static Scalar criticalVolume() {return 9.4118e-2; } +}; + +template +class ChiwomsBrine : public Opm::H2O +{ +public: + /// Chemical name + static const char* name() { return "H20-NaCl"; } + + /// Molar mass in \f$\mathrm{[kg/mol]}\f$ + static Scalar molarMass() { return 0.0180158; } + + /// Critical temperature in \f$\mathrm[K]}\f$ + static Scalar criticalTemperature() { return 647.096; } + + /// Critical pressure in \f$\mathrm[Pa]}\f$ + static Scalar criticalPressure() { return 2.21e7; } + + /// Acentric factor + static Scalar acentricFactor() { return 0.344; } + + // Critical volume [m3/kmol] + static Scalar criticalVolume() {return 5.595e-2; } +}; + +struct EOS +{ + template + static LhsEval oleic_enthalpy(LhsEval T, LhsEval p, LhsEval x) { + return 0; + } + + template + static LhsEval aqueous_enthalpy(LhsEval T, LhsEval p, LhsEval x) { + return 0; + } +}; + +} // namespace opm + +#endif // COMPONENTS_HH diff --git a/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh b/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh new file mode 100644 index 000000000..e0a91c619 --- /dev/null +++ b/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh @@ -0,0 +1,360 @@ +#ifndef TWOPHASEFLUIDSYSTEM_HH +#define TWOPHASEFLUIDSYSTEM_HH + +#include "components.hh" +#include "chiwoms.h" +#include "LBCviscosity.hpp" + +#include +#include +#include // invalid_argument +#include +#include +#include +#include // mt19937, normal_distribution +#include // epsilon +#include // boost::format + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include "ChiParameterCache.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + +namespace Opm { +/*! + * \ingroup Fluidsystems + * + * \brief A two-phase fluid system with three components + */ +template +class TwoPhaseThreeComponentFluidSystem + : public Opm::BaseFluidSystem > +{ + using ThisType = TwoPhaseThreeComponentFluidSystem; + using Base = Opm::BaseFluidSystem; + using PengRobinson = typename Opm::PengRobinson; + using PengRobinsonMixture = typename Opm::PengRobinsonMixture; + using LBCviscosity = typename Opm::LBCviscosity; + using H2O = typename Opm::H2O; + using Brine = typename Opm::Brine; + +public: + //! \copydoc BaseFluidSystem::ParameterCache + //template + //using ParameterCache = Opm::NullParameterCache; + + //! \copydoc BaseFluidSystem::ParameterCache + template + using ParameterCache = Opm::ChiParameterCache; + + /**************************************** + * Fluid phase related static parameters + ****************************************/ + + //! \copydoc BaseFluidSystem::numPhases + static const int numPhases = 2; + + //! Index of the liquid phase + static const int oilPhaseIdx = 0; + static const int gasPhaseIdx = 1; + + //! \copydoc BaseFluidSystem::phaseName + static const char* phaseName(unsigned phaseIdx) + { + static const char* name[] = {"o", // oleic phase + "g"}; // gas phase + + assert(0 <= phaseIdx && phaseIdx < numPhases); + return name[phaseIdx]; + } + + //! \copydoc BaseFluidSystem::isIdealMixture + static bool isIdealMixture(unsigned phaseIdx) + { + if (phaseIdx == oilPhaseIdx) + return false; + + // CO2 have associative effects + return true; + } + + + /**************************************** + * Component related static parameters + ****************************************/ + + //! \copydoc BaseFluidSystem::numComponents + static const int numComponents = 2; // Comp0, Comp1 and Comp2 + + //! first comp idx + static const int Comp0Idx = 0; + + //! second comp idx + static const int Comp1Idx = 1; + + // TODO: make this a loop over choises in chiwoms.hh + // using Comp0 = Opm::Methane; + using Comp0 = Opm::ChiwomsBrine; + using Comp1 = Opm::ChiwomsCO2; + + static void init(Scalar minT = 273.15, + Scalar maxT = 373.15, + Scalar minP = 1e4, + Scalar maxP = 100e6) + { + Opm::PengRobinsonParamsMixture prParams; + + // find envelopes of the 'a' and 'b' parameters for the range + // minT <= T <= maxT and minP <= p <= maxP. For + // this we take advantage of the fact that 'a' and 'b' for + // mixtures is just a convex combination of the attractive and + // repulsive parameters of the pure components + + Scalar minA = 1e30, maxA = -1e30; + Scalar minB = 1e30, maxB = -1e30; + + prParams.updatePure(minT, minP); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + minA = std::min(prParams.pureParams(compIdx).a(), minA); + maxA = std::max(prParams.pureParams(compIdx).a(), maxA); + minB = std::min(prParams.pureParams(compIdx).b(), minB); + maxB = std::max(prParams.pureParams(compIdx).b(), maxB); + }; + + prParams.updatePure(maxT, minP); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + minA = std::min(prParams.pureParams(compIdx).a(), minA); + maxA = std::max(prParams.pureParams(compIdx).a(), maxA); + minB = std::min(prParams.pureParams(compIdx).b(), minB); + maxB = std::max(prParams.pureParams(compIdx).b(), maxB); + }; + + prParams.updatePure(minT, maxP); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + minA = std::min(prParams.pureParams(compIdx).a(), minA); + maxA = std::max(prParams.pureParams(compIdx).a(), maxA); + minB = std::min(prParams.pureParams(compIdx).b(), minB); + maxB = std::max(prParams.pureParams(compIdx).b(), maxB); + }; + + prParams.updatePure(maxT, maxP); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + minA = std::min(prParams.pureParams(compIdx).a(), minA); + maxA = std::max(prParams.pureParams(compIdx).a(), maxA); + minB = std::min(prParams.pureParams(compIdx).b(), minB); + maxB = std::max(prParams.pureParams(compIdx).b(), maxB); + }; + // PengRobinson::init(/*aMin=*/minA, /*aMax=*/maxA, /*na=*/100, + // /*bMin=*/minB, /*bMax=*/maxB, /*nb=*/200); + } + + //! \copydoc BaseFluidSystem::componentName + static const char* componentName(unsigned compIdx) + { + static const char* name[] = { + Comp0::name(), + Comp1::name(), + }; + + assert(0 <= compIdx && compIdx < numComponents); + return name[compIdx]; + } + + //! \copydoc BaseFluidSystem::molarMass + static Scalar molarMass(unsigned compIdx) + { + return (compIdx == Comp0Idx) + ? Comp0::molarMass() + : (compIdx == Comp1Idx) + ? Comp1::molarMass() + : throw std::invalid_argument("Molar mass component index"); + } + + /*! + * \brief Critical temperature of a component [K]. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar criticalTemperature(unsigned compIdx) + { + return (compIdx == Comp0Idx) + ? Comp0::criticalTemperature() + : (compIdx == Comp1Idx) + ? Comp1::criticalTemperature() + : throw std::invalid_argument("Critical temperature component index"); + } + + /*! + * \brief Critical pressure of a component [Pa]. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar criticalPressure(unsigned compIdx) + { + return (compIdx == Comp0Idx) + ? Comp0::criticalPressure() + : (compIdx == Comp1Idx) + ? Comp1::criticalPressure() + : throw std::invalid_argument("Critical pressure component index"); + } + + /*! + * \brief Critical volume of a component [m3]. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar criticalVolume(unsigned compIdx) + { + return (compIdx == Comp0Idx) + ? Comp0::criticalVolume() + : (compIdx == Comp1Idx) + ? Comp1::criticalVolume() + : throw std::invalid_argument("Critical volume component index"); + } + + /*! + * \brief The acentric factor of a component []. + * + * \copydetails Doxygen::compIdxParam + */ + static Scalar acentricFactor(unsigned compIdx) + { + return (compIdx == Comp0Idx) + ? Comp0::acentricFactor() + : (compIdx == Comp1Idx) + ? Comp1::acentricFactor() + : throw std::invalid_argument("Molar mass component index"); + } + + /**************************************** + * thermodynamic relations + ****************************************/ + + /*! + * \copydoc BaseFluidSystem::density + */ + template + static LhsEval density(const FluidState& fluidState, + const ParameterCache& paramCache, + unsigned phaseIdx) + { + + LhsEval dens; + if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) { + // paramCache.updatePhase(fluidState, phaseIdx); + dens = fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx); + } + return dens; + + } + + //! \copydoc BaseFluidSystem::viscosity + template + static LhsEval viscosity(const FluidState& fluidState, + const ParameterCache& paramCache, + unsigned phaseIdx) + { + // Use LBC method to calculate viscosity + LhsEval mu; + // if (phaseIdx == gasPhaseIdx) { + mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx); + // } + // else { + // const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); + // const auto& p = Opm::decay(fluidState.pressure(0)); + // mu = Brine::liquidViscosity(T, p); + // } + return mu; + + } + + //! \copydoc BaseFluidSystem::enthalpy + template + static LhsEval enthalpy(const FluidState& fluidState, + const ParameterCache& /*paramCache*/, + unsigned phaseIdx) + { + const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); + const auto& p = Opm::decay(fluidState.pressure(phaseIdx)); + const auto& x = Opm::decay(fluidState.moleFraction(phaseIdx, Comp1Idx)); + + if(phaseIdx == oilPhaseIdx) { + return EOS::oleic_enthalpy(T, p, x); //TODO + } + else { + return EOS::aqueous_enthalpy(T, p, x); //TODO + } + } + + //! \copydoc BaseFluidSystem::fugacityCoefficient + template + static LhsEval fugacityCoefficient(const FluidState& fluidState, + const ParameterCache& paramCache, + unsigned phaseIdx, + unsigned compIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + assert(0 <= compIdx && compIdx < numComponents); + + Scalar phi = Opm::getValue( + PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); + return phi; + + + throw std::invalid_argument("crap!"); + } + + //! \copydoc BaseFluidSystem::diffusionCoefficient + template + static LhsEval diffusionCoefficient(const FluidState& /*fluidState*/, + const ParameterCache& /*paramCache*/, + unsigned /*phaseIdx*/, + unsigned /*compIdx*/) + { + return DIFFUSIVITY; + } + + /*! + * \brief Returns the interaction coefficient for two components. + *. + */ + static Scalar interactionCoefficient(unsigned comp1Idx, unsigned comp2Idx) + { + return 0.0; //-0.101;//0.1089; + } + +}; + +};//namespace opm + +#endif // TWOPHASEFLUIDSYSTEM_HH diff --git a/tests/test_chiflash.cpp b/tests/test_chiflash.cpp index 743501b0a..a8e3d99f2 100644 --- a/tests/test_chiflash.cpp +++ b/tests/test_chiflash.cpp @@ -29,6 +29,7 @@ #include "config.h" #include +#include #include #include @@ -40,471 +41,15 @@ #include -template -void createSurfaceGasFluidSystem(FluidState& gasFluidState) -{ - static const int gasPhaseIdx = FluidSystem::gasPhaseIdx; - - // temperature - gasFluidState.setTemperature(273.15 + 20); - - // gas pressure - gasFluidState.setPressure(gasPhaseIdx, 1e5); - - // gas saturation - gasFluidState.setSaturation(gasPhaseIdx, 1.0); - - // gas composition: mostly methane, a bit of propane - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::H2OIdx, 0.0); - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C1Idx, 0.94); - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C3Idx, 0.06); - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C6Idx, 0.00); - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C10Idx, 0.00); - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C15Idx, 0.00); - gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C20Idx, 0.00); - - // gas density - typename FluidSystem::template ParameterCache paramCache; - paramCache.updatePhase(gasFluidState, gasPhaseIdx); - gasFluidState.setDensity(gasPhaseIdx, - FluidSystem::density(gasFluidState, paramCache, gasPhaseIdx)); -} - -template -Scalar computeSumxg(FluidState& resultFluidState, - const FluidState& prestineFluidState, - const FluidState& gasFluidState, - Scalar additionalGas) -{ - static const int oilPhaseIdx = FluidSystem::oilPhaseIdx; - static const int gasPhaseIdx = FluidSystem::gasPhaseIdx; - static const int numComponents = FluidSystem::numComponents; - - typedef Dune::FieldVector ComponentVector; - typedef Opm::NcpFlash Flash; - - resultFluidState.assign(prestineFluidState); - - // add a bit of additional gas components - ComponentVector totalMolarities; - for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++ compIdx) - totalMolarities = - prestineFluidState.molarity(oilPhaseIdx, compIdx) - + additionalGas*gasFluidState.moleFraction(gasPhaseIdx, compIdx); - - // "flash" the modified fluid state - typename FluidSystem::ParameterCache paramCache; - Flash::solve(resultFluidState, totalMolarities); - - Scalar sumxg = 0; - for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) - sumxg += resultFluidState.moleFraction(gasPhaseIdx, compIdx); - - return sumxg; -} - -template -void makeOilSaturated(FluidState& fluidState, const FluidState& gasFluidState) -{ - static const int gasPhaseIdx = FluidSystem::gasPhaseIdx; - - FluidState prestineFluidState; - prestineFluidState.assign(fluidState); - - Scalar sumxg = 0; - for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) - sumxg += fluidState.moleFraction(gasPhaseIdx, compIdx); - - // Newton method - Scalar tol = 1e-8; - Scalar additionalGas = 0; // [mol] - for (int i = 0; std::abs(sumxg - 1) > tol; ++i) { - if (i > 50) - throw std::runtime_error("Newton method did not converge after 50 iterations"); - - Scalar eps = std::max(1e-8, additionalGas*1e-8); - - Scalar f = 1 - computeSumxg(prestineFluidState, - fluidState, - gasFluidState, - additionalGas); - Scalar fStar = 1 - computeSumxg(prestineFluidState, - fluidState, - gasFluidState, - additionalGas + eps); - Scalar fPrime = (fStar - f)/eps; - - additionalGas -= f/fPrime; - }; -} - -template -void guessInitial(FluidState& fluidState, unsigned phaseIdx) -{ - if (phaseIdx == FluidSystem::gasPhaseIdx) { - fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.74785); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.0121364); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.00606028); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.00268136); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.000204256); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 8.78291e-06); - } - else if (phaseIdx == FluidSystem::oilPhaseIdx) { - fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.50); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.03); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.07); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.20); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.15); - fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 0.05); - } - else { - assert(phaseIdx == FluidSystem::waterPhaseIdx); - } -} - -template -Scalar bringOilToSurface(FluidState& surfaceFluidState, Scalar alpha, const FluidState& reservoirFluidState, bool guessInitial) -{ - enum { - numPhases = FluidSystem::numPhases, - waterPhaseIdx = FluidSystem::waterPhaseIdx, - gasPhaseIdx = FluidSystem::gasPhaseIdx, - oilPhaseIdx = FluidSystem::oilPhaseIdx, - - numComponents = FluidSystem::numComponents - }; - - typedef Opm::NcpFlash Flash; - typedef Opm::ThreePhaseMaterialTraits MaterialTraits; - typedef Opm::LinearMaterial MaterialLaw; - typedef typename MaterialLaw::Params MaterialLawParams; - typedef Dune::FieldVector ComponentVector; - - const Scalar refPressure = 1.0135e5; // [Pa] - - // set the parameters for the capillary pressure law - MaterialLawParams matParams; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - matParams.setPcMinSat(phaseIdx, 0.0); - matParams.setPcMaxSat(phaseIdx, 0.0); - } - matParams.finalize(); - - // retieve the global volumetric component molarities - surfaceFluidState.setTemperature(273.15 + 20); - - ComponentVector molarities; - for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) - molarities[compIdx] = reservoirFluidState.molarity(oilPhaseIdx, compIdx); - - if (guessInitial) { - // we start at a fluid state with reservoir oil. - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) { - surfaceFluidState.setMoleFraction(phaseIdx, - compIdx, - reservoirFluidState.moleFraction(phaseIdx, compIdx)); - } - surfaceFluidState.setDensity(phaseIdx, reservoirFluidState.density(phaseIdx)); - surfaceFluidState.setPressure(phaseIdx, reservoirFluidState.pressure(phaseIdx)); - surfaceFluidState.setSaturation(phaseIdx, 0.0); - } - surfaceFluidState.setSaturation(oilPhaseIdx, 1.0); - surfaceFluidState.setSaturation(gasPhaseIdx, 1.0 - surfaceFluidState.saturation(oilPhaseIdx)); - } - - typename FluidSystem::template ParameterCache paramCache; - paramCache.updateAll(surfaceFluidState); - - // increase volume until we are at surface pressure. use the - // newton method for this - ComponentVector tmpMolarities; - for (int i = 0;; ++i) { - if (i >= 20) - throw Opm::NumericalIssue("Newton method did not converge after 20 iterations"); - - // calculate the deviation from the standard pressure - tmpMolarities = molarities; - tmpMolarities /= alpha; - Flash::template solve(surfaceFluidState, matParams, paramCache, tmpMolarities); - Scalar f = surfaceFluidState.pressure(gasPhaseIdx) - refPressure; - - // calculate the derivative of the deviation from the standard - // pressure - Scalar eps = alpha*1e-10; - tmpMolarities = molarities; - tmpMolarities /= alpha + eps; - Flash::template solve(surfaceFluidState, matParams, paramCache, tmpMolarities); - Scalar fStar = surfaceFluidState.pressure(gasPhaseIdx) - refPressure; - Scalar fPrime = (fStar - f)/eps; - - // newton update - Scalar delta = f/fPrime; - alpha -= delta; - if (std::abs(delta) < std::abs(alpha)*1e-9) { - break; - } - } - - // calculate the final result - tmpMolarities = molarities; - tmpMolarities /= alpha; - Flash::template solve(surfaceFluidState, matParams, paramCache, tmpMolarities); - return alpha; -} - -template -void printResult(const RawTable& rawTable, - const std::string& fieldName, - size_t firstIdx, - size_t secondIdx, - double hiresThres) -{ - std::cout << "std::vector > "< hiresThres; ++numRawHires) - {} - - for (; sampleIdx < numSamples; ++sampleIdx) { - size_t rawIdx = sampleIdx*numRawHires/numSamples; - std::cout << "{ " << rawTable[rawIdx][firstIdx] << ", " - << rawTable[rawIdx][secondIdx] << " }" - << ",\n"; - } - - numSamples = 15; - for (sampleIdx = 0; sampleIdx < numSamples; ++sampleIdx) { - size_t rawIdx = sampleIdx*(rawTable.size() - numRawHires)/numSamples + numRawHires; - std::cout << "{ " << rawTable[rawIdx][firstIdx] << ", " - << rawTable[rawIdx][secondIdx] << " }"; - if (sampleIdx < numSamples - 1) - std::cout << ",\n"; - else - std::cout << "\n"; - } - - std::cout << "};\n"; -} - -template -inline void testAll() -{ - typedef Opm::Spe5FluidSystem FluidSystem; - - enum { - numPhases = FluidSystem::numPhases, - waterPhaseIdx = FluidSystem::waterPhaseIdx, - gasPhaseIdx = FluidSystem::gasPhaseIdx, - oilPhaseIdx = FluidSystem::oilPhaseIdx, - - numComponents = FluidSystem::numComponents, - H2OIdx = FluidSystem::H2OIdx, - C1Idx = FluidSystem::C1Idx, - C3Idx = FluidSystem::C3Idx, - C6Idx = FluidSystem::C6Idx, - C10Idx = FluidSystem::C10Idx, - C15Idx = FluidSystem::C15Idx, - C20Idx = FluidSystem::C20Idx - }; - - typedef Opm::NcpFlash Flash; - typedef Dune::FieldVector ComponentVector; - typedef Opm::CompositionalFluidState FluidState; - - typedef Opm::ThreePhaseMaterialTraits MaterialTraits; - typedef Opm::LinearMaterial MaterialLaw; - typedef typename MaterialLaw::Params MaterialLawParams; - - typedef typename FluidSystem::template ParameterCache ParameterCache; - - //////////// - // Initialize the fluid system and create the capillary pressure - // parameters - //////////// - Scalar T = 273.15 + 20; // 20 deg Celsius - FluidSystem::init(/*minTemperature=*/T - 1, - /*maxTemperature=*/T + 1, - /*minPressure=*/1.0e4, - /*maxTemperature=*/40.0e6); - - // set the parameters for the capillary pressure law - MaterialLawParams matParams; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - matParams.setPcMinSat(phaseIdx, 0.0); - matParams.setPcMaxSat(phaseIdx, 0.0); - } - matParams.finalize(); - - //////////// - // Create a fluid state - //////////// - FluidState gasFluidState; - createSurfaceGasFluidSystem(gasFluidState); - - FluidState fluidState; - ParameterCache paramCache; - - // temperature - fluidState.setTemperature(T); - - // oil pressure - fluidState.setPressure(oilPhaseIdx, 4000 * 6894.7573); // 4000 PSI - - // oil saturation - fluidState.setSaturation(oilPhaseIdx, 1.0); - fluidState.setSaturation(gasPhaseIdx, 1.0 - fluidState.saturation(oilPhaseIdx)); - - // oil composition: SPE-5 reservoir oil - fluidState.setMoleFraction(oilPhaseIdx, H2OIdx, 0.0); - fluidState.setMoleFraction(oilPhaseIdx, C1Idx, 0.50); - fluidState.setMoleFraction(oilPhaseIdx, C3Idx, 0.03); - fluidState.setMoleFraction(oilPhaseIdx, C6Idx, 0.07); - fluidState.setMoleFraction(oilPhaseIdx, C10Idx, 0.20); - fluidState.setMoleFraction(oilPhaseIdx, C15Idx, 0.15); - fluidState.setMoleFraction(oilPhaseIdx, C20Idx, 0.05); - - //makeOilSaturated(fluidState, gasFluidState); - - // set the saturations and pressures of the other phases - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (phaseIdx != oilPhaseIdx) { - fluidState.setSaturation(phaseIdx, 0.0); - fluidState.setPressure(phaseIdx, fluidState.pressure(oilPhaseIdx)); - } - - // initial guess for the composition (needed by the ComputeFromReferencePhase - // constraint solver. TODO: bug in ComputeFromReferencePhase?) - guessInitial(fluidState, phaseIdx); - } - - typedef Opm::ComputeFromReferencePhase CFRP; - CFRP::solve(fluidState, - paramCache, - /*refPhaseIdx=*/oilPhaseIdx, - /*setViscosity=*/false, - /*setEnthalpy=*/false); - - //////////// - // Calculate the total molarities of the components - //////////// - ComponentVector totalMolarities; - for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) - totalMolarities[compIdx] = fluidState.saturation(oilPhaseIdx)*fluidState.molarity(oilPhaseIdx, compIdx); - - //////////// - // Gradually increase the volume for and calculate the gas - // formation factor, oil formation volume factor and gas formation - // volume factor. - //////////// - - FluidState flashFluidState, surfaceFluidState; - flashFluidState.assign(fluidState); - //Flash::guessInitial(flashFluidState, totalMolarities); - Flash::template solve(flashFluidState, matParams, paramCache, totalMolarities); - - Scalar surfaceAlpha = 1; - surfaceAlpha = bringOilToSurface(surfaceFluidState, surfaceAlpha, flashFluidState, /*guessInitial=*/true); - Scalar rho_gRef = surfaceFluidState.density(gasPhaseIdx); - Scalar rho_oRef = surfaceFluidState.density(oilPhaseIdx); - - std::vector > resultTable; - - Scalar minAlpha = 0.98; - Scalar maxAlpha = surfaceAlpha; - - std::cout << "alpha[-] p[Pa] S_g[-] rho_o[kg/m^3] rho_g[kg/m^3] [kg/mol] [kg/mol] R_s[m^3/m^3] B_g[-] B_o[-]\n"; - int n = 300; - for (int i = 0; i < n; ++i) { - // ratio between the original and the current volume - Scalar alpha = minAlpha + (maxAlpha - minAlpha)*i/(n - 1); - - // increasing the volume means decreasing the molartity - ComponentVector curTotalMolarities = totalMolarities; - curTotalMolarities /= alpha; - - // "flash" the modified reservoir oil - Flash::template solve(flashFluidState, matParams, paramCache, curTotalMolarities); - - surfaceAlpha = bringOilToSurface(surfaceFluidState, - surfaceAlpha, - flashFluidState, - /*guessInitial=*/false); - Scalar Rs = - surfaceFluidState.saturation(gasPhaseIdx) - / surfaceFluidState.saturation(oilPhaseIdx); - std::cout << alpha << " " - << flashFluidState.pressure(oilPhaseIdx) << " " - << flashFluidState.saturation(gasPhaseIdx) << " " - << flashFluidState.density(oilPhaseIdx) << " " - << flashFluidState.density(gasPhaseIdx) << " " - << flashFluidState.averageMolarMass(oilPhaseIdx) << " " - << flashFluidState.averageMolarMass(gasPhaseIdx) << " " - << Rs << " " - << rho_gRef/flashFluidState.density(gasPhaseIdx) << " " - << rho_oRef/flashFluidState.density(oilPhaseIdx) << " " - << "\n"; - - std::array tmp; - tmp[0] = alpha; - tmp[1] = flashFluidState.pressure(oilPhaseIdx); - tmp[2] = flashFluidState.saturation(gasPhaseIdx); - tmp[3] = flashFluidState.density(oilPhaseIdx); - tmp[4] = flashFluidState.density(gasPhaseIdx); - tmp[5] = flashFluidState.averageMolarMass(oilPhaseIdx); - tmp[6] = flashFluidState.averageMolarMass(gasPhaseIdx); - tmp[7] = Rs; - tmp[8] = rho_gRef/flashFluidState.density(gasPhaseIdx); - tmp[9] = rho_oRef/flashFluidState.density(oilPhaseIdx); - - resultTable.push_back(tmp); - } - - std::cout << "reference density oil [kg/m^3]: " << rho_oRef << "\n"; - std::cout << "reference density gas [kg/m^3]: " << rho_gRef << "\n"; - - Scalar hiresThresholdPressure = resultTable[20][1]; - printResult(resultTable, - "Bg", /*firstIdx=*/1, /*secondIdx=*/8, - /*hiresThreshold=*/hiresThresholdPressure); - printResult(resultTable, - "Bo", /*firstIdx=*/1, /*secondIdx=*/9, - /*hiresThreshold=*/hiresThresholdPressure); - printResult(resultTable, - "Rs", /*firstIdx=*/1, /*secondIdx=*/7, - /*hiresThreshold=*/hiresThresholdPressure); -} - void testChiFlash() { using Scalar = double; - typedef Opm::Spe5FluidSystem FluidSystem; + using FluidSystem = Opm::TwoPhaseThreeComponentFluidSystem; - enum { - numPhases = FluidSystem::numPhases, - waterPhaseIdx = FluidSystem::waterPhaseIdx, - gasPhaseIdx = FluidSystem::gasPhaseIdx, - oilPhaseIdx = FluidSystem::oilPhaseIdx, - - numComponents = FluidSystem::numComponents, - H2OIdx = FluidSystem::H2OIdx, - C1Idx = FluidSystem::C1Idx, - C3Idx = FluidSystem::C3Idx, - C6Idx = FluidSystem::C6Idx, - C10Idx = FluidSystem::C10Idx, - C15Idx = FluidSystem::C15Idx, - C20Idx = FluidSystem::C20Idx - }; - - //typedef Opm::NcpFlash Flash; + constexpr auto numComponents = FluidSystem::numComponents; typedef Dune::FieldVector ComponentVector; typedef Opm::CompositionalFluidState FluidState; - typedef Opm::ThreePhaseMaterialTraits MaterialTraits; typedef Opm::LinearMaterial MaterialLaw; typedef typename MaterialLaw::Params MaterialLawParams;