My Project
Loading...
Searching...
No Matches
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 <opm/input/eclipse/EclipseState/EclipseState.hpp>
36#include <cassert>
37#include <cmath>
38#include <memory>
39
40namespace Opm {
47template <class EffLawT>
49{
50 using EffLawParams = typename EffLawT::Params;
51 using Scalar = typename EffLawParams::Traits::Scalar;
52
53public:
54 using Traits = typename EffLawParams::Traits;
55
57 {
58 // These are initialized to two (even though they represent saturations)
59 // to signify that the values are larger than physically possible, and force
60 // using the drainage curve before the first saturation update.
61 pcSwMdc_ = 2.0;
62 krnSwMdc_ = 2.0;
63 krnSwDrainRevert_ = 2.0;
64 krnSwDrainStart_ = -2.0;
65 krnSwWAG_ = 2.0;
66
67 // krwSwMdc_ = 2.0;
68
69 pcSwMic_ = -1.0;
70 initialImb_ = false;
71 oilWaterSystem_ = false;
72 gasOilSystem_ = false;
73 pcmaxd_ = 0.0;
74 pcmaxi_ = 0.0;
75
76 deltaSwImbKrn_ = 0.0;
77 // deltaSwImbKrw_ = 0.0;
78
79 Swco_ = 0.0;
80 swatImbStart_ = 0.0;
81 isDrain_ = true;
82 cTransf_ = 0.0;
83 tolWAG_ = 0.001;
84 nState_ = 0;
85 }
86
87 static EclHysteresisTwoPhaseLawParams serializationTestObject()
88 {
90 result.deltaSwImbKrn_ = 1.0;
91 result.Sncrt_ = 2.0;
92 result.initialImb_ = true;
93 result.pcSwMic_ = 3.0;
94 result.krnSwMdc_ = 4.0;
95 result.KrndHy_ = 5.0;
96
97 return result;
98 }
99
104 void finalize()
105 {
106 if (config().enableHysteresis()) {
107 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
108 C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
109 curvatureCapPrs_ = config().curvatureCapPrs();
110 }
111
112 updateDynamicParams_();
113 }
114
115 EnsureFinalized :: finalize();
116 }
117
121 void setConfig(std::shared_ptr<EclHysteresisConfig> value)
122 { config_ = *value; }
123
128 { return config_; }
129
133 void setWagConfig(std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> value)
134 {
135 wagConfig_ = value;
136 cTransf_ = wagConfig().wagLandsParam();
137 }
138
143 { return *wagConfig_; }
144
148 void setDrainageParams(const EffLawParams& value,
150 EclTwoPhaseSystemType twoPhaseSystem)
151 {
152 drainageParams_ = value;
153
154 oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
155 gasOilSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::GasOil);
156
157 if (!config().enableHysteresis())
158 return;
159 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0 || gasOilHysteresisWAG()) {
160 Swco_ = info.Swl;
161 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
162 Sncrd_ = info.Sgcr+info.Swl;
163 Snmaxd_ = info.Sgu+info.Swl;
164 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
165 }
166 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
167 Sncrd_ = info.Sgcr;
168 Snmaxd_ = info.Sgu;
169 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
170 }
171 else {
172 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
173 Sncrd_ = info.Sowcr;
174 Snmaxd_ = 1.0 - info.Swl - info.Sgl;
175 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
176 }
177 }
178
179 // Additional Killough hysteresis model for pc
180 if (config().pcHysteresisModel() == 0) {
181 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
182 Swcrd_ = info.Sogcr;
183 pcmaxd_ = info.maxPcgo;
184 } else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
185 Swcrd_ = info.Swcr;
186 pcmaxd_ = info.maxPcgo + info.maxPcow;
187 }
188 else {
189 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
190 Swcrd_ = info.Swcr;
191 pcmaxd_ = -17.0; // At this point 'info.maxPcow' holds pre-swatinit value ...;
192 }
193 }
194
195 // For WAG hysteresis, assume initial state along primary drainage curve.
196 if (gasOilHysteresisWAG()) {
197 swatImbStart_ = Swco_;
198 swatImbStartNxt_ = -1.0; // Trigger check for saturation gt Swco at first update ...
199 cTransf_ = wagConfig().wagLandsParam();
200 krnSwDrainStart_ = Sncrd_;
201 krnSwDrainStartNxt_ = Sncrd_;
202 krnImbStart_ = 0.0;
203 krnImbStartNxt_ = 0.0;
204 krnDrainStart_ = 0.0;
205 krnDrainStartNxt_ = 0.0;
206 isDrain_ = true;
207 wasDrain_ = true;
208 krnSwImbStart_ = Sncrd_;
209 SncrtWAG_ = Sncrd_;
210 nState_ = 1;
211 }
212 }
213
217 const EffLawParams& drainageParams() const
218 { return drainageParams_; }
219
220 EffLawParams& drainageParams()
221 { return drainageParams_; }
222
226 void setImbibitionParams(const EffLawParams& value,
228 EclTwoPhaseSystemType twoPhaseSystem)
229 {
230 imbibitionParams_ = value;
231
232 if (!config().enableHysteresis())
233 return;
234
235 // Store critical nonwetting saturation
236 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
237 Sncri_ = info.Sgcr+info.Swl;
238 }
239 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
240 Sncri_ = info.Sgcr;
241 }
242 else {
243 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
244 Sncri_ = info.Sowcr;
245 }
246
247
248 // Killough hysteresis model for pc
249 if (config().pcHysteresisModel() == 0) {
250 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
251 Swcri_ = info.Sogcr;
252 Swmaxi_ = 1.0 - info.Sgl - info.Swl;
253 pcmaxi_ = info.maxPcgo;
254 } else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
255 Swcri_ = info.Swcr;
256 Swmaxi_ = 1.0 - info.Sgl;
257 pcmaxi_ = info.maxPcgo + info.maxPcow;
258 }
259 else {
260 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
261 Swcri_ = info.Swcr;
262 Swmaxi_ = info.Swu;
263 pcmaxi_ = info.maxPcow;
264 }
265 }
266 }
267
271 const EffLawParams& imbibitionParams() const
272 { return imbibitionParams_; }
273
274 EffLawParams& imbibitionParams()
275 { return imbibitionParams_; }
276
281 Scalar pcSwMdc() const
282 { return pcSwMdc_; }
283
284 Scalar pcSwMic() const
285 { return pcSwMic_; }
286
290 bool initialImb() const
291 { return initialImb_; }
292
298 void setKrwSwMdc(Scalar /* value */)
299 {}
300 // { krwSwMdc_ = value; };
301
307 Scalar krwSwMdc() const
308 { return 2.0; }
309 // { return krwSwMdc_; };
310
316 void setKrnSwMdc(Scalar value)
317 { krnSwMdc_ = value; }
318
324 Scalar krnSwMdc() const
325 { return krnSwMdc_; }
326
334 void setDeltaSwImbKrw(Scalar /* value */)
335 {}
336 // { deltaSwImbKrw_ = value; }
337
345 Scalar deltaSwImbKrw() const
346 { return 0.0; }
347// { return deltaSwImbKrw_; }
348
356 void setDeltaSwImbKrn(Scalar value)
357 { deltaSwImbKrn_ = value; }
358
366 Scalar deltaSwImbKrn() const
367 { return deltaSwImbKrn_; }
368
369
370 Scalar Swcri() const
371 { return Swcri_; }
372
373 Scalar Swcrd() const
374 { return Swcrd_; }
375
376 Scalar Swmaxi() const
377 { return Swmaxi_; }
378
379 Scalar Sncri() const
380 { return Sncri_; }
381
382 Scalar Sncrd() const
383 { return Sncrd_; }
384
385 Scalar Sncrt() const
386 { return Sncrt_; }
387
388 Scalar SnTrapped() const
389 {
390 // For Killough the trapped saturation is already computed
391 if( config().krHysteresisModel() > 1 )
392 return Sncrt_;
393 else // For Carlson we use the shift to compute it from the critial saturation
394 return Sncri_ + deltaSwImbKrn_;
395 }
396 Scalar SncrtWAG() const
397 { return SncrtWAG_; }
398
399 Scalar Snmaxd() const
400 { return Snmaxd_; }
401
402 Scalar Snhy() const
403 { return 1.0 - krnSwMdc_; }
404
405 Scalar Swco() const
406 { return Swco_; }
407
408 Scalar krnWght() const
409 { return KrndHy_/KrndMax_; }
410
411 Scalar pcWght() const // Aligning pci and pcd at Swir
412 {
413 if (pcmaxd_ < 0.0)
414 return EffLawT::twoPhaseSatPcnw(drainageParams(), 0.0)/(pcmaxi_+1e-6);
415 else
416 return pcmaxd_/(pcmaxi_+1e-6);
417 }
418
419 Scalar curvatureCapPrs() const
420 { return curvatureCapPrs_;}
421
422 bool gasOilHysteresisWAG() const
423 { return (config().enableWagHysteresis() && gasOilSystem_ && wagConfig().wagGasFlag()) ; }
424
425 Scalar reductionDrain() const
426 { return std::pow(Swco_/(swatImbStart_+tolWAG_*wagConfig().wagWaterThresholdSaturation()), wagConfig().wagSecondaryDrainageReduction());}
427
428 Scalar reductionDrainNxt() const
429 { return std::pow(Swco_/(swatImbStartNxt_+tolWAG_*wagConfig().wagWaterThresholdSaturation()), wagConfig().wagSecondaryDrainageReduction());}
430
431 bool threePhaseState() const
432 { return (swatImbStart_ > (Swco_ + wagConfig().wagWaterThresholdSaturation()) ); }
433
434 Scalar nState() const
435 { return nState_;}
436
437 Scalar krnSwDrainRevert() const
438 { return krnSwDrainRevert_;}
439
440 Scalar krnDrainStart() const
441 { return krnDrainStart_;}
442
443 Scalar krnDrainStartNxt() const
444 { return krnDrainStartNxt_;}
445
446 Scalar krnImbStart() const
447 { return krnImbStart_;}
448
449 Scalar krnImbStartNxt() const
450 { return krnImbStartNxt_;}
451
452 Scalar krnSwWAG() const
453 { return krnSwWAG_;}
454
455 Scalar krnSwDrainStart() const
456 { return krnSwDrainStart_;}
457
458 Scalar krnSwDrainStartNxt() const
459 { return krnSwDrainStartNxt_;}
460
461 Scalar krnSwImbStart() const
462 { return krnSwImbStart_;}
463
464 Scalar tolWAG() const
465 { return tolWAG_;}
466
467 template <class Evaluation>
468 Evaluation computeSwf(const Evaluation& Sw) const
469 {
470 Evaluation SgT = 1.0 - Sw - SncrtWAG(); // Sg-Sg_crit_trapped
471 Scalar SgCut = wagConfig().wagImbCurveLinearFraction()*(Snhy()- SncrtWAG());
472 Evaluation Swf = 1.0;
473 //Scalar C = wagConfig().wagLandsParam();
474 Scalar C = cTransf_;
475
476 if (SgT > SgCut) {
477 Swf -= (Sncrd() + 0.5*( SgT + Opm::sqrt( SgT*SgT + 4.0/C*SgT))); // 1-Sgf
478 }
479 else {
480 SgCut = std::max(Scalar(0.000001), SgCut);
481 Scalar SgCutValue = 0.5*( SgCut + Opm::sqrt( SgCut*SgCut + 4.0/C*SgCut));
482 Scalar SgCutSlope = SgCutValue/SgCut;
483 SgT *= SgCutSlope;
484 Swf -= (Sncrd() + SgT);
485 }
486
487 return Swf;
488 }
489
490 template <class Evaluation>
491 Evaluation computeKrImbWAG(const Evaluation& Sw) const
492 {
493 Evaluation Swf = Sw;
494 if (nState_ <= 2) // Skipping for "higher order" curves seems consistent with benchmark, further investigations needed ...
495 Swf = computeSwf(Sw);
496 if (Swf <= krnSwDrainStart_) { // Use secondary drainage curve
497 Evaluation Krg = EffLawT::twoPhaseSatKrn(drainageParams_, Swf);
498 Evaluation KrgImb2 = (Krg-krnDrainStart_)*reductionDrain() + krnImbStart_;
499 return KrgImb2;
500 }
501 else { // Fallback to primary drainage curve
502 Evaluation Sn = Sncrd_;
503 if (Swf < 1.0-SncrtWAG_) {
504 // Notation: Sn.. = Sg.. + Swco
505 Evaluation dd = (1.0-krnSwImbStart_ - Sncrd_) / (1.0-krnSwDrainStart_ - SncrtWAG_);
506 Sn += (1.0-Swf-SncrtWAG_)*dd;
507 }
508 Evaluation KrgDrn1 = EffLawT::twoPhaseSatKrn(drainageParams_, 1.0 - Sn);
509 return KrgDrn1;
510 }
511 }
512
519 bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
520 {
521 bool updateParams = false;
522
523 if (config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
524 if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
525 initialImb_ = true;
526 }
527 pcSwMdc_ = pcSw;
528 updateParams = true;
529 }
530
531 if (initialImb_ && pcSw > pcSwMic_) {
532 pcSwMic_ = pcSw;
533 updateParams = true;
534 }
535
536/*
537 // This is quite hacky: Eclipse says that it only uses relperm hysteresis for the
538 // wetting phase (indicated by '0' for the second item of the EHYSTER keyword),
539 // even though this makes about zero sense: one would expect that hysteresis can
540 // be limited to the oil phase, but the oil phase is the wetting phase of the
541 // gas-oil twophase system whilst it is non-wetting for water-oil.
542 if (krwSw < krwSwMdc_)
543 {
544 krwSwMdc_ = krwSw;
545 updateParams = true;
546 }
547*/
548
549 if (krnSw < krnSwMdc_) {
550 krnSwMdc_ = krnSw;
551 KrndHy_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
552 updateParams = true;
553 }
554
555 if (gasOilHysteresisWAG()) {
556
557 wasDrain_ = isDrain_;
558
559 if (swatImbStartNxt_ < 0.0) { // Initial check ...
560 swatImbStartNxt_ = std::max(Swco_, Swco_ + krnSw - krwSw);
561 // check if we are in threephase state sw > swco + tolWag and so > tolWag
562 // (sw = swco + krnSw - krwSw and so = krwSw for oil/gas params)
563 if ( (swatImbStartNxt_ > Swco_ + tolWAG_) && krwSw > tolWAG_) {
564 swatImbStart_ = swatImbStartNxt_;
565 krnSwWAG_ = krnSw;
566 krnSwDrainStartNxt_ = krnSwWAG_;
567 krnSwDrainStart_ = krnSwDrainStartNxt_;
568 wasDrain_ = false; // Signal start from threephase state ...
569 }
570 }
571
572 if (isDrain_) {
573 if (krnSw <= krnSwWAG_+tolWAG_) { // continue along drainage curve
574 krnSwWAG_ = std::min(krnSw, krnSwWAG_);
575 krnSwDrainRevert_ = krnSwWAG_;
576 updateParams = true;
577 }
578 else { // start new imbibition curve
579 isDrain_ = false;
580 krnSwWAG_ = krnSw;
581 updateParams = true;
582 }
583 }
584 else {
585 if (krnSw >= krnSwWAG_-tolWAG_) { // continue along imbibition curve
586 krnSwWAG_ = std::max(krnSw, krnSwWAG_);
587 krnSwDrainStartNxt_ = krnSwWAG_;
588 swatImbStartNxt_ = std::max(swatImbStartNxt_, Swco_ + krnSw - krwSw);
589 updateParams = true;
590 }
591 else { // start new drainage curve
592 isDrain_ = true;
593 krnSwDrainStart_ = krnSwDrainStartNxt_;
594 swatImbStart_ = swatImbStartNxt_;
595 krnSwWAG_ = krnSw;
596 updateParams = true;
597 }
598 }
599
600 }
601
602 if (updateParams)
603 updateDynamicParams_();
604
605 return updateParams;
606 }
607
608 template<class Serializer>
609 void serializeOp(Serializer& serializer)
610 {
611 // only serializes dynamic state - see update() and updateDynamic_()
612 serializer(deltaSwImbKrn_);
613 serializer(Sncrt_);
614 serializer(initialImb_);
615 serializer(pcSwMic_);
616 serializer(krnSwMdc_);
617 serializer(KrndHy_);
618 }
619
620 bool operator==(const EclHysteresisTwoPhaseLawParams& rhs) const
621 {
622 return this->deltaSwImbKrn_ == rhs.deltaSwImbKrn_ &&
623 this->Sncrt_ == rhs.Sncrt_ &&
624 this->initialImb_ == rhs.initialImb_ &&
625 this->pcSwMic_ == rhs.pcSwMic_ &&
626 this->krnSwMdc_ == rhs.krnSwMdc_ &&
627 this->KrndHy_ == rhs.KrndHy_;
628 }
629
630private:
631 void updateDynamicParams_()
632 {
633 // HACK: Eclipse seems to disable the wetting-phase relperm even though this is
634 // quite pointless from the physical POV. (see comment above)
635/*
636 // calculate the saturation deltas for the relative permeabilities
637 Scalar krwMdcDrainage = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
638 Scalar SwKrwMdcImbibition = EffLawT::twoPhaseSatKrwInv(imbibitionParams(), krwMdcDrainage);
639 deltaSwImbKrw_ = SwKrwMdcImbibition - krwSwMdc_;
640*/
641 if (config().krHysteresisModel() == 0 || config().krHysteresisModel() == 1) {
642 Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
643 Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage);
644 deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
645 assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
646 - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
647 }
648
649 // Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_);
650 // Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage);
651 // deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_;
652
653// assert(std::abs(EffLawT::twoPhaseSatPcnw(imbibitionParams(), pcSwMdc_ + deltaSwImbPc_)
654// - EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_)) < 1e-8);
655// assert(std::abs(EffLawT::twoPhaseSatKrn(imbibitionParams(), krnSwMdc_ + deltaSwImbKrn_)
656// - EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_)) < 1e-8);
657// assert(std::abs(EffLawT::twoPhaseSatKrw(imbibitionParams(), krwSwMdc_ + deltaSwImbKrw_)
658// - EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_)) < 1e-8);
659
660 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().pcHysteresisModel() == 0) {
661 Scalar Snhy = 1.0 - krnSwMdc_;
662 if (Snhy > Sncrd_)
663 Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/((1.0+config().modParamTrapped()*(Snmaxd_-Snhy)) + C_*(Snhy - Sncrd_));
664 else
665 {
666 Sncrt_ = Sncrd_;
667 }
668 }
669
670
671 if (gasOilHysteresisWAG()) {
672 if (isDrain_ && krnSwMdc_ == krnSwWAG_) {
673 Scalar Snhy = 1.0 - krnSwMdc_;
674 SncrtWAG_ = Sncrd_;
675 if (Snhy > Sncrd_) {
676 SncrtWAG_ += (Snhy - Sncrd_)/(1.0+config().modParamTrapped()*(Snmaxd_-Snhy) + wagConfig().wagLandsParam()*(Snhy - Sncrd_));
677 }
678 }
679
680 if (isDrain_ && (1.0-krnSwDrainRevert_) > SncrtWAG_) { //Reversal from drain to imb
681 cTransf_ = 1.0/(SncrtWAG_-Sncrd_ + 1.0e-12) - 1.0/(1.0-krnSwDrainRevert_-Sncrd_);
682 }
683
684 if (!wasDrain_ && isDrain_) { // Start of new drainage cycle
685 if (threePhaseState() || nState_>1) { // Never return to primary (two-phase) state after leaving
686 nState_ += 1;
687 krnDrainStart_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwDrainStart_);
688 krnImbStart_ = krnImbStartNxt_;
689 // Scanning shift for primary drainage
690 krnSwImbStart_ = EffLawT::twoPhaseSatKrnInv(drainageParams(), krnImbStart_);
691 }
692 }
693
694 if (!wasDrain_ && !isDrain_) { //Moving along current imb curve
695 krnDrainStartNxt_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwWAG_);
696 if (threePhaseState()) {
697 krnImbStartNxt_ = computeKrImbWAG(krnSwWAG_);
698 }
699 else {
700 Scalar swf = computeSwf(krnSwWAG_);
701 krnImbStartNxt_ = EffLawT::twoPhaseSatKrn(drainageParams(), swf);
702 }
703 }
704
705 }
706
707 }
708
709 EclHysteresisConfig config_;
710 std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> wagConfig_;
711 EffLawParams imbibitionParams_;
712 EffLawParams drainageParams_;
713
714 // largest wettinging phase saturation which is on the main-drainage curve. These are
715 // three different values because the sourounding code can choose to use different
716 // definitions for the saturations for different quantities
717 Scalar krwSwMdc_;
718 Scalar krnSwMdc_;
719 Scalar pcSwMdc_;
720
721 // largest wettinging phase saturation along main imbibition curve
722 Scalar pcSwMic_;
723 // Initial process is imbibition (for initial saturations at or below critical drainage saturation)
724 bool initialImb_;
725
726 bool oilWaterSystem_;
727 bool gasOilSystem_;
728
729
730 // offsets added to wetting phase saturation uf using the imbibition curves need to
731 // be used to calculate the wetting phase relperm, the non-wetting phase relperm and
732 // the capillary pressure
733 //Scalar deltaSwImbKrw_;
734 Scalar deltaSwImbKrn_;
735 //Scalar deltaSwImbPc_;
736
737 // the following uses the conventions of the Eclipse technical description:
738 //
739 // Sncrd_: critical non-wetting phase saturation for the drainage curve
740 // Sncri_: critical non-wetting phase saturation for the imbibition curve
741 // Swcri_: critical wetting phase saturation for the imbibition curve
742 // Swcrd_: critical wetting phase saturation for the drainage curve
743 // Swmaxi_; maximum wetting phase saturation for the imbibition curve
744 // Snmaxd_: non-wetting phase saturation where the non-wetting relperm reaches its
745 // maximum on the drainage curve
746 // C_: factor required to calculate the trapped non-wetting phase saturation using
747 // the Killough approach
748 // Swcod_: connate water saturation value used for wag hysteresis (2. drainage)
749 Scalar Sncrd_;
750 Scalar Sncri_;
751 Scalar Swcri_;
752 Scalar Swcrd_;
753 Scalar Swmaxi_;
754 Scalar Snmaxd_;
755 Scalar C_;
756
757 Scalar KrndMax_; // Krn_drain(Snmaxd_)
758 Scalar KrndHy_; // Krn_drain(1-krnSwMdc_)
759
760 Scalar pcmaxd_; // max pc for drain
761 Scalar pcmaxi_; // max pc for imb
762
763 Scalar curvatureCapPrs_; // curvature parameter used for capillary pressure hysteresis
764
765 Scalar Sncrt_; // trapped non-wetting phase saturation
766
767 // Used for WAG hysteresis
768 Scalar Swco_; // Connate water.
769 Scalar swatImbStart_; // Water saturation at start of current drainage curve (end of previous imb curve).
770 Scalar swatImbStartNxt_; // Water saturation at start of next drainage curve (end of current imb curve).
771 Scalar krnSwWAG_; // Saturation value after latest completed timestep.
772 Scalar krnSwDrainRevert_; // Saturation value at end of current drainage curve.
773 Scalar cTransf_; // Modified Lands constant used for free gas calculations to obtain consistent scanning curve
774 // when reversion to imb occurs above historical maximum gas saturation (i.e. Sw > krwSwMdc_).
775 Scalar krnSwDrainStart_; // Saturation value at start of current drainage curve (end of previous imb curve).
776 Scalar krnSwDrainStartNxt_; // Saturation value at start of current drainage curve (end of previous imb curve).
777 Scalar krnImbStart_; // Relperm at start of current drainage curve (end of previous imb curve).
778 Scalar krnImbStartNxt_; // Relperm at start of next drainage curve (end of current imb curve).
779 Scalar krnDrainStart_; // Primary (input) relperm evaluated at start of current drainage curve.
780 Scalar krnDrainStartNxt_; // Primary (input) relperm evaluated at start of next drainage curve.
781 bool isDrain_; // Status is either drainage or imbibition
782 bool wasDrain_; // Previous status.
783 Scalar krnSwImbStart_; // Saturation value where primary drainage relperm equals krnImbStart_
784
785 int nState_; // Number of cycles. Primary cycle is nState_=1.
786
787 Scalar SncrtWAG_;
788 Scalar tolWAG_;
789};
790
791} // namespace Opm
792
793#endif
Specifies the configuration used by the endpoint scaling code.
Specifies the configuration used by the ECL kr/pC hysteresis code.
Default implementation for asserting finalization of parameter objects.
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
double modParamTrapped() const
Regularisation parameter used for Killough model.
Definition EclHysteresisConfig.hpp:105
A default implementation of the parameters for the material law which implements the ECL relative per...
Definition EclHysteresisTwoPhaseLawParams.hpp:49
Scalar pcSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:281
bool initialImb() const
Status of initial process.
Definition EclHysteresisTwoPhaseLawParams.hpp:290
const EclHysteresisConfig & config() const
Returns the hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:127
void setDeltaSwImbKrn(Scalar value)
Sets the saturation value which must be added if krn is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:356
void setConfig(std::shared_ptr< EclHysteresisConfig > value)
Set the hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:121
Scalar krnSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:324
bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
Notify the hysteresis law that a given wetting-phase saturation has been seen.
Definition EclHysteresisTwoPhaseLawParams.hpp:519
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:316
void setDeltaSwImbKrw(Scalar)
Sets the saturation value which must be added if krw is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:334
void setWagConfig(std::shared_ptr< WagHysteresisConfig::WagHysteresisConfigRecord > value)
Set the WAG-hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:133
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:298
void finalize()
Calculate all dependent quantities once the independent quantities of the parameter object have been ...
Definition EclHysteresisTwoPhaseLawParams.hpp:104
void setDrainageParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the drainage curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:148
const WagHysteresisConfig::WagHysteresisConfigRecord & wagConfig() const
Returns the WAG-hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:142
Scalar deltaSwImbKrn() const
Returns the saturation value which must be added if krn is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:366
void setImbibitionParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:226
Scalar deltaSwImbKrw() const
Returns the saturation value which must be added if krw is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:345
const EffLawParams & imbibitionParams() const
Returns the parameters used for the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:271
Scalar krwSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:307
const EffLawParams & drainageParams() const
Returns the parameters used for the drainage curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:217
Default implementation for asserting finalization of parameter objects.
Definition EnsureFinalized.hpp:47
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
Definition WagHysteresisConfig.hpp:30