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
30#include <opm/input/eclipse/EclipseState/WagHysteresisConfig.hpp>
31
37
38#include <cassert>
39#include <cmath>
40#include <memory>
41namespace Opm {
48template <class EffLawT>
50{
51 using EffLawParams = typename EffLawT::Params;
52 using Scalar = typename EffLawParams::Traits::Scalar;
53
54public:
55 using Traits = typename EffLawParams::Traits;
56
58 {
59 // These are initialized to two (even though they represent saturations)
60 // to signify that the values are larger than physically possible, and force
61 // using the drainage curve before the first saturation update.
62 pcSwMdc_ = 2.0;
63 krnSwMdc_ = 2.0;
64 krwSwMdc_ = -2.0;
65 krnSwDrainRevert_ = 2.0;
66 krnSwDrainStart_ = -2.0;
67 krnSwWAG_ = 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.deltaSwImbKrw_ = 1.0;
92 result.Sncrt_ = 2.0;
93 result.Swcrt_ = 2.5;
94 result.initialImb_ = true;
95 result.pcSwMic_ = 3.0;
96 result.krnSwMdc_ = 4.0;
97 result.krwSwMdc_ = 4.5;
98 result.KrndHy_ = 5.0;
99 result.KrwdHy_ = 6.0;
100
101 return result;
102 }
103
108 void finalize()
109 {
110 if (config().enableHysteresis()) {
111 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().krHysteresisModel() == 4 || config().pcHysteresisModel() == 0) {
112 C_ = 1.0/(Sncri_ - Sncrd_ + 1.0e-12) - 1.0/(Snmaxd_ - Sncrd_);
113 curvatureCapPrs_ = config().curvatureCapPrs();
114 }
115 if (config().krHysteresisModel() == 4) {
116 Cw_ = 1.0/(Swcri_ - Swcrd_ + 1.0e-12) - 1.0/(Swmaxd_ - Swcrd_);
117 }
118 updateDynamicParams_();
119 }
120
121 EnsureFinalized :: finalize();
122 }
123
127 void setConfig(std::shared_ptr<EclHysteresisConfig> value)
128 { config_ = *value; }
129
134 { return config_; }
135
139 void setWagConfig(std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> value)
140 {
141 wagConfig_ = value;
142 cTransf_ = wagConfig().wagLandsParam();
143 }
144
149 { return *wagConfig_; }
150
154 void setDrainageParams(const EffLawParams& value,
156 EclTwoPhaseSystemType twoPhaseSystem)
157 {
158 drainageParams_ = value;
159
160 oilWaterSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
161 gasOilSystem_ = (twoPhaseSystem == EclTwoPhaseSystemType::GasOil);
162
163 if (!config().enableHysteresis())
164 return;
165 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().krHysteresisModel() == 4 || config().pcHysteresisModel() == 0 || gasOilHysteresisWAG()) {
166 Swco_ = info.Swl;
167 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
168 Sncrd_ = info.Sgcr+info.Swl;
169 Swcrd_ = info.Sogcr + info.Swl;
170 Snmaxd_ = info.Sgu+info.Swl;
171 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
172 }
173 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
174 Sncrd_ = info.Sgcr;
175 Swcrd_ = info.Swcr;
176 Snmaxd_ = info.Sgu;
177 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
178 }
179 else {
180 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
181 Sncrd_ = info.Sowcr;
182 Swcrd_ = info.Swcr;
183 Snmaxd_ = 1.0 - info.Swl - info.Sgl;
184 KrndMax_ = EffLawT::twoPhaseSatKrn(drainageParams(), 1.0-Snmaxd_);
185 }
186 }
187
188 if (config().krHysteresisModel() == 4) {
189 //Swco_ = info.Swl;
190 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
191 Swmaxd_ = 1.0 - info.Sgl - info.Swl;
192 KrwdMax_ = EffLawT::twoPhaseSatKrw(drainageParams(), Swmaxd_);
193 }
194 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
195 Swmaxd_ = info.Swu;
196 KrwdMax_ = EffLawT::twoPhaseSatKrw(drainageParams(), Swmaxd_);
197 }
198 else {
199 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
200 Swmaxd_ = info.Swu;
201 KrwdMax_ = EffLawT::twoPhaseSatKrw(drainageParams(), Swmaxd_);
202 }
203 }
204
205 // Additional Killough hysteresis model for pc
206 if (config().pcHysteresisModel() == 0) {
207 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
208 pcmaxd_ = info.maxPcgo;
209 } else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
210 pcmaxd_ = info.maxPcgo + info.maxPcow;
211 }
212 else {
213 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
214 pcmaxd_ = -17.0; // At this point 'info.maxPcow' holds pre-swatinit value ...;
215 }
216 }
217
218 // For WAG hysteresis, assume initial state along primary drainage curve.
219 if (gasOilHysteresisWAG()) {
220 swatImbStart_ = Swco_;
221 swatImbStartNxt_ = -1.0; // Trigger check for saturation gt Swco at first update ...
222 cTransf_ = wagConfig().wagLandsParam();
223 krnSwDrainStart_ = Sncrd_;
224 krnSwDrainStartNxt_ = Sncrd_;
225 krnImbStart_ = 0.0;
226 krnImbStartNxt_ = 0.0;
227 krnDrainStart_ = 0.0;
228 krnDrainStartNxt_ = 0.0;
229 isDrain_ = true;
230 wasDrain_ = true;
231 krnSwImbStart_ = Sncrd_;
232 SncrtWAG_ = Sncrd_;
233 nState_ = 1;
234 }
235 }
236
240 const EffLawParams& drainageParams() const
241 { return drainageParams_; }
242
243 EffLawParams& drainageParams()
244 { return drainageParams_; }
245
249 void setImbibitionParams(const EffLawParams& value,
251 EclTwoPhaseSystemType twoPhaseSystem)
252 {
253 imbibitionParams_ = value;
254
255 if (!config().enableHysteresis())
256 return;
257
258 // Store critical nonwetting saturation
259 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
260 Sncri_ = info.Sgcr+info.Swl;
261 Swcri_ = info.Sogcr;
262 }
263 else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
264 Sncri_ = info.Sgcr;
265 Swcri_ = info.Swcr;
266 }
267 else {
268 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
269 Sncri_ = info.Sowcr;
270 Swcri_ = info.Swcr;
271 }
272
273
274 // Killough hysteresis model for pc
275 if (config().pcHysteresisModel() == 0) {
276 if (twoPhaseSystem == EclTwoPhaseSystemType::GasOil) {
277 Swmaxi_ = 1.0 - info.Sgl - info.Swl;
278 pcmaxi_ = info.maxPcgo;
279 } else if (twoPhaseSystem == EclTwoPhaseSystemType::GasWater) {
280 Swmaxi_ = 1.0 - info.Sgl;
281 pcmaxi_ = info.maxPcgo + info.maxPcow;
282 }
283 else {
284 assert(twoPhaseSystem == EclTwoPhaseSystemType::OilWater);
285 Swmaxi_ = info.Swu;
286 pcmaxi_ = info.maxPcow;
287 }
288 }
289 }
290
294 const EffLawParams& imbibitionParams() const
295 { return imbibitionParams_; }
296
297 EffLawParams& imbibitionParams()
298 { return imbibitionParams_; }
299
304 Scalar pcSwMdc() const
305 { return pcSwMdc_; }
306
307 Scalar pcSwMic() const
308 { return pcSwMic_; }
309
313 bool initialImb() const
314 { return initialImb_; }
315
321 void setKrwSwMdc(Scalar value)
322 { krwSwMdc_ = value; };
323
329 Scalar krwSwMdc() const
330 { return krwSwMdc_; };
331
337 void setKrnSwMdc(Scalar value)
338 { krnSwMdc_ = value; }
339
345 Scalar krnSwMdc() const
346 { return krnSwMdc_; }
347
355 //void setDeltaSwImbKrw(Scalar value)
356 //{ deltaSwImbKrw_ = value; }
357
365 //Scalar deltaSwImbKrw() const
366 //{ return deltaSwImbKrw_; }
367
375 void setDeltaSwImbKrn(Scalar value)
376 { deltaSwImbKrn_ = value; }
377
385 Scalar deltaSwImbKrn() const
386 { return deltaSwImbKrn_; }
387
388
389 Scalar Swcri() const
390 { return Swcri_; }
391
392 Scalar Swcrd() const
393 { return Swcrd_; }
394
395 Scalar Swmaxi() const
396 { return Swmaxi_; }
397
398 Scalar Sncri() const
399 { return Sncri_; }
400
401 Scalar Sncrd() const
402 { return Sncrd_; }
403
404 Scalar Sncrt() const
405 { return Sncrt_; }
406
407 Scalar Swcrt() const
408 { return Swcrt_; }
409
410 Scalar SnTrapped() const
411 {
412
413 if(isDrain_)
414 return 0.0;
415
416 // For Killough the trapped saturation is already computed
417 if( config().krHysteresisModel() > 1 )
418 return Sncrt_;
419 else // For Carlson we use the shift to compute it from the critial saturation
420 return Sncri_ + deltaSwImbKrn_;
421 }
422
423 Scalar SwTrapped() const
424 {
425 if( config().krHysteresisModel() == 0 || config().krHysteresisModel() == 2)
426 return Swcrd_;
427
428 if( config().krHysteresisModel() == 1 || config().krHysteresisModel() == 3)
429 return Swcri_;
430
431 // For Killough the trapped saturation is already computed
432 if( config().krHysteresisModel() == 4 )
433 return Swcrt_;
434
435 return 0.0;
436 //else // For Carlson we use the shift to compute it from the critial saturation
437 // return Swcri_ + deltaSwImbKrw_;
438 }
439
440 Scalar SncrtWAG() const
441 { return SncrtWAG_; }
442
443 Scalar Snmaxd() const
444 { return Snmaxd_; }
445
446 Scalar Swmaxd() const
447 { return Swmaxd_; }
448
449 Scalar Snhy() const
450 { return 1.0 - krnSwMdc_; }
451
452 Scalar Swhy() const
453 { return krwSwMdc_; }
454
455 Scalar Swco() const
456 { return Swco_; }
457
458 Scalar krnWght() const
459 { return KrndHy_/KrndMax_; }
460
461 Scalar krwWght() const
462 {
463 return KrwdHy_/KrwdMax_; }
464
465 Scalar krwdMax() const
466 {
467 return KrwdMax_; }
468
469 Scalar pcWght() const // Aligning pci and pcd at Swir
470 {
471 if (pcmaxd_ < 0.0)
472 return EffLawT::twoPhaseSatPcnw(drainageParams(), 0.0)/(pcmaxi_+1e-6);
473 else
474 return pcmaxd_/(pcmaxi_+1e-6);
475 }
476
477 Scalar curvatureCapPrs() const
478 { return curvatureCapPrs_;}
479
480 bool gasOilHysteresisWAG() const
481 { return (config().enableWagHysteresis() && gasOilSystem_ && wagConfig().wagGasFlag()) ; }
482
483 Scalar reductionDrain() const
484 { return std::pow(Swco_/(swatImbStart_+tolWAG_*wagConfig().wagWaterThresholdSaturation()), wagConfig().wagSecondaryDrainageReduction());}
485
486 Scalar reductionDrainNxt() const
487 { return std::pow(Swco_/(swatImbStartNxt_+tolWAG_*wagConfig().wagWaterThresholdSaturation()), wagConfig().wagSecondaryDrainageReduction());}
488
489 bool threePhaseState() const
490 { return (swatImbStart_ > (Swco_ + wagConfig().wagWaterThresholdSaturation()) ); }
491
492 Scalar nState() const
493 { return nState_;}
494
495 Scalar krnSwDrainRevert() const
496 { return krnSwDrainRevert_;}
497
498 Scalar krnDrainStart() const
499 { return krnDrainStart_;}
500
501 Scalar krnDrainStartNxt() const
502 { return krnDrainStartNxt_;}
503
504 Scalar krnImbStart() const
505 { return krnImbStart_;}
506
507 Scalar krnImbStartNxt() const
508 { return krnImbStartNxt_;}
509
510 Scalar krnSwWAG() const
511 { return krnSwWAG_;}
512
513 Scalar krnSwDrainStart() const
514 { return krnSwDrainStart_;}
515
516 Scalar krnSwDrainStartNxt() const
517 { return krnSwDrainStartNxt_;}
518
519 Scalar krnSwImbStart() const
520 { return krnSwImbStart_;}
521
522 Scalar tolWAG() const
523 { return tolWAG_;}
524
525 template <class Evaluation>
526 Evaluation computeSwf(const Evaluation& Sw) const
527 {
528 Evaluation SgT = 1.0 - Sw - SncrtWAG(); // Sg-Sg_crit_trapped
529 Scalar SgCut = wagConfig().wagImbCurveLinearFraction()*(Snhy()- SncrtWAG());
530 Evaluation Swf = 1.0;
531 //Scalar C = wagConfig().wagLandsParam();
532 Scalar C = cTransf_;
533
534 if (SgT > SgCut) {
535 Swf -= (Sncrd() + 0.5*( SgT + Opm::sqrt( SgT*SgT + 4.0/C*SgT))); // 1-Sgf
536 }
537 else {
538 SgCut = std::max(Scalar(0.000001), SgCut);
539 Scalar SgCutValue = 0.5*( SgCut + Opm::sqrt( SgCut*SgCut + 4.0/C*SgCut));
540 Scalar SgCutSlope = SgCutValue/SgCut;
541 SgT *= SgCutSlope;
542 Swf -= (Sncrd() + SgT);
543 }
544
545 return Swf;
546 }
547
548 template <class Evaluation>
549 Evaluation computeKrImbWAG(const Evaluation& Sw) const
550 {
551 Evaluation Swf = Sw;
552 if (nState_ <= 2) // Skipping for "higher order" curves seems consistent with benchmark, further investigations needed ...
553 Swf = computeSwf(Sw);
554 if (Swf <= krnSwDrainStart_) { // Use secondary drainage curve
555 Evaluation Krg = EffLawT::twoPhaseSatKrn(drainageParams_, Swf);
556 Evaluation KrgImb2 = (Krg-krnDrainStart_)*reductionDrain() + krnImbStart_;
557 return KrgImb2;
558 }
559 else { // Fallback to primary drainage curve
560 Evaluation Sn = Sncrd_;
561 if (Swf < 1.0-SncrtWAG_) {
562 // Notation: Sn.. = Sg.. + Swco
563 Evaluation dd = (1.0-krnSwImbStart_ - Sncrd_) / (1.0-krnSwDrainStart_ - SncrtWAG_);
564 Sn += (1.0-Swf-SncrtWAG_)*dd;
565 }
566 Evaluation KrgDrn1 = EffLawT::twoPhaseSatKrn(drainageParams_, 1.0 - Sn);
567 return KrgDrn1;
568 }
569 }
570
577 bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
578 {
579 bool updateParams = false;
580
581 if (config().pcHysteresisModel() == 0 && pcSw < pcSwMdc_) {
582 if (pcSwMdc_ == 2.0 && pcSw+1.0e-6 < Swcrd_ && oilWaterSystem_) {
583 initialImb_ = true;
584 }
585 pcSwMdc_ = pcSw;
586 updateParams = true;
587 }
588
589 if (initialImb_ && pcSw > pcSwMic_) {
590 pcSwMic_ = pcSw;
591 updateParams = true;
592 }
593
594 if (krnSw < krnSwMdc_) {
595 krnSwMdc_ = krnSw;
596 KrndHy_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
597 updateParams = true;
598 }
599
600 if (krwSw > krwSwMdc_) {
601 krwSwMdc_ = krwSw;
602 KrwdHy_ = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
603 updateParams = true;
604 }
605
606 // for non WAG hysteresis we still keep track of the process
607 // for output purpose.
608 if (!gasOilHysteresisWAG()) {
609 if (krnSw < krnSwMdc_) {
610 isDrain_ = true;
611 } else {
612 isDrain_ = false;
613 }
614 } else {
615 wasDrain_ = isDrain_;
616
617 if (swatImbStartNxt_ < 0.0) { // Initial check ...
618 swatImbStartNxt_ = std::max(Swco_, Swco_ + krnSw - krwSw);
619 // check if we are in threephase state sw > swco + tolWag and so > tolWag
620 // (sw = swco + krnSw - krwSw and so = krwSw for oil/gas params)
621 if ( (swatImbStartNxt_ > Swco_ + tolWAG_) && krwSw > tolWAG_) {
622 swatImbStart_ = swatImbStartNxt_;
623 krnSwWAG_ = krnSw;
624 krnSwDrainStartNxt_ = krnSwWAG_;
625 krnSwDrainStart_ = krnSwDrainStartNxt_;
626 wasDrain_ = false; // Signal start from threephase state ...
627 }
628 }
629
630 if (isDrain_) {
631 if (krnSw <= krnSwWAG_+tolWAG_) { // continue along drainage curve
632 krnSwWAG_ = std::min(krnSw, krnSwWAG_);
633 krnSwDrainRevert_ = krnSwWAG_;
634 updateParams = true;
635 }
636 else { // start new imbibition curve
637 isDrain_ = false;
638 krnSwWAG_ = krnSw;
639 updateParams = true;
640 }
641 }
642 else {
643 if (krnSw >= krnSwWAG_-tolWAG_) { // continue along imbibition curve
644 krnSwWAG_ = std::max(krnSw, krnSwWAG_);
645 krnSwDrainStartNxt_ = krnSwWAG_;
646 swatImbStartNxt_ = std::max(swatImbStartNxt_, Swco_ + krnSw - krwSw);
647 updateParams = true;
648 }
649 else { // start new drainage curve
650 isDrain_ = true;
651 krnSwDrainStart_ = krnSwDrainStartNxt_;
652 swatImbStart_ = swatImbStartNxt_;
653 krnSwWAG_ = krnSw;
654 updateParams = true;
655 }
656 }
657
658 }
659
660 if (updateParams)
661 updateDynamicParams_();
662
663 return updateParams;
664 }
665
666 template<class Serializer>
667 void serializeOp(Serializer& serializer)
668 {
669 // only serializes dynamic state - see update() and updateDynamic_()
670 serializer(deltaSwImbKrn_);
671 //serializer(deltaSwImbKrw_);
672 serializer(Sncrt_);
673 serializer(Swcrt_);
674 serializer(initialImb_);
675 serializer(pcSwMic_);
676 serializer(krnSwMdc_);
677 serializer(krwSwMdc_);
678 serializer(KrndHy_);
679 serializer(KrwdHy_);
680 }
681
682 bool operator==(const EclHysteresisTwoPhaseLawParams& rhs) const
683 {
684 return this->deltaSwImbKrn_ == rhs.deltaSwImbKrn_ &&
685 //this->deltaSwImbKrw_ == rhs.deltaSwImbKrw_ &&
686 this->Sncrt_ == rhs.Sncrt_ &&
687 this->Swcrt_ == rhs.Swcrt_ &&
688 this->initialImb_ == rhs.initialImb_ &&
689 this->pcSwMic_ == rhs.pcSwMic_ &&
690 this->krnSwMdc_ == rhs.krnSwMdc_ &&
691 this->krwSwMdc_ == rhs.krwSwMdc_ &&
692 this->KrndHy_ == rhs.KrndHy_ &&
693 this->KrwdHy_ == rhs.KrwdHy_;
694 }
695
696private:
697 void updateDynamicParams_()
698 {
699 // calculate the saturation deltas for the relative permeabilities
700 //if (false) { // we dont support Carlson for wetting phase hysteresis
701 //Scalar krwMdcDrainage = EffLawT::twoPhaseSatKrw(drainageParams(), krwSwMdc_);
702 //Scalar SwKrwMdcImbibition = EffLawT::twoPhaseSatKrwInv(imbibitionParams(), krwMdcDrainage);
703 //deltaSwImbKrw_ = SwKrwMdcImbibition - krwSwMdc_;
704 //}
705
706 if (config().krHysteresisModel() == 0 || config().krHysteresisModel() == 1) {
707 Scalar krnMdcDrainage = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwMdc_);
708 Scalar SwKrnMdcImbibition = EffLawT::twoPhaseSatKrnInv(imbibitionParams(), krnMdcDrainage);
709 deltaSwImbKrn_ = SwKrnMdcImbibition - krnSwMdc_;
710 }
711
712 // Scalar pcMdcDrainage = EffLawT::twoPhaseSatPcnw(drainageParams(), pcSwMdc_);
713 // Scalar SwPcMdcImbibition = EffLawT::twoPhaseSatPcnwInv(imbibitionParams(), pcMdcDrainage);
714 // deltaSwImbPc_ = SwPcMdcImbibition - pcSwMdc_;
715
716 if (config().krHysteresisModel() == 2 || config().krHysteresisModel() == 3 || config().krHysteresisModel() == 4 || config().pcHysteresisModel() == 0) {
717 Scalar Snhy = 1.0 - krnSwMdc_;
718 if (Snhy > Sncrd_) {
719 Sncrt_ = Sncrd_ + (Snhy - Sncrd_)/((1.0+config().modParamTrapped()*(Snmaxd_-Snhy)) + C_*(Snhy - Sncrd_));
720 }
721 else
722 {
723 Sncrt_ = Sncrd_;
724 }
725 }
726
727 if (config().krHysteresisModel() == 4) {
728 Scalar Swhy = krwSwMdc_;
729 if (Swhy >= Swcrd_) {
730 Swcrt_ = Swcrd_ + (Swhy - Swcrd_)/((1.0+config().modParamTrapped()*(Swmaxd_-Swhy)) + Cw_*(Swhy - Swcrd_));
731 } else
732 {
733 Swcrt_ = Swcrd_;
734 }
735 }
736
737
738 if (gasOilHysteresisWAG()) {
739 if (isDrain_ && krnSwMdc_ == krnSwWAG_) {
740 Scalar Snhy = 1.0 - krnSwMdc_;
741 SncrtWAG_ = Sncrd_;
742 if (Snhy > Sncrd_) {
743 SncrtWAG_ += (Snhy - Sncrd_)/(1.0+config().modParamTrapped()*(Snmaxd_-Snhy) + wagConfig().wagLandsParam()*(Snhy - Sncrd_));
744 }
745 }
746
747 if (isDrain_ && (1.0-krnSwDrainRevert_) > SncrtWAG_) { //Reversal from drain to imb
748 cTransf_ = 1.0/(SncrtWAG_-Sncrd_ + 1.0e-12) - 1.0/(1.0-krnSwDrainRevert_-Sncrd_);
749 }
750
751 if (!wasDrain_ && isDrain_) { // Start of new drainage cycle
752 if (threePhaseState() || nState_>1) { // Never return to primary (two-phase) state after leaving
753 nState_ += 1;
754 krnDrainStart_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwDrainStart_);
755 krnImbStart_ = krnImbStartNxt_;
756 // Scanning shift for primary drainage
757 krnSwImbStart_ = EffLawT::twoPhaseSatKrnInv(drainageParams(), krnImbStart_);
758 }
759 }
760
761 if (!wasDrain_ && !isDrain_) { //Moving along current imb curve
762 krnDrainStartNxt_ = EffLawT::twoPhaseSatKrn(drainageParams(), krnSwWAG_);
763 if (threePhaseState()) {
764 krnImbStartNxt_ = computeKrImbWAG(krnSwWAG_);
765 }
766 else {
767 Scalar swf = computeSwf(krnSwWAG_);
768 krnImbStartNxt_ = EffLawT::twoPhaseSatKrn(drainageParams(), swf);
769 }
770 }
771
772 }
773
774 }
775
776 EclHysteresisConfig config_{};
777 std::shared_ptr<WagHysteresisConfig::WagHysteresisConfigRecord> wagConfig_{};
778 EffLawParams imbibitionParams_{};
779 EffLawParams drainageParams_{};
780
781 // largest wettinging phase saturation which is on the main-drainage curve. These are
782 // three different values because the sourounding code can choose to use different
783 // definitions for the saturations for different quantities
784 Scalar krwSwMdc_;
785 Scalar krnSwMdc_;
786 Scalar pcSwMdc_;
787
788 // largest wettinging phase saturation along main imbibition curve
789 Scalar pcSwMic_;
790 // Initial process is imbibition (for initial saturations at or below critical drainage saturation)
791 bool initialImb_;
792
793 bool oilWaterSystem_;
794 bool gasOilSystem_;
795
796
797 // offsets added to wetting phase saturation uf using the imbibition curves need to
798 // be used to calculate the wetting phase relperm, the non-wetting phase relperm and
799 // the capillary pressure
800 //Scalar deltaSwImbKrw_;
801 Scalar deltaSwImbKrn_;
802 //Scalar deltaSwImbPc_;
803
804 // the following uses the conventions of the Eclipse technical description:
805 //
806 // Sncrd_: critical non-wetting phase saturation for the drainage curve
807 // Sncri_: critical non-wetting phase saturation for the imbibition curve
808 // Swcri_: critical wetting phase saturation for the imbibition curve
809 // Swcrd_: critical wetting phase saturation for the drainage curve
810 // Swmaxi_; maximum wetting phase saturation for the imbibition curve
811 // Snmaxd_: non-wetting phase saturation where the non-wetting relperm reaches its
812 // maximum on the drainage curve
813 // C_: factor required to calculate the trapped non-wetting phase saturation using
814 // the Killough approach
815 // Cw_: factor required to calculate the trapped wetting phase saturation using
816 // the Killough approach
817 // Swcod_: connate water saturation value used for wag hysteresis (2. drainage)
818 Scalar Sncrd_{};
819 Scalar Sncri_{};
820 Scalar Swcri_{};
821 Scalar Swcrd_{};
822 Scalar Swmaxi_{};
823 Scalar Snmaxd_{};
824 Scalar Swmaxd_;
825 Scalar C_{};
826 Scalar Cw_{};
827
828 Scalar KrndMax_{}; // Krn_drain(Snmaxd_)
829 Scalar KrwdMax_{}; // Krw_drain(Swmaxd_)
830 Scalar KrndHy_{}; // Krn_drain(1-krnSwMdc_)
831 Scalar KrwdHy_{}; // Krw_drain(1-krwSwMdc_)
832
833 Scalar pcmaxd_; // max pc for drain
834 Scalar pcmaxi_; // max pc for imb
835
836 Scalar curvatureCapPrs_{}; // curvature parameter used for capillary pressure hysteresis
837
838 Scalar Sncrt_{}; // trapped non-wetting phase saturation
839 Scalar Swcrt_{}; // trapped wetting phase saturation
840
841 // Used for WAG hysteresis
842 Scalar Swco_; // Connate water.
843 Scalar swatImbStart_; // Water saturation at start of current drainage curve (end of previous imb curve).
844 Scalar swatImbStartNxt_{}; // Water saturation at start of next drainage curve (end of current imb curve).
845 Scalar krnSwWAG_; // Saturation value after latest completed timestep.
846 Scalar krnSwDrainRevert_; // Saturation value at end of current drainage curve.
847 Scalar cTransf_; // Modified Lands constant used for free gas calculations to obtain consistent scanning curve
848 // when reversion to imb occurs above historical maximum gas saturation (i.e. Sw > krwSwMdc_).
849 Scalar krnSwDrainStart_; // Saturation value at start of current drainage curve (end of previous imb curve).
850 Scalar krnSwDrainStartNxt_{}; // Saturation value at start of current drainage curve (end of previous imb curve).
851 Scalar krnImbStart_{}; // Relperm at start of current drainage curve (end of previous imb curve).
852 Scalar krnImbStartNxt_{}; // Relperm at start of next drainage curve (end of current imb curve).
853 Scalar krnDrainStart_{}; // Primary (input) relperm evaluated at start of current drainage curve.
854 Scalar krnDrainStartNxt_{}; // Primary (input) relperm evaluated at start of next drainage curve.
855 bool isDrain_; // Status is either drainage or imbibition
856 bool wasDrain_{}; // Previous status.
857 Scalar krnSwImbStart_{}; // Saturation value where primary drainage relperm equals krnImbStart_
858
859 int nState_; // Number of cycles. Primary cycle is nState_=1.
860
861 Scalar SncrtWAG_{};
862 Scalar tolWAG_;
863};
864
865} // namespace Opm
866
867#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.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
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:50
Scalar pcSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:304
bool initialImb() const
Status of initial process.
Definition EclHysteresisTwoPhaseLawParams.hpp:313
const EclHysteresisConfig & config() const
Returns the hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:133
void setDeltaSwImbKrn(Scalar value)
Sets the saturation value which must be added if krw is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:375
void setConfig(std::shared_ptr< EclHysteresisConfig > value)
Set the hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:127
Scalar krnSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:345
bool update(Scalar pcSw, Scalar krwSw, Scalar krnSw)
Notify the hysteresis law that a given wetting-phase saturation has been seen.
Definition EclHysteresisTwoPhaseLawParams.hpp:577
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:337
void setWagConfig(std::shared_ptr< WagHysteresisConfig::WagHysteresisConfigRecord > value)
Set the WAG-hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:139
void finalize()
Calculate all dependent quantities once the independent quantities of the parameter object have been ...
Definition EclHysteresisTwoPhaseLawParams.hpp:108
void setKrwSwMdc(Scalar value)
Set the saturation of the wetting phase where the last switch from the main drainage curve (MDC) to i...
Definition EclHysteresisTwoPhaseLawParams.hpp:321
void setDrainageParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the drainage curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:154
const WagHysteresisConfig::WagHysteresisConfigRecord & wagConfig() const
Returns the WAG-hysteresis configuration object.
Definition EclHysteresisTwoPhaseLawParams.hpp:148
Scalar deltaSwImbKrn() const
Returns the saturation value which must be added if krn is calculated using the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:385
void setImbibitionParams(const EffLawParams &value, const EclEpsScalingPointsInfo< Scalar > &info, EclTwoPhaseSystemType twoPhaseSystem)
Sets the parameters used for the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:249
const EffLawParams & imbibitionParams() const
Returns the parameters used for the imbibition curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:294
Scalar krwSwMdc() const
Get the saturation of the wetting phase where the last switch from the main drainage curve to imbibit...
Definition EclHysteresisTwoPhaseLawParams.hpp:329
const EffLawParams & drainageParams() const
Returns the parameters used for the drainage curve.
Definition EclHysteresisTwoPhaseLawParams.hpp:240
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:33