My Project
LBC.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 LBC_HPP
31 #define LBC_HPP
32 
33 #include <cmath>
34 #include <vector>
35 
36 namespace Opm
37 {
38 template <class Scalar, class FluidSystem>
40 {
41 
42 public:
43 
44  // Standard LBC model. (Lohrenz, Bray & Clark: "Calculating Viscosities of Reservoir
45  // fluids from Their Compositions", JPT 16.10 (1964).
46  template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
47  static LhsEval LBC(const FluidState& fluidState,
48  const Params& /*paramCache*/,
49  unsigned phaseIdx)
50  {
51  const Scalar MPa_atm = 0.101325;
52  const Scalar R = 8.3144598e-3;//Mj/kmol*K
53  const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
54  const auto& rho = Opm::decay<LhsEval>(fluidState.density(phaseIdx));
55 
56  LhsEval sumMm = 0.0;
57  LhsEval sumVolume = 0.0;
58  for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
59  const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
60  const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
61  const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
62  const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
63  const Scalar v_c = FluidSystem::criticalVolume(compIdx); // in m3/kmol
64  sumMm += x*Mm;
65  sumVolume += x*v_c;
66  }
67 
68  LhsEval rho_pc = sumMm/sumVolume; //mixture pseudocritical density
69  LhsEval rho_r = rho/rho_pc;
70 
71 
72  LhsEval xsum_T_c = 0.0; //mixture pseudocritical temperature
73  LhsEval xsum_Mm = 0.0; //mixture molar mass
74  LhsEval xsum_p_ca = 0.0; //mixture pseudocritical pressure
75  for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
76  const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
77  const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
78  const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
79  const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
80  Scalar p_ca = p_c / MPa_atm;
81  xsum_T_c += x*T_c;
82  xsum_Mm += x*Mm;
83  xsum_p_ca += x*p_ca;
84  }
85  LhsEval zeta_tot = Opm::pow(xsum_T_c / (Opm::pow(xsum_Mm,3.0) * Opm::pow(xsum_p_ca,4.0)),1./6);
86 
87  LhsEval my0 = 0.0;
88  LhsEval sumxrM = 0.0;
89  for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
90  const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa;
91  const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
92  const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol;
93  const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
94  Scalar p_ca = p_c / MPa_atm;
95  Scalar zeta = std::pow(T_c / (std::pow(Mm,3.0) * std::pow(p_ca,4.0)),1./6);
96  LhsEval T_r = T/T_c;
97  LhsEval xrM = x * std::pow(Mm,0.5);
98  LhsEval mys = 0.0;
99  if (T_r <=1.5) {
100  mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta;
101  } else {
102  mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta;
103  }
104  my0 += xrM*mys;
105  sumxrM += xrM;
106  }
107  my0 /= sumxrM;
108 
109  std::vector<Scalar> LBC = {0.10230,
110  0.023364,
111  0.058533,
112  -0.040758, // typo in 1964-paper: -0.40758
113  0.0093324};
114 
115  LhsEval sumLBC = 0.0;
116  for (int i = 0; i < 5; ++i) {
117  sumLBC += Opm::pow(rho_r,i)*LBC[i];
118  }
119 
120  return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot)/1e3; // mPas-> Pas
121  }
122 
123 };
124 
125 }; // namespace Opm
126 
127 #endif // LBC_HPP
Definition: LBC.hpp:40