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_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP
00028 #define OPM_THREE_PHASE_PARKER_VAN_GENUCHTEN_HPP
00029
00030 #include "ThreePhaseParkerVanGenuchtenParams.hpp"
00031
00032 #include <opm/material/common/MathToolbox.hpp>
00033
00034 #include <algorithm>
00035
00036 namespace Opm {
00048 template <class TraitsT,
00049 class ParamsT = ThreePhaseParkerVanGenuchtenParams<TraitsT> >
00050 class ThreePhaseParkerVanGenuchten
00051 {
00052 public:
00053 static_assert(TraitsT::numPhases == 3,
00054 "The number of phases considered by this capillary pressure "
00055 "law is always three!");
00056
00057 typedef TraitsT Traits;
00058 typedef ParamsT Params;
00059 typedef typename Traits::Scalar Scalar;
00060
00061 static const int numPhases = 3;
00062 static const int wettingPhaseIdx = Traits::wettingPhaseIdx;
00063 static const int nonWettingPhaseIdx = Traits::nonWettingPhaseIdx;
00064 static const int gasPhaseIdx = Traits::gasPhaseIdx;
00065
00068 static const bool implementsTwoPhaseApi = false;
00069
00072 static const bool implementsTwoPhaseSatApi = false;
00073
00076 static const bool isSaturationDependent = true;
00077
00080 static const bool isPressureDependent = false;
00081
00084 static const bool isTemperatureDependent = false;
00085
00088 static const bool isCompositionDependent = false;
00089
00101 template <class ContainerT, class FluidState>
00102 static void capillaryPressures(ContainerT& values,
00103 const Params& params,
00104 const FluidState& fluidState)
00105 {
00106 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00107
00108 values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, fluidState);
00109 values[nonWettingPhaseIdx] = 0;
00110 values[wettingPhaseIdx] = - pcnw<FluidState, Evaluation>(params, fluidState);
00111 }
00112
00122 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00123 static Evaluation pcgn(const Params& params, const FluidState& fluidState)
00124 {
00125 Scalar PC_VG_REG = 0.01;
00126
00127
00128 const auto& St =
00129 Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx))
00130 + Opm::decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
00131
00132 Evaluation Se = (St - params.Swrx())/(1. - params.Swrx());
00133
00134
00135 if (Se < 0.0)
00136 Se=0.0;
00137 if (Se > 1.0)
00138 Se=1.0;
00139
00140 if (Se>PC_VG_REG && Se<1-PC_VG_REG)
00141 {
00142 const Evaluation& x = Opm::pow(Se,-1/params.vgM()) - 1;
00143 return Opm::pow(x, 1.0 - params.vgM())/params.vgAlpha();
00144 }
00145
00146
00147 Scalar Se_regu;
00148 if (Se<=PC_VG_REG)
00149 Se_regu = PC_VG_REG;
00150 else
00151 Se_regu = 1-PC_VG_REG;
00152 const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
00153 const Evaluation& pc = Opm::pow(x, 1.0/params.vgN())/params.vgAlpha();
00154 const Evaluation& pc_prime =
00155 Opm::pow(x, 1/params.vgN()-1)
00156 * std::pow(Se_regu,-1/params.vgM()-1)
00157 / (-params.vgM())
00158 / params.vgAlpha()
00159 / (1 - params.Sgr() - params.Swrx())
00160 / params.vgN();
00161
00162
00163 return ((Se-Se_regu)*pc_prime + pc)/params.betaGN();
00164 }
00165
00175 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00176 static Evaluation pcnw(const Params& params, const FluidState& fluidState)
00177 {
00178 const Evaluation& Sw =
00179 Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
00180 Evaluation Se = (Sw-params.Swr())/(1.-params.Snr());
00181
00182 Scalar PC_VG_REG = 0.01;
00183
00184
00185 if (Se<0.0)
00186 Se=0.0;
00187 if (Se>1.0)
00188 Se=1.0;
00189
00190 if (Se>PC_VG_REG && Se<1-PC_VG_REG) {
00191 Evaluation x = Opm::pow(Se,-1/params.vgM()) - 1.0;
00192 x = Opm::pow(x, 1 - params.vgM());
00193 return x/params.vgAlpha();
00194 }
00195
00196
00197 Scalar Se_regu;
00198 if (Se<=PC_VG_REG)
00199 Se_regu = PC_VG_REG;
00200 else
00201 Se_regu = 1.0 - PC_VG_REG;
00202
00203 const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
00204 const Evaluation& pc = Opm::pow(x, 1/params.vgN())/params.vgAlpha();
00205 const Evaluation& pc_prime =
00206 Opm::pow(x,1/params.vgN()-1)
00207 * std::pow(Se_regu, -1.0/params.vgM() - 1)
00208 / (-params.vgM())
00209 / params.vgAlpha()
00210 / (1-params.Snr()-params.Swr())
00211 / params.vgN();
00212
00213
00214 return ((Se-Se_regu)*pc_prime + pc)/params.betaNW();
00215 }
00216
00221 template <class ContainerT, class FluidState>
00222 static void saturations(ContainerT& ,
00223 const Params& ,
00224 const FluidState& )
00225 { OPM_THROW(std::logic_error, "Not implemented: inverse capillary pressures"); }
00226
00230 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00231 static Evaluation Sg(const Params& , const FluidState& )
00232 { OPM_THROW(std::logic_error, "Not implemented: Sg()"); }
00233
00237 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00238 static Evaluation Sn(const Params& , const FluidState& )
00239 { OPM_THROW(std::logic_error, "Not implemented: Sn()"); }
00240
00244 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00245 static Evaluation Sw(const Params& , const FluidState& )
00246 { OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
00247
00251 template <class ContainerT, class FluidState>
00252 static void relativePermeabilities(ContainerT& values,
00253 const Params& params,
00254 const FluidState& fluidState)
00255 {
00256 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00257
00258 values[wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
00259 values[nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
00260 values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
00261 }
00262
00271 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00272 static Evaluation krw(const Params& params, const FluidState& fluidState)
00273 {
00274 const Evaluation& Sw =
00275 Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
00276
00277 const Evaluation& Se = (Sw - params.Swr()) / (1-params.Swr());
00278
00279
00280 if(Se > 1.0) return 1.;
00281 if(Se < 0.0) return 0.;
00282
00283 const Evaluation& r = 1. - Opm::pow(1 - Opm::pow(Se, 1/params.vgM()), params.vgM());
00284 return Opm::sqrt(Se)*r*r;
00285 }
00286
00299 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00300 static Evaluation krn(const Params& params, const FluidState& fluidState)
00301 {
00302 const Evaluation& Sn =
00303 Opm::decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
00304 const Evaluation& Sw =
00305 Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
00306 Evaluation Swe = Opm::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
00307 Evaluation Ste = Opm::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
00308
00309
00310 if(Swe <= 0.0) Swe = 0.;
00311 if(Ste <= 0.0) Ste = 0.;
00312 if(Ste - Swe <= 0.0) return 0.;
00313
00314 Evaluation krn_;
00315 krn_ = Opm::pow(1 - Opm::pow(Swe, 1/params.vgM()), params.vgM());
00316 krn_ -= Opm::pow(1 - Opm::pow(Ste, 1/params.vgM()), params.vgM());
00317 krn_ *= krn_;
00318
00319 if (params.krRegardsSnr())
00320 {
00321
00322
00323 const Evaluation& resIncluded =
00324 Opm::max(Opm::min(Sw - params.Snr() / (1-params.Swr()), 1.0), 0.0);
00325 krn_ *= Opm::sqrt(resIncluded );
00326 }
00327 else
00328 krn_ *= Opm::sqrt(Sn / (1 - params.Swr()));
00329
00330 return krn_;
00331 }
00332
00333
00344 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00345 static Evaluation krg(const Params& params, const FluidState& fluidState)
00346 {
00347 const Evaluation& Sg =
00348 Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
00349 const Evaluation& Se = Opm::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
00350
00351
00352 if(Se > 1.0)
00353 return 0.0;
00354 if(Se < 0.0)
00355 return 1.0;
00356
00357 Evaluation scaleFactor = 1.;
00358 if (Sg<=0.1) {
00359 scaleFactor = (Sg - params.Sgr())/(0.1 - params.Sgr());
00360 if (scaleFactor < 0.)
00361 return 0.0;
00362 }
00363
00364 return scaleFactor
00365 * Opm::pow(1 - Se, 1.0/3.)
00366 * Opm::pow(1 - Opm::pow(Se, 1/params.vgM()), 2*params.vgM());
00367 }
00368 };
00369 }
00370
00371 #endif