30 #ifndef OPM_LBC_CO2RICH_HPP
31 #define OPM_LBC_CO2RICH_HPP
38 template <
class Scalar,
class Flu
idSystem>
45 template <
class Flu
idState,
class Params,
class LhsEval =
typename Flu
idState::Scalar>
46 static LhsEval LBCco2rich(
const FluidState& fluidState,
50 const Scalar MPa_atm = 0.101325;
51 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
52 const auto& rho = Opm::decay<LhsEval>(fluidState.density(phaseIdx));
55 LhsEval sumVolume = 0.0;
56 for (
unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
57 const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6;
58 const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
59 const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000;
60 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
61 const Scalar v_c = FluidSystem::criticalVolume(compIdx);
66 LhsEval rho_pc = sumMm/sumVolume;
67 LhsEval rho_r = rho/rho_pc;
71 for (
unsigned i_compIdx = 0; i_compIdx < FluidSystem::numComponents; ++i_compIdx) {
72 const Scalar& T_c_i = FluidSystem::criticalTemperature(i_compIdx);
73 const Scalar& p_c_i = FluidSystem::criticalPressure(i_compIdx)/1e6;
74 const auto& x_i = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, i_compIdx));
75 for (
unsigned j_compIdx = 0; j_compIdx < FluidSystem::numComponents; ++j_compIdx) {
76 const Scalar& T_c_j = FluidSystem::criticalTemperature(j_compIdx);
77 const Scalar& p_c_j = FluidSystem::criticalPressure(j_compIdx)/1e6;
78 const auto& x_j = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, j_compIdx));
80 const Scalar T_c_ij = std::sqrt(T_c_i*T_c_j);
81 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);
83 xxT_p += x_i*x_j*T_c_ij/p_c_ij;
84 xxT2_p += x_i*x_j*T_c_ij*T_c_ij/p_c_ij;
88 const LhsEval T_pc = xxT2_p/xxT_p;
89 const LhsEval p_pc = T_pc/xxT_p;
91 LhsEval p_pca = p_pc / MPa_atm;
92 LhsEval zeta_tot = Opm::pow(T_pc / (Opm::pow(sumMm,3.0) * Opm::pow(p_pca,4.0)),1./6);
96 for (
unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
97 const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6;
98 const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
99 const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000;
100 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
101 Scalar p_ca = p_c / MPa_atm;
102 Scalar zeta = std::pow(T_c / (std::pow(Mm,3.0) * std::pow(p_ca,4.0)),1./6);
104 LhsEval xrM = x * std::pow(Mm,0.5);
107 mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta;
109 mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta;
116 std::vector<Scalar> LBC = {0.10230,
122 LhsEval sumLBC = 0.0;
123 for (
int i = 0; i < 5; ++i) {
124 sumLBC += Opm::pow(rho_r,i)*LBC[i];
127 return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot -1.8366e-8*Opm::pow(rho_r,13.992))/1e3;