00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_PARKER_LENHARD_HPP
00028 #define OPM_PARKER_LENHARD_HPP
00029
00030 #include "ParkerLenhardParams.hpp"
00031
00032 #include <opm/material/fluidmatrixinteractions/VanGenuchten.hpp>
00033
00034 #include <algorithm>
00035 #include <cassert>
00036
00037 namespace Opm {
00038
00046 template <class ScalarT>
00047 class PLScanningCurve
00048 {
00049 public:
00050 typedef ScalarT Scalar;
00051
00058 PLScanningCurve(Scalar Swr)
00059 {
00060 loopNum_ = 0;
00061 prev_ = new PLScanningCurve(NULL,
00062 this,
00063 -1,
00064 Swr,
00065 1e12,
00066 Swr,
00067 Swr);
00068 next_ = NULL;
00069
00070 Sw_ = 1.0;
00071 pcnw_ = 0.0;
00072 SwMic_ = 1.0;
00073 SwMdc_ = 1.0;
00074 }
00075
00076 protected:
00077 PLScanningCurve(PLScanningCurve* prevSC,
00078 PLScanningCurve* nextSC,
00079 int loopN,
00080 Scalar SwReversal,
00081 Scalar pcnwReversal,
00082 Scalar SwMiCurve,
00083 Scalar SwMdCurve)
00084 {
00085 prev_ = prevSC;
00086 next_ = nextSC;
00087 loopNum_ = loopN;
00088 Sw_ = SwReversal;
00089 pcnw_ = pcnwReversal;
00090 SwMic_ = SwMiCurve;
00091 SwMdc_ = SwMdCurve;
00092 }
00093
00094 public:
00100 ~PLScanningCurve()
00101 {
00102 if (loopNum_ == 0)
00103 delete prev_;
00104 if (loopNum_ >= 0)
00105 delete next_;
00106 }
00107
00112 PLScanningCurve* prev() const
00113 { return prev_; }
00114
00119 PLScanningCurve* next() const
00120 { return next_; }
00121
00130 void setNext(Scalar SwReversal,
00131 Scalar pcnwReversal,
00132 Scalar SwMiCurve,
00133 Scalar SwMdCurve)
00134 {
00135
00136
00137 delete next_;
00138
00139 next_ = new PLScanningCurve(this,
00140 NULL,
00141 loopNum() + 1,
00142 SwReversal,
00143 pcnwReversal,
00144 SwMiCurve,
00145 SwMdCurve);
00146 }
00147
00154 bool isValidAt_Sw(Scalar SwReversal)
00155 {
00156 if (isImbib())
00157
00158
00159
00160
00161 return this->Sw() < SwReversal && SwReversal < prev_->Sw();
00162 else
00163
00164
00165
00166
00167 return prev_->Sw() < SwReversal && SwReversal < this->Sw();
00168 }
00169
00174 bool isImbib()
00175 { return loopNum()%2 == 1; }
00176
00181 bool isDrain()
00182 { return !isImbib(); }
00183
00189 int loopNum()
00190 { return loopNum_; }
00191
00196 Scalar Sw() const
00197 { return Sw_; }
00198
00202 Scalar pcnw() const
00203 { return pcnw_; }
00204
00209 Scalar SwMic()
00210 { return SwMic_; }
00211
00216 Scalar SwMdc()
00217 { return SwMdc_; }
00218
00219 private:
00220 PLScanningCurve* prev_;
00221 PLScanningCurve* next_;
00222
00223 int loopNum_;
00224
00225 Scalar Sw_;
00226 Scalar pcnw_;
00227
00228 Scalar SwMdc_;
00229 Scalar SwMic_;
00230 };
00231
00238 template <class TraitsT, class ParamsT = ParkerLenhardParams<TraitsT> >
00239 class ParkerLenhard : public TraitsT
00240 {
00241 public:
00242 typedef TraitsT Traits;
00243 typedef ParamsT Params;
00244 typedef typename Traits::Scalar Scalar;
00245
00247 static const int numPhases = Traits::numPhases;
00248 static_assert(numPhases == 2,
00249 "The Parker-Lenhard capillary pressure law only "
00250 "applies to the case of two fluid phases");
00251
00254 static const bool implementsTwoPhaseApi = true;
00255
00258 static const bool implementsTwoPhaseSatApi = true;
00259
00262 static const bool isSaturationDependent = true;
00263
00266 static const bool isPressureDependent = false;
00267
00270 static const bool isTemperatureDependent = false;
00271
00274 static const bool isCompositionDependent = false;
00275
00276 static_assert(Traits::numPhases == 2,
00277 "The number of fluid phases must be two if you want to use "
00278 "this material law!");
00279
00280 private:
00281 typedef typename ParamsT::VanGenuchten VanGenuchten;
00282 typedef Opm::PLScanningCurve<Scalar> ScanningCurve;
00283
00284 public:
00289 static void reset(Params& params)
00290 {
00291 delete params.mdc();
00292 params.setMdc(new ScanningCurve(params.SwrPc()));
00293 params.setCsc(params.mdc());
00294 params.setPisc(NULL);
00295 params.setCurrentSnr(0.0);
00296 }
00297
00302 template <class FluidState>
00303 static void update(Params& params, const FluidState& fs)
00304 {
00305 Scalar Sw = Opm::scalarValue(fs.saturation(Traits::wettingPhaseIdx));
00306
00307 if (Sw > 1 - 1e-5) {
00308
00309
00310 reset(params);
00311 return;
00312 }
00313
00314
00315
00316 ScanningCurve* curve = findScanningCurve_(params, Sw);
00317
00318
00319
00320
00321 Scalar pc = pcnw<FluidState, Scalar>(params, fs);
00322 Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
00323 Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
00324
00325 curve->setNext(Sw, pc, Sw_mic, Sw_mdc);
00326 if (!curve->next())
00327 return;
00328
00329 params.setCsc(curve);
00330
00331
00332 if (params.csc() == params.mdc()) {
00333 params.setPisc(params.mdc()->next());
00334 params.setCurrentSnr(computeCurrentSnr_(params, Sw));
00335 }
00336 }
00337
00342 template <class Container, class FluidState>
00343 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
00344 {
00345 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00346
00347 values[Traits::wettingPhaseIdx] = 0.0;
00348 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
00349 }
00350
00355 template <class Container, class FluidState>
00356 static void saturations(Container& , const Params& , const FluidState& )
00357 { OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::saturations()"); }
00358
00363 template <class Container, class FluidState>
00364 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
00365 {
00366 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00367
00368 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
00369 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
00370 }
00371
00376 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00377 static Evaluation pcnw(const Params& params, const FluidState& fs)
00378 {
00379 const Evaluation& Sw =
00380 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00381
00382 return twoPhaseSatPcnw(params, Sw);
00383 }
00384
00385 template <class Evaluation>
00386 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00387 {
00388
00389 ScanningCurve* sc = findScanningCurve_(params, Opm::scalarValue(Sw));
00390
00391
00392 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
00393
00394
00395
00396 if (Sw_app > 1) {
00397 return 0.0;
00398 }
00399
00400
00401
00402 Scalar SwAppCurSC = absoluteToApparentSw_(params, sc->Sw());
00403 Scalar SwAppPrevSC = absoluteToApparentSw_(params, sc->prev()->Sw());
00404 const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
00405 if (sc->isImbib()) {
00406 const Evaluation& SwMic =
00407 pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
00408
00409 return VanGenuchten::twoPhaseSatPcnw(params.micParams(), SwMic);
00410 }
00411 else {
00412 const Evaluation& SwMdc =
00413 pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
00414
00415 return VanGenuchten::twoPhaseSatPcnw(params.mdcParams(), SwMdc);
00416 }
00417 }
00418
00423 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00424 static Evaluation Sw(const Params& , const FluidState& )
00425 { OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::Sw()"); }
00426
00427 template <class Evaluation>
00428 static Evaluation twoPhaseSatSw(const Params& , const Evaluation& )
00429 { OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
00430
00435 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00436 static Evaluation Sn(const Params& params, const FluidState& fs)
00437 { return 1 - Sw<FluidState, Evaluation>(params, fs); }
00438
00439 template <class Evaluation>
00440 static Evaluation twoPhaseSatSn(const Params& , const Evaluation& )
00441 { OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
00442
00447 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00448 static Evaluation krw(const Params& params, const FluidState& fs)
00449 {
00450 const Evaluation& Sw =
00451 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00452
00453 return twoPhaseSatKrw(params, Sw);
00454 }
00455
00456 template <class Evaluation>
00457 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00458 {
00459
00460
00461 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
00462 return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
00463 }
00464
00469 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00470 static Evaluation krn(const Params& params, const FluidState& fs)
00471 {
00472 const Evaluation& Sw =
00473 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00474
00475 return twoPhaseSatKrn(params, Sw);
00476 }
00477
00478 template <class Evaluation>
00479 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
00480 {
00481
00482
00483 const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
00484 return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
00485 }
00486
00490 template <class Evaluation>
00491 static Evaluation absoluteToApparentSw_(const Params& params, const Evaluation& Sw)
00492 {
00493 return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params, Sw));
00494 }
00495
00496 private:
00506 template <class Evaluation>
00507 static Evaluation absoluteToEffectiveSw_(const Params& params, const Evaluation& Sw)
00508 { return (Sw - params.SwrPc())/(1 - params.SwrPc()); }
00509
00519 template <class Evaluation>
00520 static Evaluation effectiveToAbsoluteSw_(const Params& params, const Evaluation& Swe)
00521 { return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
00522
00523
00524
00525 template <class Evaluation>
00526 static Evaluation computeCurrentSnr_(const Params& params, const Evaluation& Sw)
00527 {
00528
00529 if (Sw > 1 - params.Snr())
00530 return 0.0;
00531 if (Sw < params.SwrPc())
00532 return params.Snr();
00533
00534 if (params.Snr() == 0.0)
00535 return 0.0;
00536
00537
00538 Scalar R = 1.0/params.Snr() - 1;
00539 const Evaluation& curSnr = (1 - Sw)/(1 + R*(1 - Sw));
00540
00541
00542
00543 assert(curSnr <= params.Snr());
00544
00545 return curSnr;
00546 }
00547
00548
00549
00550 template <class Evaluation>
00551 static Evaluation trappedEffectiveSn_(const Params& params, const Evaluation& Sw)
00552 {
00553 const Evaluation& Swe = absoluteToEffectiveSw_(params, Sw);
00554 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
00555
00556 Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
00557 return Snre*(Swe - SwePisc) / (1 - Snre - SwePisc);
00558 }
00559
00560
00561
00562 template <class Evaluation>
00563 static Evaluation effectiveToApparentSw_(const Params& params, const Evaluation& Swe)
00564 {
00565 if (params.pisc() == NULL ||
00566 Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
00567 {
00568
00569
00570
00571 return Swe;
00572 }
00573
00574
00575
00576
00577 return Swe + trappedEffectiveSn_(params, Swe);
00578 }
00579
00580
00581 template <class Evaluation>
00582 static Evaluation apparentToEffectiveSw_(const Params& params, const Evaluation& Swapp)
00583 {
00584 Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
00585 if (params.pisc() == NULL || Swapp <= SwePisc) {
00586
00587
00588
00589 return Swapp;
00590 }
00591
00592 Scalar Snre = absoluteToEffectiveSw_(params.currentSnr());
00593 return
00594 (Swapp*(1 - Snre - SwePisc) + Snre*SwePisc)
00595 /(1 - SwePisc);
00596 }
00597
00598
00599
00600
00601 static ScanningCurve* findScanningCurve_(const Params& params, Scalar Sw)
00602 {
00603 if (params.pisc() == NULL || Sw <= params.pisc()->Sw()) {
00604
00605
00606
00607 return params.mdc();
00608 }
00609
00610
00611
00612
00613
00614 if (params.pisc()->next() == NULL ||
00615 params.pisc()->next()->Sw() < Sw)
00616 {
00617 return params.pisc();
00618 }
00619
00620 ScanningCurve* curve = params.csc()->next();
00621 while (true) {
00622 assert(curve != params.mdc()->prev());
00623 if (curve->isValidAt_Sw(Sw)) {
00624 return curve;
00625 }
00626 curve = curve->prev();
00627 }
00628 }
00629 };
00630
00631 }
00632
00633 #endif // PARKER_LENHARD_HPP