OGS
RelPermBrooksCoreyNonwettingPhase.cpp
Go to the documentation of this file.
1 
16 
17 #include <algorithm>
18 #include <cmath>
19 
20 #include "MaterialLib/MPL/Medium.h"
21 
22 namespace MaterialPropertyLib
23 {
25  std::string name,
26  const double residual_liquid_saturation,
27  const double residual_gas_saturation,
28  const double min_relative_permeability,
29  const double exponent)
30  : residual_liquid_saturation_(residual_liquid_saturation),
31  residual_gas_saturation_(residual_gas_saturation),
32  min_relative_permeability_(min_relative_permeability),
33  exponent_(exponent)
34 {
35  name_ = std::move(name);
36 };
37 
39  VariableArray const& variable_array,
40  ParameterLib::SpatialPosition const& pos, double const t,
41  double const dt) const
42 {
47  auto const s_L = std::visit(
48  [&variable_array, &pos, t, dt](auto&& scale) -> double
49  {
50  return scale->property(PropertyType::saturation)
51  .template value<double>(variable_array, pos, t, dt);
52  },
53  scale_);
54 
55  auto const s_L_res = residual_liquid_saturation_;
56  auto const s_L_max = 1. - residual_gas_saturation_;
57 
58  auto const lambda = exponent_;
59 
60  auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
61 
62  if (s_eff >= 1.0)
63  {
64  // fully saturated medium
66  }
67  if (s_eff <= 0.0)
68  {
69  // dry medium
70  return 1.0;
71  }
72 
73  auto const k_rel_GR = (1. - s_eff) * (1. - s_eff) *
74  (1. - std::pow(s_eff, (2. + lambda) / lambda));
75 
76  return std::max(k_rel_GR, min_relative_permeability_);
77 }
79  VariableArray const& variable_array, Variable const variable,
80  ParameterLib::SpatialPosition const& pos, double const t,
81  double const dt) const
82 {
83  if (variable != Variable::liquid_saturation)
84  {
85  OGS_FATAL(
86  "RelPermBrooksCoreyNonwettingPhase::dValue is implemented for "
87  "derivatives with respect to liquid saturation only.");
88  }
93  auto const s_L = std::visit(
94  [&variable_array, &pos, t, dt](auto&& scale) -> double
95  {
96  return scale->property(PropertyType::saturation)
97  .template value<double>(variable_array, pos, t, dt);
98  },
99  scale_);
100 
101  auto const s_L_res = residual_liquid_saturation_;
102  auto const s_L_max = 1. - residual_gas_saturation_;
103  auto const lambda = exponent_;
104 
105  auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
106  if ((s_eff < 0.) || (s_eff > 1.))
107  {
108  return 0.0;
109  }
110 
111  auto const twoL_L = (2. + lambda) / lambda;
112  auto const d_se_d_sL = 1. / (s_L_max - s_L_res);
113  auto const dk_rel_GRdse =
114  -2. * (1 - s_eff) * (1. - std::pow(s_eff, twoL_L)) -
115  twoL_L * std::pow(s_eff, twoL_L - 1.) * (1. - s_eff) * (1. - s_eff);
116 
117  return dk_rel_GRdse * d_se_d_sL;
118 }
119 
120 } // 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
RelPermBrooksCoreyNonwettingPhase(std::string name, const double residual_liquid_saturation, const double residual_gas_saturation, const double min_relative_permeability, const double exponent)
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