00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef REGULARIZED_BROOKS_COREY_HPP
00028 #define REGULARIZED_BROOKS_COREY_HPP
00029
00030 #include "BrooksCorey.hpp"
00031 #include "RegularizedBrooksCoreyParams.hpp"
00032
00033 #include <opm/material/common/Spline.hpp>
00034
00035 namespace Opm {
00064 template <class TraitsT, class ParamsT = RegularizedBrooksCoreyParams<TraitsT> >
00065 class RegularizedBrooksCorey : public TraitsT
00066 {
00067 typedef Opm::BrooksCorey<TraitsT, ParamsT> BrooksCorey;
00068
00069 public:
00070 typedef TraitsT Traits;
00071 typedef ParamsT Params;
00072 typedef typename Traits::Scalar Scalar;
00073
00075 static const int numPhases = Traits::numPhases;
00076 static_assert(numPhases == 2,
00077 "The regularized Brooks-Corey capillary pressure law only "
00078 "applies to the case of two fluid phases");
00079
00082 static const bool implementsTwoPhaseApi = true;
00083
00086 static const bool implementsTwoPhaseSatApi = true;
00087
00090 static const bool isSaturationDependent = true;
00091
00094 static const bool isPressureDependent = false;
00095
00098 static const bool isTemperatureDependent = false;
00099
00102 static const bool isCompositionDependent = false;
00103
00104 static_assert(Traits::numPhases == 2,
00105 "The number of fluid phases must be two if you want to use "
00106 "this material law!");
00107
00118 template <class Container, class FluidState>
00119 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
00120 {
00121 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00122
00123 values[Traits::wettingPhaseIdx] = 0.0;
00124 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
00125 }
00126
00131 template <class Container, class FluidState>
00132 static void saturations(Container& values, const Params& params, const FluidState& fs)
00133 {
00134 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00135
00136 values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
00137 values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
00138 }
00139
00150 template <class Container, class FluidState>
00151 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
00152 {
00153 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00154
00155 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
00156 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
00157 }
00158
00173 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00174 static Evaluation pcnw(const Params& params, const FluidState& fs)
00175 {
00176 const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00177 return twoPhaseSatPcnw(params, Sw);
00178 }
00179
00180 template <class Evaluation>
00181 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00182 {
00183 const Scalar Sthres = params.pcnwLowSw();
00184
00185 if (Sw <= Sthres) {
00186 Scalar m = params.pcnwSlopeLow();
00187 Scalar pcnw_SwLow = params.pcnwLow();
00188 return pcnw_SwLow + m*(Sw - Sthres);
00189 }
00190 else if (Sw >= 1.0) {
00191 Scalar m = params.pcnwSlopeHigh();
00192 Scalar pcnw_SwHigh = params.pcnwHigh();
00193 return pcnw_SwHigh + m*(Sw - 1.0);
00194 }
00195
00196
00197
00198 return BrooksCorey::twoPhaseSatPcnw(params, Sw);
00199 }
00200
00207 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00208 static Evaluation Sw(const Params& params, const FluidState& fs)
00209 {
00210 const Evaluation& pC =
00211 Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
00212 - Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
00213 return twoPhaseSatSw(params, pC);
00214 }
00215
00216 template <class Evaluation>
00217 static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pcnw)
00218 {
00219 const Scalar Sthres = params.pcnwLowSw();
00220
00221
00222
00223
00224
00225
00226 Evaluation Sw = 1.5;
00227 if (pcnw >= params.entryPressure())
00228 Sw = BrooksCorey::twoPhaseSatSw(params, pcnw);
00229
00230
00231
00232
00233
00234
00235
00236 if (Sw <= Sthres) {
00237
00238 Scalar m = params.pcnwSlopeLow();
00239 Scalar pcnw_SwLow = params.pcnwLow();
00240 return Sthres + (pcnw - pcnw_SwLow)/m;
00241 }
00242 else if (Sw > 1.0) {
00243 Scalar m = params.pcnwSlopeHigh();
00244 Scalar pcnw_SwHigh = params.pcnwHigh();
00245 return 1.0 + (pcnw - pcnw_SwHigh)/m;;
00246 }
00247
00248 return BrooksCorey::twoPhaseSatSw(params, pcnw);
00249 }
00250
00255 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00256 static Evaluation Sn(const Params& params, const FluidState& fs)
00257 { return 1 - Sw<FluidState, Evaluation>(params, fs); }
00258
00259 template <class Evaluation>
00260 static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pcnw)
00261 { return 1 - twoPhaseSatSw(params, pcnw); }
00262
00277 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00278 static Evaluation krw(const Params& params, const FluidState& fs)
00279 {
00280 const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00281 return twoPhaseSatKrw(params, Sw);
00282 }
00283
00284 template <class Evaluation>
00285 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00286 {
00287 if (Sw <= 0.0)
00288 return 0.0;
00289 else if (Sw >= 1.0)
00290 return 1.0;
00291
00292 return BrooksCorey::twoPhaseSatKrw(params, Sw);
00293 }
00294
00295 template <class Evaluation>
00296 static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
00297 {
00298 if (krw <= 0.0)
00299 return 0.0;
00300 else if (krw >= 1.0)
00301 return 1.0;
00302
00303 return BrooksCorey::twoPhaseSatKrwInv(params, krw);
00304 }
00305
00320 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00321 static Evaluation krn(const Params& params, const FluidState& fs)
00322 {
00323 const Evaluation& Sw =
00324 1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00325 return twoPhaseSatKrn(params, Sw);
00326 }
00327
00328 template <class Evaluation>
00329 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
00330 {
00331 if (Sw >= 1.0)
00332 return 0.0;
00333 else if (Sw <= 0.0)
00334 return 1.0;
00335
00336 return BrooksCorey::twoPhaseSatKrn(params, Sw);
00337 }
00338
00339 template <class Evaluation>
00340 static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
00341 {
00342 if (krn <= 0.0)
00343 return 1.0;
00344 else if (krn >= 1.0)
00345 return 0.0;
00346
00347 return BrooksCorey::twoPhaseSatKrnInv(params, krn);
00348 }
00349 };
00350 }
00351
00352 #endif