OGS
CapillaryPressureRegularizedVanGenuchten.cpp
Go to the documentation of this file.
1
13
14#include <algorithm>
15#include <cmath>
16
17#include "BaseLib/Error.h"
21
22namespace MaterialPropertyLib
23{
24void 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 = variable_array.liquid_saturation;
54
56
57 double const Sg = 1 - Sl;
58 if (!(Sg < Sg_r_ || Sg > Sg_max_))
59 {
60 return getPcBarvGSg(Sg);
61 }
62 if (Sg < Sg_r_)
63 {
64 return 0.0;
65 }
66
68}
69
71 VariableArray const& variable_array, Variable const variable,
72 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
73 double const /*dt*/) const
74{
75 if (variable != Variable::liquid_saturation)
76 {
78 "CapillaryPressureRegularizedVanGenuchten::dValue is implemented "
79 "for derivatives with respect to liquid saturation only.");
80 }
81
82 const double Sl = variable_array.liquid_saturation;
83
85
86 double const Sg = 1 - Sl;
87 if (!(Sg < Sg_r_ || Sg > Sg_max_))
88 {
89 return -getdPcdSvGBar(Sg);
90 }
91 if (Sg < Sg_r_)
92 {
93 return 0.0;
94 }
95
96 return -dPcdSvGBarSg_max_;
97}
98
100 double const Sg) const
101{
102 double const S_bar = getSBar(Sg);
103 return getPcvGSg(S_bar) - getPcvGSg(Sg_r_ + (Sg_max_ - Sg_r_) * xi_ / 2);
104}
105
107{
108 return Sg_r_ + (1 - xi_) * (Sg - Sg_r_) + 0.5 * xi_ * (Sg_max_ - Sg_r_);
109}
110
112 double const Sg) const
113{
114 double const Se = (Sg_max_ - Sg) / (Sg_max_ - Sg_r_);
115 return p_b_ * std::pow(std::pow(Se, (-1.0 / m_)) - 1.0, 1.0 - m_);
116}
117
119 double const Sg) const
120{
121 double S_bar = getSBar(Sg);
122 return getdPcdSvG(S_bar) * (1 - xi_);
123}
124
126 const double Sg) const
127{
128 double const n = 1 / (1 - m_);
129 double const Se = (Sg_max_ - Sg) / (Sg_max_ - Sg_r_);
130 auto const temp = std::pow(Se, (-1 / m_));
131 return p_b_ * (1 / (m_ * n)) * (1 / (Sg_max_ - Sg_r_)) *
132 std::pow(temp - 1, (1 / n) - 1) * temp / Se;
133}
134
135} // 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:76
void checkVanGenuchtenExponentRange(const double m)
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 >, Eigen::MatrixXd > PropertyDataType
Definition Property.h:31