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_REGULARIZED_VAN_GENUCHTEN_HPP
00028 #define OPM_REGULARIZED_VAN_GENUCHTEN_HPP
00029
00030 #include "VanGenuchten.hpp"
00031 #include "RegularizedVanGenuchtenParams.hpp"
00032
00033 #include <opm/material/common/Spline.hpp>
00034
00035 #include <algorithm>
00036
00037 namespace Opm {
00038
00070 template <class TraitsT, class ParamsT = RegularizedVanGenuchtenParams<TraitsT> >
00071 class RegularizedVanGenuchten : public TraitsT
00072 {
00073 typedef Opm::VanGenuchten<TraitsT, ParamsT> VanGenuchten;
00074
00075 public:
00076 typedef TraitsT Traits;
00077 typedef ParamsT Params;
00078
00079 typedef typename Traits::Scalar Scalar;
00080
00082 static const int numPhases = Traits::numPhases;
00083 static_assert(numPhases == 2,
00084 "The regularized van Genuchten capillary pressure law only "
00085 "applies to the case of two fluid phases");
00086
00089 static const bool implementsTwoPhaseApi = true;
00090
00093 static const bool implementsTwoPhaseSatApi = true;
00094
00097 static const bool isSaturationDependent = true;
00098
00101 static const bool isPressureDependent = false;
00102
00105 static const bool isTemperatureDependent = false;
00106
00109 static const bool isCompositionDependent = false;
00110
00115 template <class Container, class FluidState>
00116 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
00117 {
00118 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00119
00120 values[Traits::wettingPhaseIdx] = 0.0;
00121 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
00122 }
00123
00128 template <class Container, class FluidState>
00129 static void saturations(Container& values, const Params& params, const FluidState& fs)
00130 {
00131 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00132
00133 values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
00134 values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
00135 }
00136
00141 template <class Container, class FluidState>
00142 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
00143 {
00144 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00145
00146 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
00147 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
00148 }
00149
00162 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00163 static Evaluation pcnw(const Params& params, const FluidState& fs)
00164 {
00165 const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00166 return twoPhaseSatPcnw(params, Sw);
00167 }
00168
00169 template <class Evaluation>
00170 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00171 {
00172
00173
00174 const Scalar SwThLow = params.pcnwLowSw();
00175 const Scalar SwThHigh = params.pcnwHighSw();
00176
00177
00178
00179
00180
00181
00182 if (Sw < SwThLow) {
00183 return params.pcnwLow() + params.pcnwSlopeLow()*(Sw - SwThLow);
00184 }
00185 else if (Sw > SwThHigh)
00186 {
00187 Scalar yTh = params.pcnwHigh();
00188 Scalar m1 = (0.0 - yTh)/(1.0 - SwThHigh)*2;
00189
00190 if (Sw < 1.0) {
00191
00192 const Spline<Scalar>& sp = params.pcnwHighSpline();
00193
00194 return sp.eval(Sw);
00195 }
00196 else {
00197
00198 return m1*(Sw - 1.0) + 0.0;
00199 }
00200 }
00201
00202
00203
00204 return VanGenuchten::twoPhaseSatPcnw(params, Sw);
00205 }
00206
00221 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00222 static Evaluation Sw(const Params& params, const FluidState& fs)
00223 {
00224 const Evaluation& pC =
00225 Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
00226 - Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
00227 return twoPhaseSatSw(params, pC);
00228 }
00229
00230 template <class Evaluation>
00231 static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pC)
00232 {
00233
00234
00235 const Scalar SwThLow = params.pcnwLowSw();
00236 const Scalar SwThHigh = params.pcnwHighSw();
00237
00238
00239
00240
00241 if (pC <= 0) {
00242
00243 Scalar m1 = params.pcnwSlopeHigh();
00244 return pC/m1 + 1.0;
00245 }
00246
00247 Evaluation Sw = VanGenuchten::twoPhaseSatSw(params, pC);
00248
00249
00250 if (Sw <= SwThLow) {
00251
00252 Scalar pC_SwLow = VanGenuchten::twoPhaseSatPcnw(params, SwThLow);
00253 return (pC - pC_SwLow)/params.pcnwSlopeLow() + SwThLow;
00254 }
00255 else if (SwThHigh < Sw )
00256 {
00257
00258 const Spline<Scalar>& spline = params.pcnwHighSpline();
00259
00260 return spline.template intersectInterval<Evaluation>(SwThHigh, 1.0,
00261 0, 0, 0, pC);
00262 }
00263
00264 return Sw;
00265 }
00266
00271 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00272 static Evaluation Sn(const Params& params, const FluidState& fs)
00273 { return 1 - Sw<FluidState, Evaluation>(params, fs); }
00274
00275 template <class Evaluation>
00276 static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pc)
00277 { return 1 - twoPhaseSatSw(params, pc); }
00278
00293 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00294 static Evaluation krw(const Params& params, const FluidState& fs)
00295 {
00296 const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00297 return twoPhaseSatKrw(params, Sw);
00298 }
00299
00300 template <class Evaluation>
00301 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00302 {
00303
00304 if (Sw <= 0)
00305 return 0.0;
00306 else if (Sw >= 1)
00307 return 1.0;
00308
00309 return VanGenuchten::twoPhaseSatKrw(params, Sw);
00310 }
00311
00326 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00327 static Evaluation krn(const Params& params, const FluidState& fs)
00328 {
00329 const Evaluation& Sw =
00330 1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00331 return twoPhaseSatKrn(params, Sw);
00332 }
00333
00334 template <class Evaluation>
00335 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
00336 {
00337
00338 if (Sw <= 0)
00339 return 1.0;
00340 else if (Sw >= 1)
00341 return 0.0;
00342
00343 return VanGenuchten::twoPhaseSatKrn(params, Sw);
00344 }
00345 };
00346
00347 }
00348
00349 #endif