diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md new file mode 100644 index 00000000000..5637495c26a --- /dev/null +++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md @@ -0,0 +1 @@ +The relative permeability function proposed by Brooks&Corey (1964) \cite BrooksCorey1964. diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md new file mode 100644 index 00000000000..47e4b2a9521 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md @@ -0,0 +1 @@ +The minimal relative permeability of the gas (or non-wetting) phase. diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md new file mode 100644 index 00000000000..9cb5411f550 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md @@ -0,0 +1 @@ +The minimal relative permeability of the liquid (or wetting) phase. diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md new file mode 100644 index 00000000000..5015ac5ec53 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md @@ -0,0 +1 @@ +Exponent of the Brooks&Corey relative permeability function. diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md new file mode 100644 index 00000000000..d8631eb27f4 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md @@ -0,0 +1 @@ +The smallest degree of saturation of the gas (or non-wetting) phase in the medium. diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md new file mode 100644 index 00000000000..95ce0ed4b31 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md @@ -0,0 +1 @@ +The smallest degree of saturation of the liquid (or wetting) phase in the medium. diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp index 02646eddf90..0c165ece888 100644 --- a/MaterialLib/MPL/CreateProperty.cpp +++ b/MaterialLib/MPL/CreateProperty.cpp @@ -65,6 +65,11 @@ std::unique_ptr createProperty( return createSaturationBrooksCorey(config); } + if (boost::iequals(property_type, "RelPermBrooksCorey")) + { + return createRelPermBrooksCorey(config); + } + // If none of the above property types are found, OGS throws an error. OGS_FATAL("The specified component property type '%s' was not recognized", property_type.c_str()); diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h index 646a6e2213a..386d7be0313 100644 --- a/MaterialLib/MPL/Properties/CreateProperties.h +++ b/MaterialLib/MPL/Properties/CreateProperties.h @@ -17,4 +17,5 @@ #include "CreateIdealGasLaw.h" #include "CreateLinearProperty.h" #include "CreateParameterProperty.h" +#include "CreateRelPermBrooksCorey.h" #include "CreateSaturationBrooksCorey.h" \ No newline at end of file diff --git a/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp new file mode 100644 index 00000000000..b9b0f89539e --- /dev/null +++ b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp @@ -0,0 +1,51 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "BaseLib/ConfigTree.h" +#include "RelPermBrooksCorey.h" + +namespace MaterialPropertyLib +{ +std::unique_ptr createRelPermBrooksCorey( + BaseLib::ConfigTree const& config) +{ + // check is reading the parameter, not peeking it... + //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey} + config.checkConfigParameter("type", "RelPermBrooksCorey"); + DBUG("Create RelPermBrooksCorey medium property"); + + auto const residual_liquid_saturation = + //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__residual_liquid_saturation} + config.getConfigParameter("residual_liquid_saturation"); + auto const residual_gas_saturation = + //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__residual_gas_saturation} + config.getConfigParameter("residual_gas_saturation"); + auto const min_relative_permeability_liquid = + //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__k_rel_min_liquid} + config.getConfigParameter("min_relative_permeability_liquid"); + auto const min_relative_permeability_gas = + //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__k_rel_min_gas} + config.getConfigParameter("min_relative_permeability_gas"); + auto const exponent = + //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__lambda} + config.getConfigParameter("lambda"); + if (exponent <= 0.) + { + OGS_FATAL("Exponent 'lambda' must be positive."); + } + + return std::make_unique( + residual_liquid_saturation, + residual_gas_saturation, + min_relative_permeability_liquid, + min_relative_permeability_gas, + exponent); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h new file mode 100644 index 00000000000..691059e25b1 --- /dev/null +++ b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h @@ -0,0 +1,27 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#pragma once + +namespace BaseLib +{ +class ConfigTree; +} + +namespace MaterialPropertyLib +{ +class RelPermBrooksCorey; +} + +namespace MaterialPropertyLib +{ +std::unique_ptr createRelPermBrooksCorey( + BaseLib::ConfigTree const& config); +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h index bd4b40595de..24797c11709 100644 --- a/MaterialLib/MPL/Properties/Properties.h +++ b/MaterialLib/MPL/Properties/Properties.h @@ -17,4 +17,5 @@ #include "IdealGasLaw.h" #include "LinearProperty.h" #include "ParameterProperty.h" +#include "RelPermBrooksCorey.h" #include "SaturationBrooksCorey.h" diff --git a/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp b/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp new file mode 100644 index 00000000000..6ff47331871 --- /dev/null +++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp @@ -0,0 +1,111 @@ +/** + * \file + * \author Norbert Grunwald + * \date 02.07.2018 + * \brief + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "MaterialLib/MPL/Properties/RelPermBrooksCorey.h" +#include "MaterialLib/MPL/Medium.h" + +#include +#include + +namespace MaterialPropertyLib +{ +RelPermBrooksCorey::RelPermBrooksCorey( + const double residual_liquid_saturation, + const double residual_gas_saturation, + const double min_relative_permeability_liquid, + const double min_relative_permeability_gas, + const double exponent) + : _residual_liquid_saturation(residual_liquid_saturation), + _residual_gas_saturation(residual_gas_saturation), + _min_relative_permeability_liquid(min_relative_permeability_liquid), + _min_relative_permeability_gas(min_relative_permeability_gas), + _exponent(exponent){}; + +PropertyDataType RelPermBrooksCorey::value( + VariableArray const& variable_array, + ParameterLib::SpatialPosition const& pos, + double const t) const +{ + /// here, an extra computation of saturation is forced, guaranteeing a + /// correct value. In order to speed up the computing time, saturation could + /// be insertred into the primary variable array after it is computed in the + /// FEM assembly. + auto const s_L = _medium->property(PropertyType::saturation) + .template value(variable_array, pos, t); + + auto const s_L_res = _residual_liquid_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 lambda = _exponent; + + auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res); + + if (s_eff >= 1.0) + { + // fully saturated medium + return Pair{{1.0, k_rel_min_GR}}; + } + if (s_eff <= 0.0) + { + // dry medium + return Pair{{k_rel_min_LR, 1.0}}; + } + + auto const k_rel_LR = std::pow(s_eff, (2. + 3. * lambda) / lambda); + auto const k_rel_GR = (1. - s_eff) * (1. - s_eff) * + (1. - std::pow(s_eff, (2. + lambda) / lambda)); + + const Pair kRel = {std::max(k_rel_LR, k_rel_min_LR), + std::max(k_rel_GR, k_rel_min_GR)}; + return kRel; +} +PropertyDataType RelPermBrooksCorey::dValue( + VariableArray const& variable_array, Variable const primary_variable, + ParameterLib::SpatialPosition const& pos, double const t) const +{ + (void)primary_variable; + assert((primary_variable == Variable::liquid_saturation) && + "RelPermBrooksCorey::dValue is implemented for " + " derivatives with respect to liquid saturation only."); + auto const s_L = _medium->property(PropertyType::saturation) + .template value(variable_array, pos, t); + + auto const s_L_res = _residual_liquid_saturation; + auto const s_L_max = 1. - _residual_gas_saturation; + auto const lambda = _exponent; + + auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res); + if ((s_eff < 0.) || (s_eff > 1.)) + return Pair{0., 0.}; + + auto const d_se_d_sL = 1. / (s_L_max - s_L_res); + auto const dk_rel_LRdse = + (3 * lambda + 2.) / lambda * std::pow(s_eff, 2. / lambda + 2.); + + auto const dk_rel_LRdsL = dk_rel_LRdse * d_se_d_sL; + + auto const _2L_L = (2. + lambda) / lambda; + auto const dk_rel_GRdse = + -2. * (1 - s_eff) * (1. - std::pow(s_eff, _2L_L)) - + _2L_L * std::pow(s_eff, _2L_L - 1.) * (1. - s_eff) * (1. - s_eff); + + auto const dk_rel_GRdsL = dk_rel_GRdse * d_se_d_sL; + const Pair dkReldsL = {{dk_rel_LRdsL, dk_rel_GRdsL}}; + + return dkReldsL; +} + +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/RelPermBrooksCorey.h b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h new file mode 100644 index 00000000000..a387ac129a5 --- /dev/null +++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h @@ -0,0 +1,76 @@ +/** + * \file + * \author Norbert Grunwald + * \date 27.06.2018 + * \brief + * + * \copyright + * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "BaseLib/ConfigTree.h" +#include "MaterialLib/MPL/Property.h" + +namespace MaterialPropertyLib +{ +class Medium; +class Phase; +class Component; +/** + * \class RelPermBrooksCorey + * \brief Relative permeability function proposed by Brooks&Corey + * \details This property must be a medium property, it + * computes the permeability reduction due to saturation as function + * of capillary pressure. + */ +class RelPermBrooksCorey final : public Property +{ +private: + Medium* _medium; + const double _residual_liquid_saturation; + const double _residual_gas_saturation; + const double _min_relative_permeability_liquid; + const double _min_relative_permeability_gas; + const double _exponent; + +public: + RelPermBrooksCorey(const double /*residual_liquid_saturation*/, + const double /*residual_gas_saturation*/, + const double /*_min_relative_permeability_liquid*/, + const double /*_min_relative_permeability_gas*/, + const double /*exponent*/ + ); + /// This method assigns a pointer to the material object that is the owner + /// of this property + void setScale( + std::variant scale_pointer) override + { + if (std::holds_alternative(scale_pointer)) + { + _medium = std::get(scale_pointer); + } + else + { + OGS_FATAL( + "The property 'RelPermBrooksCorey' is implemented on the " + "'media' scale only."); + } + } + + /// Those methods override the base class implementations and + /// actually compute and set the property _values and _dValues. + PropertyDataType value(VariableArray const& variable_array, + ParameterLib::SpatialPosition const& pos, + double const t) const override; + PropertyDataType dValue(VariableArray const& variable_array, + Variable const variable, + ParameterLib::SpatialPosition const& pos, + double const t) const override; +}; + +} // namespace MaterialPropertyLib diff --git a/Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp b/Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp new file mode 100644 index 00000000000..a20582c2fc3 --- /dev/null +++ b/Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp @@ -0,0 +1,125 @@ +/** + * \file + * + * \copyright + * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#include +#include + +#include "TestMPL.h" +#include "Tests/TestTools.h" + +#include "MaterialLib/MPL/Medium.h" +#include "MaterialLib/MPL/Properties/RelPermBrooksCorey.h" + +TEST(MaterialPropertyLib, RelPermBrooksCorey) +{ + const double ref_lambda = 2.5; + const double ref_residual_liquid_saturation = 0.01; + const double ref_residual_gas_saturation = 0.01; + const double ref_k_rel_L_min = 1.e-9; + const double ref_k_rel_G_min = 1.e-7; + + std::stringstream m_beg; + std::stringstream m_end; + + m_beg << "\n"; + m_beg << "\n"; + m_beg << "\n"; + m_beg << " \n"; + m_beg << " saturation\n"; + m_beg << " Constant\n"; + // constant saturation value will be inserted here... + m_end << " \n"; + m_end << " \n"; + m_end << " relative_permeability\n"; + m_end << " RelPermBrooksCorey\n"; + m_end << " " + << ref_residual_liquid_saturation + << "\n"; + m_end << " " << ref_residual_gas_saturation + << "\n"; + m_end << " " << ref_lambda << "\n"; + m_end << " " << ref_k_rel_L_min + << "\n"; + m_end << " " << ref_k_rel_G_min + << "\n"; + m_end << " \n"; + m_end << "\n"; + m_end << "\n"; + + std::array ref_saturation = { + -0.01, 0.00, 0.01, 0.02, 0.04, 0.10, 0.20, 0.40, + 0.60, 0.80, 0.90, 0.96, 0.98, 0.99, 1.00, 1.01}; + + std::array ref_k_rel_L = { + 1.0000000000E-09, 1.0000000000E-09, 1.0000000000E-09, 2.7123199126E-08, + 1.7636064574E-06, 1.1467333635E-04, 1.9615735463E-03, 3.0156880778E-02, + 1.4540473747E-01, 4.4088351375E-01, 6.9346243084E-01, 8.8856788446E-01, + 9.6177504114E-01, 1.0000000000E+00, 1.0000000000E+00, 1.0000000000E+00}; + + std::array ref_k_rel_G = { + 1.0000000000E+00, 1.0000000000E+00, 1.0000000000E+00, 9.7944075784E-01, + 9.3794411438E-01, 8.1354659673E-01, 6.1592154540E-01, 2.9343532599E-01, + 9.4837870477E-02, 1.2086349932E-02, 1.3426518034E-03, 5.1003060118E-05, + 1.9046571241E-06, 1.0000000000E-07, 1.0000000000E-07, 1.0000000000E-07}; + + std::array ref_dk_rel_L_ds_L = { + 0.0000000000E+00, 0.0000000000E+00, 0.0000000000E+00, 1.0306815668E-05, + 2.2339015126E-04, 4.8417630905E-03, 3.9231470925E-02, 2.9383627425E-01, + 9.3650508877E-01, 2.1207055092E+00, 2.9608508283E+00, 3.5542715378E+00, + 3.7677785117E+00, 3.8775510204E+00, 0.0000000000E+00, 0.0000000000E+00}; + + std::array ref_dk_rel_G_ds_L = { + 0.0000000000E+00, 0.0000000000E+00, -2.0408163265E+00, + -2.0654018726E+00, -2.0807295100E+00, -2.0524729938E+00, + -1.8805652791E+00, -1.3132397982E+00, -6.8017950205E-01, + -1.8533091175E-01, -4.4178730635E-02, -5.0791425970E-03, + -5.7061547092E-04, -0.0000000000E+00, 0.0000000000E+00, + 0.0000000000E+00}; + + for (size_t idx = 0; idx < ref_saturation.size(); idx++) + { + std::stringstream m_sat; + m_sat << " " << ref_saturation[idx] << "\n"; + std::stringstream m; + m << m_beg.str() << m_sat.str() << m_end.str(); + + auto const& medium = createTestMaterial(m.str()); + MaterialPropertyLib::VariableArray variable_array; + ParameterLib::SpatialPosition const pos; + double const time = std::numeric_limits::quiet_NaN(); + + auto k_rel = + medium + ->property( + MaterialPropertyLib::PropertyType::relative_permeability) + .template value(variable_array, pos, + time); + + auto dk_rel_ds_L = + medium + ->property( + MaterialPropertyLib::PropertyType::relative_permeability) + .template dValue( + variable_array, + MaterialPropertyLib::Variable::liquid_saturation, pos, + time); + + auto k_rel_L = k_rel[0]; + auto k_rel_G = k_rel[1]; + + auto dk_rel_L_ds_L = dk_rel_ds_L[0]; + auto dk_rel_G_ds_L = dk_rel_ds_L[1]; + + ASSERT_NEAR(k_rel_L, ref_k_rel_L[idx], 1.0e-10); + ASSERT_NEAR(k_rel_G, ref_k_rel_G[idx], 1.0e-10); + ASSERT_NEAR(dk_rel_L_ds_L, ref_dk_rel_L_ds_L[idx], 1.0e-10); + ASSERT_NEAR(dk_rel_G_ds_L, ref_dk_rel_G_ds_L[idx], 1.0e-10); + } +} \ No newline at end of file