Skip to content
Snippets Groups Projects
Commit 3b60014c authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Merge branch 'GeneralizeVanGenuchtensExponents' into 'master'

Generalize van-Genuchten's exponents

See merge request ogs/ogs!5126
parents cd9e9ad1 19a55e4b
No related branches found
No related tags found
No related merge requests found
Showing
with 83 additions and 45 deletions
\copydoc MaterialPropertyLib::SaturationVanGenuchten \copydoc MaterialPropertyLib::SaturationVanGenuchten
Either the `exponent` (which is the pressure exponent \f$m\f$ in the equations)
must be set and then the saturation exponent is \f$n = 1 / (1 - m)\f$.
Or both values the `pressure_exponent` and the `saturation_exponent` must be set
independently.
The exponent of the van Genuchten saturation function. The "pressure" exponent \f$m\f$ of the van Genuchten saturation function.
The "saturation" exponent is then set to \f$n = 1 / (1 - m)\f$.
See `pressure_exponent` and `saturation_exponent` to set the values
independently.
The "pressure" exponent \f$m\f$ of the van Genuchten saturation function.
The "saturation" exponent \f$n\f$ must be set as well.
The "saturation" exponent \f$n\f$ of the van Genuchten saturation function.
The "pressure" exponent \f$m\f$ must be set as well.
...@@ -31,14 +31,31 @@ std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten( ...@@ -31,14 +31,31 @@ std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
auto const residual_gas_saturation = auto const residual_gas_saturation =
//! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_gas_saturation} //! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_gas_saturation}
config.getConfigParameter<double>("residual_gas_saturation"); config.getConfigParameter<double>("residual_gas_saturation");
auto const exponent =
//! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent} double pressure_exponent;
config.getConfigParameter<double>("exponent"); double saturation_exponent;
if (auto const optional_exponent =
//! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent}
config.getConfigParameterOptional<double>("exponent"))
{
pressure_exponent = *optional_exponent;
saturation_exponent = 1.0 / (1.0 - pressure_exponent);
}
else
{
pressure_exponent =
//! \ogs_file_param{properties__property__SaturationVanGenuchten__pressure_exponent}
config.getConfigParameter<double>("pressure_exponent");
saturation_exponent =
//! \ogs_file_param{properties__property__SaturationVanGenuchten__saturation_exponent}
config.getConfigParameter<double>("saturation_exponent");
}
//! \ogs_file_param{properties__property__SaturationVanGenuchten__p_b} //! \ogs_file_param{properties__property__SaturationVanGenuchten__p_b}
auto const p_b = config.getConfigParameter<double>("p_b"); auto const p_b = config.getConfigParameter<double>("p_b");
return std::make_unique<SaturationVanGenuchten>( return std::make_unique<SaturationVanGenuchten>(
std::move(property_name), residual_liquid_saturation, std::move(property_name), residual_liquid_saturation,
residual_gas_saturation, exponent, p_b); residual_gas_saturation, pressure_exponent, saturation_exponent, p_b);
} }
} // namespace MaterialPropertyLib } // namespace MaterialPropertyLib
...@@ -21,11 +21,13 @@ SaturationVanGenuchten::SaturationVanGenuchten( ...@@ -21,11 +21,13 @@ SaturationVanGenuchten::SaturationVanGenuchten(
std::string name, std::string name,
double const residual_liquid_saturation, double const residual_liquid_saturation,
double const residual_gas_saturation, double const residual_gas_saturation,
double const exponent, double const pressure_exponent,
double const saturation_exponent,
double const p_b) double const p_b)
: S_L_res_(residual_liquid_saturation), : S_L_res_(residual_liquid_saturation),
S_L_max_(1. - residual_gas_saturation), S_L_max_(1. - residual_gas_saturation),
m_(exponent), m_(pressure_exponent),
n_(saturation_exponent),
p_b_(p_b) p_b_(p_b)
{ {
name_ = std::move(name); name_ = std::move(name);
...@@ -33,10 +35,18 @@ SaturationVanGenuchten::SaturationVanGenuchten( ...@@ -33,10 +35,18 @@ SaturationVanGenuchten::SaturationVanGenuchten(
if (!(m_ > 0 && m_ < 1)) if (!(m_ > 0 && m_ < 1))
{ {
OGS_FATAL( OGS_FATAL(
"The exponent value m = {:g} of van Genuchten saturation model, is " "The pressure exponent value m = {:g} of van Genuchten saturation "
"out of its range of (0, 1)", "model, is out of its range of (0, 1)",
m_); m_);
} }
if (n_ < 1)
{
OGS_FATAL(
"The saturation exponent value n = {:g} of van Genuchten "
"saturation model, is out of its range of [1, +inf)",
n_);
}
} }
PropertyDataType SaturationVanGenuchten::value( PropertyDataType SaturationVanGenuchten::value(
...@@ -52,8 +62,7 @@ PropertyDataType SaturationVanGenuchten::value( ...@@ -52,8 +62,7 @@ PropertyDataType SaturationVanGenuchten::value(
} }
double const p = p_cap / p_b_; double const p = p_cap / p_b_;
double const n = 1. / (1. - m_); double const p_to_n = std::pow(p, n_);
double const p_to_n = std::pow(p, n);
double const S_eff = std::pow(p_to_n + 1., -m_); double const S_eff = std::pow(p_to_n + 1., -m_);
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_; double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
...@@ -80,8 +89,7 @@ PropertyDataType SaturationVanGenuchten::dValue( ...@@ -80,8 +89,7 @@ PropertyDataType SaturationVanGenuchten::dValue(
} }
double const p = p_cap / p_b_; double const p = p_cap / p_b_;
double const n = 1. / (1. - m_); double const p_to_n = std::pow(p, n_);
double const p_to_n = std::pow(p, n);
double const S_eff = std::pow(p_to_n + 1., -m_); double const S_eff = std::pow(p_to_n + 1., -m_);
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_; double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
...@@ -91,9 +99,8 @@ PropertyDataType SaturationVanGenuchten::dValue( ...@@ -91,9 +99,8 @@ PropertyDataType SaturationVanGenuchten::dValue(
return 0.; return 0.;
} }
double const dS_eff_dp_cap = -m_ * std::pow(p, n - 1) * double const dS_eff_dp_cap =
std::pow(1 + p_to_n, -1. - m_) / -m_ * n_ * p_to_n * S_eff / (p_cap * (p_to_n + 1.));
(p_b_ * (1. - m_));
return dS_eff_dp_cap * (S_L_max_ - S_L_res_); return dS_eff_dp_cap * (S_L_max_ - S_L_res_);
} }
...@@ -117,8 +124,7 @@ PropertyDataType SaturationVanGenuchten::d2Value( ...@@ -117,8 +124,7 @@ PropertyDataType SaturationVanGenuchten::d2Value(
} }
double const p = p_cap / p_b_; double const p = p_cap / p_b_;
double const n = 1. / (1. - m_); double const p_to_n = std::pow(p, n_);
double const p_to_n = std::pow(p, n);
double const S_eff = std::pow(p_to_n + 1., -m_); double const S_eff = std::pow(p_to_n + 1., -m_);
double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_; double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
...@@ -129,8 +135,8 @@ PropertyDataType SaturationVanGenuchten::d2Value( ...@@ -129,8 +135,8 @@ PropertyDataType SaturationVanGenuchten::d2Value(
} }
double const d2S_eff_dp_cap2 = double const d2S_eff_dp_cap2 =
m_ * p_to_n * std::pow(p_to_n + 1., -m_ - 2.) * (p_to_n - m_) / m_ * n_ * n_ * p_to_n * S_eff * (p_to_n - m_) /
(p_cap * p_cap * (m_ - 1.) * (m_ - 1.)); ((p_cap * p_to_n + p_cap) * (p_cap * p_to_n + p_cap));
return d2S_eff_dp_cap2 * (S_L_max_ - S_L_res_); return d2S_eff_dp_cap2 * (S_L_max_ - S_L_res_);
} }
} // namespace MaterialPropertyLib } // namespace MaterialPropertyLib
...@@ -24,13 +24,12 @@ class Component; ...@@ -24,13 +24,12 @@ class Component;
* \f[S_\text{eff}=\frac{S-S_r}{S_{\text{max}}-S_r}.\f] * \f[S_\text{eff}=\frac{S-S_r}{S_{\text{max}}-S_r}.\f]
* Above, \f$S_r\f$ and \f$S_{\text{max}}\f$ are the residual and the maximum * Above, \f$S_r\f$ and \f$S_{\text{max}}\f$ are the residual and the maximum
* saturations. * saturations.
* The exponent \f$m \in (0,1)\f$ and the pressure scaling parameter \f$p_b\f$ * The (pressure) exponent \f$m \in (0,1)\f$ and the pressure scaling parameter
* (it is equal to \f$\rho g/\alpha\f$ in original publication) are given by the * \f$p_b\f$ (it is equal to \f$\rho g/\alpha\f$ in original publication) are
* user. * given by the user. The scaling parameter \f$p_b\f$ is given in same units as
* The scaling parameter \f$p_b\f$ is given in same units as pressure. * pressure.
* *
* In the original work another exponent \f$n\f$ is used, but usually set to * Another (saturaton) exponent \f$n\f$ is usually set to \f$n = 1 / (1 - m)\f$.
* \f$n = 1 / (1 - m)\f$, and also in this implementation.
* *
* The saturation is computed from the capillary pressure as follows: * The saturation is computed from the capillary pressure as follows:
* \f[S(p_c)= * \f[S(p_c)=
...@@ -49,7 +48,8 @@ public: ...@@ -49,7 +48,8 @@ public:
SaturationVanGenuchten(std::string name, SaturationVanGenuchten(std::string name,
double const residual_liquid_saturation, double const residual_liquid_saturation,
double const residual_gas_saturation, double const residual_gas_saturation,
double const exponent, double const pressure_exponent,
double const saturation_exponent,
double const p_b); double const p_b);
void checkScale() const override void checkScale() const override
...@@ -83,6 +83,7 @@ private: ...@@ -83,6 +83,7 @@ private:
double const S_L_res_; double const S_L_res_;
double const S_L_max_; double const S_L_max_;
double const m_; double const m_;
double const n_;
double const p_b_; double const p_b_;
}; };
} // namespace MaterialPropertyLib } // namespace MaterialPropertyLib
...@@ -21,21 +21,20 @@ TEST(MaterialPropertyLib, CapillaryPressureVanGenuchten) ...@@ -21,21 +21,20 @@ TEST(MaterialPropertyLib, CapillaryPressureVanGenuchten)
{ {
double const residual_liquid_saturation = 0.1; double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05; double const residual_gas_saturation = 0.05;
double const exponent = 0.79; double const pressure_exponent = 0.79;
double const saturation_exponent = 1. / (1. - pressure_exponent);
double const p_b = 10000; double const p_b = 10000;
double const maximum_capillary_pressure = 20000; double const maximum_capillary_pressure = 20000;
MPL::Property const& pressure = MPL::Property const& pressure = MPL::CapillaryPressureVanGenuchten{
MPL::CapillaryPressureVanGenuchten{"capillary_pressure", "capillary_pressure", residual_liquid_saturation,
residual_liquid_saturation, residual_gas_saturation, pressure_exponent,
residual_gas_saturation, // TODO? saturation_exponent,
exponent, p_b, maximum_capillary_pressure};
p_b,
maximum_capillary_pressure}; MPL::Property const& saturation = MPL::SaturationVanGenuchten{
"saturation", residual_liquid_saturation, residual_gas_saturation,
MPL::Property const& saturation = pressure_exponent, saturation_exponent, p_b};
MPL::SaturationVanGenuchten{"saturation", residual_liquid_saturation,
residual_gas_saturation, exponent, p_b};
MPL::VariableArray variable_array; MPL::VariableArray variable_array;
ParameterLib::SpatialPosition const pos; ParameterLib::SpatialPosition const pos;
......
...@@ -20,12 +20,13 @@ TEST(MaterialPropertyLib, SaturationVanGenuchten) ...@@ -20,12 +20,13 @@ TEST(MaterialPropertyLib, SaturationVanGenuchten)
{ {
double const residual_liquid_saturation = 0.1; double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05; double const residual_gas_saturation = 0.05;
double const exponent = 0.79; double const pressure_exponent = 0.79;
double const saturation_exponent = 1. / (1. - pressure_exponent);
double const p_b = 5000; double const p_b = 5000;
MPL::Property const& pressure_saturation = MPL::Property const& pressure_saturation = MPL::SaturationVanGenuchten{
MPL::SaturationVanGenuchten{"saturation", residual_liquid_saturation, "saturation", residual_liquid_saturation, residual_gas_saturation,
residual_gas_saturation, exponent, p_b}; pressure_exponent, saturation_exponent, p_b};
MPL::VariableArray variable_array; MPL::VariableArray variable_array;
ParameterLib::SpatialPosition const pos; ParameterLib::SpatialPosition const pos;
......
...@@ -78,11 +78,12 @@ TEST(RichardsMechanics, computeMicroPorosity) ...@@ -78,11 +78,12 @@ TEST(RichardsMechanics, computeMicroPorosity)
// Saturation // Saturation
constexpr double residual_liquid_saturation = 0; constexpr double residual_liquid_saturation = 0;
constexpr double residual_gas_saturation = 0; constexpr double residual_gas_saturation = 0;
constexpr double exponent = 0.4; constexpr double pressure_exponent = 0.4;
constexpr double saturation_exponent = 1. / (1. - pressure_exponent);
constexpr double p_b = 300e6; constexpr double p_b = 300e6;
MPL::Property const& saturation_micro = MPL::SaturationVanGenuchten{ MPL::Property const& saturation_micro = MPL::SaturationVanGenuchten{
"saturation_micro", residual_liquid_saturation, residual_gas_saturation, "saturation_micro", residual_liquid_saturation, residual_gas_saturation,
exponent, p_b}; pressure_exponent, saturation_exponent, p_b};
// Swelling // Swelling
constexpr std::array swelling_pressures{1e7, 1e7, 1e7}; constexpr std::array swelling_pressures{1e7, 1e7, 1e7};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment