ProcessLib::RichardsMechanics Namespace Reference

Namespaces
	detail

Classes
struct	MicroPorosityStateSpace
struct	MicroPorosityParameters
struct	IntegrationPointData
struct	LocalAssemblerInterface
class	LocalDataInitializer
struct	SecondaryData
class	RichardsMechanicsLocalAssembler
class	RichardsMechanicsProcess
struct	RichardsMechanicsProcessData

Functions
template<int DisplacementDim>
std::ostream &	operator<< (std::ostream &os, MicroPorosityStateSpace< DisplacementDim > const &state)
template<int DisplacementDim>
MicroPorosityStateSpace< DisplacementDim >	computeMicroPorosity (MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &I_2_C_el_inverse, double const rho_LR_m, double const mu_LR, MicroPorosityParameters const &micro_porosity_parameters, double const alpha_B, double const phi, double const p_L, double const p_L_m_prev, MaterialPropertyLib::VariableArray const &, double const S_L_m_prev, double const phi_m_prev, ParameterLib::SpatialPosition const pos, double const t, double const dt, MaterialPropertyLib::Property const &saturation_micro, MaterialPropertyLib::Property const &swelling_stress_rate)
template<int GlobalDim, template< typename, typename, typename, int > class LocalAssemblerImplementation, typename LocalAssemblerInterface , typename... ExtraCtorArgs>
void	createLocalAssemblers (const unsigned, std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, const unsigned shapefunction_order, std::vector< std::unique_ptr< LocalAssemblerInterface >> &local_assemblers, ExtraCtorArgs &&... extra_ctor_args)
void	checkMPLProperties (std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &media)
template<int DisplacementDim>
std::unique_ptr< Process >	createRichardsMechanicsProcess (std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &media)
template std::unique_ptr< Process >	createRichardsMechanicsProcess< 2 > (std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &media)
template std::unique_ptr< Process >	createRichardsMechanicsProcess< 3 > (std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &media)
template<int DisplacementDim, typename IPData >
void	updateSwellingStressAndVolumetricStrain (IPData &ip_data, MaterialPropertyLib::Medium const &medium, MaterialPropertyLib::Phase const &solid_phase, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_el, double const rho_LR, double const mu, std::optional< MicroPorosityParameters > micro_porosity_parameters, double const alpha, double const phi, double const p_cap_ip, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt)

Function Documentation

checkMPLProperties()

void ProcessLib::RichardsMechanics::checkMPLProperties

(

std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &

media

)

Definition at line 31 of file CreateRichardsMechanicsProcess.cpp.

{
    std::array const required_medium_properties = {
        MaterialPropertyLib::reference_temperature,
        MaterialPropertyLib::bishops_effective_stress,
        MaterialPropertyLib::relative_permeability,
        MaterialPropertyLib::saturation,
        MaterialPropertyLib::porosity,
        MaterialPropertyLib::biot_coefficient};
    std::array const required_liquid_properties = {
        MaterialPropertyLib::viscosity, MaterialPropertyLib::density};
    std::array const required_solid_properties = {MaterialPropertyLib::density};
 
    for (auto const& m : media)
    {
        checkRequiredProperties(*m.second, required_medium_properties);
        checkRequiredProperties(m.second->phase("AqueousLiquid"),
                                required_liquid_properties);
        checkRequiredProperties(m.second->phase("Solid"),
                                required_solid_properties);
    }
}

References MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::checkRequiredProperties(), MaterialPropertyLib::density, MaterialPropertyLib::porosity, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, and MaterialPropertyLib::viscosity.

Referenced by createRichardsMechanicsProcess().

computeMicroPorosity()

template<int DisplacementDim>

MicroPorosityStateSpace<DisplacementDim> ProcessLib::RichardsMechanics::computeMicroPorosity	(	MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &	I_2_C_el_inverse,
		double const	rho_LR_m,
		double const	mu_LR,
		MicroPorosityParameters const &	micro_porosity_parameters,
		double const	alpha_B,
		double const	phi,
		double const	p_L,
		double const	p_L_m_prev,
		MaterialPropertyLib::VariableArray const &	,
		double const	S_L_m_prev,
		double const	phi_m_prev,
		ParameterLib::SpatialPosition const	pos,
		double const	t,
		double const	dt,
		MaterialPropertyLib::Property const &	saturation_micro,
		MaterialPropertyLib::Property const &	swelling_stress_rate
	)

