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_DEFAULT_MATERIAL_HPP
00028 #define OPM_ECL_DEFAULT_MATERIAL_HPP
00029
00030 #include "EclDefaultMaterialParams.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
00040 namespace Opm {
00041
00055 template <class TraitsT,
00056 class GasOilMaterialLawT,
00057 class OilWaterMaterialLawT,
00058 class ParamsT = EclDefaultMaterialParams<TraitsT,
00059 typename GasOilMaterialLawT::Params,
00060 typename OilWaterMaterialLawT::Params> >
00061 class EclDefaultMaterial : public TraitsT
00062 {
00063 public:
00064 typedef GasOilMaterialLawT GasOilMaterialLaw;
00065 typedef OilWaterMaterialLawT OilWaterMaterialLaw;
00066
00067
00068 static_assert(TraitsT::numPhases == 3,
00069 "The number of phases considered by this capillary pressure "
00070 "law is always three!");
00071 static_assert(GasOilMaterialLaw::numPhases == 2,
00072 "The number of phases considered by the gas-oil capillary "
00073 "pressure law must be two!");
00074 static_assert(OilWaterMaterialLaw::numPhases == 2,
00075 "The number of phases considered by the oil-water capillary "
00076 "pressure law must be two!");
00077 static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
00078 typename OilWaterMaterialLaw::Scalar>::value,
00079 "The two two-phase capillary pressure laws must use the same "
00080 "type of floating point values.");
00081
00082 static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
00083 "The gas-oil material law must implement the two-phase saturation "
00084 "only API to for the default Ecl capillary pressure law!");
00085 static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
00086 "The oil-water material law must implement the two-phase saturation "
00087 "only API to for the default Ecl capillary pressure law!");
00088
00089 typedef TraitsT Traits;
00090 typedef ParamsT Params;
00091 typedef typename Traits::Scalar Scalar;
00092
00093 static const int numPhases = 3;
00094 static const int waterPhaseIdx = Traits::wettingPhaseIdx;
00095 static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
00096 static const int gasPhaseIdx = Traits::gasPhaseIdx;
00097
00100 static const bool implementsTwoPhaseApi = false;
00101
00104 static const bool implementsTwoPhaseSatApi = false;
00105
00108 static const bool isSaturationDependent = true;
00109
00112 static const bool isPressureDependent = false;
00113
00116 static const bool isTemperatureDependent = false;
00117
00120 static const bool isCompositionDependent = false;
00121
00136 template <class ContainerT, class FluidState>
00137 static void capillaryPressures(ContainerT& values,
00138 const Params& params,
00139 const FluidState& state)
00140 {
00141 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00142 values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
00143 values[oilPhaseIdx] = 0;
00144 values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
00145
00146 Valgrind::CheckDefined(values[gasPhaseIdx]);
00147 Valgrind::CheckDefined(values[oilPhaseIdx]);
00148 Valgrind::CheckDefined(values[waterPhaseIdx]);
00149 }
00150
00151
00152
00153
00154
00155
00156
00157 static void oilWaterHysteresisParams(Scalar& pcSwMdc,
00158 Scalar& krnSwMdc,
00159 const Params& params)
00160 {
00161 pcSwMdc = params.oilWaterParams().pcSwMdc();
00162 krnSwMdc = params.oilWaterParams().krnSwMdc();
00163
00164 Valgrind::CheckDefined(pcSwMdc);
00165 Valgrind::CheckDefined(krnSwMdc);
00166 }
00167
00168
00169
00170
00171
00172
00173
00174 static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
00175 const Scalar& krnSwMdc,
00176 Params& params)
00177 {
00178 const double krwSw = 2.0;
00179 params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
00180 }
00181
00182
00183
00184
00185
00186
00187
00188 static void gasOilHysteresisParams(Scalar& pcSwMdc,
00189 Scalar& krnSwMdc,
00190 const Params& params)
00191 {
00192 pcSwMdc = params.gasOilParams().pcSwMdc();
00193 krnSwMdc = params.gasOilParams().krnSwMdc();
00194
00195 Valgrind::CheckDefined(pcSwMdc);
00196 Valgrind::CheckDefined(krnSwMdc);
00197 }
00198
00199
00200
00201
00202
00203
00204
00205 static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
00206 const Scalar& krnSwMdc,
00207 Params& params)
00208 {
00209 const double krwSw = 2.0;
00210 params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc);
00211 }
00212
00222 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00223 static Evaluation pcgn(const Params& params,
00224 const FluidState& fs)
00225 {
00226 const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
00227 return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
00228 }
00229
00239 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00240 static Evaluation pcnw(const Params& params,
00241 const FluidState& fs)
00242 {
00243 const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
00244 return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
00245 }
00246
00250 template <class ContainerT, class FluidState>
00251 static void saturations(ContainerT& ,
00252 const Params& ,
00253 const FluidState& )
00254 {
00255 OPM_THROW(std::logic_error, "Not implemented: saturations()");
00256 }
00257
00261 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00262 static Evaluation Sg(const Params& ,
00263 const FluidState& )
00264 {
00265 OPM_THROW(std::logic_error, "Not implemented: Sg()");
00266 }
00267
00271 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00272 static Evaluation Sn(const Params& ,
00273 const FluidState& )
00274 {
00275 OPM_THROW(std::logic_error, "Not implemented: Sn()");
00276 }
00277
00281 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00282 static Evaluation Sw(const Params& ,
00283 const FluidState& )
00284 {
00285 OPM_THROW(std::logic_error, "Not implemented: Sw()");
00286 }
00287
00303 template <class ContainerT, class FluidState>
00304 static void relativePermeabilities(ContainerT& values,
00305 const Params& params,
00306 const FluidState& fluidState)
00307 {
00308 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00309
00310 values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
00311 values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
00312 values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
00313 }
00314
00318 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00319 static Evaluation krg(const Params& params,
00320 const FluidState& fluidState)
00321 {
00322 const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
00323 return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
00324 }
00325
00329 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00330 static Evaluation krw(const Params& params,
00331 const FluidState& fluidState)
00332 {
00333 const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
00334 return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
00335 }
00336
00340 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00341 static Evaluation krn(const Params& params,
00342 const FluidState& fluidState)
00343 {
00344 Scalar Swco = params.Swl();
00345
00346 Evaluation Sw =
00347 Opm::max(Evaluation(Swco),
00348 Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
00349 Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
00350
00351 Evaluation Sw_ow = Sg + Sw;
00352 Evaluation So_go = 1.0 - Sw_ow;
00353 const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
00354 const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
00355
00356
00357
00358
00359 const Scalar epsilon = 1e-5;
00360 if (Opm::scalarValue(Sw_ow) - Swco < epsilon) {
00361 Evaluation kro2 = (kro_ow + kro_go)/2;;
00362 if (Opm::scalarValue(Sw_ow) - Swco > epsilon/2) {
00363 Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
00364 Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
00365 return kro2*alpha + kro1*(1 - alpha);
00366 }
00367
00368 return kro2;
00369 }
00370 else
00371 return (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
00372 }
00373
00381 template <class FluidState>
00382 static void updateHysteresis(Params& params, const FluidState& fluidState)
00383 {
00384 Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
00385 Scalar So = Opm::scalarValue(fluidState.saturation(oilPhaseIdx));
00386 Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
00387
00388 if (params.inconsistentHysteresisUpdate()) {
00389 Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399 params.oilWaterParams().update(1 - So, 1 - So, 1 - So);
00400 params.gasOilParams().update(1 - Sg, 1 - Sg, 1 - Sg);
00401 }
00402 else {
00403 Scalar Swco = params.Swl();
00404 Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw));
00405 Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
00406
00407 Scalar Sw_ow = Sg + std::max(Swco, Sw);
00408 Scalar So_go = 1 + Sw_ow;
00409
00410 params.oilWaterParams().update(Sw, 1 - Sg, Sw_ow);
00411 params.gasOilParams().update(1 - Sg, So_go, 1 - Sg);
00412 }
00413 }
00414 };
00415 }
00416
00417 #endif