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_BROOKS_COREY_HPP
00028 #define OPM_BROOKS_COREY_HPP
00029
00030 #include "BrooksCoreyParams.hpp"
00031
00032 #include <opm/material/common/MathToolbox.hpp>
00033 #include <opm/common/ErrorMacros.hpp>
00034 #include <opm/common/Exceptions.hpp>
00035
00036 #include <algorithm>
00037 #include <cassert>
00038 #include <cmath>
00039
00040 namespace Opm {
00053 template <class TraitsT, class ParamsT = BrooksCoreyParams<TraitsT> >
00054 class BrooksCorey : public TraitsT
00055 {
00056 public:
00057 typedef TraitsT Traits;
00058 typedef ParamsT Params;
00059 typedef typename Traits::Scalar Scalar;
00060
00062 static const int numPhases = Traits::numPhases;
00063 static_assert(numPhases == 2,
00064 "The Brooks-Corey capillary pressure law only applies "
00065 "to the case of two fluid phases");
00066
00069 static const bool implementsTwoPhaseApi = true;
00070
00073 static const bool implementsTwoPhaseSatApi = true;
00074
00077 static const bool isSaturationDependent = true;
00078
00081 static const bool isPressureDependent = false;
00082
00085 static const bool isTemperatureDependent = false;
00086
00089 static const bool isCompositionDependent = false;
00090
00091 static_assert(Traits::numPhases == 2,
00092 "The number of fluid phases must be two if you want to use "
00093 "this material law!");
00094
00098 template <class Container, class FluidState>
00099 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
00100 {
00101 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00102
00103 values[Traits::wettingPhaseIdx] = 0.0;
00104 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
00105 }
00106
00111 template <class Container, class FluidState>
00112 static void saturations(Container& values, const Params& params, const FluidState& fs)
00113 {
00114 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00115
00116 values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
00117 values[Traits::nonWettingPhaseIdx] = 1.0 - values[Traits::wettingPhaseIdx];
00118 }
00119
00130 template <class Container, class FluidState>
00131 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
00132 {
00133 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00134
00135 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
00136 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
00137 }
00138
00152 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00153 static Evaluation pcnw(const Params& params, const FluidState& fs)
00154 {
00155 const Evaluation& Sw =
00156 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00157
00158 assert(0.0 <= Sw && Sw <= 1.0);
00159
00160 return twoPhaseSatPcnw(params, Sw);
00161 }
00162
00163 template <class Evaluation>
00164 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00165 {
00166 assert(0.0 <= Sw && Sw <= 1.0);
00167
00168 return params.entryPressure()*Opm::pow(Sw, -1/params.lambda());
00169 }
00170
00171 template <class Evaluation>
00172 static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
00173 {
00174 assert(pcnw > 0.0);
00175
00176 return Opm::pow(params.entryPressure()/pcnw, -params.lambda());
00177 }
00178
00191 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00192 static Evaluation Sw(const Params& params, const FluidState& fs)
00193 {
00194 Evaluation pC =
00195 Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
00196 - Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
00197 return twoPhaseSatSw(params, pC);
00198 }
00199
00200 template <class Evaluation>
00201 static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pc)
00202 {
00203 assert(pc > 0.0);
00204
00205 return Opm::pow(pc/params.entryPressure(), -params.lambda());
00206 }
00207
00212 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00213 static Evaluation Sn(const Params& params, const FluidState& fs)
00214 { return 1.0 - Sw<FluidState, Evaluation>(params, fs); }
00215
00216 template <class Evaluation>
00217 static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pc)
00218 { return 1.0 - twoPhaseSatSw(params, pc); }
00227 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00228 static Evaluation krw(const Params& params, const FluidState& fs)
00229 {
00230 const auto& Sw =
00231 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00232
00233 return twoPhaseSatKrw(params, Sw);
00234 }
00235
00236 template <class Evaluation>
00237 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00238 {
00239 assert(0.0 <= Sw && Sw <= 1.0);
00240
00241 return Opm::pow(Sw, 2.0/params.lambda() + 3.0);
00242 }
00243
00244 template <class Evaluation>
00245 static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
00246 {
00247 return Opm::pow(krw, 1.0/(2.0/params.lambda() + 3.0));
00248 }
00249
00258 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00259 static Evaluation krn(const Params& params, const FluidState& fs)
00260 {
00261 const Evaluation& Sw =
00262 1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00263
00264 return twoPhaseSatKrn(params, Sw);
00265 }
00266
00267 template <class Evaluation>
00268 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
00269 {
00270 assert(0.0 <= Sw && Sw <= 1.0);
00271
00272 Scalar exponent = 2.0/params.lambda() + 1.0;
00273 const Evaluation Sn = 1.0 - Sw;
00274 return Sn*Sn*(1. - Opm::pow(Sw, exponent));
00275 }
00276
00277 template <class Evaluation>
00278 static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
00279 {
00280
00281
00282 Evaluation Sw = 0.5;
00283 Scalar eps = 1e-10;
00284 for (int i = 0; i < 20; ++i) {
00285 Evaluation f = krn - twoPhaseSatKrn(params, Sw);
00286 Evaluation fStar = krn - twoPhaseSatKrn(params, Sw + eps);
00287 Evaluation fPrime = (fStar - f)/eps;
00288
00289 Evaluation delta = f/fPrime;
00290 Sw -= delta;
00291 if (Sw < 0)
00292 Sw = 0.0;
00293 if (Opm::abs(delta) < 1e-10)
00294 return Sw;
00295 }
00296
00297 OPM_THROW(NumericalProblem,
00298 "Couldn't invert the Brooks-Corey non-wetting phase"
00299 " relperm within 20 iterations");
00300 }
00301 };
00302 }
00303
00304 #endif // BROOKS_COREY_HPP