Definition at line 63 of file ComputeMicroPorosity.h.

{
    namespace MPL = MaterialPropertyLib;
    static constexpr int kelvin_vector_size =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    // phi_m, e_sw, p_L_m, sigma_sw
    static constexpr int nls_size = 1 + 1 + 1 + kelvin_vector_size;
 
    static constexpr int i_phi_m = 0;
    static constexpr int i_e_sw = 1;
    static constexpr int i_p_L_m = 2;
    static constexpr int i_sigma_sw = 3;
 
    using ResidualVectorType = Eigen::Matrix<double, nls_size, 1>;
    using JacobianMatrix =
        Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>;
 
    Eigen::FullPivLU<Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>>
        linear_solver;
 
    JacobianMatrix jacobian;
 
    // Agglomerated solution vector construction.
    ResidualVectorType solution = ResidualVectorType::Zero();
 
    if (p_L >= 0 && p_L_m_prev >= 0)
    {
        return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
                solution.template segment<kelvin_vector_size>(i_sigma_sw)};
    }
 
    double const alpha_bar =
        micro_porosity_parameters.mass_exchange_coefficient;
 
    auto const update_residual = [&](ResidualVectorType& residual)
    {
        double const delta_phi_m = solution[i_phi_m];
        double const delta_e_sw = solution[i_e_sw];
        auto const& delta_sigma_sw =
            solution.template segment<kelvin_vector_size>(i_sigma_sw);
        double const delta_p_L_m = solution[i_p_L_m];
 
        double const phi_m = phi_m_prev + delta_phi_m;
        double const p_L_m = p_L_m_prev + delta_p_L_m;
        MPL::VariableArray variables_prev;
        variables_prev[static_cast<int>(MPL::Variable::capillary_pressure)] =
            -p_L_m_prev;
        variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
            S_L_m_prev;
 
        MPL::VariableArray variables;
        variables[static_cast<int>(MPL::Variable::capillary_pressure)] = -p_L_m;
 
        double const S_L_m =
            saturation_micro.template value<double>(variables, pos, t, dt);
        variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L_m;
        double const delta_S_L_m = S_L_m - S_L_m_prev;
 
        auto const sigma_sw_dot =
            MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                swelling_stress_rate.template value<Eigen::Matrix3d>(
                    variables, variables_prev, pos, t, dt));
 
        residual[i_phi_m] = delta_phi_m - (alpha_B - phi) * delta_e_sw;
        residual[i_e_sw] = delta_e_sw + I_2_C_el_inverse.dot(delta_sigma_sw);
        residual.template segment<kelvin_vector_size>(i_sigma_sw).noalias() =
            delta_sigma_sw - sigma_sw_dot * dt;
 
        residual[i_p_L_m] =
            rho_LR_m *
                (phi_m * delta_S_L_m - (alpha_B - phi) * S_L_m * delta_e_sw) +
            phi_m * S_L_m * rho_LR_m * delta_e_sw -
            micro_porosity_parameters.mass_exchange_coefficient / mu_LR *
                (p_L - p_L_m) * dt;
    };
 
    auto const update_jacobian = [&](JacobianMatrix& jacobian)
    {
        jacobian = JacobianMatrix::Identity();
 
        double const delta_phi_m = solution[i_phi_m];
        double const delta_e_sw = solution[i_e_sw];
        double const delta_p_L_m = solution[i_p_L_m];
 
        double const phi_m = phi_m_prev + delta_phi_m;
        double const p_L_m = p_L_m_prev + delta_p_L_m;
        MPL::VariableArray variables_prev;
        variables_prev[static_cast<int>(MPL::Variable::capillary_pressure)] =
            -p_L_m_prev;
        MPL::VariableArray variables;
        variables[static_cast<int>(MPL::Variable::capillary_pressure)] = -p_L_m;
 
        double const S_L_m =
            saturation_micro.template value<double>(variables, pos, t, dt);
        variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
            S_L_m_prev;
        variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L_m;
        double const delta_S_L_m = S_L_m - S_L_m_prev;
 
        double const dS_L_m_dp_cap_m = saturation_micro.template dValue<double>(
            variables, MPL::Variable::capillary_pressure, pos, t, dt);
        auto const dsigma_sw_dS_L_m =
            MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                swelling_stress_rate.template dValue<Eigen::Matrix3d>(
                    variables, variables_prev, MPL::Variable::liquid_saturation,
                    pos, t, dt));
 
        jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi);
 
        jacobian.template block<1, kelvin_vector_size>(i_e_sw, i_sigma_sw) =
            I_2_C_el_inverse.transpose();
 
        jacobian.template block<kelvin_vector_size, 1>(i_sigma_sw, i_p_L_m) =
            -dsigma_sw_dS_L_m * dS_L_m_dp_cap_m;
 
        jacobian(i_p_L_m, i_phi_m) =
            rho_LR_m * (delta_S_L_m + S_L_m * delta_e_sw);
 
        jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi - phi_m);
 
        jacobian(i_p_L_m, i_p_L_m) =
            alpha_bar / mu_LR * dt -
            rho_LR_m * (phi_m - (alpha_B - phi - phi_m) * delta_e_sw) *
                dS_L_m_dp_cap_m;
    };
 
    auto const update_solution = [&](ResidualVectorType const& increment)
    { solution += increment; };
 
    auto newton_solver =
        NumLib::NewtonRaphson<decltype(linear_solver), JacobianMatrix,
                              decltype(update_jacobian), ResidualVectorType,
                              decltype(update_residual),
                              decltype(update_solution)>(
            linear_solver, update_jacobian, update_residual, update_solution,
            micro_porosity_parameters.nonlinear_solver_parameters);
 
    auto const success_iterations = newton_solver.solve(jacobian);
 
    if (!success_iterations)
    {
        OGS_FATAL(
            "Could not find solution for local double structure nonlinear "
            "problem.");
    }
 
    return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
            solution.template segment<kelvin_vector_size>(i_sigma_sw)};
}

