My Project
Loading...
Searching...
No Matches
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
36namespace Opm
37{
38template <class Scalar, class FluidSystem>
40{
41
42public:
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 = Opm::Constants<Scalar>::R;
53 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
54 const auto& P = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
55 const auto& Z = Opm::decay<LhsEval>(fluidState.compressFactor(phaseIdx));
56
57 LhsEval sumVolume = 0.0;
58 for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
59 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
60 const Scalar v_c = FluidSystem::criticalVolume(compIdx) / 1000; // converting to m3/mol from m3/kmol
61 sumVolume += x*v_c;
62 }
63
64 LhsEval rho_pc = 1.0 / sumVolume;
65 LhsEval V = (R * T * Z)/P;
66 LhsEval rho = 1.0 / V;
67 LhsEval rho_r = rho / rho_pc;
68
69 LhsEval xsum_T_c = 0.0; // mixture pseudocritical temperature
70 LhsEval xsum_Mm = 0.0; // mixture molar mass
71 LhsEval xsum_p_ca = 0.0; // mixture pseudocritical pressure
72 for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
73 const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6; // converting to Mpa from pascal
74 const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
75 const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; // converting to kg/kmol from kg/mol;
76 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
77 Scalar p_ca = p_c / MPa_atm;
78 xsum_T_c += x * T_c;
79 xsum_Mm += x * Mm;
80 xsum_p_ca += x * p_ca;
81 }
82 LhsEval zeta_tot = Opm::pow(xsum_T_c / (Opm::pow(xsum_Mm,3.0) * Opm::pow(xsum_p_ca,4.0)),1./6);
83
84 LhsEval my0 = 0.0;
85 LhsEval sumxrM = 0.0;
86 for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
87 const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6; // converting to Mpa from pa;
88 const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
89 const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; // converting to kg/kmol from kg/mol;
90 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
91 Scalar p_ca = p_c / MPa_atm;
92 Scalar zeta = std::pow(T_c / (std::pow(Mm,3.0) * std::pow(p_ca,4.0)),1./6);
93 LhsEval T_r = T/T_c;
94 LhsEval xrM = x * std::pow(Mm,0.5);
95 LhsEval mys = 0.0;
96 if (T_r <= 1.5) {
97 mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta;
98 } else {
99 mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta;
100 }
101 my0 += xrM*mys;
102 sumxrM += xrM;
103 }
104 my0 /= sumxrM;
105
106 std::vector<Scalar> LBC = {0.10230,
107 0.023364,
108 0.058533,
109 -0.040758, // typo in 1964-paper: -0.40758
110 0.0093324};
111
112 LhsEval sumLBC = 0.0;
113 for (int i = 0; i < 5; ++i) {
114 sumLBC += Opm::pow(rho_r,i)*LBC[i];
115 }
116
117 return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot)/1e3; // mPas-> Pas
118 }
119
120};
121
122} // namespace Opm
123
124#endif // LBC_HPP
A central place for various physical constants occuring in some equations.
Definition Constants.hpp:41
Definition LBC.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30