My Project
Loading...
Searching...
No Matches
EclDefaultMaterial.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_DEFAULT_MATERIAL_HPP
28#define OPM_ECL_DEFAULT_MATERIAL_HPP
29
30#include <opm/common/TimingMacros.hpp>
31
35
36#include <algorithm>
37#include <stdexcept>
38#include <type_traits>
39
40namespace Opm {
41
55template <class TraitsT,
56 class GasOilMaterialLawT,
57 class OilWaterMaterialLawT,
58 class ParamsT = EclDefaultMaterialParams<TraitsT,
59 typename GasOilMaterialLawT::Params,
60 typename OilWaterMaterialLawT::Params> >
61class EclDefaultMaterial : public TraitsT
62{
63public:
64 using GasOilMaterialLaw = GasOilMaterialLawT;
65 using OilWaterMaterialLaw = OilWaterMaterialLawT;
66
67 // some safety checks
68 static_assert(TraitsT::numPhases == 3,
69 "The number of phases considered by this capillary pressure "
70 "law is always three!");
71 static_assert(GasOilMaterialLaw::numPhases == 2,
72 "The number of phases considered by the gas-oil capillary "
73 "pressure law must be two!");
74 static_assert(OilWaterMaterialLaw::numPhases == 2,
75 "The number of phases considered by the oil-water capillary "
76 "pressure law must be two!");
77 static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
78 typename OilWaterMaterialLaw::Scalar>::value,
79 "The two two-phase capillary pressure laws must use the same "
80 "type of floating point values.");
81
82 static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
83 "The gas-oil material law must implement the two-phase saturation "
84 "only API to for the default Ecl capillary pressure law!");
85 static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
86 "The oil-water material law must implement the two-phase saturation "
87 "only API to for the default Ecl capillary pressure law!");
88
89 using Traits = TraitsT;
90 using Params = ParamsT;
91 using Scalar = typename Traits::Scalar;
92
93 static constexpr int numPhases = 3;
94 static constexpr int waterPhaseIdx = Traits::wettingPhaseIdx;
95 static constexpr int oilPhaseIdx = Traits::nonWettingPhaseIdx;
96 static constexpr int gasPhaseIdx = Traits::gasPhaseIdx;
97
100 static constexpr bool implementsTwoPhaseApi = false;
101
104 static constexpr bool implementsTwoPhaseSatApi = false;
105
108 static constexpr bool isSaturationDependent = true;
109
112 static constexpr bool isPressureDependent = false;
113
116 static constexpr bool isTemperatureDependent = false;
117
120 static constexpr bool isCompositionDependent = false;
121
136 template <class ContainerT, class FluidState>
137 static void capillaryPressures(ContainerT& values,
138 const Params& params,
139 const FluidState& state)
140 {
141 OPM_TIMEFUNCTION_LOCAL();
142 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
143 values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
144 values[oilPhaseIdx] = 0;
145 values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
146
147 Valgrind::CheckDefined(values[gasPhaseIdx]);
148 Valgrind::CheckDefined(values[oilPhaseIdx]);
149 Valgrind::CheckDefined(values[waterPhaseIdx]);
150 }
151
152 /*
153 * Hysteresis parameters for oil-water
154 * @see EclHysteresisTwoPhaseLawParams::soMax(...)
155 * @see EclHysteresisTwoPhaseLawParams::swMax(...)
156 * @see EclHysteresisTwoPhaseLawParams::swMin(...)
157 * \param params Parameters
158 */
159 static void oilWaterHysteresisParams(Scalar& soMax,
160 Scalar& swMax,
161 Scalar& swMin,
162 const Params& params)
163 {
164 soMax = 1.0 - params.oilWaterParams().krnSwMdc();
165 swMax = params.oilWaterParams().krwSwMdc();
166 swMin = params.oilWaterParams().pcSwMdc();
167 Valgrind::CheckDefined(soMax);
168 Valgrind::CheckDefined(swMax);
169 Valgrind::CheckDefined(swMin);
170 }
171
172 /*
173 * Hysteresis parameters for oil-water
174 * @see EclHysteresisTwoPhaseLawParams::soMax(...)
175 * @see EclHysteresisTwoPhaseLawParams::swMax(...)
176 * @see EclHysteresisTwoPhaseLawParams::swMin(...)
177 * \param params Parameters
178 */
179 static void setOilWaterHysteresisParams(const Scalar& soMax,
180 const Scalar& swMax,
181 const Scalar& swMin,
182 Params& params)
183 {
184 params.oilWaterParams().update(swMin, swMax, 1.0 - soMax);
185 }
186
187 /*
188 * Hysteresis parameters for gas-oil
189 * @see EclHysteresisTwoPhaseLawParams::sgMax(...)
190 * @see EclHysteresisTwoPhaseLawParams::shMax(...)
191 * @see EclHysteresisTwoPhaseLawParams::soMin(...)
192 * \param params Parameters
193 */
194 static void gasOilHysteresisParams(Scalar& sgMax,
195 Scalar& shMax,
196 Scalar& soMin,
197 const Params& params)
198 {
199 const auto Swco = params.Swl();
200 sgMax = 1.0 - params.gasOilParams().krnSwMdc() - Swco;
201 shMax = params.gasOilParams().krwSwMdc();
202 soMin = params.gasOilParams().pcSwMdc();
203
204 Valgrind::CheckDefined(sgMax);
205 Valgrind::CheckDefined(shMax);
206 Valgrind::CheckDefined(soMin);
207 }
208
209 /*
210 * Hysteresis parameters for gas-oil
211 * @see EclHysteresisTwoPhaseLawParams::sgMax(...)
212 * @see EclHysteresisTwoPhaseLawParams::shMax(...)
213 * @see EclHysteresisTwoPhaseLawParams::soMin(...)
214 * \param params Parameters
215 */
216 static void setGasOilHysteresisParams(const Scalar& sgMax,
217 const Scalar& shMax,
218 const Scalar& soMin,
219 Params& params)
220 {
221 const auto Swco = params.Swl();
222 params.gasOilParams().update(soMin, shMax, 1.0 - sgMax - Swco);
223 }
224
225 static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping)
226 {
227 const auto Swco = params.Swl();
228 return params.gasOilParams().SnTrapped(maximumTrapping) - Swco;
229 }
230
231 static Scalar trappedOilSaturation(const Params& params, bool maximumTrapping)
232 {
233 return params.oilWaterParams().SnTrapped(maximumTrapping) + params.gasOilParams().SwTrapped();
234 }
235
236 static Scalar trappedWaterSaturation(const Params& params)
237 {
238 return params.oilWaterParams().SwTrapped();
239 }
240
241 static Scalar strandedGasSaturation(const Params& params, Scalar Sg, Scalar Kg)
242 {
243 const auto Swco = params.Swl();
244 return params.gasOilParams().SnStranded(Sg, Kg) - Swco;
245 }
246
256 template <class FluidState, class Evaluation = typename FluidState::Scalar>
257 static Evaluation pcgn(const Params& params,
258 const FluidState& fs)
259 {
260 OPM_TIMEFUNCTION_LOCAL();
261 // Maximum attainable oil saturation is 1-SWL.
262 const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
263 return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
264 }
265
275 template <class FluidState, class Evaluation = typename FluidState::Scalar>
276 static Evaluation pcnw(const Params& params,
277 const FluidState& fs)
278 {
279 OPM_TIMEFUNCTION_LOCAL();
280 const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
281 return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
282 }
283
287 template <class ContainerT, class FluidState>
288 static void saturations(ContainerT& /*values*/,
289 const Params& /*params*/,
290 const FluidState& /*fluidState*/)
291 {
292 throw std::logic_error("Not implemented: saturations()");
293 }
294
298 template <class FluidState, class Evaluation = typename FluidState::Scalar>
299 static Evaluation Sg(const Params& /*params*/,
300 const FluidState& /*fluidState*/)
301 {
302 throw std::logic_error("Not implemented: Sg()");
303 }
304
308 template <class FluidState, class Evaluation = typename FluidState::Scalar>
309 static Evaluation Sn(const Params& /*params*/,
310 const FluidState& /*fluidState*/)
311 {
312 throw std::logic_error("Not implemented: Sn()");
313 }
314
318 template <class FluidState, class Evaluation = typename FluidState::Scalar>
319 static Evaluation Sw(const Params& /*params*/,
320 const FluidState& /*fluidState*/)
321 {
322 throw std::logic_error("Not implemented: Sw()");
323 }
324
340 template <class ContainerT, class FluidState>
341 static void relativePermeabilities(ContainerT& values,
342 const Params& params,
343 const FluidState& fluidState)
344 {
345 OPM_TIMEFUNCTION_LOCAL();
346 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
347
348 values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
349 values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
350 values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
351 }
352
356 template <class FluidState, class Evaluation = typename FluidState::Scalar>
357 static Evaluation krg(const Params& params,
358 const FluidState& fluidState)
359 {
360 OPM_TIMEFUNCTION_LOCAL();
361 // Maximum attainable oil saturation is 1-SWL.
362 const Evaluation sw = 1.0 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
363 return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), sw);
364 }
365
369 template <class FluidState, class Evaluation = typename FluidState::Scalar>
370 static Evaluation krw(const Params& params,
371 const FluidState& fluidState)
372 {
373 OPM_TIMEFUNCTION_LOCAL();
374 const Evaluation sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
375 return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), sw);
376 }
377
381 template <class FluidState, class Evaluation = typename FluidState::Scalar>
382 static Evaluation krn(const Params& params,
383 const FluidState& fluidState)
384 {
385 OPM_TIMEFUNCTION_LOCAL();
386 const Scalar Swco = params.Swl();
387
388 const Evaluation sw =
389 max(Evaluation(Swco),
390 decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
391
392 const Evaluation sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
393
394 const Evaluation Sw_ow = sg + sw;
395 const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
396 const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
397
398 // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
399 // < epsilon/2 and interpolate between the oridinary and the regularized kro between
400 // epsilon and epsilon/2
401 constexpr const Scalar epsilon = 1e-5;
402 if (scalarValue(Sw_ow) - Swco < epsilon) {
403 const Evaluation kro2 = (kro_ow + kro_go)/2;
404 if (scalarValue(Sw_ow) - Swco > epsilon/2) {
405 const Evaluation kro1 = (sg * kro_go + (sw - Swco) * kro_ow) / (Sw_ow - Swco);
406 const Evaluation alpha = (epsilon - (Sw_ow - Swco)) / (epsilon / 2);
407
408 return kro2 * alpha + kro1 * (1 - alpha);
409 }
410
411 return kro2;
412 }
413
414 return (sg * kro_go + (sw - Swco) * kro_ow) / (Sw_ow - Swco);
415 }
416
420 template <class Evaluation, class FluidState>
421 static Evaluation relpermOilInOilGasSystem(const Params& params,
422 const FluidState& fluidState)
423 {
424 OPM_TIMEFUNCTION_LOCAL();
425 const Evaluation sw =
426 max(Evaluation{ params.Swl() },
427 decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
428
429 const Evaluation sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
430 const Evaluation So_go = 1.0 - (sg + sw);
431
432 return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
433 }
434
438 template <class Evaluation, class FluidState>
439 static Evaluation relpermOilInOilWaterSystem(const Params& params,
440 const FluidState& fluidState)
441 {
442 OPM_TIMEFUNCTION_LOCAL();
443 const Evaluation sw =
444 max(Evaluation{ params.Swl() },
445 decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
446
447 const Evaluation sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
448 const Evaluation Sw_ow = sg + sw;
449
450 return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
451 }
452
460 template <class FluidState>
461 static bool updateHysteresis(Params& params, const FluidState& fluidState)
462 {
463 const Scalar Swco = params.Swl();
464 const Scalar sw = clampSaturation(fluidState, waterPhaseIdx);
465 const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
466 const Scalar sg = clampSaturation(fluidState, gasPhaseIdx);
467 bool owChanged = params.oilWaterParams().update(/*pcSw=*/sw, /*krwSw=*/sw, /*krnSw=*/1 - So);
468 bool gochanged = params.gasOilParams().update(/*pcSw=*/ So,
469 /*krwSw=*/ So,
470 /*krnSw=*/ 1.0 - Swco - sg);
471 return owChanged || gochanged;
472 }
473
474 template <class FluidState>
475 static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
476 {
477 OPM_TIMEFUNCTION_LOCAL();
478 const auto sat = scalarValue(fluidState.saturation(phaseIndex));
479 return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
480 }
481};
482
483} // namespace Opm
484
485#endif
Default implementation for the parameters required by the default three-phase capillary pressure mode...
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 default three phase capillary pressure law used by the ECLipse simulator.
Definition EclDefaultMaterial.hpp:62
static Evaluation relpermOilInOilWaterSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/water system.
Definition EclDefaultMaterial.hpp:439
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition EclDefaultMaterial.hpp:357
static Evaluation pcgn(const Params &params, const FluidState &fs)
Capillary pressure between the gas and the non-wetting liquid (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:257
static constexpr bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition EclDefaultMaterial.hpp:120
static constexpr bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition EclDefaultMaterial.hpp:116
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition EclDefaultMaterial.hpp:299
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/gas system.
Definition EclDefaultMaterial.hpp:421
static constexpr bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition EclDefaultMaterial.hpp:100
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 EclDefaultMaterial.hpp:137
static void saturations(ContainerT &, const Params &, const FluidState &)
The inverse of the capillary pressure.
Definition EclDefaultMaterial.hpp:288
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:309
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:382
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 EclDefaultMaterial.hpp:276
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition EclDefaultMaterial.hpp:108
static bool updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition EclDefaultMaterial.hpp:461
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition EclDefaultMaterial.hpp:341
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability of the wetting phase.
Definition EclDefaultMaterial.hpp:370
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition EclDefaultMaterial.hpp:319
static constexpr bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition EclDefaultMaterial.hpp:112
static constexpr bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition EclDefaultMaterial.hpp:104
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30