References MaterialPropertyLib::capillary_pressure, MathLib::KelvinVector::kelvin_vector_dimensions(), ProcessLib::RichardsMechanics::MicroPorosityParameters::mass_exchange_coefficient, ProcessLib::RichardsMechanics::MicroPorosityParameters::nonlinear_solver_parameters, OGS_FATAL, MaterialPropertyLib::saturation_micro, NumLib::NewtonRaphson< LinearSolver, JacobianMatrix, JacobianMatrixUpdate, ResidualVector, ResidualUpdate, SolutionUpdate >::solve(), and MaterialPropertyLib::swelling_stress_rate.

createLocalAssemblers()

template<int GlobalDim, template< typename, typename, typename, int > class LocalAssemblerImplementation, typename LocalAssemblerInterface , typename... ExtraCtorArgs>

void ProcessLib::RichardsMechanics::createLocalAssemblers	(	const unsigned	,
		std::vector< MeshLib::Element * > const &	mesh_elements,
		NumLib::LocalToGlobalIndexMap const &	dof_table,
		const unsigned	shapefunction_order,
		std::vector< std::unique_ptr< LocalAssemblerInterface >> &	local_assemblers,
		ExtraCtorArgs &&...	extra_ctor_args
	)

Creates local assemblers for each element of the given mesh.

Template Parameters

LocalAssemblerImplementation	the individual local assembler type
LocalAssemblerInterface	the general local assembler interface
ExtraCtorArgs	types of additional constructor arguments. Those arguments will be passed to the constructor of LocalAssemblerImplementation.

The first two template parameters cannot be deduced from the arguments. Therefore they always have to be provided manually.

Definition at line 70 of file CreateLocalAssemblers.h.

{
    DBUG("Create local assemblers.");
 
    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
        dof_table, shapefunction_order, mesh_elements, local_assemblers,
        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
}

References DBUG().

Referenced by ProcessLib::RichardsMechanics::RichardsMechanicsProcess< DisplacementDim >::initializeConcreteProcess().

createRichardsMechanicsProcess()

template<int DisplacementDim>

