OGS
CapillaryPressureVanGenuchten.cpp
Go to the documentation of this file.
1
13
14#include <algorithm>
15#include <cmath>
16
18
19namespace MaterialPropertyLib
20{
22 std::string name,
23 double const residual_liquid_saturation,
24 double const residual_gas_saturation,
25 double const exponent,
26 double const p_b,
27 double const maximum_capillary_pressure)
29 S_L_max_(1. - residual_gas_saturation),
30 m_(exponent),
31 p_b_(p_b),
32 p_cap_max_(maximum_capillary_pressure)
33{
34 name_ = std::move(name);
35
36 if (S_L_res_ < 0 || S_L_res_ > 1)
37 {
39 "Van Genuchten capillary pressure model: The residual liquid "
40 "saturation value S_L_res = {:g} is out of admissible range [0, "
41 "1].",
42 S_L_res_);
43 }
44 if (S_L_max_ < 0 || S_L_max_ > 1)
45 {
47 "Van Genuchten capillary pressure model: The maximum liquid "
48 "saturation value S_L_max = {:g} is out of admissible range [0, "
49 "1].",
50 S_L_max_);
51 }
52 if (S_L_res_ >= S_L_max_)
53 {
55 "Van Genuchten capillary pressure model: The maximum liquid "
56 "saturation S_L_max = {:g} must not be less or equal to the "
57 "residual liquid saturion S_L_res = {:g}.",
59 }
60 if (!(m_ > 0 && m_ < 1))
61 {
63 "Van Genuchten capillary pressure model: The exponent value m = "
64 "{:g} is out of of admissible range (0, 1).",
65 m_);
66 }
67 if (p_b_ <= 0)
68 {
70 "Van Genuchten capillary pressure model: The pressure scaling "
71 "value p_b = {:g} must be positive.",
72 p_b_);
73 }
74 if (p_cap_max_ < 0)
75 {
77 "Van Genuchten capillary pressure model: The maximum capillary "
78 "pressure value p_cap_max = {:g} must be non-negative.",
80 }
81}
82
84 VariableArray const& variable_array,
85 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
86 double const /*dt*/) const
87{
88 double const S_L = variable_array.liquid_saturation;
89
90 if (S_L <= S_L_res_)
91 {
92 return p_cap_max_;
93 }
94
95 if (S_L >= S_L_max_)
96 {
97 return 0.;
98 }
99
100 double const S_eff = (S_L - S_L_res_) / (S_L_max_ - S_L_res_);
101 assert(0 <= S_eff && S_eff <= 1);
102
103 double const p_cap =
104 p_b_ * std::pow(std::pow(S_eff, -1.0 / m_) - 1.0, 1.0 - m_);
105 assert(p_cap > 0);
106 return std::min(p_cap, p_cap_max_);
107}
108
110 VariableArray const& variable_array, Variable const variable,
111 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
112 double const /*dt*/) const
113{
114 if (variable != Variable::liquid_saturation)
115 {
116 OGS_FATAL(
117 "CapillaryPressureVanGenuchten::dValue is implemented for "
118 "derivatives with respect to liquid saturation only.");
119 }
120
121 double const S_L = variable_array.liquid_saturation;
122
123 if (S_L <= S_L_res_)
124 {
125 return 0.;
126 }
127
128 if (S_L >= S_L_max_)
129 {
130 return 0.;
131 }
132
133 double const S_eff = (S_L - S_L_res_) / (S_L_max_ - S_L_res_);
134
135 assert(0 < S_eff && S_eff < 1.0);
136
137 double const val1 = std::pow(S_eff, -1.0 / m_);
138 double const p_cap = p_b_ * std::pow(val1 - 1.0, 1.0 - m_);
139 if (p_cap >= p_cap_max_)
140 {
141 return 0.;
142 }
143
144 double const val2 = std::pow(val1 - 1.0, -m_);
145 return p_b_ * (m_ - 1.0) * val1 * val2 / (m_ * (S_L - S_L_res_));
146}
147} // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition Error.h:26
double const S_L_res_
Residual saturation of liquid phase.
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
double const S_L_max_
Maximum saturation of liquid phase.
CapillaryPressureVanGenuchten(std::string name, double const residual_liquid_saturation, double const residual_gas_saturation, double const exponent, double const p_b, double const maximum_capillary_pressure)
virtual PropertyDataType value() const
Definition Property.cpp:76
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