My Project
EclStone1Material.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_ECL_STONE1_MATERIAL_HPP
28 #define OPM_ECL_STONE1_MATERIAL_HPP
29 
31 
34 
35 #include <algorithm>
36 #include <cmath>
37 #include <stdexcept>
38 #include <type_traits>
39 
40 namespace Opm {
41 
55 template <class TraitsT,
56  class GasOilMaterialLawT,
57  class OilWaterMaterialLawT,
58  class ParamsT = EclStone1MaterialParams<TraitsT, GasOilMaterialLawT, OilWaterMaterialLawT> >
59 class EclStone1Material : public TraitsT
60 {
61 public:
62  using GasOilMaterialLaw = GasOilMaterialLawT;
63  using OilWaterMaterialLaw = OilWaterMaterialLawT;
64 
65  // some safety checks
66  static_assert(TraitsT::numPhases == 3,
67  "The number of phases considered by this capillary pressure "
68  "law is always three!");
69  static_assert(GasOilMaterialLaw::numPhases == 2,
70  "The number of phases considered by the gas-oil capillary "
71  "pressure law must be two!");
72  static_assert(OilWaterMaterialLaw::numPhases == 2,
73  "The number of phases considered by the oil-water capillary "
74  "pressure law must be two!");
75  static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
76  typename OilWaterMaterialLaw::Scalar>::value,
77  "The two two-phase capillary pressure laws must use the same "
78  "type of floating point values.");
79 
80  static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
81  "The gas-oil material law must implement the two-phase saturation "
82  "only API to for the default Ecl capillary pressure law!");
83  static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
84  "The oil-water material law must implement the two-phase saturation "
85  "only API to for the default Ecl capillary pressure law!");
86 
87  using Traits = TraitsT;
88  using Params = ParamsT;
89  using Scalar = typename Traits::Scalar;
90 
91  static constexpr int numPhases = 3;
92  static constexpr int waterPhaseIdx = Traits::wettingPhaseIdx;
93  static constexpr int oilPhaseIdx = Traits::nonWettingPhaseIdx;
94  static constexpr int gasPhaseIdx = Traits::gasPhaseIdx;
95 
98  static constexpr bool implementsTwoPhaseApi = false;
99 
102  static constexpr bool implementsTwoPhaseSatApi = false;
103 
106  static constexpr bool isSaturationDependent = true;
107 
110  static constexpr bool isPressureDependent = false;
111 
114  static constexpr bool isTemperatureDependent = false;
115 
118  static constexpr bool isCompositionDependent = false;
119 
134  template <class ContainerT, class FluidState>
135  static void capillaryPressures(ContainerT& values,
136  const Params& params,
137  const FluidState& state)
138  {
139  using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
140  values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
141  values[oilPhaseIdx] = 0;
142  values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
143  Valgrind::CheckDefined(values[gasPhaseIdx]);
144  Valgrind::CheckDefined(values[oilPhaseIdx]);
145  Valgrind::CheckDefined(values[waterPhaseIdx]);
146  }
147 
148  /*
149  * Hysteresis parameters for oil-water
150  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
151  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
152  * \param params Parameters
153  */
154  static void oilWaterHysteresisParams(Scalar& pcSwMdc,
155  Scalar& krnSwMdc,
156  const Params& params)
157  {
158  pcSwMdc = params.oilWaterParams().pcSwMdc();
159  krnSwMdc = params.oilWaterParams().krnSwMdc();
160 
161  Valgrind::CheckDefined(pcSwMdc);
162  Valgrind::CheckDefined(krnSwMdc);
163  }
164 
165  /*
166  * Hysteresis parameters for oil-water
167  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
168  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
169  * \param params Parameters
170  */
171  static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
172  const Scalar& krnSwMdc,
173  Params& params)
174  {
175  constexpr const double krwSw = 2.0; //Should not be used
176  params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
177  }
178 
179  /*
180  * Hysteresis parameters for gas-oil
181  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
182  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
183  * \param params Parameters
184  */
185  static void gasOilHysteresisParams(Scalar& pcSwMdc,
186  Scalar& krnSwMdc,
187  const Params& params)
188  {
189  const auto Swco = params.Swl();
190 
191  // Pretend oil saturation is 'Swco' larger than it really is in
192  // order to infer correct maximum Sg values in output layer.
193  pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
194  krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
195 
196  Valgrind::CheckDefined(pcSwMdc);
197  Valgrind::CheckDefined(krnSwMdc);
198  }
199 
200  /*
201  * Hysteresis parameters for gas-oil
202  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
203  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
204  * \param params Parameters
205  */
206  static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
207  const Scalar& krnSwMdc,
208  Params& params)
209  {
210  // Maximum attainable oil saturation is 1-SWL
211  const auto Swco = params.Swl();
212  constexpr const double krwSw = 2.0; //Should not be used
213  params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
214  }
215 
225  template <class FluidState, class Evaluation = typename FluidState::Scalar>
226  static Evaluation pcgn(const Params& params,
227  const FluidState& fs)
228  {
229  // Maximum attainable oil saturation is 1-SWL
230  const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
231  return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
232  }
233 
243  template <class FluidState, class Evaluation = typename FluidState::Scalar>
244  static Evaluation pcnw(const Params& params,
245  const FluidState& fs)
246  {
247  const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
248  Valgrind::CheckDefined(Sw);
249 
250  const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
251  Valgrind::CheckDefined(result);
252 
253  return result;
254  }
255 
259  template <class ContainerT, class FluidState>
260  static void saturations(ContainerT& /* values */,
261  const Params& /* params */,
262  const FluidState& /* fluidState */)
263  {
264  throw std::logic_error("Not implemented: saturations()");
265  }
266 
270  template <class FluidState, class Evaluation = typename FluidState::Scalar>
271  static Evaluation Sg(const Params& /* params */,
272  const FluidState& /* fluidState */)
273  {
274  throw std::logic_error("Not implemented: Sg()");
275  }
276 
280  template <class FluidState, class Evaluation = typename FluidState::Scalar>
281  static Evaluation Sn(const Params& /* params */,
282  const FluidState& /* fluidState */)
283  {
284  throw std::logic_error("Not implemented: Sn()");
285  }
286 
290  template <class FluidState, class Evaluation = typename FluidState::Scalar>
291  static Evaluation Sw(const Params& /* params */,
292  const FluidState& /* fluidState */)
293  {
294  throw std::logic_error("Not implemented: Sw()");
295  }
296 
312  template <class ContainerT, class FluidState>
313  static void relativePermeabilities(ContainerT& values,
314  const Params& params,
315  const FluidState& fluidState)
316  {
317  using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
318 
319  values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
320  values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
321  values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
322  }
323 
327  template <class FluidState, class Evaluation = typename FluidState::Scalar>
328  static Evaluation krg(const Params& params,
329  const FluidState& fluidState)
330  {
331  // Maximum attainable oil saturation is 1-SWL,
332  const Evaluation Sw = 1 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
333  return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
334  }
335 
339  template <class FluidState, class Evaluation = typename FluidState::Scalar>
340  static Evaluation krw(const Params& params,
341  const FluidState& fluidState)
342  {
343  const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
344  return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
345  }
346 
350  template <class FluidState, class Evaluation = typename FluidState::Scalar>
351  static Evaluation krn(const Params& params,
352  const FluidState& fluidState)
353  {
354  // the Eclipse docu is inconsistent with naming the variable of connate water: In
355  // some places the connate water saturation is represented by "Swl", in others
356  // "Swco" is used.
357  const Scalar Swco = params.Swl();
358 
359  // oil relperm at connate water saturations (with Sg=0)
360  const Scalar krocw = params.krocw();
361 
362  const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
363  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
364 
365  const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
366  const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
367 
368  Evaluation beta;
369  if (Sw <= Swco)
370  beta = 1.0;
371  else {
372  // there seems to be an error in the ECL documentation: using the approach to
373  // the scaled saturations as described there leads to significant deviations
374  // from the results produced by Eclipse 100.
375  const Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
376  const Evaluation SSg = Sg/(1.0 - Swco);
377  const Evaluation SSo = 1.0 - SSw - SSg;
378 
379  if (SSw >= 1.0 || SSg >= 1.0)
380  beta = 1.0;
381  else
382  beta = pow( SSo/((1 - SSw)*(1 - SSg)), params.eta());
383  }
384 
385  return max(0.0, min(1.0, beta*kro_ow*kro_go/krocw));
386  }
387 
391  template <class Evaluation, class FluidState>
392  static Evaluation relpermOilInOilGasSystem(const Params& params,
393  const FluidState& fluidState)
394  {
395  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
396 
397  return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - params.Swl());
398  }
399 
403  template <class Evaluation, class FluidState>
404  static Evaluation relpermOilInOilWaterSystem(const Params& params,
405  const FluidState& fluidState)
406  {
407  const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
408 
409  return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
410  }
411 
419  template <class FluidState>
420  static void updateHysteresis(Params& params, const FluidState& fluidState)
421  {
422  const Scalar Swco = params.Swl();
423  const Scalar Sw = scalarValue(fluidState.saturation(waterPhaseIdx));
424  const Scalar Sg = scalarValue(fluidState.saturation(gasPhaseIdx));
425 
426  params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
427  params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
428  /*krwSw=*/ 1.0 - Swco - Sg,
429  /*krnSw=*/ 1.0 - Swco - Sg);
430  }
431 };
432 
433 } // namespace Opm
434 
435 #endif
Default implementation for the parameters required by the three-phase capillary pressure/relperm Ston...
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Some templates to wrap the valgrind client request macros.
Implements the second phase capillary pressure/relperm law suggested by Stone as used by the ECLipse ...
Definition: EclStone1Material.hpp:60
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: EclStone1Material.hpp:106
static constexpr bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: EclStone1Material.hpp:102
static constexpr bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: EclStone1Material.hpp:118
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability of the wetting phase.
Definition: EclStone1Material.hpp:340
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition: EclStone1Material.hpp:351
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition: EclStone1Material.hpp:313
static void capillaryPressures(ContainerT &values, const Params &params, const FluidState &state)
Implements the default three phase capillary pressure law used by the ECLipse simulator.
Definition: EclStone1Material.hpp:135
static void saturations(ContainerT &, const Params &, const FluidState &)
The inverse of the capillary pressure.
Definition: EclStone1Material.hpp:260
static constexpr bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: EclStone1Material.hpp:114
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition: EclStone1Material.hpp:271
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition: EclStone1Material.hpp:328
static constexpr bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: EclStone1Material.hpp:98
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition: EclStone1Material.hpp:291
static Evaluation pcgn(const Params &params, const FluidState &fs)
Capillary pressure between the gas and the non-wetting liquid (i.e., oil) phase.
Definition: EclStone1Material.hpp:226
static constexpr bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: EclStone1Material.hpp:110
static Evaluation relpermOilInOilWaterSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/water system.
Definition: EclStone1Material.hpp:404
static Evaluation pcnw(const Params &params, const FluidState &fs)
Capillary pressure between the non-wetting liquid (i.e., oil) and the wetting liquid (i....
Definition: EclStone1Material.hpp:244
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/gas system.
Definition: EclStone1Material.hpp:392
static void updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition: EclStone1Material.hpp:420
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition: EclStone1Material.hpp:281