std::unique_ptr< Process > ProcessLib::RichardsMechanics::createRichardsMechanicsProcess	(	std::string	name,
		MeshLib::Mesh &	mesh,
		std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&	jacobian_assembler,
		std::vector< ProcessVariable > const &	variables,
		std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &	parameters,
		std::optional< ParameterLib::CoordinateSystem > const &	local_coordinate_system,
		unsigned const	integration_order,
		BaseLib::ConfigTree const &	config,
		std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &	media
	)

Input File Parameter:

prj__processes__process__type

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__coupling_scheme

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__process_variables

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__process_variables__pressure

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__process_variables__displacement

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__specific_body_force

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__initial_stress

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__micro_porosity

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__micro_porosity__nonlinear_solver

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__micro_porosity__mass_exchange_coefficient

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__mass_lumping

Input File Parameter:

prj__processes__process__RICHARDS_MECHANICS__explicit_hm_coupling_in_unsaturated_zone

Definition at line 56 of file CreateRichardsMechanicsProcess.cpp.

{
    config.checkConfigParameter("type", "RICHARDS_MECHANICS");
    DBUG("Create RichardsMechanicsProcess.");
 
    auto const coupling_scheme =
        config.getConfigParameterOptional<std::string>("coupling_scheme");
    const bool use_monolithic_scheme =
        !(coupling_scheme && (*coupling_scheme == "staggered"));
 
    // Process variable.
 
    auto const pv_config = config.getConfigSubtree("process_variables");
 
    ProcessVariable* variable_p;
    ProcessVariable* variable_u;
    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
        process_variables;
    if (use_monolithic_scheme)  // monolithic scheme.
    {
        auto per_process_variables = findProcessVariables(
            variables, pv_config,
            {
             "pressure",
             "displacement"});
        variable_p = &per_process_variables[0].get();
        variable_u = &per_process_variables[1].get();
        process_variables.push_back(std::move(per_process_variables));
    }
    else  // staggered scheme.
    {
        using namespace std::string_literals;
        for (auto const& variable_name : {"pressure"s, "displacement"s})
        {
            auto per_process_variables =
                findProcessVariables(variables, pv_config, {variable_name});
            process_variables.push_back(std::move(per_process_variables));
        }
        variable_p = &process_variables[0][0].get();
        variable_u = &process_variables[1][0].get();
    }
 
    DBUG("Associate displacement with process variable '{:s}'.",
         variable_u->getName());
 
    if (variable_u->getNumberOfGlobalComponents() != DisplacementDim)
    {
        OGS_FATAL(
            "Number of components of the process variable '{:s}' is different "
            "from the displacement dimension: got {:d}, expected {:d}",
            variable_u->getName(),
            variable_u->getNumberOfGlobalComponents(),
            DisplacementDim);
    }
 
    DBUG("Associate pressure with process variable '{:s}'.",
         variable_p->getName());
    if (variable_p->getNumberOfGlobalComponents() != 1)
    {
        OGS_FATAL(
            "Pressure process variable '{:s}' is not a scalar variable but has "
            "{:d} components.",
            variable_p->getName(),
            variable_p->getNumberOfGlobalComponents());
    }
 
    auto solid_constitutive_relations =
        MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
            parameters, local_coordinate_system, config);
 
    // Specific body force
    Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
    {
        std::vector<double> const b =
            config.getConfigParameter<std::vector<double>>(
                "specific_body_force");
        if (b.size() != DisplacementDim)
        {
            OGS_FATAL(
                "The size of the specific body force vector does not match the "
                "displacement dimension. Vector size is {:d}, displacement "
                "dimension is {:d}",
                b.size(), DisplacementDim);
        }
 
        std::copy_n(b.data(), b.size(), specific_body_force.data());
    }
 
    auto media_map =
        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
    DBUG("Check the media properties of RichardsMechanics process ...");
    checkMPLProperties(media);
    DBUG("Media properties verified.");
 
    // Initial stress conditions
    auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
        config, "initial_stress", parameters,
        // Symmetric tensor size, 4 or 6, not a Kelvin vector.
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
        &mesh);
 
    std::optional<MicroPorosityParameters> micro_porosity_parameters;
    if (auto const micro_porosity_config =
        config.getConfigSubtreeOptional("micro_porosity"))
    {
        micro_porosity_parameters = MicroPorosityParameters{
            NumLib::createNewtonRaphsonSolverParameters(
                micro_porosity_config->getConfigSubtree("nonlinear_solver")),
            micro_porosity_config->getConfigParameter<double>(
                "mass_exchange_coefficient")};
    }
 
    auto const mass_lumping =
        config.getConfigParameter<bool>("mass_lumping", false);
 
    auto const explicit_hm_coupling_in_unsaturated_zone =
        config.getConfigParameter<bool>(
            "explicit_hm_coupling_in_unsaturated_zone", false);
 
    RichardsMechanicsProcessData<DisplacementDim> process_data{
        materialIDs(mesh),
        std::move(media_map),
        std::move(solid_constitutive_relations),
        initial_stress,
        specific_body_force,
        micro_porosity_parameters,
        mass_lumping,
        explicit_hm_coupling_in_unsaturated_zone};
 
    SecondaryVariableCollection secondary_variables;
 
    ProcessLib::createSecondaryVariables(config, secondary_variables);
 
    return std::make_unique<RichardsMechanicsProcess<DisplacementDim>>(
        std::move(name), mesh, std::move(jacobian_assembler), parameters,
        integration_order, std::move(process_variables),
        std::move(process_data), std::move(secondary_variables),
        use_monolithic_scheme);
}

