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