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