My Project
Loading...
Searching...
No Matches
CubicEOSParams.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_PARAMS_HPP
24#define CUBIC_EOS_PARAMS_HPP
25
27
28#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
29
30#include <opm/material/eos/PRParams.hpp>
31#include <opm/material/eos/RKParams.hpp>
32#include <opm/material/eos/SRKParams.hpp>
33
34namespace Opm
35{
36
37template <class Scalar, class FluidSystem, unsigned phaseIdx>
39{
40 enum { numComponents = FluidSystem::numComponents };
41
42 static constexpr Scalar R = Constants<Scalar>::R;
43
44 using EOSType = CompositionalConfig::EOSType;
48
49public:
50
51 void setEOSType(const EOSType eos_type)
52 {
53 EosType_= eos_type;
54 }
55
56 void updatePure(Scalar temperature, Scalar pressure)
57 {
58 Valgrind::CheckDefined(temperature);
59 Valgrind::CheckDefined(pressure);
60
61 // calculate Ai and Bi
62 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
63 Scalar pr = pressure / FluidSystem::criticalPressure(compIdx);
64 Scalar Tr = temperature / FluidSystem::criticalTemperature(compIdx);
65 Scalar OmegaA = OmegaA_(temperature, compIdx);
66 Scalar OmegaB = OmegaB_();
67
68 Scalar newA = OmegaA * pr / (Tr * Tr);
69 Scalar newB = OmegaB * pr / Tr;
70 assert(std::isfinite(scalarValue(newA)));
71 assert(std::isfinite(scalarValue(newB)));
72
73 setAi(newA, compIdx);
74 setBi(newB, compIdx);
75 Valgrind::CheckDefined(Ai(compIdx));
76 Valgrind::CheckDefined(Bi(compIdx));
77 }
78
79 // Update Aij
80 updateACache_();
81
82 }
83
84 template <class FluidState>
85 void updateMix(const FluidState& fs)
86 {
87 using FlashEval = typename FluidState::Scalar;
88
89 FlashEval newA = 0;
90 FlashEval newB = 0;
91 for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
92 const FlashEval moleFracI = fs.moleFraction(phaseIdx, compIIdx);
93 FlashEval xi = max(0.0, min(1.0, moleFracI));
94 Valgrind::CheckDefined(xi);
95
96 for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
97 const FlashEval moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
98 FlashEval xj = max(0.0, min(1.0, moleFracJ));
99 Valgrind::CheckDefined(xj);
100
101 // Calculate A
102 newA += xi * xj * aCache_[compIIdx][compJIdx];
103 assert(std::isfinite(scalarValue(newA)));
104 }
105
106 // Calculate B
107 newB += max(0.0, xi) * Bi(compIIdx);
108 assert(std::isfinite(scalarValue(newB)));
109 }
110
111 // assign A and B
112 setA(decay<Scalar>(newA));
113 setB(decay<Scalar>(newB));
114 Valgrind::CheckDefined(A());
115 Valgrind::CheckDefined(B());
116 }
117
118 template <class FluidState>
119 void updateSingleMoleFraction(const FluidState& fs,
120 unsigned /*compIdx*/)
121 {
122 updateMix(fs);
123 }
124
125 Scalar aCache(unsigned compIIdx, unsigned compJIdx ) const
126 {
127 return aCache_[compIIdx][compJIdx];
128 }
129
130 void setAi(Scalar value, unsigned compIdx)
131 {
132 Ai_[compIdx] = value;
133 }
134
135 void setBi(Scalar value, unsigned compIdx)
136 {
137 Bi_[compIdx] = value;
138 }
139
140 Scalar Ai(unsigned compIdx) const
141 {
142 return Ai_[compIdx];
143 }
144
145 Scalar Bi(unsigned compIdx) const
146 {
147 return Bi_[compIdx];
148 }
149
150 void setA(Scalar value)
151 {
152 A_ = value;
153 }
154
155 void setB(Scalar value)
156 {
157 B_ = value;
158 }
159
160 Scalar A() const
161 {
162 return A_;
163 }
164
165 Scalar B() const
166 {
167 return B_;
168 }
169
170 Scalar m1() const
171 {
172 switch (EosType_) {
173 case EOSType::PRCORR:
174 case EOSType::PR:
175 return PR::calcm1();
176 case EOSType::RK:
177 return RK::calcm1();
178 case EOSType::SRK:
179 return SRK::calcm1();
180 default:
181 throw std::runtime_error("EOS type not implemented!");
182 }
183 }
184
185 Scalar m2() const
186 {
187 switch (EosType_) {
188 case EOSType::PRCORR:
189 case EOSType::PR:
190 return PR::calcm2();
191 case EOSType::RK:
192 return RK::calcm2();
193 case EOSType::SRK:
194 return SRK::calcm2();
195 default:
196 throw std::runtime_error("EOS type not implemented!");
197 }
198 }
199
200protected:
201 std::array<Scalar, numComponents> Ai_;
202 std::array<Scalar, numComponents> Bi_;
203 Scalar A_;
204 Scalar B_;
205 std::array<std::array<Scalar, numComponents>, numComponents> aCache_;
206
207 EOSType EosType_;
208
209private:
210 void updateACache_()
211 {
212 for (unsigned compIIdx = 0; compIIdx < numComponents; ++ compIIdx) {
213 for (unsigned compJIdx = 0; compJIdx < numComponents; ++ compJIdx) {
214 // interaction coefficient as given in SPE5
215 Scalar Psi = FluidSystem::interactionCoefficient(compIIdx, compJIdx);
216
217 aCache_[compIIdx][compJIdx] = sqrt(Ai(compIIdx) * Ai(compJIdx)) * (1 - Psi);
218 }
219 }
220 }
221
222 Scalar OmegaA_(Scalar temperature, unsigned compIdx)
223 {
224 switch (EosType_) {
225 case EOSType::PRCORR:
226 return PR::calcOmegaA(temperature, compIdx, /*modified=*/true);
227 case EOSType::PR:
228 return PR::calcOmegaA(temperature, compIdx, /*modified=*/false);
229 case EOSType::RK:
230 return RK::calcOmegaA(temperature, compIdx);
231 case EOSType::SRK:
232 return SRK::calcOmegaA(temperature, compIdx);
233 default:
234 throw std::runtime_error("EOS type not implemented!");
235 }
236 }
237
238 Scalar OmegaB_()
239 {
240 switch (EosType_) {
241 case EOSType::PRCORR:
242 case EOSType::PR:
243 return PR::calcOmegaB();
244 case EOSType::RK:
245 return RK::calcOmegaB();
246 case EOSType::SRK:
247 return SRK::calcOmegaB();
248 default:
249 throw std::runtime_error("EOS type not implemented!");
250 }
251 }
252
253}; // class CubicEOSParams
254
255} // namespace Opm
256
257#endif
Definition Constants.hpp:40
Definition CubicEOSParams.hpp:39
Definition PRParams.hpp:35
Definition RKParams.hpp:33
Definition SRKParams.hpp:33
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30