OGS
RelPermBrooksCorey.cpp
Go to the documentation of this file.
1 
15 #include "RelPermBrooksCorey.h"
16 
17 #include <algorithm>
18 #include <cmath>
19 
20 #include "MaterialLib/MPL/Medium.h"
21 
22 namespace MaterialPropertyLib
23 {
25  const double residual_liquid_saturation,
26  const double residual_gas_saturation,
27  const double min_relative_permeability,
28  const double exponent)
29  : residual_liquid_saturation_(residual_liquid_saturation),
30  residual_gas_saturation_(residual_gas_saturation),
31  min_relative_permeability_(min_relative_permeability),
32  exponent_(exponent)
33 {
34  name_ = std::move(name);
35 };
36 
38  VariableArray const& variable_array,
39  ParameterLib::SpatialPosition const& pos, double const t,
40  double const dt) const
41 {
46  auto const s_L = std::visit(
47  [&variable_array, &pos, t, dt](auto&& scale) -> double
48  {
49  return scale->property(PropertyType::saturation)
50  .template value<double>(variable_array, pos, t, dt);
51  },
52  scale_);
53 
54  auto const s_L_res = residual_liquid_saturation_;
55  auto const s_L_max = 1. - residual_gas_saturation_;
56 
57  auto const lambda = exponent_;
58 
59  auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
60 
61  if (s_eff >= 1.0)
62  {
63  // fully saturated medium
64  return 1.0;
65  }
66  if (s_eff <= 0.0)
67  {
68  // dry medium
70  }
71 
72  auto const k_rel_LR = std::pow(s_eff, (2. + 3. * lambda) / lambda);
73 
74  return std::max(k_rel_LR, min_relative_permeability_);
75 }
77  VariableArray const& variable_array, Variable const variable,
78  ParameterLib::SpatialPosition const& pos, double const t,
79  double const dt) const
80 {
81  if (variable != Variable::liquid_saturation)
82  {
83  OGS_FATAL(
84  "RelPermBrooksCorey::dValue is implemented for derivatives with "
85  "respect to liquid saturation only.");
86  }
87 
92  auto const s_L = std::visit(
93  [&variable_array, &pos, t, dt](auto&& scale) -> double
94  {
95  return scale->property(PropertyType::saturation)
96  .template value<double>(variable_array, pos, t, dt);
97  },
98  scale_);
99 
100  auto const s_L_res = residual_liquid_saturation_;
101  auto const s_L_max = 1. - residual_gas_saturation_;
102  auto const lambda = exponent_;
103 
104  auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
105  if ((s_eff < 0.) || (s_eff > 1.))
106  {
107  return 0.;
108  }
109 
110  auto const d_se_d_sL = 1. / (s_L_max - s_L_res);
111  auto const dk_rel_LRdse =
112  (3 * lambda + 2.) / lambda * std::pow(s_eff, 2. / lambda + 2.);
113 
114  return dk_rel_LRdse * d_se_d_sL;
115 }
116 
117 } // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition: Error.h:26
virtual PropertyDataType value() const
Definition: Property.cpp:72
std::variant< Medium *, Phase *, Component * > scale_
Definition: Property.h:287
RelPermBrooksCorey(std::string name, const double, const double, const double, const double)
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
std::variant< double, Eigen::Matrix< double, 2, 1 >, Eigen::Matrix< double, 3, 1 >, Eigen::Matrix< double, 2, 2 >, Eigen::Matrix< double, 3, 3 >, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > PropertyDataType
Definition: Property.h:35
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:108
void scale(PETScVector &x, double const a)
Definition: LinAlg.cpp:44