References BaseLib::ConfigTree::checkConfigParameter(), checkMPLProperties(), MaterialPropertyLib::createMaterialSpatialDistributionMap(), NumLib::createNewtonRaphsonSolverParameters(), ProcessLib::createSecondaryVariables(), DBUG(), ProcessLib::findProcessVariables(), BaseLib::ConfigTree::getConfigParameter(), BaseLib::ConfigTree::getConfigParameterOptional(), BaseLib::ConfigTree::getConfigSubtree(), BaseLib::ConfigTree::getConfigSubtreeOptional(), ProcessLib::ProcessVariable::getName(), ProcessLib::ProcessVariable::getNumberOfGlobalComponents(), MathLib::KelvinVector::kelvin_vector_dimensions(), MeshLib::materialIDs(), MaterialPropertyLib::name, and OGS_FATAL.

createRichardsMechanicsProcess< 2 >()

template std::unique_ptr< Process > ProcessLib::RichardsMechanics::createRichardsMechanicsProcess< 2 >	(	std::string	name,
		MeshLib::Mesh &	mesh,
		std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&	jacobian_assembler,
		std::vector< ProcessVariable > const &	variables,
		std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &	parameters,
		std::optional< ParameterLib::CoordinateSystem > const &	local_coordinate_system,
		unsigned const	integration_order,
		BaseLib::ConfigTree const &	config,
		std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &	media
	)

Referenced by ProjectData::parseProcesses().

createRichardsMechanicsProcess< 3 >()

template std::unique_ptr< Process > ProcessLib::RichardsMechanics::createRichardsMechanicsProcess< 3 >	(	std::string	name,
		MeshLib::Mesh &	mesh,
		std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&	jacobian_assembler,
		std::vector< ProcessVariable > const &	variables,
		std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &	parameters,
		std::optional< ParameterLib::CoordinateSystem > const &	local_coordinate_system,
		unsigned const	integration_order,
		BaseLib::ConfigTree const &	config,
		std::map< int, std::shared_ptr< MaterialPropertyLib::Medium >> const &	media
	)

Referenced by ProjectData::parseProcesses().

operator<<()

template<int DisplacementDim>

std::ostream& ProcessLib::RichardsMechanics::operator<<	(	std::ostream &	os,
		MicroPorosityStateSpace< DisplacementDim > const &	state
	)

Definition at line 43 of file ComputeMicroPorosity.h.

{
    return os << "phi_m: " << state.phi_m << ", "
              << "e_sw: " << state.e_sw << ", "
              << "p_L_m: " << state.p_L_m << ", "
              << "sigma_sw: " << state.sigma_sw.transpose();
}

References ProcessLib::RichardsMechanics::MicroPorosityStateSpace< DisplacementDim >::e_sw, ProcessLib::RichardsMechanics::MicroPorosityStateSpace< DisplacementDim >::p_L_m, ProcessLib::RichardsMechanics::MicroPorosityStateSpace< DisplacementDim >::phi_m, and ProcessLib::RichardsMechanics::MicroPorosityStateSpace< DisplacementDim >::sigma_sw.

