My Project
Loading...
Searching...
No Matches
EclHysteresisTwoPhaseLaw.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 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*/
27#ifndef OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_HPP
28#define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_HPP
29#include <opm/common/TimingMacros.hpp>
31
32#include <stdexcept>
33
34namespace Opm {
35
41template <class EffectiveLawT,
42 class ParamsT = EclHysteresisTwoPhaseLawParams<EffectiveLawT> >
43class EclHysteresisTwoPhaseLaw : public EffectiveLawT::Traits
44{
45public:
46 using EffectiveLaw = EffectiveLawT;
47 using EffectiveLawParams = typename EffectiveLaw::Params;
48
49 using Traits = typename EffectiveLaw::Traits;
50 using Params = ParamsT;
51 using Scalar = typename EffectiveLaw::Scalar;
52
53 enum { wettingPhaseIdx = Traits::wettingPhaseIdx };
54 enum { nonWettingPhaseIdx = Traits::nonWettingPhaseIdx };
55
57 static constexpr int numPhases = EffectiveLaw::numPhases;
58 static_assert(numPhases == 2,
59 "The endpoint scaling applies to the nested twophase laws, not to "
60 "the threephase one!");
61
64 static constexpr bool implementsTwoPhaseApi = true;
65
66 static_assert(EffectiveLaw::implementsTwoPhaseApi,
67 "The material laws put into EclEpsTwoPhaseLaw must implement the "
68 "two-phase material law API!");
69
72 static constexpr bool implementsTwoPhaseSatApi = true;
73
74 static_assert(EffectiveLaw::implementsTwoPhaseSatApi,
75 "The material laws put into EclEpsTwoPhaseLaw must implement the "
76 "two-phase material law saturation API!");
77
80 static constexpr bool isSaturationDependent = true;
81
84 static constexpr bool isPressureDependent = false;
85
88 static constexpr bool isTemperatureDependent = false;
89
92 static constexpr bool isCompositionDependent = false;
93
104 template <class Container, class FluidState>
105 static void capillaryPressures(Container& /* values */,
106 const Params& /* params */,
107 const FluidState& /* fs */)
108 {
109 throw std::invalid_argument("The capillaryPressures(fs) method is not yet implemented");
110 }
111
122 template <class Container, class FluidState>
123 static void relativePermeabilities(Container& /* values */,
124 const Params& /* params */,
125 const FluidState& /* fs */)
126 {
127 throw std::invalid_argument("The pcnw(fs) method is not yet implemented");
128 }
129
141 template <class FluidState, class Evaluation = typename FluidState::Scalar>
142 static Evaluation pcnw(const Params& /* params */,
143 const FluidState& /* fs */)
144 {
145 throw std::invalid_argument("The pcnw(fs) method is not yet implemented");
146 }
147
148 template <class Evaluation>
149 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
150 {
151 OPM_TIMEFUNCTION_LOCAL();
152 // if no pc hysteresis is enabled, use the drainage curve
153 if (!params.config().enableHysteresis() || params.config().pcHysteresisModel() < 0)
154 return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
155
156 // Initial imbibition process
157 if (params.initialImb()) {
158 if (Sw >= params.pcSwMic()) {
159 return EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), Sw);
160 }
161 else { // Reversal
162 const Evaluation& F = (1.0/(params.pcSwMic()-Sw+params.curvatureCapPrs())-1.0/params.curvatureCapPrs())
163 / (1.0/(params.pcSwMic()-params.Swcrd()+params.curvatureCapPrs())-1.0/params.curvatureCapPrs());
164
165 const Evaluation& Pcd = EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
166 const Evaluation& Pci = EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), Sw);
167 const Evaluation& pc_Killough = Pci+F*(Pcd-Pci);
168
169 return pc_Killough;
170 }
171 }
172
173 // Initial drainage process
174 if (Sw <= params.pcSwMdc())
175 return EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
176
177 // Reversal
178 Scalar Swma = 1.0-params.Sncrt();
179 if (Sw >= Swma) {
180 const Evaluation& Pci = EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), Sw);
181 return Pci;
182 }
183 else {
184 Scalar pciwght = params.pcWght(); // Align pci and pcd at Swir
185 //const Evaluation& SwEff = params.Swcri()+(Sw-params.pcSwMdc())/(Swma-params.pcSwMdc())*(Swma-params.Swcri());
186 const Evaluation& SwEff = Sw; // This is Killough 1976, Gives significantly better fit compared to benchmark then the above "scaling"
187 const Evaluation& Pci = pciwght*EffectiveLaw::twoPhaseSatPcnw(params.imbibitionParams(), SwEff);
188
189 const Evaluation& Pcd = EffectiveLaw::twoPhaseSatPcnw(params.drainageParams(), Sw);
190
191 if (Pci == Pcd)
192 return Pcd;
193
194 const Evaluation& F = (1.0/(Sw-params.pcSwMdc()+params.curvatureCapPrs())-1.0/params.curvatureCapPrs())
195 / (1.0/(Swma-params.pcSwMdc()+params.curvatureCapPrs())-1.0/params.curvatureCapPrs());
196
197 const Evaluation& pc_Killough = Pcd+F*(Pci-Pcd);
198
199 return pc_Killough;
200 }
201
202 return 0.0;
203 }
204
208 template <class Container, class FluidState>
209 static void saturations(Container& /* values */,
210 const Params& /* params */,
211 const FluidState& /* fs */)
212 {
213 throw std::invalid_argument("The saturations(fs) method is not yet implemented");
214 }
215
220 template <class FluidState, class Evaluation = typename FluidState::Scalar>
221 static Evaluation Sw(const Params& /* params */,
222 const FluidState& /* fs */)
223 {
224 throw std::invalid_argument("The Sw(fs) method is not yet implemented");
225 }
226
227 template <class Evaluation>
228 static Evaluation twoPhaseSatSw(const Params& /* params */,
229 const Evaluation& /* pc */)
230 {
231 throw std::invalid_argument("The twoPhaseSatSw(pc) method is not yet implemented");
232 }
233
238 template <class FluidState, class Evaluation = typename FluidState::Scalar>
239 static Evaluation Sn(const Params& /* params */,
240 const FluidState& /* fs */)
241 {
242 throw std::invalid_argument("The Sn(pc) method is not yet implemented");
243 }
244
245 template <class Evaluation>
246 static Evaluation twoPhaseSatSn(const Params& /* params */,
247 const Evaluation& /* pc */)
248 {
249 throw std::invalid_argument("The twoPhaseSatSn(pc) method is not yet implemented");
250 }
251
261 template <class FluidState, class Evaluation = typename FluidState::Scalar>
262 static Evaluation krw(const Params& /* params */,
263 const FluidState& /* fs */)
264 {
265 throw std::invalid_argument("The krw(fs) method is not yet implemented");
266 }
267
268 template <class Evaluation>
269 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
270 {
271 OPM_TIMEFUNCTION_LOCAL();
272 // if no relperm hysteresis is enabled, use the drainage curve
273 if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
274 return EffectiveLaw::twoPhaseSatKrw(params.drainageParams(), Sw);
275
276 if (params.config().krHysteresisModel() == 0 || params.config().krHysteresisModel() == 2)
277 // use drainage curve for wetting phase
278 return EffectiveLaw::twoPhaseSatKrw(params.drainageParams(), Sw);
279
280 // use imbibition curve for wetting phase
281 assert(params.config().krHysteresisModel() == 1 || params.config().krHysteresisModel() == 3);
282 return EffectiveLaw::twoPhaseSatKrw(params.imbibitionParams(), Sw);
283 }
284
288 template <class FluidState, class Evaluation = typename FluidState::Scalar>
289 static Evaluation krn(const Params& /* params */,
290 const FluidState& /* fs */)
291 {
292 throw std::invalid_argument("The krn(fs) method is not yet implemented");
293 }
294
295 template <class Evaluation>
296 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
297 {
298 OPM_TIMEFUNCTION_LOCAL();
299 // If WAG hysteresis is enabled, the convential hysteresis model is ignored.
300 // (Two-phase model, non-wetting: only gas in oil.)
301 if (params.gasOilHysteresisWAG()) {
302
303 // Primary drainage
304 if (Sw <= params.krnSwMdc()+params.tolWAG() && params.nState()==1) {
305 Evaluation krn = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
306 return krn;
307 }
308
309 // Imbibition or reversion to two-phase drainage retracing imb curve
310 // (Shift along primary drainage curve.)
311 if (params.nState()==1) {
312 Evaluation Swf = params.computeSwf(Sw);
313 Evaluation krn = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Swf);
314 return krn;
315 }
316
317 // Three-phase drainage along current secondary drainage curve
318 if (Sw <= params.krnSwDrainRevert()+params.tolWAG() /*&& params.nState()>=1 */) {
319 Evaluation Krg = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
320 Evaluation KrgDrain2 = (Krg-params.krnDrainStart())*params.reductionDrain() + params.krnImbStart();
321 return KrgDrain2;
322 }
323
324 // Subsequent imbibition: Scanning curve derived from previous secondary drainage
325 if (Sw >= params.krnSwWAG()-params.tolWAG() /*&& Sw > params.krnSwDrainRevert() && params.nState()>=1 */) {
326 Evaluation KrgImb2 = params.computeKrImbWAG(Sw);
327 return KrgImb2;
328 }
329 else {/* Sw < params.krnSwWAG() */ // Reversion along "next" drainage curve
330 Evaluation Krg = EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
331 Evaluation KrgDrainNxt = (Krg-params.krnDrainStartNxt())*params.reductionDrainNxt() + params.krnImbStartNxt();
332 return KrgDrainNxt;
333 }
334 }
335
336 // if no relperm hysteresis is enabled, use the drainage curve
337 if (!params.config().enableHysteresis() || params.config().krHysteresisModel() < 0)
338 return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
339
340
341 // if it is enabled, use either the drainage or the imbibition curve. if the
342 // imbibition curve is used, the saturation must be shifted.
343 if (Sw <= params.krnSwMdc())
344 return EffectiveLaw::twoPhaseSatKrn(params.drainageParams(), Sw);
345
346 if (params.config().krHysteresisModel() <= 1) { //Carlson
347 return EffectiveLaw::twoPhaseSatKrn(params.imbibitionParams(),
348 Sw + params.deltaSwImbKrn());
349 }
350
351 // Killough
352 assert(params.config().krHysteresisModel() == 2 || params.config().krHysteresisModel() == 3);
353 Evaluation Snorm = params.Sncri()+(1.0-Sw-params.Sncrt())*(params.Snmaxd()-params.Sncri())/(params.Snhy()-params.Sncrt());
354 return params.krnWght()*EffectiveLaw::twoPhaseSatKrn(params.imbibitionParams(),1.0-Snorm);
355 }
356};
357
358} // namespace Opm
359
360#endif
A default implementation of the parameters for the material law which implements the ECL relative per...
This material law implements the hysteresis model of the ECL file format.
Definition EclHysteresisTwoPhaseLaw.hpp:44
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition EclHysteresisTwoPhaseLaw.hpp:80
static Evaluation krn(const Params &, const FluidState &)
The relative permeability of the non-wetting phase.
Definition EclHysteresisTwoPhaseLaw.hpp:289
static void capillaryPressures(Container &, const Params &, const FluidState &)
The capillary pressure-saturation curves depending on absolute saturations.
Definition EclHysteresisTwoPhaseLaw.hpp:105
static constexpr bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition EclHysteresisTwoPhaseLaw.hpp:64
static constexpr bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition EclHysteresisTwoPhaseLaw.hpp:92
static Evaluation Sw(const Params &, const FluidState &)
Calculate wetting liquid phase saturation given that the rest of the fluid state has been initialized...
Definition EclHysteresisTwoPhaseLaw.hpp:221
static Evaluation Sn(const Params &, const FluidState &)
Calculate non-wetting liquid phase saturation given that the rest of the fluid state has been initial...
Definition EclHysteresisTwoPhaseLaw.hpp:239
static void saturations(Container &, const Params &, const FluidState &)
The saturation-capillary pressure curves.
Definition EclHysteresisTwoPhaseLaw.hpp:209
static Evaluation krw(const Params &, const FluidState &)
The relative permeability for the wetting phase.
Definition EclHysteresisTwoPhaseLaw.hpp:262
static constexpr int numPhases
The number of fluid phases.
Definition EclHysteresisTwoPhaseLaw.hpp:57
static constexpr bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition EclHysteresisTwoPhaseLaw.hpp:84
static constexpr bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition EclHysteresisTwoPhaseLaw.hpp:72
static constexpr bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition EclHysteresisTwoPhaseLaw.hpp:88
static Evaluation pcnw(const Params &, const FluidState &)
The capillary pressure-saturation curve.
Definition EclHysteresisTwoPhaseLaw.hpp:142
static void relativePermeabilities(Container &, const Params &, const FluidState &)
The relative permeability-saturation curves depending on absolute saturations.
Definition EclHysteresisTwoPhaseLaw.hpp:123
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30