My Project
Loading...
Searching...
No Matches
CubicEOS.hpp
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 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
23#ifndef CUBIC_EOS_HPP
24#define CUBIC_EOS_HPP
25
29
30namespace Opm
31{
32template<class Scalar, class FluidSystem>
34{
35 enum { numComponents = FluidSystem::numComponents };
36
37 static constexpr Scalar R = Constants<Scalar>::R;
38
39public:
40 template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
41 static LhsEval computeFugacityCoefficient(const FluidState& fs,
42 const Params& params,
43 unsigned phaseIdx,
44 unsigned compIdx)
45 {
46 // note that we normalize the component mole fractions, so
47 // that their sum is 100%. This increases numerical stability
48 // considerably if the fluid state is not physical.
49
50 // extract variables
51 LhsEval Vm = params.molarVolume(phaseIdx);
52 LhsEval T = fs.temperature(phaseIdx);
53 LhsEval p = fs.pressure(phaseIdx); // molar volume in [bar]
54 LhsEval A = params.A(phaseIdx);
55 LhsEval B = params.B(phaseIdx);
56 LhsEval Bi = params.Bi(phaseIdx, compIdx);
57 LhsEval m1 = params.m1(phaseIdx);
58 LhsEval m2 = params.m2(phaseIdx);
59
60 // Calculate Bi / B
61 LhsEval Bi_B = Bi / B;
62
63 // Calculate the compressibility factor
64 LhsEval RT = R * T;
65 LhsEval Z = p * Vm / RT; // compressibility factor
66
67 // calculate sum(A_ij * x_j)
68 LhsEval A_s = 0.0;
69 for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
70 A_s += params.aCache(phaseIdx, compIdx, compJIdx) * fs.moleFraction(phaseIdx, compJIdx);
71 }
72
73 // calculate fugacity coefficient
74 LhsEval alpha;
75 LhsEval beta;
76 LhsEval gamma;
77 LhsEval ln_phi;
78 LhsEval fugCoeff;
79
80 alpha = -log(Z - B) + Bi_B * (Z - 1);
81 beta = log((Z + m2 * B) / (Z + m1 * B)) * A / ((m1 - m2) * B);
82 gamma = (2 / A) * A_s - Bi_B;
83 ln_phi = alpha + (beta * gamma);
84
85 fugCoeff = exp(ln_phi);
86
88 // limit the fugacity coefficient to a reasonable range:
89 //
90 // on one side, we want the mole fraction to be at
91 // least 10^-3 if the fugacity is at the current pressure
92 //
93 fugCoeff = min(1e10, fugCoeff);
94 //
95 // on the other hand, if the mole fraction of the component is 100%, we want the
96 // fugacity to be at least 10^-3 Pa
97 //
98 fugCoeff = max(1e-10, fugCoeff);
100
101 return fugCoeff;
102 }
103
104 template <class FluidState, class Params>
105 static typename FluidState::Scalar computeMolarVolume(const FluidState& fs,
106 Params& params,
107 unsigned phaseIdx,
108 bool isGasPhase)
109 {
110 Valgrind::CheckDefined(fs.temperature(phaseIdx));
111 Valgrind::CheckDefined(fs.pressure(phaseIdx));
112
113 using Evaluation = typename FluidState::Scalar;
114
115 // extract variables
116 const Evaluation& T = fs.temperature(phaseIdx);
117 const Evaluation& p = fs.pressure(phaseIdx);
118 const Evaluation& A = params.A(phaseIdx);
119 const Evaluation& B = params.B(phaseIdx);
120 const Evaluation& m1 = params.m1(phaseIdx);
121 const Evaluation& m2 = params.m2(phaseIdx);
122
123 // coefficients in cubic solver
124 const Evaluation a1 = 1.0; // cubic term
125 const Evaluation a2 = (m1 + m2 - 1) * B - 1; // quadratic term
126 const Evaluation a3 = A + m1 * m2 * B * B - (m1 + m2) * B * (B + 1); // linear term
127 const Evaluation a4 = -A * B - m1 * m2 * B * B * (B + 1); // constant term
128 Valgrind::CheckDefined(a1);
129 Valgrind::CheckDefined(a2);
130 Valgrind::CheckDefined(a3);
131 Valgrind::CheckDefined(a4);
132
133 // root of cubic equation
134 Evaluation Vm = 0;
135 Valgrind::SetUndefined(Vm);
136 Evaluation Z[3] = {0.0, 0.0, 0.0};
137 int numSol = cubicRoots(Z, a1, a2, a3, a4);
138
139 // pick correct root
140 const Evaluation RT_p = R * T / p;
141 if (numSol == 3) {
142 // the EOS has three intersections with the pressure,
143 // i.e. the molar volume of gas is the largest one and the
144 // molar volume of liquid is the smallest one
145 if (isGasPhase)
146 Vm = max(1e-7, Z[2] * RT_p);
147 else
148 Vm = max(1e-7, Z[0] * RT_p);
149 }
150 else if (numSol == 1) {
151 // the EOS only has one intersection with the pressure,
152 // for the other phase, we take the extremum of the EOS
153 // with the largest distance from the intersection.
154 Vm = max(1e-7, Z[0] * RT_p);
155 }
156
157 Valgrind::CheckDefined(Vm);
158 assert(std::isfinite(scalarValue(Vm)));
159 assert(Vm > 0);
160 return Vm;
161
162 }
163
164}; // class CubicEOS
165
166} // namespace Opm
167
168#endif
Provides free functions to invert polynomials of degree 1, 2 and 3.
Some templates to wrap the valgrind client request macros.
Definition Constants.hpp:40
Definition CubicEOS.hpp:34
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition PolynomialUtils.hpp:331