OGS
CapillaryPressureRegularizedVanGenuchten.cpp
Go to the documentation of this file.
1 
13 
14 #include <algorithm>
15 #include <cmath>
16 
17 #include "BaseLib/Error.h"
18 #include "MaterialLib/MPL/Medium.h"
21 
22 namespace MaterialPropertyLib
23 {
24 void checkSaturationRange(const double Sl)
25 {
26  if (Sl < 0 || Sl > 1)
27  {
28  OGS_FATAL("The saturation of {:e} is out of its range of [0, 1]", Sl);
29  }
30 }
31 
34  double const residual_liquid_saturation,
35  double const maximum_liquid_saturation,
36  double const exponent,
37  double const p_b)
38  : Sg_r_(1.0 - maximum_liquid_saturation),
39  Sg_max_(1.0 - residual_liquid_saturation),
40  m_(exponent),
41  p_b_(p_b),
42  PcBarvGSg_Sg_max_(getPcBarvGSg(Sg_max_)),
43  dPcdSvGBarSg_max_(getdPcdSvGBar(Sg_max_))
44 {
46 }
47 
49  VariableArray const& variable_array,
50  ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
51  double const /*dt*/) const
52 {
53  const double Sl = std::get<double>(
54  variable_array[static_cast<int>(Variable::liquid_saturation)]);
55 
57 
58  double const Sg = 1 - Sl;
59  if (!(Sg < Sg_r_ || Sg > Sg_max_))
60  {
61  return getPcBarvGSg(Sg);
62  }
63  if (Sg < Sg_r_)
64  {
65  return 0.0;
66  }
67 
69 }
70 
72  VariableArray const& variable_array, Variable const primary_variable,
73  ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
74  double const /*dt*/) const
75 {
76  if (primary_variable != Variable::liquid_saturation)
77  {
78  OGS_FATAL(
79  "CapillaryPressureRegularizedVanGenuchten::dValue is implemented "
80  "for derivatives with respect to liquid saturation only.");
81  }
82 
83  const double Sl = std::get<double>(
84  variable_array[static_cast<int>(Variable::liquid_saturation)]);
85 
87 
88  double const Sg = 1 - Sl;
89  if (!(Sg < Sg_r_ || Sg > Sg_max_))
90  {
91  return -getdPcdSvGBar(Sg);
92  }
93  if (Sg < Sg_r_)
94  {
95  return 0.0;
96  }
97 
98  return -dPcdSvGBarSg_max_;
99 }
100 
102  double const Sg) const
103 {
104  double const S_bar = getSBar(Sg);
105  return getPcvGSg(S_bar) - getPcvGSg(Sg_r_ + (Sg_max_ - Sg_r_) * xi_ / 2);
106 }
107 
109 {
110  return Sg_r_ + (1 - xi_) * (Sg - Sg_r_) + 0.5 * xi_ * (Sg_max_ - Sg_r_);
111 }
112 
114  double const Sg) const
115 {
116  double const Se = (Sg_max_ - Sg) / (Sg_max_ - Sg_r_);
117  return p_b_ * std::pow(std::pow(Se, (-1.0 / m_)) - 1.0, 1.0 - m_);
118 }
119 
121  double const Sg) const
122 {
123  double S_bar = getSBar(Sg);
124  return getdPcdSvG(S_bar) * (1 - xi_);
125 }
126 
128  const double Sg) const
129 {
130  double const n = 1 / (1 - m_);
131  double const Se = (Sg_max_ - Sg) / (Sg_max_ - Sg_r_);
132  auto const temp = std::pow(Se, (-1 / m_));
133  return p_b_ * (1 / (m_ * n)) * (1 / (Sg_max_ - Sg_r_)) *
134  std::pow(temp - 1, (1 / n) - 1) * temp / Se;
135 }
136 
137 } // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition: Error.h:26
static constexpr double xi_
parameter in regularized van Genuchten model
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
CapillaryPressureRegularizedVanGenuchten(double const residual_liquid_saturation, double const maximum_liquid_saturation, double const exponent, double const p_b)
virtual PropertyDataType value() const
Definition: Property.cpp:72
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 checkVanGenuchtenExponentRange(const double m)