Detailed Description

The Lu - Mc Cartney water retention behavior.

This property must be a medium property, it computes the saturation of the wetting phase as function of capillary pressure.

Original source: Lu, Y., McCartney, J.S. Temperature effects on adsorption and capillarity water retention mechanisms in constrained unsaturated soils. Acta Geotech. 19, 6467-6482 (2024). https://doi.org/10.1007/s11440-024-02341-9 and Guo, G., Zheng, L., Lu, Y., Behbehani, F., & McCartney, J. (2026). Coupled THM modeling of bentonite heating and hydration in tank TemperatureTypes with a new temperature-dependent water retention model. Computers and Geotechnics, 190, 107736. https://doi.org/10.1016/j.compgeo.2025.107736 Get rights and content

Definition at line 41 of file SaturationLuMcCartney.h.

#include <SaturationLuMcCartney.h>

Inheritance diagram for MaterialPropertyLib::SaturationLuMcCartney:

[legend]

Collaboration diagram for MaterialPropertyLib::SaturationLuMcCartney:

[legend]

Public Member Functions
	SaturationLuMcCartney (std::string name, std::string const &material)
void	checkScale () const override
PropertyDataType	value (VariableArray const &variable_array, ParameterLib::SpatialPosition const &, double const, double const) const override
PropertyDataType	dValue (VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &, double const, double const) const override
PropertyDataType	d2Value (VariableArray const &variable_array, Variable const variable1, Variable const variable2, ParameterLib::SpatialPosition const &, double const, double const) const override
	Default implementation: 2nd derivative of any constant property is zero.
