My Project
LBCco2rich.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  Copyright 2022 NORCE.
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 2 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 
21  Consult the COPYING file in the top-level source directory of this
22  module for the precise wording of the license and the list of
23  copyright holders.
24 */
30 #ifndef OPM_LBC_CO2RICH_HPP
31 #define OPM_LBC_CO2RICH_HPP
32 
33 #include <cmath>
34 #include <vector>
35 
36 namespace Opm
37 {
38 template <class Scalar, class FluidSystem>
39 class ViscosityModels
40 {
41 
42 public:
43 
44  // Improved LBC model for CO2 rich mixtures. (Lansangan, Taylor, Smith & Kovarik - 1993)
45  template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
46  static LhsEval LBCco2rich(const FluidState& fluidState,
47  const Params& /*paramCache*/,
48  unsigned phaseIdx)
49  {
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));
53 
54  LhsEval sumMm = 0.0;
55  LhsEval sumVolume = 0.0;
56  for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
57  const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
58  const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
59  const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
60  const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
61  const Scalar v_c = FluidSystem::criticalVolume(compIdx); // in m3/kmol
62  sumMm += x*Mm;
63  sumVolume += x*v_c;
64  }
65 
66  LhsEval rho_pc = sumMm/sumVolume; //mixture pseudocritical density
67  LhsEval rho_r = rho/rho_pc;
68 
69  LhsEval xxT_p = 0.0; // x*x*T_c/p_c
70  LhsEval xxT2_p = 0.0; // x*x*T^2_c/p_c
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; // in Mpa;
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; // in Mpa;
78  const auto& x_j = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, j_compIdx));
79 
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);
82 
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;
85  }
86  }
87 
88  const LhsEval T_pc = xxT2_p/xxT_p; //mixture pseudocritical temperature
89  const LhsEval p_pc = T_pc/xxT_p; //mixture pseudocritical pressure
90 
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);
93 
94  LhsEval my0 = 0.0;
95  LhsEval sumxrM = 0.0;
96  for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
97  const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
98  const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
99  const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
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);
103  LhsEval T_r = T/T_c;
104  LhsEval xrM = x * std::pow(Mm,0.5);
105  LhsEval mys = 0.0;
106  if (T_r <=1.5) {
107  mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta;
108  } else {
109  mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta;
110  }
111  my0 += xrM*mys;
112  sumxrM += xrM;
113  }
114  my0 /= sumxrM;
115 
116  std::vector<Scalar> LBC = {0.10230,
117  0.023364,
118  0.058533,
119  -0.040758, // typo in 1964-paper: -0.40758
120  0.0093324};
121 
122  LhsEval sumLBC = 0.0;
123  for (int i = 0; i < 5; ++i) {
124  sumLBC += Opm::pow(rho_r,i)*LBC[i];
125  }
126 
127  return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot -1.8366e-8*Opm::pow(rho_r,13.992))/1e3; // mPas-> Pas
128  }
129 };
130 
131 }; // namespace Opm
132 
133 #endif // OPM_LBC_CO2RICH_HPP