diff --git a/opm/material/fluidsystems/Co2BrineFluidSystem.hh b/opm/material/fluidsystems/Co2BrineFluidSystem.hh index dceb0818c..f17152c53 100644 --- a/opm/material/fluidsystems/Co2BrineFluidSystem.hh +++ b/opm/material/fluidsystems/Co2BrineFluidSystem.hh @@ -7,7 +7,8 @@ #include #include -#include +#include +//#include namespace Opm { /*! @@ -36,9 +37,12 @@ namespace Opm { template using ParameterCache = Opm::PTFlashParameterCache>; - using LBCviscosity = typename Opm::LBCviscosity>; + using ViscosityModel = typename Opm::ViscosityModels>; + //using ViscosityModel = typename Opm::ViscosityModels>; + using PengRobinsonMixture = typename Opm::PengRobinsonMixture>; + /*! * \brief The acentric factor of a component []. * @@ -159,9 +163,7 @@ namespace Opm { { // Use LBC method to calculate viscosity LhsEval mu; - // mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx); - // mu = LBCviscosity::LBC(fluidState, paramCache, phaseIdx); - mu = LBCviscosity::LBCJulia(fluidState, paramCache, phaseIdx); + mu = ViscosityModel::LBC(fluidState, paramCache, phaseIdx); return mu; diff --git a/opm/material/fluidsystems/ThreeComponentFluidSystem.hh b/opm/material/fluidsystems/ThreeComponentFluidSystem.hh index df6143c23..952573974 100644 --- a/opm/material/fluidsystems/ThreeComponentFluidSystem.hh +++ b/opm/material/fluidsystems/ThreeComponentFluidSystem.hh @@ -9,7 +9,7 @@ // TODO: this is something else need to check #include -#include +#include namespace Opm { /*! @@ -43,7 +43,7 @@ namespace Opm { template using ParameterCache = Opm::PTFlashParameterCache>; - using LBCviscosity = typename Opm::LBCviscosity>; + using ViscosityModel = typename Opm::ViscosityModels>; using PengRobinsonMixture = typename Opm::PengRobinsonMixture>; /*! @@ -171,9 +171,7 @@ namespace Opm { { // Use LBC method to calculate viscosity LhsEval mu; - // mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx); - //mu = LBCviscosity::LBC(fluidState, paramCache, phaseIdx); - mu = LBCviscosity::LBCJulia(fluidState, paramCache, phaseIdx); + mu = ViscosityModel::LBC(fluidState, paramCache, phaseIdx); } diff --git a/opm/material/viscositymodels/LBCviscosity.hpp b/opm/material/viscositymodels/LBC.hpp similarity index 71% rename from opm/material/viscositymodels/LBCviscosity.hpp rename to opm/material/viscositymodels/LBC.hpp index 23f15eb58..7fe28ef8a 100644 --- a/opm/material/viscositymodels/LBCviscosity.hpp +++ b/opm/material/viscositymodels/LBC.hpp @@ -22,11 +22,11 @@ */ /*! * \file - * \copydoc Opm::LBCviscosity + * \copydoc Opm::LBC */ -#ifndef LBC_VISCOSITY_HPP -#define LBC_VISCOSITY_HPP +#ifndef LBC_HPP +#define LBC_HPP #include #include @@ -34,7 +34,7 @@ namespace Opm { template -class LBCviscosity +class ViscosityModels { public: @@ -119,12 +119,12 @@ public: } - // Improved LBC model for CO2 rich mixtures. (Lansangan, Taylor, Smith & Kovarik - 1993) - template - static LhsEval LBCmod(const FluidState& fluidState, + // Improved LBC model for CO2 rich mixtures. (Lansangan, Taylor, Smith & Kovarik - 1993) + template + static LhsEval LBCmodified(const FluidState& fluidState, const Params& /*paramCache*/, unsigned phaseIdx) - { + { const Scalar MPa_atm = 0.101325; const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); const auto& rho = Opm::decay(fluidState.density(phaseIdx)); @@ -188,87 +188,24 @@ public: } my0 += xrM*mys; sumxrM += xrM; - } - my0 /= sumxrM; + } + my0 /= sumxrM; - std::vector LBC = {0.10230, + 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 - } - - - // translation of the viscosity code from the Julia code - template - static LhsEval LBCJulia(const FluidState& fluidState, - const Params& /*paramCache*/, - unsigned phaseIdx) { - constexpr Scalar mol_factor = 1000.; - constexpr Scalar rankine = 5. / 9.; - constexpr Scalar psia = 6.894757293168360e+03; - constexpr Scalar R = 8.3144598; - const LhsEval T = Opm::decay(fluidState.temperature(phaseIdx)); - const LhsEval P = Opm::decay(fluidState.pressure(phaseIdx)); - const LhsEval Z = Opm::decay(fluidState.compressFactor(phaseIdx)); - const LhsEval rho = P / (R * T * Z); - - LhsEval P_pc = 0.; - LhsEval T_pc = 0.; - LhsEval Vc = 0.; - LhsEval mwc = 0.; - LhsEval a = 0.; - LhsEval b = 0.; - for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { - const LhsEval mol_frac = Opm::decay(fluidState.moleFraction(phaseIdx, compIdx)); - const LhsEval mol_weight = FluidSystem::molarMass(compIdx); // TODO: check values - const Scalar p_c = FluidSystem::criticalPressure(compIdx); - const Scalar T_c = FluidSystem::criticalTemperature(compIdx); - const Scalar v_c = FluidSystem::criticalVolume(compIdx); - mwc += mol_frac * mol_weight; - P_pc += mol_frac * p_c; - T_pc += mol_frac * T_c; - Vc += mol_frac * v_c; - - const LhsEval tr = T / T_c; - const Scalar Tc = T_c / rankine; - const Scalar Pc = p_c / psia; - - const LhsEval mwi = Opm::sqrt(mol_factor * mol_weight); - const LhsEval e_i = 5.4402 * Opm::pow(Tc, 1./6.) / (mwi * Opm::pow(Pc, 2./3.) * 1.e-3); - - LhsEval mu_i; - if (tr > 1.5) { - mu_i = 17.78e-5 * Opm::pow(4.58*tr - 1.67, 0.625) / e_i; - } else { - mu_i = 34.e-5 * Opm::pow(tr, 0.94) / e_i; + LhsEval sumLBC = 0.0; + for (int i = 0; i < 5; ++i) { + sumLBC += Opm::pow(rho_r,i)*LBC[i]; } - a += mol_frac * mu_i * mwi; - b += mol_frac * mwi; - } - const LhsEval mu_atm = a / b; - const LhsEval e_mix = 5.4402 * Opm::pow(T_pc/rankine, 1./6.) / - (Opm::sqrt(mol_factor * mwc) * Opm::pow(P_pc/psia, 2./3.) * (1e-3)); - const LhsEval rhor = Vc * rho; - LhsEval corr = 0.; - const std::vector LBC{0.10230, 0.023364, 0.058533, -0.040758, 0.0093324}; - const Scalar shift = -1.e-4; - for (unsigned i = 0; i < 5; ++i) { - corr += LBC[i] * Opm::pow(rhor, i); + return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot -1.8366e-8*Opm::pow(rho_r,13.992))/1e3; // mPas-> Pas } - LhsEval mu = mu_atm + (corr * corr * corr * corr + shift)/e_mix; - return mu; - } }; -} // namespace Opm +}; // namespace Opm -#endif // LBC_VISCOSITY_HPP +#endif // LBC_HPP diff --git a/opm/material/viscositymodels/LBCmodified.hpp b/opm/material/viscositymodels/LBCmodified.hpp new file mode 100644 index 000000000..e008dae4c --- /dev/null +++ b/opm/material/viscositymodels/LBCmodified.hpp @@ -0,0 +1,131 @@ +// -*- 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::ViscosityModels + */ + +#ifndef LBC_MODIFIED_HPP +#define LBC_MODIFIED_HPP + +#include +#include + +namespace Opm +{ +template +class ViscosityModels +{ + +public: + + // Improved LBC model for CO2 rich mixtures. (Lansangan, Taylor, Smith & Kovarik - 1993) + template + static LhsEval LBCmodified(const FluidState& fluidState, + const Params& /*paramCache*/, + unsigned phaseIdx) + { + const Scalar MPa_atm = 0.101325; + 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_MODIFIED_HPP