My Project
Loading...
Searching...
No Matches
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
40namespace Opm {
41
55template <class TraitsT,
56 class GasOilMaterialLawT,
57 class OilWaterMaterialLawT,
58 class ParamsT = EclStone1MaterialParams<TraitsT, GasOilMaterialLawT, OilWaterMaterialLawT> >
59class EclStone1Material : public TraitsT
60{
61public:
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::soMax(...)
151 * @see EclHysteresisTwoPhaseLawParams::swMax(...)
152 * @see EclHysteresisTwoPhaseLawParams::swMin(...)
153 * \param params Parameters
154 */
155 static void oilWaterHysteresisParams(Scalar& soMax,
156 Scalar& swMax,
157 Scalar& swMin,
158 const Params& params)
159 {
160 soMax = 1.0 - params.oilWaterParams().krnSwMdc();
161 swMax = params.oilWaterParams().krwSwMdc();
162 swMin = params.oilWaterParams().pcSwMdc();
163 Valgrind::CheckDefined(soMax);
164 Valgrind::CheckDefined(swMax);
165 Valgrind::CheckDefined(swMin);
166 }
167
168 /*
169 * Hysteresis parameters for oil-water
170 * @see EclHysteresisTwoPhaseLawParams::soMax(...)
171 * @see EclHysteresisTwoPhaseLawParams::swMax(...)
172 * @see EclHysteresisTwoPhaseLawParams::swMin(...)
173 * \param params Parameters
174 */
175 static void setOilWaterHysteresisParams(const Scalar& soMax,
176 const Scalar& swMax,
177 const Scalar& swMin,
178 Params& params)
179 {
180 params.oilWaterParams().update(swMin, swMax, 1.0 - soMax);
181 }
182
183 /*
184 * Hysteresis parameters for gas-oil
185 * @see EclHysteresisTwoPhaseLawParams::sgMax(...)
186 * @see EclHysteresisTwoPhaseLawParams::shMax(...)
187 * @see EclHysteresisTwoPhaseLawParams::soMin(...)
188 * \param params Parameters
189 */
190 static void gasOilHysteresisParams(Scalar& sgMax,
191 Scalar& shMax,
192 Scalar& soMin,
193 const Params& params)
194 {
195 const auto Swco = params.Swl();
196 sgMax = 1.0 - params.gasOilParams().krnSwMdc() - Swco;
197 shMax = params.gasOilParams().krwSwMdc();
198 soMin = params.gasOilParams().pcSwMdc();
199
200 Valgrind::CheckDefined(sgMax);
201 Valgrind::CheckDefined(shMax);
202 Valgrind::CheckDefined(soMin);
203 }
204
205 /*
206 * Hysteresis parameters for gas-oil
207 * @see EclHysteresisTwoPhaseLawParams::sgMax(...)
208 * @see EclHysteresisTwoPhaseLawParams::shMax(...)
209 * @see EclHysteresisTwoPhaseLawParams::soMin(...)
210 * \param params Parameters
211 */
212 static void setGasOilHysteresisParams(const Scalar& sgMax,
213 const Scalar& shMax,
214 const Scalar& soMin,
215 Params& params)
216 {
217 const auto Swco = params.Swl();
218 params.gasOilParams().update(soMin, shMax, 1.0 - sgMax - Swco);
219 }
220
221 static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping)
222 {
223 const auto Swco = params.Swl();
224 return params.gasOilParams().SnTrapped(maximumTrapping) - Swco;
225 }
226
227 static Scalar trappedOilSaturation(const Params& params, bool maximumTrapping)
228 {
229 return params.oilWaterParams().SnTrapped(maximumTrapping) + params.gasOilParams().SwTrapped();
230 }
231
232 static Scalar trappedWaterSaturation(const Params& params)
233 {
234 return params.oilWaterParams().SwTrapped();
235 }
236
237 static Scalar strandedGasSaturation(const Params& params, Scalar Sg, Scalar Kg)
238 {
239 const auto Swco = params.Swl();
240 return params.gasOilParams().SnStranded(Sg, Kg) - Swco;
241 }
251 template <class FluidState, class Evaluation = typename FluidState::Scalar>
252 static Evaluation pcgn(const Params& params,
253 const FluidState& fs)
254 {
255 // Maximum attainable oil saturation is 1-SWL
256 const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
257 return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
258 }
259
269 template <class FluidState, class Evaluation = typename FluidState::Scalar>
270 static Evaluation pcnw(const Params& params,
271 const FluidState& fs)
272 {
273 const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
274 Valgrind::CheckDefined(Sw);
275
276 const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
277 Valgrind::CheckDefined(result);
278
279 return result;
280 }
281
285 template <class ContainerT, class FluidState>
286 static void saturations(ContainerT& /* values */,
287 const Params& /* params */,
288 const FluidState& /* fluidState */)
289 {
290 throw std::logic_error("Not implemented: saturations()");
291 }
292
296 template <class FluidState, class Evaluation = typename FluidState::Scalar>
297 static Evaluation Sg(const Params& /* params */,
298 const FluidState& /* fluidState */)
299 {
300 throw std::logic_error("Not implemented: Sg()");
301 }
302
306 template <class FluidState, class Evaluation = typename FluidState::Scalar>
307 static Evaluation Sn(const Params& /* params */,
308 const FluidState& /* fluidState */)
309 {
310 throw std::logic_error("Not implemented: Sn()");
311 }
312
316 template <class FluidState, class Evaluation = typename FluidState::Scalar>
317 static Evaluation Sw(const Params& /* params */,
318 const FluidState& /* fluidState */)
319 {
320 throw std::logic_error("Not implemented: Sw()");
321 }
322
338 template <class ContainerT, class FluidState>
339 static void relativePermeabilities(ContainerT& values,
340 const Params& params,
341 const FluidState& fluidState)
342 {
343 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
344
345 values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
346 values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
347 values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
348 }
349
353 template <class FluidState, class Evaluation = typename FluidState::Scalar>
354 static Evaluation krg(const Params& params,
355 const FluidState& fluidState)
356 {
357 // Maximum attainable oil saturation is 1-SWL,
358 const Evaluation Sw = 1 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
359 return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
360 }
361
365 template <class FluidState, class Evaluation = typename FluidState::Scalar>
366 static Evaluation krw(const Params& params,
367 const FluidState& fluidState)
368 {
369 const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
370 return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
371 }
372
376 template <class FluidState, class Evaluation = typename FluidState::Scalar>
377 static Evaluation krn(const Params& params,
378 const FluidState& fluidState)
379 {
380 // the Eclipse docu is inconsistent with naming the variable of connate water: In
381 // some places the connate water saturation is represented by "Swl", in others
382 // "Swco" is used.
383 const Scalar Swco = params.Swl();
384
385 // oil relperm at connate water saturations (with Sg=0)
386 const Scalar krocw = params.krocw();
387
388 const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
389 const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
390
391 const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
392 const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
393
394 Evaluation beta;
395 if (Sw <= Swco)
396 beta = 1.0;
397 else {
398 // there seems to be an error in the ECL documentation: using the approach to
399 // the scaled saturations as described there leads to significant deviations
400 // from the results produced by Eclipse 100.
401 const Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
402 const Evaluation SSg = Sg/(1.0 - Swco);
403 const Evaluation SSo = 1.0 - SSw - SSg;
404
405 if (SSw >= 1.0 || SSg >= 1.0)
406 beta = 1.0;
407 else
408 beta = pow( SSo/((1 - SSw)*(1 - SSg)), params.eta());
409 }
410
411 return max(0.0, min(1.0, beta*kro_ow*kro_go/krocw));
412 }
413
417 template <class Evaluation, class FluidState>
418 static Evaluation relpermOilInOilGasSystem(const Params& params,
419 const FluidState& fluidState)
420 {
421 const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
422
423 return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - params.Swl());
424 }
425
429 template <class Evaluation, class FluidState>
430 static Evaluation relpermOilInOilWaterSystem(const Params& params,
431 const FluidState& fluidState)
432 {
433 const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
434
435 return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
436 }
437
445 template <class FluidState>
446 static bool updateHysteresis(Params& params, const FluidState& fluidState)
447 {
448 const Scalar Swco = params.Swl();
449 const Scalar Sw = clampSaturation(fluidState, waterPhaseIdx);
450 const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
451 const Scalar Sg = clampSaturation(fluidState, gasPhaseIdx);
452 bool owChanged = params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/1 - So);
453 bool gochanged = params.gasOilParams().update(/*pcSw=*/ So,
454 /*krwSw=*/ So,
455 /*krnSw=*/ 1.0 - Swco - Sg);
456 return owChanged || gochanged;
457 }
458
459 template <class FluidState>
460 static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
461 {
462 OPM_TIMEFUNCTION_LOCAL();
463 const auto sat = scalarValue(fluidState.saturation(phaseIndex));
464 return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
465 }
466};
467
468} // namespace Opm
469
470#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:366
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition EclStone1Material.hpp:377
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition EclStone1Material.hpp:339
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:286
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:297
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition EclStone1Material.hpp:354
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:317
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:252
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:430
static bool updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition EclStone1Material.hpp:446
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:270
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/gas system.
Definition EclStone1Material.hpp:418
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition EclStone1Material.hpp:307
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30