My Project
Loading...
Searching...
No Matches
ParkerLenhard.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_PARKER_LENHARD_HPP
28#define OPM_PARKER_LENHARD_HPP
29
31
33
34#include <algorithm>
35#include <cassert>
36
37namespace Opm {
38
46template <class ScalarT>
48{
49public:
50 typedef ScalarT Scalar;
51
58 explicit PLScanningCurve(Scalar Swr)
59 {
60 loopNum_ = 0;
61 prev_ = new PLScanningCurve(nullptr, // prev
62 this, // next
63 -1, // loop number
64 Swr, // Sw
65 1e12, // pcnw
66 Swr, // SwMic
67 Swr); // SwMdc
68 next_ = nullptr;
69
70 Sw_ = 1.0;
71 pcnw_ = 0.0;
72 SwMic_ = 1.0;
73 SwMdc_ = 1.0;
74 }
75
76 PLScanningCurve& operator=(const PLScanningCurve&) = delete;
77
78protected:
80 PLScanningCurve* nextSC,
81 int loopN,
82 Scalar SwReversal,
83 Scalar pcnwReversal,
84 Scalar SwMiCurve,
85 Scalar SwMdCurve)
86 : prev_(prevSC)
87 , next_(nextSC)
88 , loopNum_(loopN)
89 , Sw_(SwReversal)
90 , pcnw_(pcnwReversal)
91 , SwMdc_(SwMdCurve)
92 , SwMic_(SwMiCurve)
93 {
94 }
95
96public:
103 {
104 if (loopNum_ == 0)
105 delete prev_;
106 if (loopNum_ >= 0)
107 delete next_;
108 }
109
115 { return prev_; }
116
122 { return next_; }
123
132 void setNext(Scalar SwReversal,
133 Scalar pcnwReversal,
134 Scalar SwMiCurve,
135 Scalar SwMdCurve)
136 {
137 // if next_ is nullptr, delete does nothing, so
138 // this is valid!!
139 delete next_;
140
141 next_ = new PLScanningCurve(this, // prev
142 nullptr, // next
143 loopNum() + 1,
144 SwReversal,
145 pcnwReversal,
146 SwMiCurve,
147 SwMdCurve);
148 }
149
156 bool isValidAt_Sw(Scalar SwReversal)
157 {
158 if (isImbib())
159 // for inbibition the given saturation
160 // must be between the start of the
161 // current imbibition and the the start
162 // of the last drainage
163 return this->Sw() < SwReversal && SwReversal < prev_->Sw();
164 else
165 // for drainage the given saturation
166 // must be between the start of the
167 // last imbibition and the start
168 // of the current drainage
169 return prev_->Sw() < SwReversal && SwReversal < this->Sw();
170 }
171
176 bool isImbib()
177 { return loopNum()%2 == 1; }
178
183 bool isDrain()
184 { return !isImbib(); }
185
192 { return loopNum_; }
193
198 Scalar Sw() const
199 { return Sw_; }
200
204 Scalar pcnw() const
205 { return pcnw_; }
206
211 Scalar SwMic()
212 { return SwMic_; }
213
218 Scalar SwMdc()
219 { return SwMdc_; }
220
221private:
222 PLScanningCurve* prev_;
223 PLScanningCurve* next_;
224
225 int loopNum_;
226
227 Scalar Sw_;
228 Scalar pcnw_;
229
230 Scalar SwMdc_;
231 Scalar SwMic_;
232};
233
240template <class TraitsT, class ParamsT = ParkerLenhardParams<TraitsT> >
241class ParkerLenhard : public TraitsT
242{
243public:
244 typedef TraitsT Traits;
245 typedef ParamsT Params;
246 typedef typename Traits::Scalar Scalar;
247
249 static const int numPhases = Traits::numPhases;
250 static_assert(numPhases == 2,
251 "The Parker-Lenhard capillary pressure law only "
252 "applies to the case of two fluid phases");
253
256 static const bool implementsTwoPhaseApi = true;
257
260 static const bool implementsTwoPhaseSatApi = true;
261
264 static const bool isSaturationDependent = true;
265
268 static const bool isPressureDependent = false;
269
272 static const bool isTemperatureDependent = false;
273
276 static const bool isCompositionDependent = false;
277
278 static_assert(Traits::numPhases == 2,
279 "The number of fluid phases must be two if you want to use "
280 "this material law!");
281
282private:
283 typedef typename ParamsT::VanGenuchten VanGenuchten;
285
286public:
291 static void reset(Params& params)
292 {
293 delete params.mdc(); // this will work even if mdc_ == nullptr!
294 params.setMdc(new ScanningCurve(params.SwrPc()));
295 params.setCsc(params.mdc());
296 params.setPisc(nullptr);
297 params.setCurrentSnr(0.0);
298 }
299
304 template <class FluidState>
305 static void update(Params& params, const FluidState& fs)
306 {
307 Scalar Sw = scalarValue(fs.saturation(Traits::wettingPhaseIdx));
308
309 if (Sw > 1 - 1e-5) {
310 // if the absolute saturation is almost 1,
311 // it means that we're back to the beginning
312 reset(params);
313 return;
314 }
315
316 // find the loop number which corrosponds to the
317 // given effective saturation
318 ScanningCurve* curve = findScanningCurve_(params, Sw);
319
320 // calculate the apparent saturation on the MIC and MDC
321 // which yield the same capillary pressure as the
322 // Sw at the current scanning curve
323 Scalar pc = pcnw<FluidState, Scalar>(params, fs);
324 Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
325 Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
326
327 curve->setNext(Sw, pc, Sw_mic, Sw_mdc);
328 if (!curve->next())
329 return;
330
331 params.setCsc(curve);
332
333 // if we're back on the MDC, we also have a new PISC!
334 if (params.csc() == params.mdc()) {
335 params.setPisc(params.mdc()->next());
336 params.setCurrentSnr(computeCurrentSnr_(params, Sw));
337 }
338 }
339
344 template <class Container, class FluidState>
345 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
346 {
347 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
348
349 values[Traits::wettingPhaseIdx] = 0.0; // reference phase
350 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
351 }
352
357 template <class Container, class FluidState>
358 static void saturations(Container& /*values*/, const Params& /*params*/, const FluidState& /*fs*/)
359 { throw std::logic_error("Not implemented: ParkerLenhard::saturations()"); }
360
365 template <class Container, class FluidState>
366 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
367 {
368 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
369
370 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
371 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
372 }
373
378 template <class FluidState, class Evaluation = typename FluidState::Scalar>
379 static Evaluation pcnw(const Params& params, const FluidState& fs)
380 {
381 const Evaluation& Sw =
382 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
383
384 return twoPhaseSatPcnw(params, Sw);
385 }
386
387 template <class Evaluation>
388 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
389 {
390 // calculate the current apparent saturation
391 ScanningCurve* sc = findScanningCurve_(params, scalarValue(Sw));
392
393 // calculate the apparant saturation
394 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
395
396 // if the apparent saturation exceeds the 'legal' limits,
397 // we also the underlying material law decide what to do.
398 if (Sw_app > 1) {
399 return 0.0; // VanGenuchten::pcnw(params.mdcParams(), Sw_app);
400 }
401
402 // put the apparent saturation into the main imbibition or
403 // drainage curve
404 Scalar SwAppCurSC = absoluteToApparentSw_(params, sc->Sw());
405 Scalar SwAppPrevSC = absoluteToApparentSw_(params, sc->prev()->Sw());
406 const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
407 if (sc->isImbib()) {
408 const Evaluation& SwMic =
409 pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
410
411 return VanGenuchten::twoPhaseSatPcnw(params.micParams(), SwMic);
412 }
413 else { // sc->isDrain()
414 const Evaluation& SwMdc =
415 pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
416
417 return VanGenuchten::twoPhaseSatPcnw(params.mdcParams(), SwMdc);
418 }
419 }
420
425 template <class FluidState, class Evaluation = typename FluidState::Scalar>
426 static Evaluation Sw(const Params& /*params*/, const FluidState& /*fs*/)
427 { throw std::logic_error("Not implemented: ParkerLenhard::Sw()"); }
428
429 template <class Evaluation>
430 static Evaluation twoPhaseSatSw(const Params& /*params*/, const Evaluation& /*pc*/)
431 { throw std::logic_error("Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
432
437 template <class FluidState, class Evaluation = typename FluidState::Scalar>
438 static Evaluation Sn(const Params& params, const FluidState& fs)
439 { return 1 - Sw<FluidState, Evaluation>(params, fs); }
440
441 template <class Evaluation>
442 static Evaluation twoPhaseSatSn(const Params& /*params*/, const Evaluation& /*pc*/)
443 { throw std::logic_error("Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
444
449 template <class FluidState, class Evaluation = typename FluidState::Scalar>
450 static Evaluation krw(const Params& params, const FluidState& fs)
451 {
452 const Evaluation& Sw =
453 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
454
455 return twoPhaseSatKrw(params, Sw);
456 }
457
458 template <class Evaluation>
459 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
460 {
461 // for the effective permeability we only use Land's law and
462 // the relative permeability of the main drainage curve.
463 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
464 return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
465 }
466
471 template <class FluidState, class Evaluation = typename FluidState::Scalar>
472 static Evaluation krn(const Params& params, const FluidState& fs)
473 {
474 const Evaluation& Sw =
475 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
476
477 return twoPhaseSatKrn(params, Sw);
478 }
479
480 template <class Evaluation>
481 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
482 {
483 // for the effective permeability we only use Land's law and
484 // the relative permeability of the main drainage curve.
485 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
486 return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
487 }
488
492 template <class Evaluation>
493 static Evaluation absoluteToApparentSw_(const Params& params, const Evaluation& Sw)
494 {
495 return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params, Sw));
496 }
497
498private:
508 template <class Evaluation>
509 static Evaluation absoluteToEffectiveSw_(const Params& params, const Evaluation& Sw)
510 { return (Sw - params.SwrPc())/(1 - params.SwrPc()); }
511
521 template <class Evaluation>
522 static Evaluation effectiveToAbsoluteSw_(const Params& params, const Evaluation& Swe)
523 { return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
524
525 // return the effctive residual non-wetting saturation, given an
526 // effective wetting saturation
527 template <class Evaluation>
528 static Evaluation computeCurrentSnr_(const Params& params, const Evaluation& Sw)
529 {
530 // regularize
531 if (Sw > 1 - params.Snr())
532 return 0.0;
533 if (Sw < params.SwrPc())
534 return params.Snr();
535
536 if (params.Snr() == 0.0)
537 return 0.0;
538
539 // use Land's law
540 Scalar R = 1.0/params.Snr() - 1;
541 const Evaluation& curSnr = (1 - Sw)/(1 + R*(1 - Sw));
542
543 // the current effective residual non-wetting saturation must
544 // be smaller than the residual non-wetting saturation
545 assert(curSnr <= params.Snr());
546
547 return curSnr;
548 }
549
550 // returns the trapped effective non-wetting saturation for a
551 // given wetting phase saturation
552 template <class Evaluation>
553 static Evaluation trappedEffectiveSn_(const Params& params, const Evaluation& Sw)
554 {
555 const Evaluation& Swe = absoluteToEffectiveSw_(params, Sw);
556 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
557
558 Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
559 return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
560 }
561
562 // returns the apparent saturation of the wetting phase depending
563 // on the effective saturation
564 template <class Evaluation>
565 static Evaluation effectiveToApparentSw_(const Params& params, const Evaluation& Swe)
566 {
567 if (params.pisc() == nullptr ||
568 Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
569 {
570 // we are on the main drainage curve, i.e. no non-wetting
571 // fluid is trapped -> apparent saturation == effective
572 // saturation
573 return Swe;
574 }
575
576 // we are on a imbibition or drainage curve which is not the
577 // main drainage curve -> apparent saturation == effective
578 // saturation + trapped effective saturation
579 return Swe + trappedEffectiveSn_(params, Swe);
580 }
581
582 // Returns the effective saturation to a given apparent one
583 template <class Evaluation>
584 static Evaluation apparentToEffectiveSw_(const Params& params, const Evaluation& Swapp)
585 {
586 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
587 if (params.pisc() == nullptr || Swapp <= SwePisc) {
588 // we are on the main drainage curve, i.e.
589 // no non-wetting fluid is trapped
590 // -> apparent saturation == effective saturation
591 return Swapp;
592 }
593
594 Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
595 return
596 (Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
597 /(1 - SwePisc);
598 }
599
600
601 // find the loop on which the an effective
602 // saturation has to be
603 static ScanningCurve* findScanningCurve_(const Params& params, Scalar Sw)
604 {
605 if (params.pisc() == nullptr || Sw <= params.pisc()->Sw()) {
606 // we don't have a PISC yet, or the effective
607 // saturation is smaller than the saturation where the
608 // PISC begins. In this case are on the MDC
609 return params.mdc();
610 }
611
612 // if we have a primary imbibition curve, and our current
613 // effective saturation is higher than the beginning of
614 // the secondary drainage curve. this means we are on the
615 // PISC again.
616 if (params.pisc()->next() == nullptr ||
617 params.pisc()->next()->Sw() < Sw)
618 {
619 return params.pisc();
620 }
621
622 ScanningCurve* curve = params.csc()->next();
623 while (true) {
624 assert(curve != params.mdc()->prev());
625 if (curve->isValidAt_Sw(Sw)) {
626 return curve;
627 }
628 curve = curve->prev();
629 }
630 }
631};
632
633} // namespace Opm
634
635#endif // PARKER_LENHARD_HPP
Default parameter class for the Parker-Lenhard hysteresis model.
Implementation of the van Genuchten capillary pressure - saturation relation.
Represents a scanning curve in the Parker-Lenhard hysteresis model.
Definition ParkerLenhard.hpp:48
PLScanningCurve * next() const
Return the next scanning curve, i.e.
Definition ParkerLenhard.hpp:121
~PLScanningCurve()
Destructor.
Definition ParkerLenhard.hpp:102
PLScanningCurve(Scalar Swr)
Constructs main imbibition curve.
Definition ParkerLenhard.hpp:58
bool isDrain()
Returns true iff the scanning curve is a drainage curve.
Definition ParkerLenhard.hpp:183
void setNext(Scalar SwReversal, Scalar pcnwReversal, Scalar SwMiCurve, Scalar SwMdCurve)
Set the next scanning curve.
Definition ParkerLenhard.hpp:132
bool isImbib()
Returns true iff the scanning curve is a imbibition curve.
Definition ParkerLenhard.hpp:176
Scalar Sw() const
Absolute wetting-phase saturation at the scanning curve's reversal point.
Definition ParkerLenhard.hpp:198
bool isValidAt_Sw(Scalar SwReversal)
Returns true iff the given effective saturation Swei is within the scope of the curve,...
Definition ParkerLenhard.hpp:156
Scalar SwMdc()
Apparent saturation of the last reversal point on the pressure MDC.
Definition ParkerLenhard.hpp:218
Scalar pcnw() const
Capillary pressure at the last reversal point.
Definition ParkerLenhard.hpp:204
PLScanningCurve * prev() const
Return the previous scanning curve, i.e.
Definition ParkerLenhard.hpp:114
int loopNum()
The loop number of the scanning curve.
Definition ParkerLenhard.hpp:191
Scalar SwMic()
Apparent saturation of the last reversal point on the pressure MIC.
Definition ParkerLenhard.hpp:211
Implements the Parker-Lenhard twophase p_c-Sw hysteresis model.
Definition ParkerLenhard.hpp:242
static void reset(Params &params)
Resets the hysteresis model to the initial parameters on the main drainage curve.
Definition ParkerLenhard.hpp:291
static Evaluation Sw(const Params &, const FluidState &)
Calculate the wetting phase saturations depending on the phase pressures.
Definition ParkerLenhard.hpp:426
static Evaluation krn(const Params &params, const FluidState &fs)
The relative permeability for the non-wetting phase of the params.
Definition ParkerLenhard.hpp:472
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition ParkerLenhard.hpp:260
static Evaluation krw(const Params &params, const FluidState &fs)
The relative permeability for the wetting phase of the medium.
Definition ParkerLenhard.hpp:450
static void update(Params &params, const FluidState &fs)
Set the current absolute saturation for the current timestep.
Definition ParkerLenhard.hpp:305
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
Returns the relative permeabilities of the phases dependening on the phase saturations.
Definition ParkerLenhard.hpp:366
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition ParkerLenhard.hpp:268
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition ParkerLenhard.hpp:272
static Evaluation absoluteToApparentSw_(const Params &params, const Evaluation &Sw)
Convert an absolute wetting saturation to an apparent one.
Definition ParkerLenhard.hpp:493
static Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition ParkerLenhard.hpp:438
static Evaluation pcnw(const Params &params, const FluidState &fs)
Returns the capillary pressure dependend on the phase saturations.
Definition ParkerLenhard.hpp:379
static const int numPhases
The number of fluid phases.
Definition ParkerLenhard.hpp:249
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition ParkerLenhard.hpp:264
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition ParkerLenhard.hpp:256
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
Returns the capillary pressure dependening on the phase saturations.
Definition ParkerLenhard.hpp:345
static void saturations(Container &, const Params &, const FluidState &)
Returns the capillary pressure dependening on the phase saturations.
Definition ParkerLenhard.hpp:358
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition ParkerLenhard.hpp:276
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation &Sw)
The saturation-capillary pressure curve according to van Genuchten using a material law specific API.
Definition VanGenuchten.hpp:194
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30