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_ECL_STONE1_MATERIAL_HPP
00028 #define OPM_ECL_STONE1_MATERIAL_HPP
00029
00030 #include "EclStone1MaterialParams.hpp"
00031
00032 #include <opm/common/Valgrind.hpp>
00033 #include <opm/material/common/MathToolbox.hpp>
00034
00035 #include <opm/common/Exceptions.hpp>
00036 #include <opm/common/ErrorMacros.hpp>
00037
00038 #include <algorithm>
00039 #include <cmath>
00040
00041 namespace Opm {
00042
00056 template <class TraitsT,
00057 class GasOilMaterialLawT,
00058 class OilWaterMaterialLawT,
00059 class ParamsT = EclStone1MaterialParams<TraitsT, GasOilMaterialLawT, OilWaterMaterialLawT> >
00060 class EclStone1Material : public TraitsT
00061 {
00062 public:
00063 typedef GasOilMaterialLawT GasOilMaterialLaw;
00064 typedef OilWaterMaterialLawT OilWaterMaterialLaw;
00065
00066
00067 static_assert(TraitsT::numPhases == 3,
00068 "The number of phases considered by this capillary pressure "
00069 "law is always three!");
00070 static_assert(GasOilMaterialLaw::numPhases == 2,
00071 "The number of phases considered by the gas-oil capillary "
00072 "pressure law must be two!");
00073 static_assert(OilWaterMaterialLaw::numPhases == 2,
00074 "The number of phases considered by the oil-water capillary "
00075 "pressure law must be two!");
00076 static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
00077 typename OilWaterMaterialLaw::Scalar>::value,
00078 "The two two-phase capillary pressure laws must use the same "
00079 "type of floating point values.");
00080
00081 static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
00082 "The gas-oil material law must implement the two-phase saturation "
00083 "only API to for the default Ecl capillary pressure law!");
00084 static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
00085 "The oil-water material law must implement the two-phase saturation "
00086 "only API to for the default Ecl capillary pressure law!");
00087
00088 typedef TraitsT Traits;
00089 typedef ParamsT Params;
00090 typedef typename Traits::Scalar Scalar;
00091
00092 static const int numPhases = 3;
00093 static const int waterPhaseIdx = Traits::wettingPhaseIdx;
00094 static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
00095 static const int gasPhaseIdx = Traits::gasPhaseIdx;
00096
00099 static const bool implementsTwoPhaseApi = false;
00100
00103 static const bool implementsTwoPhaseSatApi = false;
00104
00107 static const bool isSaturationDependent = true;
00108
00111 static const bool isPressureDependent = false;
00112
00115 static const bool isTemperatureDependent = false;
00116
00119 static const bool isCompositionDependent = false;
00120
00135 template <class ContainerT, class FluidState>
00136 static void capillaryPressures(ContainerT& values,
00137 const Params& params,
00138 const FluidState& state)
00139 {
00140 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00141 values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
00142 values[oilPhaseIdx] = 0;
00143 values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
00144 Valgrind::CheckDefined(values[gasPhaseIdx]);
00145 Valgrind::CheckDefined(values[oilPhaseIdx]);
00146 Valgrind::CheckDefined(values[waterPhaseIdx]);
00147 }
00148
00149
00150
00151
00152
00153
00154
00155 static void oilWaterHysteresisParams(Scalar& pcSwMdc,
00156 Scalar& krnSwMdc,
00157 const Params& params)
00158 {
00159 pcSwMdc = params.oilWaterParams().pcSwMdc();
00160 krnSwMdc = params.oilWaterParams().krnSwMdc();
00161
00162 Valgrind::CheckDefined(pcSwMdc);
00163 Valgrind::CheckDefined(krnSwMdc);
00164 }
00165
00166
00167
00168
00169
00170
00171
00172 static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
00173 const Scalar& krnSwMdc,
00174 Params& params)
00175 {
00176 const double krwSw = 2.0;
00177 params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
00178 }
00179
00180
00181
00182
00183
00184
00185
00186 static void gasOilHysteresisParams(Scalar& pcSwMdc,
00187 Scalar& krnSwMdc,
00188 const Params& params)
00189 {
00190 pcSwMdc = params.gasOilParams().pcSwMdc();
00191 krnSwMdc = params.gasOilParams().krnSwMdc();
00192
00193 Valgrind::CheckDefined(pcSwMdc);
00194 Valgrind::CheckDefined(krnSwMdc);
00195 }
00196
00197
00198
00199
00200
00201
00202
00203 static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
00204 const Scalar& krnSwMdc,
00205 Params& params)
00206 {
00207 const double krwSw = 2.0;
00208 params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc);
00209 }
00210
00220 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00221 static Evaluation pcgn(const Params& params,
00222 const FluidState& fs)
00223 {
00224 const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
00225 return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
00226 }
00227
00237 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00238 static Evaluation pcnw(const Params& params,
00239 const FluidState& fs)
00240 {
00241 const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
00242 Valgrind::CheckDefined(Sw);
00243 const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
00244 Valgrind::CheckDefined(result);
00245 return result;
00246 }
00247
00251 template <class ContainerT, class FluidState>
00252 static void saturations(ContainerT& ,
00253 const Params& ,
00254 const FluidState& )
00255 {
00256 OPM_THROW(std::logic_error, "Not implemented: saturations()");
00257 }
00258
00262 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00263 static Evaluation Sg(const Params& ,
00264 const FluidState& )
00265 {
00266 OPM_THROW(std::logic_error, "Not implemented: Sg()");
00267 }
00268
00272 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00273 static Evaluation Sn(const Params& ,
00274 const FluidState& )
00275 {
00276 OPM_THROW(std::logic_error, "Not implemented: Sn()");
00277 }
00278
00282 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00283 static Evaluation Sw(const Params& ,
00284 const FluidState& )
00285 {
00286 OPM_THROW(std::logic_error, "Not implemented: Sw()");
00287 }
00288
00304 template <class ContainerT, class FluidState>
00305 static void relativePermeabilities(ContainerT& values,
00306 const Params& params,
00307 const FluidState& fluidState)
00308 {
00309 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00310
00311 values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
00312 values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
00313 values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
00314 }
00315
00319 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00320 static Evaluation krg(const Params& params,
00321 const FluidState& fluidState)
00322 {
00323 const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
00324 return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
00325 }
00326
00330 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00331 static Evaluation krw(const Params& params,
00332 const FluidState& fluidState)
00333 {
00334 const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
00335 return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
00336 }
00337
00341 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00342 static Evaluation krn(const Params& params,
00343 const FluidState& fluidState)
00344 {
00345
00346
00347
00348 Scalar Swco = params.Swl();
00349
00350
00351 Scalar krocw = params.krocw();
00352
00353 const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
00354 const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
00355
00356 Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
00357 Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Swco);
00358
00359 Evaluation beta;
00360 if (Sw <= Swco)
00361 beta = 1.0;
00362 else {
00363
00364
00365
00366 Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
00367 Evaluation SSg = Sg/(1.0 - Swco);
00368 Evaluation SSo = 1.0 - SSw - SSg;
00369
00370 if (SSw >= 1.0 || SSg >= 1.0)
00371 beta = 1.0;
00372 else
00373 beta = Opm::pow( SSo/((1 - SSw)*(1 - SSg)), params.eta());
00374 }
00375
00376 return Opm::max(0.0, Opm::min(1.0, beta*kro_ow*kro_go/krocw));
00377 }
00378
00386 template <class FluidState>
00387 static void updateHysteresis(Params& params, const FluidState& fluidState)
00388 {
00389 Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
00390 Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
00391
00392 params.oilWaterParams().update(Sw, Sw, Sw);
00393 params.gasOilParams().update(1 - Sg, 1 - Sg, 1 - Sg);
00394 }
00395 };
00396 }
00397
00398 #endif