Skip to content
Snippets Groups Projects
Commit b14255b7 authored by wenqing's avatar wenqing
Browse files

[MPL/RelPermUdell] Removed the non-wetting phase part

parent 94f5e4b0
No related branches found
No related tags found
No related merge requests found
...@@ -28,27 +28,22 @@ std::unique_ptr<RelPermUdell> createRelPermUdell( ...@@ -28,27 +28,22 @@ std::unique_ptr<RelPermUdell> createRelPermUdell(
DBUG("Create RelPermUdell medium property {:s}.", property_name); DBUG("Create RelPermUdell medium property {:s}.", property_name);
auto const residual_liquid_saturation = auto const residual_liquid_saturation =
//! \ogs_file_param{properties__property__RelPermUdell__residual_liquid_saturation} //! \ogs_file_param{properties__property__RelativePermeabilityUdell__residual_liquid_saturation}
config.getConfigParameter<double>("residual_liquid_saturation"); config.getConfigParameter<double>("residual_liquid_saturation");
auto const residual_gas_saturation = auto const residual_gas_saturation =
//! \ogs_file_param{properties__property__RelPermUdell__residual_gas_saturation} //! \ogs_file_param{properties__property__RelativePermeabilityUdell__residual_gas_saturation}
config.getConfigParameter<double>("residual_gas_saturation"); config.getConfigParameter<double>("residual_gas_saturation");
auto const min_relative_permeability_liquid = auto const min_relative_permeability =
//! \ogs_file_param{properties__property__RelPermUdell__min_relative_permeability_liquid} //! \ogs_file_param{properties__property__RelativePermeabilityUdell__min_relative_permeability}
config.getConfigParameter<double>("min_relative_permeability_liquid"); config.getConfigParameter<double>("min_relative_permeability");
auto const min_relative_permeability_gas =
//! \ogs_file_param{properties__property__RelPermUdell__min_relative_permeability_gas}
config.getConfigParameter<double>("min_relative_permeability_gas");
if ((min_relative_permeability_liquid < 0) || if (min_relative_permeability < 0)
(min_relative_permeability_gas < 0))
{ {
OGS_FATAL("Minimal relative permeabilities must be non-negative."); OGS_FATAL("Minimal relative permeability must be non-negative.");
} }
return std::make_unique<RelPermUdell>( return std::make_unique<RelPermUdell>(
std::move(property_name), residual_liquid_saturation, std::move(property_name), residual_liquid_saturation,
residual_gas_saturation, min_relative_permeability_liquid, residual_gas_saturation, min_relative_permeability);
min_relative_permeability_gas);
} }
} // namespace MaterialPropertyLib } // namespace MaterialPropertyLib
...@@ -23,12 +23,10 @@ namespace MaterialPropertyLib ...@@ -23,12 +23,10 @@ namespace MaterialPropertyLib
RelPermUdell::RelPermUdell(std::string name, RelPermUdell::RelPermUdell(std::string name,
const double residual_liquid_saturation, const double residual_liquid_saturation,
const double residual_gas_saturation, const double residual_gas_saturation,
const double min_relative_permeability_liquid, const double min_relative_permeability)
const double min_relative_permeability_gas)
: residual_liquid_saturation_(residual_liquid_saturation), : residual_liquid_saturation_(residual_liquid_saturation),
residual_gas_saturation_(residual_gas_saturation), residual_gas_saturation_(residual_gas_saturation),
min_relative_permeability_liquid_(min_relative_permeability_liquid), min_relative_permeability_(min_relative_permeability)
min_relative_permeability_gas_(min_relative_permeability_gas)
{ {
name_ = std::move(name); name_ = std::move(name);
} }
...@@ -38,69 +36,58 @@ PropertyDataType RelPermUdell::value( ...@@ -38,69 +36,58 @@ PropertyDataType RelPermUdell::value(
ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/, ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
double const /*dt*/) const double const /*dt*/) const
{ {
const double s_L = std::get<double>( const double S_L = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation)]); variable_array[static_cast<int>(Variable::liquid_saturation)]);
if (std::isnan(s_L)) if (std::isnan(S_L))
{ {
OGS_FATAL("Liquid saturation not set in RelPermUdell::value()."); OGS_FATAL("Liquid saturation not set in RelPermUdell::value().");
} }
auto const s_L_res = residual_liquid_saturation_; auto const S_L_res = residual_liquid_saturation_;
auto const s_L_max = 1. - residual_gas_saturation_; auto const S_L_max = 1. - residual_gas_saturation_;
auto const k_rel_min_LR = min_relative_permeability_liquid_;
auto const k_rel_min_GR = min_relative_permeability_gas_;
auto const s = (s_L - s_L_res) / (s_L_max - s_L_res); auto const S_e = (S_L - S_L_res) / (S_L_max - S_L_res);
if (s >= 1.0) if (S_e >= 1.0)
{ {
// fully saturated medium // fully saturated medium
return Eigen::Vector2d{1.0, k_rel_min_GR}; return 1.0;
} }
if (s <= 0.0) if (S_e <= 0.0)
{ {
// dry medium // dry medium
return Eigen::Vector2d{k_rel_min_LR, 1.0}; return min_relative_permeability_;
} }
auto const k_rel_LR = s * s * s; return std::max(min_relative_permeability_, S_e * S_e * S_e);
auto const k_rel_GR = (1. - s) * (1. - s) * (1. - s);
return Eigen::Vector2d{std::max(k_rel_LR, k_rel_min_LR),
std::max(k_rel_GR, k_rel_min_GR)};
} }
PropertyDataType RelPermUdell::dValue( PropertyDataType RelPermUdell::dValue(
VariableArray const& variable_array, Variable const primary_variable, VariableArray const& variable_array, Variable const variable,
ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/, ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
double const /*dt*/) const double const /*dt*/) const
{ {
(void)primary_variable; (void)variable;
assert((primary_variable == Variable::liquid_saturation) && assert((variable == Variable::liquid_saturation) &&
"RelPermUdell::dValue is implemented for " "RelPermUdell::dValue is implemented for "
" derivatives with respect to liquid saturation only."); " derivatives with respect to liquid saturation only.");
const double s_L = std::get<double>( const double S_L = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation)]); variable_array[static_cast<int>(Variable::liquid_saturation)]);
auto const s_L_res = residual_liquid_saturation_; auto const S_L_res = residual_liquid_saturation_;
auto const s_L_max = 1. - residual_gas_saturation_; auto const S_L_max = 1. - residual_gas_saturation_;
auto const s = (s_L - s_L_res) / (s_L_max - s_L_res); auto const S_e = (S_L - S_L_res) / (S_L_max - S_L_res);
if ((s < 0.) || (s > 1.)) if ((S_e < 0.) || (S_e > 1.))
{ {
return Eigen::Vector2d{0., 0.}; return 0.;
} }
auto const d_se_d_sL = 1. / (s_L_max - s_L_res); auto const dS_e_dS_L = 1. / (S_L_max - S_L_res);
auto const dk_rel_LRdse = 3. * s * s;
auto const dk_rel_LRdsL = dk_rel_LRdse * d_se_d_sL;
auto const dk_rel_GRdse = -3. * (1. - s) * (1. - s);
auto const dk_rel_GRdsL = dk_rel_GRdse * d_se_d_sL;
return Eigen::Vector2d{dk_rel_LRdsL, dk_rel_GRdsL}; auto const dk_rel_LR_dS_e = 3. * S_e * S_e;
return dk_rel_LR_dS_e * dS_e_dS_L;
} }
} // namespace MaterialPropertyLib } // namespace MaterialPropertyLib
...@@ -31,6 +31,10 @@ class Component; ...@@ -31,6 +31,10 @@ class Component;
* - \f$S^{\text{eff}}_{\alpha}\f$ is the effective saturation of * - \f$S^{\text{eff}}_{\alpha}\f$ is the effective saturation of
* phase \f$\alpha\f$ * phase \f$\alpha\f$
* *
* This class handles the wetting (liquid) phase portion of this relative
* permeability property, i,e with \f$\alpha=L\f$ for the relative permeability
* function.
*
* \details This property must be a medium property, it * \details This property must be a medium property, it
* computes the permeability reduction due to saturation as function of * computes the permeability reduction due to saturation as function of
* capillary pressure. * capillary pressure.
...@@ -40,22 +44,20 @@ class RelPermUdell final : public Property ...@@ -40,22 +44,20 @@ class RelPermUdell final : public Property
private: private:
const double residual_liquid_saturation_; const double residual_liquid_saturation_;
const double residual_gas_saturation_; const double residual_gas_saturation_;
const double min_relative_permeability_liquid_; const double min_relative_permeability_;
const double min_relative_permeability_gas_;
public: public:
RelPermUdell(std::string name, const double residual_liquid_saturation, RelPermUdell(std::string name, const double residual_liquid_saturation,
const double residual_gas_saturation, const double residual_gas_saturation,
const double min_relative_permeability_liquid, const double min_relative_permeability);
const double min_relative_permeability_gas);
void checkScale() const override void checkScale() const override
{ {
if (!std::holds_alternative<Medium*>(scale_)) if (!std::holds_alternative<Medium*>(scale_))
{ {
OGS_FATAL( OGS_FATAL(
"The property 'RelPermUdell' is implemented on the " "The property 'RelativePermeabilityUdell' is implemented on "
"'media' scale only."); "the 'media' scale only.");
} }
} }
...@@ -63,7 +65,7 @@ public: ...@@ -63,7 +65,7 @@ public:
ParameterLib::SpatialPosition const& pos, ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override; double const t, double const dt) const override;
PropertyDataType dValue(VariableArray const& variable_array, PropertyDataType dValue(VariableArray const& variable_array,
Variable const primary_variable, Variable const variable,
ParameterLib::SpatialPosition const& pos, ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override; double const t, double const dt) const override;
}; };
......
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