updateSwellingStressAndVolumetricStrain()

template<int DisplacementDim, typename IPData >

void ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain	(	IPData &	ip_data,
		MaterialPropertyLib::Medium const &	medium,
		MaterialPropertyLib::Phase const &	solid_phase,
		MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &	C_el,
		double const	rho_LR,
		double const	mu,
		std::optional< MicroPorosityParameters >	micro_porosity_parameters,
		double const	alpha,
		double const	phi,
		double const	p_cap_ip,
		MPL::VariableArray &	variables,
		MPL::VariableArray &	variables_prev,
		ParameterLib::SpatialPosition const &	x_position,
		double const	t,
		double const	dt
	)

Definition at line 32 of file RichardsMechanicsFEM-impl.h.

{
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    auto& sigma_sw = ip_data.sigma_sw;
    auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
    if (!medium.hasProperty(MPL::PropertyType::saturation_micro))
    {
        // If there is swelling, compute it. Update volumetric strain rate,
        // s.t. it corresponds to the mechanical part only.
        sigma_sw = sigma_sw_prev;
        if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
        {
            using DimMatrix = Eigen::Matrix<double, 3, 3>;
            auto const sigma_sw_dot =
                MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                    solid_phase
                        .property(MPL::PropertyType::swelling_stress_rate)
                        .template value<DimMatrix>(variables, variables_prev,
                                                   x_position, t, dt));
            sigma_sw += sigma_sw_dot * dt;
 
            // !!! Misusing volumetric strain for mechanical volumetric
            // strain just to update the transport porosity !!!
            std::get<double>(variables[static_cast<int>(
                MPL::Variable::volumetric_strain)]) +=
                identity2.transpose() * C_el.inverse() * sigma_sw;
            std::get<double>(variables_prev[static_cast<int>(
                MPL::Variable::volumetric_strain)]) +=
                identity2.transpose() * C_el.inverse() * sigma_sw_prev;
        }
    }
 
    // TODO (naumov) saturation_micro must be always defined together with
    // the micro_porosity_parameters.
    if (medium.hasProperty(MPL::PropertyType::saturation_micro))
    {
        double const phi_M_prev = ip_data.transport_porosity_prev;
        double const phi_prev = ip_data.porosity_prev;
        double const phi_m_prev = phi_prev - phi_M_prev;
        double const p_L_m_prev = ip_data.liquid_pressure_m_prev;
 
        auto const S_L_m_prev = ip_data.saturation_m_prev;
 
        auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] =
            computeMicroPorosity<DisplacementDim>(
                identity2.transpose() * C_el.inverse(), rho_LR, mu,
                *micro_porosity_parameters, alpha, phi, -p_cap_ip, p_L_m_prev,
                variables_prev, S_L_m_prev, phi_m_prev, x_position, t, dt,
                medium.property(MPL::PropertyType::saturation_micro),
                solid_phase.property(MPL::PropertyType::swelling_stress_rate));
 
        auto& phi_M = ip_data.transport_porosity;
        phi_M = phi - (phi_m_prev + delta_phi_m);
        variables_prev[static_cast<int>(MPL::Variable::transport_porosity)] =
            phi_M_prev;
        variables[static_cast<int>(MPL::Variable::transport_porosity)] = phi_M;
 
        auto& p_L_m = ip_data.liquid_pressure_m;
        p_L_m = p_L_m_prev + delta_p_L_m;
        {  // Update micro saturation.
            MPL::VariableArray variables_prev;
            variables_prev[static_cast<int>(
                MPL::Variable::capillary_pressure)] = -p_L_m_prev;
            MPL::VariableArray variables;
            variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
                -p_L_m;
 
            ip_data.saturation_m =
                medium.property(MPL::PropertyType::saturation_micro)
                    .template value<double>(variables, x_position, t, dt);
        }
        sigma_sw = sigma_sw_prev + delta_sigma_sw;
    }
}

References MaterialPropertyLib::alpha, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::Medium::hasProperty(), MaterialPropertyLib::Phase::hasProperty(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::Medium::property(), MaterialPropertyLib::Phase::property(), MaterialPropertyLib::saturation_micro, MaterialPropertyLib::swelling_stress_rate, and MaterialPropertyLib::transport_porosity.