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_VAN_GENUCHTEN_HPP
00028 #define OPM_VAN_GENUCHTEN_HPP
00029
00030 #include "VanGenuchtenParams.hpp"
00031
00032 #include <opm/material/common/MathToolbox.hpp>
00033
00034 #include <algorithm>
00035 #include <cmath>
00036 #include <cassert>
00037
00038 namespace Opm {
00054 template <class TraitsT, class ParamsT = VanGenuchtenParams<TraitsT> >
00055 class VanGenuchten : public TraitsT
00056 {
00057 public:
00059 typedef TraitsT Traits;
00060
00062 typedef ParamsT Params;
00063
00065 typedef typename Traits::Scalar Scalar;
00066
00068 static const int numPhases = Traits::numPhases;
00069 static_assert(numPhases == 2,
00070 "The van Genuchten capillary pressure law only "
00071 "applies to the case of two fluid phases");
00072
00075 static const bool implementsTwoPhaseApi = true;
00076
00079 static const bool implementsTwoPhaseSatApi = true;
00080
00083 static const bool isSaturationDependent = true;
00084
00087 static const bool isPressureDependent = false;
00088
00091 static const bool isTemperatureDependent = false;
00092
00095 static const bool isCompositionDependent = false;
00096
00113 template <class Container, class FluidState>
00114 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
00115 {
00116 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00117
00118 values[Traits::wettingPhaseIdx] = 0.0;
00119 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
00120 }
00121
00126 template <class Container, class FluidState>
00127 static void saturations(Container& values, const Params& params, const FluidState& fs)
00128 {
00129 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00130
00131 values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
00132 values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
00133 }
00134
00145 template <class Container, class FluidState>
00146 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
00147 {
00148 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00149
00150 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
00151 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
00152 }
00153
00168 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00169 static Evaluation pcnw(const Params& params, const FluidState& fs)
00170 {
00171 const Evaluation& Sw =
00172 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00173
00174 assert(0 <= Sw && Sw <= 1);
00175
00176 return twoPhaseSatPcnw(params, Sw);
00177 }
00178
00193 template <class Evaluation>
00194 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00195 {
00196 return Opm::pow(Opm::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha();
00197 }
00198
00211 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00212 static Evaluation Sw(const Params& params, const FluidState& fs)
00213 {
00214 Evaluation pC =
00215 Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
00216 - Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
00217 return twoPhaseSatSw(params, pC);
00218 }
00219
00220 template <class Evaluation>
00221 static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pC)
00222 {
00223 assert(pC >= 0);
00224
00225 return Opm::pow(Opm::pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM());
00226 }
00227
00232 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00233 static Evaluation Sn(const Params& params, const FluidState& fs)
00234 { return 1 - Sw<FluidState, Evaluation>(params, fs); }
00235
00236 template <class Evaluation>
00237 static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pC)
00238 { return 1 - twoPhaseSatSw(params, pC); }
00239
00250 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00251 static Evaluation krw(const Params& params, const FluidState& fs)
00252 {
00253 const Evaluation& Sw =
00254 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00255
00256 return twoPhaseSatKrw(params, Sw);
00257 }
00258
00259 template <class Evaluation>
00260 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00261 {
00262 assert(0.0 <= Sw && Sw <= 1.0);
00263
00264 Evaluation r = 1.0 - Opm::pow(1.0 - Opm::pow(Sw, 1/params.vgM()), params.vgM());
00265 return Opm::sqrt(Sw)*r*r;
00266 }
00267
00277 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00278 static Evaluation krn(const Params& params, const FluidState& fs)
00279 {
00280 const Evaluation& Sw =
00281 1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00282
00283 return twoPhaseSatKrn(params, Sw);
00284 }
00285
00286 template <class Evaluation>
00287 static Evaluation twoPhaseSatKrn(const Params& params, Evaluation Sw)
00288 {
00289 assert(0 <= Sw && Sw <= 1);
00290
00291 return
00292 Opm::pow(1 - Sw, 1.0/3) *
00293 Opm::pow(1 - Opm::pow(Sw, 1/params.vgM()), 2*params.vgM());
00294 }
00295 };
00296 }
00297
00298 #endif