/
RelPermBrooksCorey.cpp
111 lines (92 loc) · 3.82 KB
/
RelPermBrooksCorey.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
/**
* \author Norbert Grunwald
* \date 02.07.2018
* \brief
*
* \copyright
* Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#include "MaterialLib/MPL/Properties/RelPermBrooksCorey.h"
#include "MaterialLib/MPL/Medium.h"
#include <algorithm>
#include <cmath>
namespace MaterialPropertyLib
{
RelPermBrooksCorey::RelPermBrooksCorey(
const double residual_liquid_saturation,
const double residual_gas_saturation,
const double min_relative_permeability_liquid,
const double min_relative_permeability_gas,
const double exponent)
: _residual_liquid_saturation(residual_liquid_saturation),
_residual_gas_saturation(residual_gas_saturation),
_min_relative_permeability_liquid(min_relative_permeability_liquid),
_min_relative_permeability_gas(min_relative_permeability_gas),
_exponent(exponent){};
PropertyDataType RelPermBrooksCorey::value(
VariableArray const& variable_array,
ParameterLib::SpatialPosition const& pos,
double const t) const
{
/// here, an extra computation of saturation is forced, guaranteeing a
/// correct value. In order to speed up the computing time, saturation could
/// be insertred into the primary variable array after it is computed in the
/// FEM assembly.
auto const s_L = _medium->property(PropertyType::saturation)
.template value<double>(variable_array, pos, t);
auto const s_L_res = _residual_liquid_saturation;
auto const s_L_max = 1. - _residual_gas_saturation;
auto const k_rel_min_LR = _min_relative_permeability_liquid;
auto const k_rel_min_GR = _min_relative_permeability_gas;
auto const lambda = _exponent;
auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
if (s_eff >= 1.0)
{
// fully saturated medium
return Pair{{1.0, k_rel_min_GR}};
}
if (s_eff <= 0.0)
{
// dry medium
return Pair{{k_rel_min_LR, 1.0}};
}
auto const k_rel_LR = std::pow(s_eff, (2. + 3. * lambda) / lambda);
auto const k_rel_GR = (1. - s_eff) * (1. - s_eff) *
(1. - std::pow(s_eff, (2. + lambda) / lambda));
const Pair kRel = {std::max(k_rel_LR, k_rel_min_LR),
std::max(k_rel_GR, k_rel_min_GR)};
return kRel;
}
PropertyDataType RelPermBrooksCorey::dValue(
VariableArray const& variable_array, Variable const primary_variable,
ParameterLib::SpatialPosition const& pos, double const t) const
{
(void)primary_variable;
assert((primary_variable == Variable::liquid_saturation) &&
"RelPermBrooksCorey::dValue is implemented for "
" derivatives with respect to liquid saturation only.");
auto const s_L = _medium->property(PropertyType::saturation)
.template value<double>(variable_array, pos, t);
auto const s_L_res = _residual_liquid_saturation;
auto const s_L_max = 1. - _residual_gas_saturation;
auto const lambda = _exponent;
auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
if ((s_eff < 0.) || (s_eff > 1.))
return Pair{0., 0.};
auto const d_se_d_sL = 1. / (s_L_max - s_L_res);
auto const dk_rel_LRdse =
(3 * lambda + 2.) / lambda * std::pow(s_eff, 2. / lambda + 2.);
auto const dk_rel_LRdsL = dk_rel_LRdse * d_se_d_sL;
auto const _2L_L = (2. + lambda) / lambda;
auto const dk_rel_GRdse =
-2. * (1 - s_eff) * (1. - std::pow(s_eff, _2L_L)) -
_2L_L * std::pow(s_eff, _2L_L - 1.) * (1. - s_eff) * (1. - s_eff);
auto const dk_rel_GRdsL = dk_rel_GRdse * d_se_d_sL;
const Pair dkReldsL = {{dk_rel_LRdsL, dk_rel_GRdsL}};
return dkReldsL;
}
} // namespace MaterialPropertyLib