My Project
EclHysteresisTwoPhaseLawParams.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_PARAMS_HPP
28#define OPM_ECL_HYSTERESIS_TWO_PHASE_LAW_PARAMS_HPP
29
34
35#include <cassert>
36#include <cmath>
37#include <memory>
38
39namespace Opm {
46template <class EffLawT>
48{
49 using EffLawParams = typename EffLawT::Params;
50 using Scalar = typename EffLawParams::Traits::Scalar;
51
52public:
53 using Traits = typename EffLawParams::Traits;
54
56 {
57 // These are initialized to two (even though they represent saturations)
58 // to signify that the values are larger than physically possible, and force
59 // using the drainage curve before the first saturation update.
60 pcSwMdc_ = 2.0;
61 krnSwMdc_ = 2.0;
62 // krwSwMdc_ = 2.0;
63
64 pcSwMic_ = -1.0;
65 initialImb_ = false;
66 oilWaterSystem_ = false;
67 pcmaxd_ = 0.0;
68 pcmaxi_ = 0.0;
69
70 deltaSwImbKrn_ = 0.0;
71 // deltaSwImbKrw_ = 0.0;
72 }
73
74 static EclHysteresisTwoPhaseLawParams serializationTestObject()
75 {
77 result.deltaSwImbKrn_ = 1.0;
78 result.Sncrt_ = 2.0;
79 result.initialImb_ = true;
80 result.pcSwMic_ = 3.0;
81 result.krnSwMdc_ = 4.0;
82 result.KrndHy_ = 5.0;
83
84 return result;
85 }
86
91 void finalize()
92 {
93 if (config().enableHysteresis()) {
94 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
95 C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
96 curvatureCapPrs_ = config().curvatureCapPrs();
97 }
98
99 updateDynamicParams_();
100 }
101
103 }
104
108 void setConfig(std::shared_ptr<EclHysteresisConfig> value)
109 { config_ = *value; }
110
115 { return config_; }
116
120 void setDrainageParams(const EffLawParams& value,
122 EclTwoPhaseSystemType twoPhaseSystem)
123 {
124 drainageParams_ = value;
125
126 oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
127
128 if (!config().enableHysteresis())
129 return;
130
131 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
132 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
133 Sncrd_ = info.Sgcr+info.Swl;
134 Snmaxd_ = info.Sgu+info.Swl;
135 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
136 }
137 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
138 Sncrd_ = info.Sgcr;
139 Snmaxd_ = info.Sgu;
140 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
141 }
142 else {
143 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
144 Sncrd_ = info.Sowcr;
145 Snmaxd_ = 1.0 - info.Swl - info.Sgl;
146 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
147 }
148 }
149
150 // Additional Killough hysteresis model for pc
151 if (config().pcHysteresisModel() == 0) {
152 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
153 Swcrd_ = info.Sogcr;
154 pcmaxd_ = info.maxPcgo;
155 } else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
156 Swcrd_ = info.Swcr;
157 pcmaxd_ = info.maxPcgo + info.maxPcow;
158 }
159 else {
160 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
161 Swcrd_ = info.Swcr;
162 pcmaxd_ = -17.0; // At this point 'info.maxPcow' holds pre-swatinit value ...;
163 }
164 }
165 }
166
170 const EffLawParams& drainageParams() const
171 { return drainageParams_; }
172
173 EffLawParams& drainageParams()
174 { return drainageParams_; }
175
179 void setImbibitionParams(const EffLawParams& value,
181 EclTwoPhaseSystemType twoPhaseSystem)
182 {
183 imbibitionParams_ = value;
184
185 if (!config().enableHysteresis())
186 return;
187
188 // Killough hysteresis model for nonw kr
189 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
190 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
191 Sncri_ = info.Sgcr+info.Swl;
192 }
193 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
194 Sncri_ = info.Sgcr;
195 }
196 else {
197 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
198 Sncri_ = info.Sowcr;
199 }
200 }
201
202 // Killough hysteresis model for pc
203 if (config().pcHysteresisModel() == 0) {
204 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
205 Swcri_ = info.Sogcr;
206 Swmaxi_ = 1.0 - info.Sgl - info.Swl;
207 pcmaxi_ = info.maxPcgo;
208 } else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
209 Swcri_ = info.Swcr;
210 Swmaxi_ = 1.0 - info.Sgl;
211 pcmaxi_ = info.maxPcgo + info.maxPcow;
212 }
213 else {
214 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
215 Swcri_ = info.Swcr;
216 Swmaxi_ = info.Swu;
217 pcmaxi_ = info.maxPcow;
218 }
219 }
220 }
221
225 const EffLawParams& imbibitionParams() const
226 { return imbibitionParams_; }
227
228 EffLawParams& imbibitionParams()
229 { return imbibitionParams_; }
230
235 Scalar pcSwMdc() const
236 { return pcSwMdc_; }
237
238 Scalar pcSwMic() const
239 { return pcSwMic_; }
240
244 bool initialImb() const
245 { return initialImb_; }
246
252 void setKrwSwMdc(Scalar /* value */)
253 {}
254 // { krwSwMdc_ = value; };
255
261 Scalar krwSwMdc() const
262 { return 2.0; }
263 // { return krwSwMdc_; };
264
270 void setKrnSwMdc(Scalar value)
271 { krnSwMdc_ = value; }
272
278 Scalar krnSwMdc() const
279 { return krnSwMdc_; }
280
288 void setDeltaSwImbKrw(Scalar /* value */)
289 {}
290 // { deltaSwImbKrw_ = value; }
291
299 Scalar deltaSwImbKrw() const
300 { return 0.0; }
301// { return deltaSwImbKrw_; }
302
310 void setDeltaSwImbKrn(Scalar value)
311 { deltaSwImbKrn_ = value; }
312
320 Scalar deltaSwImbKrn() const
321 { return deltaSwImbKrn_; }
322
323
324 Scalar Swcri() const
325 { return Swcri_; }
326
327 Scalar Swcrd() const
328 { return Swcrd_; }
329
330 Scalar Swmaxi() const
331 { return Swmaxi_; }
332
333 Scalar Sncri() const
334 { return Sncri_; }
335
336 Scalar Sncrd() const
337 { return Sncrd_; }
338
339 Scalar Sncrt() const
340 { return Sncrt_; }
341
342 Scalar Snmaxd() const
343 { return Snmaxd_; }
344
345 Scalar Snhy() const
346 { return 1.0 - krnSwMdc_; }
347
348 Scalar krnWght() const
349 { return KrndHy_/KrndMax_; }
350
351 Scalar pcWght() const // Aligning pci and pcd at Swir
352 {
353 if (pcmaxd_ < 0.0)
354 return EffLawT::twoPhaseSatPcnw(drainageParams(), 0.0)/(pcmaxi_+1e-6);
355 else
356 return pcmaxd_/(pcmaxi_+1e-6);
357 }
358
359 Scalar curvatureCapPrs() const
360 { return curvatureCapPrs_;}
361
362
369 void update(Scalar pcSw, Scalar /* krwSw */, Scalar krnSw)
370 {
371 bool updateParams = false;
372
373 if (config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
374 if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
375 initialImb_ = true;
376 }
377 pcSwMdc_ = pcSw;
378 updateParams = true;
379 }
380
381 if (initialImb_ && pcSw > pcSwMic_) {
382 pcSwMic_ = pcSw;
383 updateParams = true;
384 }
385
386/*
387 // This is quite hacky: Eclipse says that it only uses relperm hysteresis for the
388 // wetting phase (indicated by '0' for the second item of the EHYSTER keyword),
389 // even though this makes about zero sense: one would expect that hysteresis can
390 // be limited to the oil phase, but the oil phase is the wetting phase of the
391 // gas-oil twophase system whilst it is non-wetting for water-oil.
392 if (krwSw < krwSwMdc_)
393 {
394 krwSwMdc_ = krwSw;
395 updateParams = true;
396 }
397*/
398
399 if (krnSw < krnSwMdc_) {
400 krnSwMdc_ = krnSw;
401 KrndHy_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
402 updateParams = true;
403 }
404
405 if (updateParams)
406 updateDynamicParams_();
407 }
408
409 template<class Serializer>
410 void serializeOp(Serializer& serializer)
411 {
412 // only serializes dynamic state - see update() and updateDynamic_()
413 serializer(deltaSwImbKrn_);
414 serializer(Sncrt_);
415 serializer(initialImb_);
416 serializer(pcSwMic_);
417 serializer(krnSwMdc_);
418 serializer(KrndHy_);
419 }
420
421 bool operator==(const EclHysteresisTwoPhaseLawParams& rhs) const
422 {
423 return this->deltaSwImbKrn_ == rhs.deltaSwImbKrn_ &&
424 this->Sncrt_ == rhs.Sncrt_ &&
425 this->initialImb_ == rhs.initialImb_ &&
426 this->pcSwMic_ == rhs.pcSwMic_ &&
427 this->krnSwMdc_ == rhs.krnSwMdc_ &&
428 this->KrndHy_ == rhs.KrndHy_;
429 }
430
431private:
432 void updateDynamicParams_()
433 {
434 // HACK: Eclipse seems to disable the wetting-phase relperm even though this is
435 // quite pointless from the physical POV. (see comment above)
436/*
437 // calculate the saturation deltas for the relative permeabilities
438 Scalar krwMdcDrainage = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
439 Scalar SwKrwMdcImbibition = EffLawT::twoPhaseSatKrwInv(imbibitionParams(), krwMdcDrainage);
440 deltaSwImbKrw_ = SwKrwMdcImbibition - krwSwMdc_;
441*/
442 if (config().krHysteresisModel() == 0 || config().krHysteresisModel() == 1) {
443 Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
444 Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage);
445 deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
446 assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
447 - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
448 }
449
450 // Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_);
451 // Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage);
452 // deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_;
453
454// assert(std::abs(EffLawT::twoPhaseSatPcnw(imbibitionParams(), pcSwMdc_ + deltaSwImbPc_)
455// - EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_)) < 1e-8);
456// assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
457// - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
458// assert(std::abs(EffLawT::twoPhaseSatKrw(imbibitionParams(), krwSwMdc_ + deltaSwImbKrw_)
459// - EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_)) < 1e-8);
460
461 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
462 Scalar Snhy = 1.0 - krnSwMdc_;
463 if (Snhy > Sncrd_)
464 Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/((1.0+config().modParamTrapped()*(Snmaxd_-Snhy)) + C_*(Snhy - Sncrd_));
465 else
466 {
467 Sncrt_ = Sncrd_;
468 }
469 }
470 }
471
472 EclHysteresisConfig config_;
473 EffLawParams imbibitionParams_;
474 EffLawParams drainageParams_;
475
476 // largest wettinging phase saturation which is on the main-drainage curve. These are
477 // three different values because the sourounding code can choose to use different
478 // definitions for the saturations for different quantities
479 //Scalar krwSwMdc_;
480 Scalar krnSwMdc_;
481 Scalar pcSwMdc_;
482
483 // largest wettinging phase saturation along main imbibition curve
484 Scalar pcSwMic_;
485 // Initial process is imbibition (for initial saturations at or below critical drainage saturation)
486 bool initialImb_;
487
488 bool oilWaterSystem_;
489
490 // offsets added to wetting phase saturation uf using the imbibition curves need to
491 // be used to calculate the wetting phase relperm, the non-wetting phase relperm and
492 // the capillary pressure
493 //Scalar deltaSwImbKrw_;
494 Scalar deltaSwImbKrn_;
495 //Scalar deltaSwImbPc_;
496
497 // trapped non-wetting phase saturation
498 Scalar Sncrt_;
499
500 // the following uses the conventions of the Eclipse technical description:
501 //
502 // Sncrd_: critical non-wetting phase saturation for the drainage curve
503 // Sncri_: critical non-wetting phase saturation for the imbibition curve
504 // Swcri_: critical wetting phase saturation for the imbibition curve
505 // Swcrd_: critical wetting phase saturation for the drainage curve
506 // Swmaxi_; maximum wetting phase saturation for the imbibition curve
507 // Snmaxd_: non-wetting phase saturation where the non-wetting relperm reaches its
508 // maximum on the drainage curve
509 // C_: factor required to calculate the trapped non-wetting phase saturation using
510 // the Killough approach
511 Scalar Sncrd_;
512 Scalar Sncri_;
513 Scalar Swcri_;
514 Scalar Swcrd_;
515 Scalar Swmaxi_;
516 Scalar Snmaxd_;
517 Scalar C_;
518
519 Scalar KrndMax_; // Krn_drain(Snmaxd_)
520 Scalar KrndHy_; // Krn_drain(1-krnSwMdc_)
521
522 Scalar pcmaxd_; // max pc for drain
523 Scalar pcmaxi_; // max pc for imb
524
525 Scalar curvatureCapPrs_; // curvature parameter used for capillary pressure hysteresis
526};
527
528} // namespace Opm
529
530#endif
Specifies the configuration used by the ECL kr/pC hysteresis code.
Definition: EclHysteresisConfig.hpp:42
int pcHysteresisModel() const
Return the type of the hysteresis model which is used for capillary pressure.
Definition: EclHysteresisConfig.hpp:71
int krHysteresisModel() const
Return the type of the hysteresis model which is used for relative permeability.
Definition: EclHysteresisConfig.hpp:97
double curvatureCapPrs() const
Curvature parameter used for capillary pressure hysteresis.
Definition: EclHysteresisConfig.hpp:113
bool enableHysteresis() const
Returns whether hysteresis is enabled.
Definition: EclHysteresisConfig.hpp:53
A default implementation of the parameters for the material law which implements the ECL relative per...
Definition: EclHysteresisTwoPhaseLawParams.hpp:48
Scalar pcSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:235
bool initialImb() const
Status of initial process.
Definition: EclHysteresisTwoPhaseLawParams.hpp:244
const EclHysteresisConfig & config() const
Returns the endpoint scaling configuration object.
Definition: EclHysteresisTwoPhaseLawParams.hpp:114
void setDeltaSwImbKrn(Scalar value)
Sets the saturation value which must be added if krn is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:310
void setConfig(std::shared_ptr< EclHysteresisConfig > value)
Set the endpoint scaling configuration object.
Definition: EclHysteresisTwoPhaseLawParams.hpp:108
Scalar krnSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:278
void setKrnSwMdc(Scalar value)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition: EclHysteresisTwoPhaseLawParams.hpp:270
void setDeltaSwImbKrw(Scalar)
Sets the saturation value which must be added if krw is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:288
void setKrwSwMdc(Scalar)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition: EclHysteresisTwoPhaseLawParams.hpp:252
void finalize()
Calculate all dependent quantities once the independent quantities of the parameter object have been ...
Definition: EclHysteresisTwoPhaseLawParams.hpp:91
void setDrainageParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the drainage curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:120
Scalar deltaSwImbKrn() const
Returns the saturation value which must be added if krn is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:320
void setImbibitionParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:179
Scalar deltaSwImbKrw() const
Returns the saturation value which must be added if krw is calculated using the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:299
const EffLawParams & imbibitionParams() const
Returns the parameters used for the imbibition curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:225
Scalar krwSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition: EclHysteresisTwoPhaseLawParams.hpp:261
const EffLawParams & drainageParams() const
Returns the parameters used for the drainage curve.
Definition: EclHysteresisTwoPhaseLawParams.hpp:170
void update(Scalar pcSw, Scalar, Scalar krnSw)
Notify the hysteresis law that a given wetting-phase saturation has been seen.
Definition: EclHysteresisTwoPhaseLawParams.hpp:369
Default implementation for asserting finalization of parameter objects.
Definition: EnsureFinalized.hpp:47
void finalize()
Mark the object as finalized.
Definition: EnsureFinalized.hpp:75
Class for (de-)serializing.
Definition: Serializer.hpp:84
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition: Exceptions.hpp:30
EclTwoPhaseSystemType
Specified which fluids are involved in a given twophase material law for endpoint scaling.
Definition: EclEpsConfig.hpp:42
This structure represents all values which can be possibly used as scaling points in the endpoint sca...
Definition: EclEpsScalingPoints.hpp:57