template<typename PressureType, typename TemperatureType>
promote< PressureType, TemperatureType >	water_content (PressureType p_cap, TemperatureType T, VariableArray const &variable_array) const
template<typename PressureType, typename TemperatureType>
promote< PressureType, TemperatureType >	adsorptive_water_content (PressureType p_cap, TemperatureType T, VariableArray const &variable_array) const
template<typename PressureType, typename TemperatureType>
promote< PressureType, TemperatureType >	capillary_water_content (PressureType p_cap, TemperatureType T, VariableArray const &variable_array) const
template<typename TemperatureType>
promote< TemperatureType >	theta_a_max (TemperatureType T, VariableArray const &variable_array) const
template<typename TemperatureType>
promote< TemperatureType >	psi_max (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	c (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	CEC (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	psi_c (TemperatureType T, VariableArray const &variable_array) const
template<typename TemperatureType>
promote< TemperatureType >	A_H (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	chi (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	delta_h (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	alpha (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	epsilon_w (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	zeta_s (TemperatureType T) const
template<typename TemperatureType>
promote< TemperatureType >	rho_w (TemperatureType T) const
Public Member Functions inherited from MaterialPropertyLib::Property
virtual	~Property ()
virtual PropertyDataType	initialValue (ParameterLib::SpatialPosition const &pos, double const t) const
virtual PropertyDataType	value () const
virtual PropertyDataType	value (VariableArray const &variable_array, VariableArray const &variable_array_prev, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
virtual PropertyDataType	dValue (VariableArray const &variable_array, VariableArray const &variable_array_prev, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
virtual void	setProperties (std::vector< std::unique_ptr< Phase > > const &phases)
	Default implementation:
void	setScale (std::variant< Medium , Phase , Component * > scale)
template<typename T>
T	initialValue (ParameterLib::SpatialPosition const &pos, double const t) const
template<typename T>
T	value () const
template<typename T>
T	value (VariableArray const &variable_array, VariableArray const &variable_array_prev, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
template<typename T>
T	value (VariableArray const &variable_array, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
template<typename T>
T	dValue (VariableArray const &variable_array, VariableArray const &variable_array_prev, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
template<typename T>
T	dValue (VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
template<typename T>
T	d2Value (VariableArray const &variable_array, Variable const &variable1, Variable const &variable2, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

Private Member Functions
template<typename PressureType, typename TemperatureType>
boost::math::differentiation::promote< PressureType, TemperatureType >	water_content (PressureType p_cap, TemperatureType T, VariableArray const &variable_array) const
template<typename PressureType, typename TemperatureType>
boost::math::differentiation::promote< PressureType, TemperatureType >	adsorptive_water_content (PressureType p_cap, TemperatureType T, VariableArray const &variable_array) const
template<typename PressureType, typename TemperatureType>
boost::math::differentiation::promote< PressureType, TemperatureType >	capillary_water_content (PressureType p_cap, TemperatureType T, VariableArray const &variable_array) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	theta_a_max (TemperatureType T, VariableArray const &variable_array) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	psi_max (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	c (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	CEC (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	psi_c (TemperatureType T, VariableArray const &variable_array) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	A_H (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	chi (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	delta_h (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	alpha (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	epsilon_w (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	zeta_s (TemperatureType T) const
template<typename TemperatureType>
boost::math::differentiation::promote< TemperatureType >	rho_w (TemperatureType T) const

Private Attributes
std::string const	material_
const MaterialLib::Fluid::DimensionLessGibbsFreeEnergyRegion1	gibbs_free_energy_
double const	nu_w = 1.8e-5
	<-— Lu and McCartney 2024 parameters ---->
double const	A = 0.92
	fitting parameter for cation exchange capacity
double const	B = 0.9
	fitting parameter for cation exchange capacity
double const	T1 = 373.15
	parameter for temperature effect of CEC in K
double const	T2
	parameter for temperature effect of CEC in K
double const	nu_e
	main absorption frequency of UV light in Hz
double const	Tr = 293.15
	reference temperature in K
double const	C1 = -0.00151
	parameter for temperature correction
double const	n_w = 1.33
	diffraction index
double	E1_minus_EL
	material parameters for adsorptive water content
double	nu_mr
double	CEC_max
	maximum cation exchange capacity in meq/g
double	bw
double	M
	exponent for adsorptive water content
double	SSA
	material parameters for capillary water content
double	epsilon_s
	dielectic constant of soil particles
double	n_s
	refractive index of soil particles
double	eta_alpha
	soil related coefficient
double	alpha_0_inv
	air entry suction at reference temperature Tr in Pa
double	chi_r
double	N
	exponent for capillary water content

Static Private Attributes
static constexpr double	ref_T_ = 1386
	reference temperature in K.
static constexpr double	ref_p_ = 1.653e7

Additional Inherited Members
Protected Attributes inherited from MaterialPropertyLib::Property
std::string	name_
PropertyDataType	value_
	The single value of a property.
PropertyDataType	dvalue_
std::variant< Medium , Phase , Component * >	scale_

Constructor & Destructor Documentation

◆ SaturationLuMcCartney()

MaterialPropertyLib::SaturationLuMcCartney::SaturationLuMcCartney	(	std::string	name,
		std::string const &	material )

Definition at line 157 of file SaturationLuMcCartney.cpp.

    : material_(material),
      gibbs_free_energy_(
          MaterialLib::Fluid::DimensionLessGibbsFreeEnergyRegion1())
{
    name_ = std::move(name);
    if (material_ == "MX80")
    {
        E1_minus_EL = 7533.930402019812;
        nu_mr = 2.752035e-01;
        CEC_max = 0.8004082526701221;
        bw = 0.18;
        M = 0.105;
        N = 1.15;
        SSA = 700.0e3;
        epsilon_s = 3.2089251122574;
        n_s = 1.3298392371979677;
        eta_alpha = 5.0;
        alpha_0_inv = 3.3e6;
        chi_r = delta_h(Tr) / C1;
    }
    else if (material_ == "BoomClay")
    {
        E1_minus_EL = 6803;
        nu_mr = 7.994125e-02;
        CEC_max = 0.2534798545419938;
        bw = 0.11;
        M = 0.085;
        N = 1.29;
        SSA = 260e3;
        epsilon_s = 1.3267932723717908;
        n_s = 35.435652398543425;
        eta_alpha = 0.6;
        alpha_0_inv = 36e3;
        chi_r = delta_h(Tr) / C1;
    }
    else if (material_ == "FEBEX")
    {
        E1_minus_EL = 7904.547103777231;
        nu_mr = 3.374814e-01;
        CEC_max = 1.0108839312584308;
        bw = 0.18;
        M = 0.132;
        N = 1.22;
        SSA = 860.0e3;
        epsilon_s = 0.940367803768519;
        n_s = 748.7804881357005;
        eta_alpha = 3.5;
        alpha_0_inv = 330e3;
        chi_r = delta_h(Tr) / C1;
    }
    else if (material_ == "GMZ01")
    {
        E1_minus_EL = 8369.804359014617;
        nu_mr = 1.505956e-01;
        CEC_max = 0.8004082526701221;
        bw = 0.25;
        M = 0.135;
        N = 1.3;
        SSA = 700.0e3;
        epsilon_s = 2.1590265530454893;
        n_s = 2175.885110079861;
        eta_alpha = 8.5;
        alpha_0_inv = 8.0e6;
        chi_r = delta_h(Tr) / C1;
    }
 
    else
    {
        OGS_FATAL(
            "Material '{}' not implemented in SaturationLuMcCartney model.",
            material_);
    }
}

References alpha_0_inv, bw, C1, CEC_max, chi_r, delta_h(), E1_minus_EL, epsilon_s, eta_alpha, gibbs_free_energy_, M, material_, N, n_s, MaterialPropertyLib::name, MaterialPropertyLib::Property::name_, nu_mr, OGS_FATAL, SSA, and Tr.

Member Function Documentation

◆ A_H() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::A_H ( TemperatureType T ) const

Definition at line 94 of file SaturationLuMcCartney.cpp.

{
    auto const dielectric_term =
        pow((epsilon_s - epsilon_w(T)) / (epsilon_s + epsilon_w(T)), 2);
    auto const RI_term =
        pow(n_s * n_s - n_w * n_w, 2) / pow(n_s * n_s + n_w * n_w, 3. / 2.);
    auto const prefactor_dielectic =
        (3 * MaterialLib::PhysicalConstant::BoltzmannConstant * T) / 4.;
    auto const prefactor_RI =
        (3 * MaterialLib::PhysicalConstant::PlanckConstant * nu_e) /
        (16 * sqrt(2));
    return prefactor_dielectic * dielectric_term + prefactor_RI * RI_term;
}

References MaterialLib::PhysicalConstant::BoltzmannConstant, epsilon_s, epsilon_w(), n_s, n_w, nu_e, and MaterialLib::PhysicalConstant::PlanckConstant.

◆ A_H() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::A_H ( TemperatureType T ) const

private

Referenced by psi_c().

◆ adsorptive_water_content() [1/2]

template<typename PressureType, typename TemperatureType>

promote< PressureType, TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::adsorptive_water_content	(	PressureType	p_cap,
		TemperatureType	T,
		VariableArray const &	variable_array ) const

Definition at line 29 of file SaturationLuMcCartney.cpp.

{
    return theta_a_max(T, variable_array) *
           (1 - pow(exp((p_cap - psi_max(T)) / p_cap), M));
}

References M, psi_max(), and theta_a_max().

◆ adsorptive_water_content() [2/2]

template<typename PressureType, typename TemperatureType>

boost::math::differentiation::promote< PressureType, TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::adsorptive_water_content	(	PressureType	p_cap,
		TemperatureType	T,
		VariableArray const &	variable_array ) const

private

Referenced by capillary_water_content(), and water_content().

◆ alpha() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::alpha ( TemperatureType T ) const

Definition at line 123 of file SaturationLuMcCartney.cpp.

{
    auto const psi_aev = alpha_0_inv * exp(eta_alpha * (Tr / T - 1));
    return 1.0 / psi_aev;
}

References alpha_0_inv, eta_alpha, and Tr.

◆ alpha() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::alpha ( TemperatureType T ) const

private

◆ c() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::c ( TemperatureType T ) const

Definition at line 71 of file SaturationLuMcCartney.cpp.

{
    return exp((E1_minus_EL) /
               (MaterialLib::PhysicalConstant::IdealGasConstant * T));
}

References E1_minus_EL, and MaterialLib::PhysicalConstant::IdealGasConstant.

◆ c() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::c ( TemperatureType T ) const

private

◆ capillary_water_content() [1/2]

template<typename PressureType, typename TemperatureType>

promote< PressureType, TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::capillary_water_content	(	PressureType	p_cap,
		TemperatureType	T,
		VariableArray const &	variable_array ) const

Definition at line 39 of file SaturationLuMcCartney.cpp.

{
    auto const theta_mean =
        (variable_array.porosity -
         adsorptive_water_content(p_cap, T, variable_array)) /
        2.0;
    auto const chi_frac = (chi(T) + T) / (chi_r + Tr);
    auto const erf_term =
        erf(sqrt(2) * (p_cap / psi_c(T, variable_array) * chi_frac - 1));
    auto const power_term = pow(alpha(T) * p_cap * chi_frac, N);
    auto const cap_term =
        theta_mean * (1 - erf_term) * pow(1 + power_term, ((1 / N) - 1));
    return cap_term;
}

References adsorptive_water_content(), MaterialPropertyLib::alpha, chi(), chi_r, N, MaterialPropertyLib::VariableArray::porosity, psi_c(), and Tr.

◆ capillary_water_content() [2/2]

template<typename PressureType, typename TemperatureType>

boost::math::differentiation::promote< PressureType, TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::capillary_water_content	(	PressureType	p_cap,
		TemperatureType	T,
		VariableArray const &	variable_array ) const

private

Referenced by water_content().

◆ CEC() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::CEC ( TemperatureType T ) const

Definition at line 78 of file SaturationLuMcCartney.cpp.

{
    auto const cos_arg = A * std::numbers::pi * (B * T - T1) / (T2 - T1);
    return 0.5 * CEC_max * (cos(cos_arg) + 1);
}

References A, B, CEC_max, T1, and T2.

◆ CEC() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::CEC ( TemperatureType T ) const

private

Referenced by theta_a_max(), and zeta_s().

◆ checkScale()

void MaterialPropertyLib::SaturationLuMcCartney::checkScale ( ) const

inlineoverridevirtual

Reimplemented from MaterialPropertyLib::Property.

Definition at line 46 of file SaturationLuMcCartney.h.

    {
        if (!std::holds_alternative<Medium*>(scale_))
        {
            OGS_FATAL(
                "The property 'SaturationLuMcCartney' is implemented on the "
                "'media' scale only.");
        }
    }

References OGS_FATAL, and MaterialPropertyLib::Property::scale_.

◆ chi() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::chi ( TemperatureType T ) const

Definition at line 109 of file SaturationLuMcCartney.cpp.

{
    return delta_h(T) / C1;
}

References C1, and delta_h().

◆ chi() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::chi ( TemperatureType T ) const

private

Referenced by capillary_water_content().

◆ d2Value()

PropertyDataType MaterialPropertyLib::SaturationLuMcCartney::d2Value	(	VariableArray const &	variable_array,
		Variable const	variable1,
		Variable const	variable2,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

Default implementation: 2nd derivative of any constant property is zero.

This virtual method will compute the second derivative of a property with respect to the given variables pv1 and pv2.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 332 of file SaturationLuMcCartney.cpp.

{
    (void)variable1;
    (void)variable2;
    assert(((variable1 == Variable::capillary_pressure) ||
            (variable1 == Variable::temperature)) &&
           ((variable2 == Variable::capillary_pressure) ||
            (variable2 == Variable::temperature)) &&
           "SaturationLuMcCartney::d2Value is implemented for  derivatives "
           "with respect to capillary pressure only.");
 
    const double p_cap = variable_array.capillary_pressure;
    const double T = variable_array.temperature;
 
    auto const& medium = *std::get<Medium*>(scale_);
    auto const& solid_phase = medium.phase(PhaseName::Solid);
    auto const& porosity_property = medium[PropertyType::porosity];
    auto const& solid_denisty_property = solid_phase[PropertyType::density];
 
    auto const phi =
        std::get<double>(porosity_property.value(variable_array, pos, t, dt));
    auto const rho_s = std::get<double>(
        solid_denisty_property.value(variable_array, pos, t, dt));
    VariableArray rho_d_and_phi;
    rho_d_and_phi.density = rho_s * (1.0 - phi);
    rho_d_and_phi.porosity = phi;
 
    auto const p_cap_max = psi_max(T);
    auto const p_cut = 1e-10;
 
    if (p_cap >= p_cap_max)
    {
        return 0.;
    }
    if (p_cap <= p_cut && (variable1 == Variable::capillary_pressure ||
                           variable2 == Variable::capillary_pressure))
    {
        return 0.;
    }
 
    auto const variables = make_ftuple<double, 2, 2>(p_cap, T);
    auto const& v = std::get<0>(variables);  // Up to Nw derivatives at w=11
    auto const& w = std::get<1>(variables);
    auto const y = water_content(v, w, rho_d_and_phi);
    double DDS = 0.0;
    auto returnVariableTuple = [](const Variable v1,
                                  const Variable v2) -> std::tuple<int, int>
    {
        if (v1 == Variable::capillary_pressure &&
            v2 == Variable::capillary_pressure)
        {
            return {2, 0};
        }
        if (v1 == Variable::temperature && v2 == Variable::temperature)
        {
            return {0, 2};
        }
        return {1, 1};
    };
    auto [i, j] = returnVariableTuple(variable1, variable2);
    if (j > 0)
    {
        if (p_cap <= p_cut)
        {
            auto const variables_b2 = make_ftuple<double, 2, 2>(p_cut, T);
            auto const& v_b2 = std::get<0>(variables_b2);
            auto const& w_b2 = std::get<1>(variables_b2);
            auto const y_b2 = water_content(v_b2, w_b2, rho_d_and_phi);
            return y_b2.derivative(i, j) / phi;
        }
    }
    DDS = y.derivative(i, j);
    return DDS / phi;
}

References MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, psi_max(), MaterialPropertyLib::Property::scale_, MaterialPropertyLib::Solid, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, and water_content().

◆ delta_h() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::delta_h ( TemperatureType T ) const

Definition at line 115 of file SaturationLuMcCartney.cpp.

{
    // delta_hr = -0.516; //J / m ^ 2 high plasticity clay
    auto const delta_hr = -0.5151936137548038;
    return delta_hr * pow((1 - Tr) / (1 - T), 0.38);
}

References Tr.

◆ delta_h() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::delta_h ( TemperatureType T ) const

private

Referenced by SaturationLuMcCartney(), and chi().

◆ dValue()

PropertyDataType MaterialPropertyLib::SaturationLuMcCartney::dValue	(	VariableArray const &	variable_array,
		Variable const	variable,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

This virtual method will compute the property derivative value based on the variables that are passed as arguments with the default implementation using empty variables array for the previous time step.

The default implementation of this method only returns the property value derivative without altering it.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 268 of file SaturationLuMcCartney.cpp.

{
    if (variable != Variable::capillary_pressure &&
        variable != Variable::temperature)
    {
        OGS_FATAL(
            "SaturationLuMcCartney::dValue is implemented for derivatives "
            "with respect to capillary pressure or temperature only.");
    }
 
    const double p_cap = variable_array.capillary_pressure;
    const double T = variable_array.temperature;
 
    auto const& medium = *std::get<Medium*>(scale_);
    auto const& solid_phase = medium.phase(PhaseName::Solid);
    auto const& porosity_property = medium[PropertyType::porosity];
    auto const& solid_denisty_property = solid_phase[PropertyType::density];
 
    auto const phi =
        std::get<double>(porosity_property.value(variable_array, pos, t, dt));
    auto const rho_s = std::get<double>(
        solid_denisty_property.value(variable_array, pos, t, dt));
    VariableArray rho_d_and_phi;
    rho_d_and_phi.density = rho_s * (1.0 - phi);
    rho_d_and_phi.porosity = phi;
 
    auto const p_cap_max = psi_max(T);
    auto const p_cut = 1.0e-10;
 
    auto const variables = make_ftuple<double, 1, 1>(p_cap, T);
    auto const& v = std::get<0>(variables);
    auto const& w = std::get<1>(variables);
    auto const y = water_content(v, w, rho_d_and_phi);
    double DS = 0.0;
    if (p_cap >= p_cap_max)
    {
        return 0.0;
    }
    if (variable == Variable::capillary_pressure)
    {
        if (p_cap <= 1.0e-10)
        {
            return 0.;
        }
        return y.derivative(1, 0) / phi;
    }
    if (variable == Variable::temperature)
    {
        if (p_cap <= p_cut)
        {
            auto const variables_b2 = make_ftuple<double, 1, 1>(p_cut, T);
            auto const& v_b2 = std::get<0>(variables_b2);
            auto const& w_b2 = std::get<1>(variables_b2);
            auto const y_b2 = water_content(v_b2, w_b2, rho_d_and_phi);
            return y_b2.derivative(0, 1) / phi;
        }
        DS = y.derivative(0, 1);
    }
    return DS / phi;
}

References MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, OGS_FATAL, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, psi_max(), MaterialPropertyLib::Property::scale_, MaterialPropertyLib::Solid, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, and water_content().

◆ epsilon_w() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::epsilon_w ( TemperatureType T ) const

Definition at line 130 of file SaturationLuMcCartney.cpp.

{
    auto const r_w = rho_w(T) / 1000.0;
    return pow(
        10, (0.7017 + 642.0 / T - 1.167e5 / pow(T, 2) + 9.190e6 / pow(T, 3) +
             (1.667 - 11.41 / T - 3.526e4 / pow(T, 2)) * log(r_w) / log(10)) +
                1);
}

References rho_w().

◆ epsilon_w() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::epsilon_w ( TemperatureType T ) const

private

Referenced by A_H().

◆ psi_c() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::psi_c	(	TemperatureType	T,
		VariableArray const &	variable_array ) const

Definition at line 85 of file SaturationLuMcCartney.cpp.

{
    return A_H(T) / (6 * std::numbers::pi) *
           pow(theta_a_max(T, variable_array) / (SSA * variable_array.density),
               -3.0);
}

References A_H(), MaterialPropertyLib::VariableArray::density, SSA, and theta_a_max().

◆ psi_c() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::psi_c	(	TemperatureType	T,
		VariableArray const &	variable_array ) const

private

Referenced by capillary_water_content().

◆ psi_max() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::psi_max ( TemperatureType T ) const

Definition at line 64 of file SaturationLuMcCartney.cpp.

{
    return MaterialLib::PhysicalConstant::IdealGasConstant * T * c(T) /
           (3.0 * nu_w);
}

References MaterialPropertyLib::c, MaterialLib::PhysicalConstant::IdealGasConstant, and nu_w.

◆ psi_max() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::psi_max ( TemperatureType T ) const

private

Referenced by adsorptive_water_content(), d2Value(), dValue(), and value().

◆ rho_w() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::rho_w ( TemperatureType T ) const

Definition at line 147 of file SaturationLuMcCartney.cpp.

{
    double const pi = 0.0;
    auto const tau = ref_T_ / T;
 
    return ref_p_ /
           (MaterialLib::PhysicalConstant::SpecificGasConstant::WaterVapour *
            T * gibbs_free_energy_.get_dgamma_dpi(tau, pi));
}

References gibbs_free_energy_, ref_p_, ref_T_, and MaterialLib::PhysicalConstant::SpecificGasConstant::WaterVapour.

◆ rho_w() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::rho_w ( TemperatureType T ) const

private

Referenced by epsilon_w().

◆ theta_a_max() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::theta_a_max	(	TemperatureType	T,
		VariableArray const &	variable_array ) const

Definition at line 57 of file SaturationLuMcCartney.cpp.

{
    return (1 - variable_array.porosity) * (CEC(T) / zeta_s(T) + bw);
}

References bw, CEC(), MaterialPropertyLib::VariableArray::porosity, and zeta_s().

◆ theta_a_max() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::theta_a_max	(	TemperatureType	T,
		VariableArray const &	variable_array ) const

private

Referenced by adsorptive_water_content(), and psi_c().

◆ value()

PropertyDataType MaterialPropertyLib::SaturationLuMcCartney::value	(	VariableArray const &	variable_array,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

Those methods override the base class implementations and actually compute and set the property values_ and dValues_.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 233 of file SaturationLuMcCartney.cpp.

{
    const double p_cap = variable_array.capillary_pressure;
    const double T = variable_array.temperature;
 
    auto const& medium = *std::get<Medium*>(scale_);
    auto const& solid_phase = medium.phase(PhaseName::Solid);
    auto const& porosity_property = medium[PropertyType::porosity];
    auto const& solid_denisty_property = solid_phase[PropertyType::density];
 
    auto const phi =
        std::get<double>(porosity_property.value(variable_array, pos, t, dt));
    auto const rho_s = std::get<double>(
        solid_denisty_property.value(variable_array, pos, t, dt));
    VariableArray rho_d_and_phi;
    rho_d_and_phi.density = rho_s * (1.0 - phi);
    rho_d_and_phi.porosity = phi;
    double const p_cut = 1e-10;
 
    auto const S_L_max = water_content(p_cut, T, rho_d_and_phi) / phi;
    auto const p_cap_max = psi_max(T);
    if (p_cap <= p_cut)
    {
        return S_L_max;
    }
    if (p_cap >= p_cap_max)
    {
        return 0.0;
    }
    return water_content(p_cap, T, rho_d_and_phi) / phi;
}

References MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, psi_max(), MaterialPropertyLib::Property::scale_, MaterialPropertyLib::Solid, MaterialPropertyLib::VariableArray::temperature, and water_content().

◆ water_content() [1/2]

template<typename PressureType, typename TemperatureType>

promote< PressureType, TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::water_content	(	PressureType	p_cap,
		TemperatureType	T,
		VariableArray const &	variable_array ) const

Definition at line 19 of file SaturationLuMcCartney.cpp.

{
    return adsorptive_water_content(p_cap, T, variable_array) +
           capillary_water_content(p_cap, T, variable_array);
}

References adsorptive_water_content(), and capillary_water_content().

◆ water_content() [2/2]

template<typename PressureType, typename TemperatureType>

boost::math::differentiation::promote< PressureType, TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::water_content	(	PressureType	p_cap,
		TemperatureType	T,
		VariableArray const &	variable_array ) const

private

Referenced by d2Value(), dValue(), and value().

◆ zeta_s() [1/2]

template<typename TemperatureType>

promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::zeta_s ( TemperatureType T ) const

Definition at line 141 of file SaturationLuMcCartney.cpp.

{
    return (c(T) - 1) / c(T) * CEC(T) / nu_mr;
}

References MaterialPropertyLib::c, CEC(), and nu_mr.

◆ zeta_s() [2/2]

template<typename TemperatureType>

boost::math::differentiation::promote< TemperatureType > MaterialPropertyLib::SaturationLuMcCartney::zeta_s ( TemperatureType T ) const

private

Referenced by theta_a_max().

Member Data Documentation

◆ A

double const MaterialPropertyLib::SaturationLuMcCartney::A = 0.92

private

fitting parameter for cation exchange capacity

Definition at line 133 of file SaturationLuMcCartney.h.

Referenced by CEC().

◆ alpha_0_inv

double MaterialPropertyLib::SaturationLuMcCartney::alpha_0_inv

private

air entry suction at reference temperature Tr in Pa

Definition at line 159 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and alpha().

◆ B

double const MaterialPropertyLib::SaturationLuMcCartney::B = 0.9

private

fitting parameter for cation exchange capacity

Definition at line 134 of file SaturationLuMcCartney.h.

Referenced by CEC().

◆ bw

double MaterialPropertyLib::SaturationLuMcCartney::bw

private

Definition at line 150 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and theta_a_max().

◆ C1

double const MaterialPropertyLib::SaturationLuMcCartney::C1 = -0.00151

private

parameter for temperature correction

Definition at line 141 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and chi().

◆ CEC_max

double MaterialPropertyLib::SaturationLuMcCartney::CEC_max

private

maximum cation exchange capacity in meq/g

Definition at line 149 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and CEC().

◆ chi_r

double MaterialPropertyLib::SaturationLuMcCartney::chi_r

private

Definition at line 160 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and capillary_water_content().

◆ E1_minus_EL

double MaterialPropertyLib::SaturationLuMcCartney::E1_minus_EL

private

material parameters for adsorptive water content

energy difference between heat of adorption of first and second layer in J/mol

Definition at line 145 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and c().

◆ epsilon_s

double MaterialPropertyLib::SaturationLuMcCartney::epsilon_s

private

dielectic constant of soil particles

Definition at line 155 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and A_H().

◆ eta_alpha

double MaterialPropertyLib::SaturationLuMcCartney::eta_alpha

private

soil related coefficient

Definition at line 157 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and alpha().

◆ gibbs_free_energy_

const MaterialLib::Fluid::DimensionLessGibbsFreeEnergyRegion1 MaterialPropertyLib::SaturationLuMcCartney::gibbs_free_energy_

private

Definition at line 126 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and rho_w().

◆ M

double MaterialPropertyLib::SaturationLuMcCartney::M

private

exponent for adsorptive water content

Definition at line 151 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and adsorptive_water_content().

◆ material_

std::string const MaterialPropertyLib::SaturationLuMcCartney::material_

private

Definition at line 123 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney().

◆ N

double MaterialPropertyLib::SaturationLuMcCartney::N

private

exponent for capillary water content

Definition at line 161 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and capillary_water_content().

◆ n_s

double MaterialPropertyLib::SaturationLuMcCartney::n_s

private

refractive index of soil particles

Definition at line 156 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and A_H().

◆ n_w

double const MaterialPropertyLib::SaturationLuMcCartney::n_w = 1.33

private

diffraction index

Definition at line 142 of file SaturationLuMcCartney.h.

Referenced by A_H().

◆ nu_e

double const MaterialPropertyLib::SaturationLuMcCartney::nu_e

private

Initial value:

=

2.45e9

main absorption frequency of UV light in Hz

Definition at line 138 of file SaturationLuMcCartney.h.

Referenced by A_H().

◆ nu_mr

double MaterialPropertyLib::SaturationLuMcCartney::nu_mr

private

gravimetric moisture content when first layer is fully saturated

Definition at line 147 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and zeta_s().

◆ nu_w

double const MaterialPropertyLib::SaturationLuMcCartney::nu_w = 1.8e-5

private

<-— Lu and McCartney 2024 parameters ---->

molar volume of water in m^3/mol

Definition at line 132 of file SaturationLuMcCartney.h.

Referenced by psi_max().

◆ ref_p_

double MaterialPropertyLib::SaturationLuMcCartney::ref_p_ = 1.653e7

staticconstexprprivate

Definition at line 129 of file SaturationLuMcCartney.h.

Referenced by rho_w().

◆ ref_T_

double MaterialPropertyLib::SaturationLuMcCartney::ref_T_ = 1386

staticconstexprprivate

reference temperature in K.

Definition at line 128 of file SaturationLuMcCartney.h.

Referenced by rho_w().

◆ SSA

double MaterialPropertyLib::SaturationLuMcCartney::SSA

private

material parameters for capillary water content

specific surface area in m^2/g

Definition at line 154 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), and psi_c().

◆ T1

double const MaterialPropertyLib::SaturationLuMcCartney::T1 = 373.15

private

parameter for temperature effect of CEC in K

Definition at line 135 of file SaturationLuMcCartney.h.

Referenced by CEC().

◆ T2

double const MaterialPropertyLib::SaturationLuMcCartney::T2

private

Initial value:

=

1273.15

parameter for temperature effect of CEC in K

Definition at line 136 of file SaturationLuMcCartney.h.

Referenced by CEC().

◆ Tr

double const MaterialPropertyLib::SaturationLuMcCartney::Tr = 293.15

private

reference temperature in K

Definition at line 140 of file SaturationLuMcCartney.h.

Referenced by SaturationLuMcCartney(), alpha(), capillary_water_content(), and delta_h().

The documentation for this class was generated from the following files:

Detailed Description

Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ SaturationLuMcCartney()

Member Function Documentation

◆ A_H() [1/2]

◆ A_H() [2/2]

◆ adsorptive_water_content() [1/2]

◆ adsorptive_water_content() [2/2]

◆ alpha() [1/2]

◆ alpha() [2/2]

◆ c() [1/2]

◆ c() [2/2]

◆ capillary_water_content() [1/2]

◆ capillary_water_content() [2/2]

◆ CEC() [1/2]

◆ CEC() [2/2]

◆ checkScale()

◆ chi() [1/2]

◆ chi() [2/2]

◆ d2Value()

◆ delta_h() [1/2]

◆ delta_h() [2/2]

◆ dValue()

◆ epsilon_w() [1/2]

◆ epsilon_w() [2/2]

◆ psi_c() [1/2]

◆ psi_c() [2/2]

◆ psi_max() [1/2]

◆ psi_max() [2/2]

◆ rho_w() [1/2]

◆ rho_w() [2/2]

◆ theta_a_max() [1/2]

◆ theta_a_max() [2/2]

◆ value()

◆ water_content() [1/2]

◆ water_content() [2/2]

◆ zeta_s() [1/2]

◆ zeta_s() [2/2]

Member Data Documentation

◆ A

◆ alpha_0_inv

◆ B

◆ bw

◆ C1

◆ CEC_max

◆ chi_r

◆ E1_minus_EL

◆ epsilon_s

◆ eta_alpha

◆ gibbs_free_energy_

◆ M

◆ material_

◆ N

◆ n_s

◆ n_w

◆ nu_e

◆ nu_mr

◆ nu_w

◆ ref_p_

◆ ref_T_

◆ SSA

◆ T1

◆ T2

◆ Tr