diff --git a/doc/documentation/src/_assets/mohr_coulomb_yield_curve.png b/doc/documentation/src/_assets/mohr_coulomb_yield_curve.png new file mode 100644 index 00000000000..fd214476520 Binary files /dev/null and b/doc/documentation/src/_assets/mohr_coulomb_yield_curve.png differ diff --git a/doc/documentation/src/developer_guide/development_parts.rst b/doc/documentation/src/developer_guide/development_parts.rst index 1d42c739672..d91e762fdc9 100644 --- a/doc/documentation/src/developer_guide/development_parts.rst +++ b/doc/documentation/src/developer_guide/development_parts.rst @@ -6,3 +6,4 @@ Development of specific code blocks .. toctree:: materialdevelopment + mohr_coulomb_material diff --git a/doc/documentation/src/developer_guide/mohr_coulomb_material.rst b/doc/documentation/src/developer_guide/mohr_coulomb_material.rst new file mode 100644 index 00000000000..0b689e7fc9f --- /dev/null +++ b/doc/documentation/src/developer_guide/mohr_coulomb_material.rst @@ -0,0 +1,81 @@ +Mohr-Coulomb material +===================== + +``MAT_Struct_MohrCoulomb`` is a finite-strain isotropic elastoplastic material. It uses the +principal elastic logarithmic (Hencky) strains and their work-conjugate Kirchhoff stresses. +The multiplicative history variable is the inverse plastic right Cauchy-Green tensor. + +Yield surface and flow rule +--------------------------- + +For ordered principal Kirchhoff stresses +:math:`\tau_1 \geq \tau_2 \geq \tau_3`, the active Mohr-Coulomb plane is + +.. math:: + + f = k\tau_1-\tau_3-2\sqrt{k}\,c(\bar{\varepsilon}^{p}) \leq 0, + \qquad + k=\frac{1+\sin\phi}{1-\sin\phi}. + +The plastic potential has the same form with +:math:`m=(1+\sin\psi)/(1-\sin\psi)`, where :math:`\psi` is the dilatancy angle. +Setting ``DILATANCY_ANGLE`` equal to ``FRICTION_ANGLE`` gives associative plasticity. +The return mapping works directly with active planes in principal stress space and treats a +smooth face, both triaxial edges, and the hydrostatic apex separately. + +The cohesion hardening law combines linear and Voce terms: + +.. math:: + + c(\bar{\varepsilon}^{p}) = + c_0 + H\bar{\varepsilon}^{p} + + (c_\infty-c_0)\left(1-\exp(-b\bar{\varepsilon}^{p})\right). + +Set ``SATURATION_HARDENING`` to zero to disable the Voce term. + +Input +----- + +.. code-block:: yaml + + MAT_Struct_MohrCoulomb: + YOUNG_MODULUS: 1.0e5 + POISSON_RATIO: 0.3 + DENSITY: 2200 + COHESION: 100 + FRICTION_ANGLE: 0.5235987755982 + DILATANCY_ANGLE: 0.4363323129986 + LINEAR_HARDENING: 0 + SATURATION_HARDENING: 500 + HARDENING_EXP: 10 + TOLERANCE: 1.0e-10 + MAX_ITERATIONS: 100 + +The material requires ``KINEM: nonlinear``. Angles are specified in radians. The dilatancy +angle must be strictly positive and must not exceed the friction angle; a zero dilatancy angle +does not provide an admissible hydrostatic apex flow direction. + +Output +------ + +The material registers ``plastic_strain``, ``accumulated_plastic_strain``, +``accumulated_plastic_volumetric_strain``, and ``local_dissipated_energy`` as Gauss-point +output quantities. + +Verification +------------ + +The calculated triaxial-compression and triaxial-tension meridians coincide with the analytical +Mohr-Coulomb yield surface: + +.. image:: ../_assets/mohr_coulomb_yield_curve.png + :alt: Analytical and calculated Mohr-Coulomb yield meridians + :width: 80% + +The constitutive tests also cover elastic response, associative and non-associated flow, +linear and Voce hardening, both edges, the apex, the consistent smooth-face tangent, output +variables, and serialization. + +The return algorithm follows Clausen, Damkilde, and Andersen, +*Efficient return algorithms for associated plasticity with multiple yield planes*, +Computers & Structures 85 (2007), doi:10.1016/j.compstruc.2007.04.002. diff --git a/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.cpp b/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.cpp index a342a8a6e74..6b17a2021ab 100644 --- a/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.cpp +++ b/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.cpp @@ -101,6 +101,8 @@ std::string_view Core::Materials::to_string(Core::Materials::MaterialType materi return "MAT_Struct_ThermoStVenantK"; case m_pldruckprag: return "MAT_Struct_DruckerPrager"; + case m_plmohrcoulomb: + return "MAT_Struct_MohrCoulomb"; case m_thermopllinelast: return "MAT_Struct_ThermoPlasticLinElast"; case m_superelast: diff --git a/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.hpp b/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.hpp index e3805640d7f..680a6b5c91f 100644 --- a/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.hpp +++ b/src/core/legacy_enum_definitions/4C_legacy_enum_definitions_materials.hpp @@ -133,6 +133,7 @@ namespace Core::Materials m_particle_wall_dem, ///< particle wall material for DEM m_particle_pd, ///< particle material for PD m_pldruckprag, ///< Plastic linear elastic St.Venant Kirchhoff / Drucker Prager plasticity + m_plmohrcoulomb, ///< Finite-strain logarithmic Mohr-Coulomb plasticity m_plgtn, ///< Plastic linear elastic St.Venant Kirchhoff / GTN plasticity m_plelasthyper, ///< general hyperelastic material for finite strain von-Mises plasticity ///< using a semi-smooth Newton strategy (only in combination with such diff --git a/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp b/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp index 95cf63b34fa..3e3e5b4d493 100644 --- a/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp +++ b/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp @@ -20,12 +20,14 @@ #include "4C_mat_inelastic_defgrad_factors_service.hpp" #include "4C_mat_micromaterial.hpp" #include "4C_mat_plasticdruckerprager.hpp" +#include "4C_mat_plasticmohrcoulomb.hpp" #include "4C_mat_scatra_growth_remodel.hpp" #include "4C_mat_scatra_nonlocal_stimulus.hpp" #include "4C_porofluid_pressure_based_elast_scatra_input.hpp" #include "4C_structure_new_input.hpp" #include +#include #include #include @@ -58,6 +60,55 @@ std::unordered_map Global::v {.description = "Newtonian fluid"}); } + /*----------------------------------------------------------------------*/ + // Finite-strain logarithmic Mohr-Coulomb plasticity + { + using namespace Core::IO::InputSpecBuilders::Validators; + constexpr double half_pi = std::numbers::pi / 2.0; + + known_materials[Core::Materials::m_plmohrcoulomb] = group("MAT_Struct_MohrCoulomb", + { + parameter("YOUNG_MODULUS", + {.description = "Young's modulus", .validator = positive()}), + parameter("POISSON_RATIO", + {.description = "Poisson's ratio (must be in (-1, 0.5) for stability)", + .validator = in_range(excl(-1.0), excl(0.5))}), + parameter("DENSITY", + {.description = "Mass density", .validator = positive_or_zero()}), + parameter( + "COHESION", {.description = "Initial cohesion", .validator = positive()}), + parameter( + "FRICTION_ANGLE", {.description = "Friction angle in radians", + .validator = in_range(excl(0.0), excl(half_pi))}), + parameter("DILATANCY_ANGLE", + {.description = "Dilatancy angle in radians " + "(equal to friction angle for associated plasticity)", + .validator = in_range(excl(0.0), excl(half_pi))}), + parameter( + "LINEAR_HARDENING", {.description = "Linear isotropic cohesion-hardening modulus", + .default_value = 0.0, + .validator = positive_or_zero()}), + parameter("SATURATION_HARDENING", + {.description = "Asymptotic cohesion for Voce hardening; zero disables it", + .default_value = 0.0, + .validator = positive_or_zero()}), + parameter("HARDENING_EXP", {.description = "Voce hardening exponent", + .default_value = 0.0, + .validator = positive_or_zero()}), + parameter("TOLERANCE", {.description = "Local return-mapping tolerance", + .default_value = 1.0e-10, + .validator = positive()}), + parameter( + "MAX_ITERATIONS", {.description = "Maximum local return-mapping iterations", + .default_value = 50, + .validator = positive()}), + }, + {.description = + "Finite-strain Mohr-Coulomb plasticity based on elastic logarithmic strain and " + "Kirchhoff stress, with non-associated flow and linear plus Voce isotropic " + "cohesion hardening."}); + } + /*----------------------------------------------------------------------*/ // Weakly compressible fluid according to Murnaghan-Tait { diff --git a/src/mat/4C_mat_material_factory.cpp b/src/mat/4C_mat_material_factory.cpp index c9b20a9bdf8..868717ae691 100644 --- a/src/mat/4C_mat_material_factory.cpp +++ b/src/mat/4C_mat_material_factory.cpp @@ -114,6 +114,7 @@ #include "4C_mat_plasticelasthyper.hpp" #include "4C_mat_plasticgtn.hpp" #include "4C_mat_plasticlinelast.hpp" +#include "4C_mat_plasticmohrcoulomb.hpp" #include "4C_mat_plasticnlnlogneohooke.hpp" #include "4C_mat_poro_density_law.hpp" #include "4C_mat_poro_law.hpp" @@ -249,6 +250,10 @@ std::unique_ptr Mat::make_parameter( { return make_parameter_impl(id, type, input_data); } + case Core::Materials::m_plmohrcoulomb: + { + return make_parameter_impl(id, type, input_data); + } case Core::Materials::m_plgtn: { return make_parameter_impl(id, type, input_data); diff --git a/src/mat/4C_mat_plasticmohrcoulomb.cpp b/src/mat/4C_mat_plasticmohrcoulomb.cpp new file mode 100644 index 00000000000..db27a2426a9 --- /dev/null +++ b/src/mat/4C_mat_plasticmohrcoulomb.cpp @@ -0,0 +1,723 @@ +// This file is part of 4C multiphysics licensed under the +// GNU Lesser General Public License v3.0 or later. +// +// See the LICENSE.md file in the top-level for license information. +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include "4C_mat_plasticmohrcoulomb.hpp" + +#include "4C_comm_pack_helpers.hpp" +#include "4C_global_data.hpp" +#include "4C_linalg_fixedsizematrix_tensor_products.hpp" +#include "4C_linalg_symmetric_tensor_eigen.hpp" +#include "4C_linalg_tensor.hpp" +#include "4C_linalg_tensor_generators.hpp" +#include "4C_linalg_utils_densematrix_inverse.hpp" +#include "4C_mat_par_bundle.hpp" + +#include +#include +#include +#include + +FOUR_C_NAMESPACE_OPEN + +namespace +{ + /** + * Result of the material-local active-set return mapping in ordered principal-stress space. + * + * This type and the following helper functions remain in the anonymous namespace because they + * implement details of PlasticMohrCoulomb rather than reusable material or linear-algebra + * operations. + */ + struct ActiveSetResult + { + bool valid = false; + bool numerical_failure = false; + Core::LinAlg::Tensor stress{}; + Core::LinAlg::Tensor plastic_strain_increment{}; + Core::LinAlg::Tensor tangent{}; + double accumulated_plastic_strain = 0.0; + }; + + double sum_components(const Core::LinAlg::Tensor& vector) + { + return vector(0) + vector(1) + vector(2); + } + + Core::LinAlg::Tensor deviator(const Core::LinAlg::Tensor& vector) + { + const double mean = sum_components(vector) / 3.0; + return Core::LinAlg::Tensor{{vector(0) - mean, vector(1) - mean, vector(2) - mean}}; + } + + double equivalent_plastic_strain(const Core::LinAlg::Tensor& plastic_strain_increment) + { + const auto deviatoric_increment = deviator(plastic_strain_increment); + return std::sqrt(2.0 / 3.0) * Core::LinAlg::norm2(deviatoric_increment); + } + + Core::LinAlg::Tensor yield_normal( + int positive_index, int negative_index, double factor) + { + Core::LinAlg::Tensor normal{}; + normal(positive_index) = factor; + normal(negative_index) = -1.0; + return normal; + } + + double max_yield_function( + const Core::LinAlg::Tensor& stress, double friction_stress_ratio, double q) + { + double maximum = -std::numeric_limits::infinity(); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + if (i != j) maximum = std::max(maximum, friction_stress_ratio * stress(i) - stress(j) - q); + return maximum; + } + + double cohesion(const Mat::PAR::PlasticMohrCoulomb& parameters, double plastic_strain) + { + double value = parameters.cohesion_ + parameters.linear_hardening_ * plastic_strain; + if (parameters.saturation_hardening_ > 0.0) + value += (parameters.saturation_hardening_ - parameters.cohesion_) * + (1.0 - std::exp(-parameters.hardening_exponent_ * plastic_strain)); + return value; + } + + double cohesion_derivative(const Mat::PAR::PlasticMohrCoulomb& parameters, double plastic_strain) + { + double value = parameters.linear_hardening_; + if (parameters.saturation_hardening_ > 0.0) + value += (parameters.saturation_hardening_ - parameters.cohesion_) * + parameters.hardening_exponent_ * + std::exp(-parameters.hardening_exponent_ * plastic_strain); + return value; + } + + Core::LinAlg::Tensor elastic_principal_tangent( + double shear_modulus, double lame_parameter) + { + Core::LinAlg::Tensor tangent{}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + tangent(i, j) = lame_parameter + (i == j ? 2.0 * shear_modulus : 0.0); + return tangent; + } + + bool belongs_to_apex_flow_cone(const Core::LinAlg::Tensor& plastic_strain_increment, + double dilatancy_stress_ratio, double tolerance) + { + if (Core::LinAlg::norm2(plastic_strain_increment) <= tolerance) return true; + + std::array, 6> flow_normals{}; + int normal_index = 0; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + if (i != j) flow_normals[normal_index++] = yield_normal(i, j, dilatancy_stress_ratio); + + for (int first = 0; first < 4; ++first) + for (int second = first + 1; second < 5; ++second) + for (int third = second + 1; third < 6; ++third) + { + Core::LinAlg::Tensor flow_basis{}; + for (int row = 0; row < 3; ++row) + { + flow_basis(row, 0) = flow_normals[first](row); + flow_basis(row, 1) = flow_normals[second](row); + flow_basis(row, 2) = flow_normals[third](row); + } + const double determinant = Core::LinAlg::det(flow_basis); + if (std::abs(determinant) <= std::numeric_limits::epsilon()) continue; + + const auto determinant_with_column = [&](int column) + { + auto matrix = flow_basis; + for (int row = 0; row < 3; ++row) matrix(row, column) = plastic_strain_increment(row); + return Core::LinAlg::det(matrix); + }; + + const Core::LinAlg::Tensor multipliers{ + {determinant_with_column(0) / determinant, determinant_with_column(1) / determinant, + determinant_with_column(2) / determinant}}; + if (multipliers(0) < -tolerance || multipliers(1) < -tolerance || + multipliers(2) < -tolerance) + continue; + + const auto reconstructed = multipliers(0) * flow_normals[first] + + multipliers(1) * flow_normals[second] + + multipliers(2) * flow_normals[third]; + const auto error = reconstructed - plastic_strain_increment; + if (Core::LinAlg::norm2(error) <= + tolerance * std::max(1.0, Core::LinAlg::norm2(plastic_strain_increment))) + return true; + } + return false; + } + + ActiveSetResult return_to_active_set(const Core::LinAlg::Tensor& trial_stress, + double accumulated_plastic_strain_last, + const Core::LinAlg::Tensor& elastic_tangent, double friction_stress_ratio, + double dilatancy_stress_ratio, const std::array, 2>& planes, + int number_of_planes, const Mat::PAR::PlasticMohrCoulomb& parameters) + { + std::array, 2> yield_normals{}; + std::array, 2> flow_normals{}; + std::array, 2> elastic_yield_normals{}; + std::array, 2> elastic_flow_normals{}; + for (int alpha = 0; alpha < number_of_planes; ++alpha) + { + yield_normals[alpha] = + yield_normal(planes[alpha].first, planes[alpha].second, friction_stress_ratio); + flow_normals[alpha] = + yield_normal(planes[alpha].first, planes[alpha].second, dilatancy_stress_ratio); + elastic_yield_normals[alpha] = elastic_tangent * yield_normals[alpha]; + elastic_flow_normals[alpha] = elastic_tangent * flow_normals[alpha]; + } + + std::array, 2> active_matrix{}; + for (int alpha = 0; alpha < number_of_planes; ++alpha) + for (int beta = 0; beta < number_of_planes; ++beta) + active_matrix[alpha][beta] = + Core::LinAlg::dot(yield_normals[alpha], elastic_flow_normals[beta]); + + const double determinant = + number_of_planes == 1 + ? active_matrix[0][0] + : active_matrix[0][0] * active_matrix[1][1] - active_matrix[0][1] * active_matrix[1][0]; + if (std::abs(determinant) <= std::numeric_limits::epsilon()) + return {.numerical_failure = true}; + + auto solve_active_matrix = [&](const std::array& right_hand_side) + { + std::array result{}; + if (number_of_planes == 1) + { + result[0] = right_hand_side[0] / active_matrix[0][0]; + } + else + { + result[0] = + (right_hand_side[0] * active_matrix[1][1] - active_matrix[0][1] * right_hand_side[1]) / + determinant; + result[1] = + (active_matrix[0][0] * right_hand_side[1] - right_hand_side[0] * active_matrix[1][0]) / + determinant; + } + return result; + }; + + const double q_factor = 2.0 * std::sqrt(friction_stress_ratio); + auto state_at_accumulated_strain = [&](double accumulated_plastic_strain) + { + const double q = q_factor * cohesion(parameters, accumulated_plastic_strain); + std::array right_hand_side{}; + for (int alpha = 0; alpha < number_of_planes; ++alpha) + right_hand_side[alpha] = Core::LinAlg::dot(yield_normals[alpha], trial_stress) - q; + const std::array plastic_multipliers = solve_active_matrix(right_hand_side); + + Core::LinAlg::Tensor plastic_strain_increment{}; + for (int alpha = 0; alpha < number_of_planes; ++alpha) + plastic_strain_increment += plastic_multipliers[alpha] * flow_normals[alpha]; + + return std::make_tuple(plastic_multipliers, plastic_strain_increment, + accumulated_plastic_strain - accumulated_plastic_strain_last - + equivalent_plastic_strain(plastic_strain_increment)); + }; + + const auto [multipliers_at_lower, plastic_increment_at_lower, residual_at_lower] = + state_at_accumulated_strain(accumulated_plastic_strain_last); + if (residual_at_lower > parameters.tolerance_) return {}; + + double lower = accumulated_plastic_strain_last; + double upper = lower + std::max(2.0 * equivalent_plastic_strain(plastic_increment_at_lower), + 10.0 * parameters.tolerance_); + double residual_upper = std::get<2>(state_at_accumulated_strain(upper)); + for (int iteration = 0; residual_upper < 0.0 && iteration < parameters.max_iterations_; + ++iteration) + { + upper = lower + 2.0 * (upper - lower); + residual_upper = std::get<2>(state_at_accumulated_strain(upper)); + } + if (residual_upper < 0.0) return {.numerical_failure = true}; + + double accumulated_plastic_strain = + std::clamp(lower + equivalent_plastic_strain(plastic_increment_at_lower), lower, upper); + for (int iteration = 0; iteration < parameters.max_iterations_; ++iteration) + { + const auto [multipliers, plastic_strain_increment, residual] = + state_at_accumulated_strain(accumulated_plastic_strain); + if (std::abs(residual) <= parameters.tolerance_) break; + + if (residual > 0.0) + upper = accumulated_plastic_strain; + else + lower = accumulated_plastic_strain; + + const double plastic_increment_norm = equivalent_plastic_strain(plastic_strain_increment); + double residual_derivative = 1.0; + if (plastic_increment_norm > parameters.tolerance_) + { + const auto plastic_increment_deviator = deviator(plastic_strain_increment); + const std::array inverse_times_ones = + solve_active_matrix({1.0, number_of_planes == 2 ? 1.0 : 0.0}); + Core::LinAlg::Tensor flow_correction{}; + for (int alpha = 0; alpha < number_of_planes; ++alpha) + flow_correction += inverse_times_ones[alpha] * flow_normals[alpha]; + const double q_derivative = + q_factor * cohesion_derivative(parameters, accumulated_plastic_strain); + residual_derivative += + q_derivative * 2.0 / 3.0 * + Core::LinAlg::dot(plastic_increment_deviator, deviator(flow_correction)) / + plastic_increment_norm; + } + + double candidate = accumulated_plastic_strain - residual / residual_derivative; + if (!(candidate > lower && candidate < upper) || !std::isfinite(candidate)) + candidate = 0.5 * (lower + upper); + accumulated_plastic_strain = candidate; + } + + const auto [plastic_multipliers, plastic_strain_increment, residual] = + state_at_accumulated_strain(accumulated_plastic_strain); + if (std::abs(residual) > 10.0 * parameters.tolerance_) return {.numerical_failure = true}; + for (int alpha = 0; alpha < number_of_planes; ++alpha) + if (plastic_multipliers[alpha] < -10.0 * parameters.tolerance_) return {}; + + const auto updated_stress = trial_stress - elastic_tangent * plastic_strain_increment; + const double q = q_factor * cohesion(parameters, accumulated_plastic_strain); + if (max_yield_function(updated_stress, friction_stress_ratio, q) > + 100.0 * parameters.tolerance_) + return {}; + if (updated_stress(0) + parameters.tolerance_ < updated_stress(1) || + updated_stress(1) + parameters.tolerance_ < updated_stress(2)) + return {}; + + Core::LinAlg::Tensor algorithmic_tangent = elastic_tangent; + std::array, 2> hardening_matrix = active_matrix; + const double plastic_increment_norm = equivalent_plastic_strain(plastic_strain_increment); + if (plastic_increment_norm > parameters.tolerance_) + { + const Core::LinAlg::Tensor hardening_direction = + 2.0 / (3.0 * plastic_increment_norm) * deviator(plastic_strain_increment); + const double q_derivative = + q_factor * cohesion_derivative(parameters, accumulated_plastic_strain); + for (int beta = 0; beta < number_of_planes; ++beta) + { + const double hardening_component = + Core::LinAlg::dot(hardening_direction, flow_normals[beta]); + for (int alpha = 0; alpha < number_of_planes; ++alpha) + hardening_matrix[alpha][beta] += q_derivative * hardening_component; + } + } + + const double hardening_determinant = number_of_planes == 1 + ? hardening_matrix[0][0] + : hardening_matrix[0][0] * hardening_matrix[1][1] - + hardening_matrix[0][1] * hardening_matrix[1][0]; + if (std::abs(hardening_determinant) <= std::numeric_limits::epsilon()) + return {.numerical_failure = true}; + + std::array, 2> inverse_hardening_matrix{}; + if (number_of_planes == 1) + { + inverse_hardening_matrix[0][0] = 1.0 / hardening_matrix[0][0]; + } + else + { + inverse_hardening_matrix[0][0] = hardening_matrix[1][1] / hardening_determinant; + inverse_hardening_matrix[0][1] = -hardening_matrix[0][1] / hardening_determinant; + inverse_hardening_matrix[1][0] = -hardening_matrix[1][0] / hardening_determinant; + inverse_hardening_matrix[1][1] = hardening_matrix[0][0] / hardening_determinant; + } + + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + for (int alpha = 0; alpha < number_of_planes; ++alpha) + for (int beta = 0; beta < number_of_planes; ++beta) + algorithmic_tangent(i, j) -= elastic_flow_normals[alpha](i) * + inverse_hardening_matrix[alpha][beta] * + elastic_yield_normals[beta](j); + + return {.valid = true, + .stress = updated_stress, + .plastic_strain_increment = plastic_strain_increment, + .tangent = algorithmic_tangent, + .accumulated_plastic_strain = accumulated_plastic_strain}; + } +} // namespace + +Mat::PAR::PlasticMohrCoulomb::PlasticMohrCoulomb(const Core::Mat::PAR::Parameter::Data& matdata) + : Parameter(matdata), + youngs_(matdata.parameters.get("YOUNG_MODULUS")), + poisson_ratio_(matdata.parameters.get("POISSON_RATIO")), + density_(matdata.parameters.get("DENSITY")), + cohesion_(matdata.parameters.get("COHESION")), + friction_angle_(matdata.parameters.get("FRICTION_ANGLE")), + dilatancy_angle_(matdata.parameters.get("DILATANCY_ANGLE")), + linear_hardening_(matdata.parameters.get("LINEAR_HARDENING")), + saturation_hardening_(matdata.parameters.get("SATURATION_HARDENING")), + hardening_exponent_(matdata.parameters.get("HARDENING_EXP")), + tolerance_(matdata.parameters.get("TOLERANCE")), + max_iterations_(matdata.parameters.get("MAX_ITERATIONS")) +{ + if (dilatancy_angle_ <= 0.0) + FOUR_C_THROW("DILATANCY_ANGLE must be strictly positive for an admissible apex return."); + if (dilatancy_angle_ > friction_angle_) + FOUR_C_THROW("DILATANCY_ANGLE must not exceed FRICTION_ANGLE."); + if (saturation_hardening_ > 0.0 && saturation_hardening_ < cohesion_) + FOUR_C_THROW("SATURATION_HARDENING must be zero or at least COHESION."); + if (saturation_hardening_ > 0.0 && hardening_exponent_ <= 0.0) + FOUR_C_THROW("HARDENING_EXP must be positive when saturation hardening is enabled."); +} + +std::shared_ptr Mat::PAR::PlasticMohrCoulomb::create_material() +{ + return std::make_shared(this); +} + +Mat::PlasticMohrCoulombType Mat::PlasticMohrCoulombType::instance_; + +Core::Communication::ParObject* Mat::PlasticMohrCoulombType::create( + Core::Communication::UnpackBuffer& buffer) +{ + auto* material = new Mat::PlasticMohrCoulomb(); + material->unpack(buffer); + return material; +} + +Mat::PlasticMohrCoulomb::PlasticMohrCoulomb() : params_(nullptr) {} + +Mat::PlasticMohrCoulomb::PlasticMohrCoulomb(Mat::PAR::PlasticMohrCoulomb* params) : params_(params) +{ +} + +void Mat::PlasticMohrCoulomb::pack(Core::Communication::PackBuffer& data) const +{ + add_to_pack(data, unique_par_object_id()); + add_to_pack(data, params_ != nullptr ? params_->id() : -1); + add_to_pack(data, inv_plastic_rcg_last_); + add_to_pack(data, inv_plastic_rcg_current_); + add_to_pack(data, accumulated_plastic_strain_last_); + add_to_pack(data, accumulated_plastic_strain_current_); + add_to_pack(data, accumulated_plastic_volumetric_strain_last_); + add_to_pack(data, accumulated_plastic_volumetric_strain_current_); + add_to_pack(data, dissipated_energy_last_); + add_to_pack(data, dissipated_energy_current_); +} + +void Mat::PlasticMohrCoulomb::unpack(Core::Communication::UnpackBuffer& buffer) +{ + Core::Communication::extract_and_assert_id(buffer, unique_par_object_id()); + int material_id; + extract_from_pack(buffer, material_id); + params_ = nullptr; + if (Global::Problem::instance()->materials() != nullptr && + Global::Problem::instance()->materials()->num() != 0) + { + const int problem_instance = Global::Problem::instance()->materials()->get_read_from_problem(); + Core::Mat::PAR::Parameter* material = + Global::Problem::instance(problem_instance)->materials()->parameter_by_id(material_id); + if (material->type() != material_type()) + FOUR_C_THROW("Packed material type {} does not match Mohr-Coulomb type {}.", material->type(), + material_type()); + params_ = static_cast(material); + } + + extract_from_pack(buffer, inv_plastic_rcg_last_); + extract_from_pack(buffer, inv_plastic_rcg_current_); + extract_from_pack(buffer, accumulated_plastic_strain_last_); + extract_from_pack(buffer, accumulated_plastic_strain_current_); + extract_from_pack(buffer, accumulated_plastic_volumetric_strain_last_); + extract_from_pack(buffer, accumulated_plastic_volumetric_strain_current_); + extract_from_pack(buffer, dissipated_energy_last_); + extract_from_pack(buffer, dissipated_energy_current_); +} + +void Mat::PlasticMohrCoulomb::setup(int numgp, const Discret::Elements::Fibers& fibers, + const std::optional& coord_system) +{ + const auto identity = Core::LinAlg::TensorGenerators::identity; + inv_plastic_rcg_last_.assign(numgp, identity); + inv_plastic_rcg_current_.assign(numgp, identity); + accumulated_plastic_strain_last_.assign(numgp, 0.0); + accumulated_plastic_strain_current_.assign(numgp, 0.0); + accumulated_plastic_volumetric_strain_last_.assign(numgp, 0.0); + accumulated_plastic_volumetric_strain_current_.assign(numgp, 0.0); + dissipated_energy_last_.assign(numgp, 0.0); + dissipated_energy_current_.assign(numgp, 0.0); +} + +void Mat::PlasticMohrCoulomb::update() +{ + inv_plastic_rcg_last_ = inv_plastic_rcg_current_; + accumulated_plastic_strain_last_ = accumulated_plastic_strain_current_; + accumulated_plastic_volumetric_strain_last_ = accumulated_plastic_volumetric_strain_current_; + dissipated_energy_last_ = dissipated_energy_current_; +} + +void Mat::PlasticMohrCoulomb::evaluate(const Core::LinAlg::Tensor* defgrad, + const Core::LinAlg::SymmetricTensor& glstrain, + const Teuchos::ParameterList& params, const EvaluationContext<3>& context, + Core::LinAlg::SymmetricTensor& stress, + Core::LinAlg::SymmetricTensor& cmat, int gp, int eleGID) +{ + FOUR_C_ASSERT_ALWAYS(defgrad != nullptr, "Mohr-Coulomb requires the deformation gradient."); + const double determinant_deformation_gradient = Core::LinAlg::det(*defgrad); + FOUR_C_ASSERT_ALWAYS(determinant_deformation_gradient > 0.0, + "Mohr-Coulomb requires a positive deformation-gradient determinant, got {}.", + determinant_deformation_gradient); + + const double shear_modulus = params_->youngs_ / (2.0 * (1.0 + params_->poisson_ratio_)); + const double lame_parameter = + params_->youngs_ * params_->poisson_ratio_ / + ((1.0 + params_->poisson_ratio_) * (1.0 - 2.0 * params_->poisson_ratio_)); + const double bulk_modulus = params_->youngs_ / (3.0 * (1.0 - 2.0 * params_->poisson_ratio_)); + const auto elastic_tangent = elastic_principal_tangent(shear_modulus, lame_parameter); + + const double sin_friction = std::sin(params_->friction_angle_); + const double sin_dilatancy = std::sin(params_->dilatancy_angle_); + const double friction_stress_ratio = (1.0 + sin_friction) / (1.0 - sin_friction); + const double dilatancy_stress_ratio = (1.0 + sin_dilatancy) / (1.0 - sin_dilatancy); + + const auto inverse_deformation_gradient = Core::LinAlg::inv(*defgrad); + const auto elastic_lcg_trial = Core::LinAlg::assume_symmetry( + *defgrad * inv_plastic_rcg_last_.at(gp) * Core::LinAlg::transpose(*defgrad)); + const auto& [ascending_stretch_squares, ascending_eigenvectors] = + Core::LinAlg::eig(elastic_lcg_trial); + + std::array stretch_squares{}; + std::array, 3> spatial_principal_directions{}; + std::array, 3> material_principal_directions{}; + Core::LinAlg::Tensor trial_logarithmic_strain{}; + for (int i = 0; i < 3; ++i) + { + const int ascending_index = 2 - i; + stretch_squares[i] = ascending_stretch_squares[ascending_index]; + FOUR_C_ASSERT_ALWAYS(stretch_squares[i] > 0.0, + "Trial elastic left Cauchy-Green tensor is not positive definite."); + trial_logarithmic_strain(i) = 0.5 * std::log(stretch_squares[i]); + for (int j = 0; j < 3; ++j) + spatial_principal_directions[i](j) = ascending_eigenvectors(j, ascending_index); + material_principal_directions[i] = + inverse_deformation_gradient * spatial_principal_directions[i]; + } + + const auto trial_stress = elastic_tangent * trial_logarithmic_strain; + const double q_last = 2.0 * std::sqrt(friction_stress_ratio) * + cohesion(*params_, accumulated_plastic_strain_last_.at(gp)); + + ActiveSetResult return_result{}; + // checking for plasticity + if (max_yield_function(trial_stress, friction_stress_ratio, q_last) <= params_->tolerance_) + { + // elastic regime + return_result.valid = true; + return_result.stress = trial_stress; + return_result.tangent = elastic_tangent; + return_result.accumulated_plastic_strain = accumulated_plastic_strain_last_.at(gp); + } + else + { + // plastic regime + constexpr std::array, 2> face_planes{{{0, 2}, {0, 0}}}; + return_result = return_to_active_set(trial_stress, accumulated_plastic_strain_last_.at(gp), + elastic_tangent, friction_stress_ratio, dilatancy_stress_ratio, face_planes, 1, *params_); + + // An inadmissible smooth-face return requires testing the adjacent edge returns. + if (!return_result.valid) + { + constexpr std::array, 2> compression_edge_planes{{{0, 2}, {1, 2}}}; + const ActiveSetResult compression_edge = return_to_active_set(trial_stress, + accumulated_plastic_strain_last_.at(gp), elastic_tangent, friction_stress_ratio, + dilatancy_stress_ratio, compression_edge_planes, 2, *params_); + + constexpr std::array, 2> tension_edge_planes{{{0, 2}, {0, 1}}}; + const ActiveSetResult tension_edge = return_to_active_set(trial_stress, + accumulated_plastic_strain_last_.at(gp), elastic_tangent, friction_stress_ratio, + dilatancy_stress_ratio, tension_edge_planes, 2, *params_); + + if (compression_edge.valid && tension_edge.valid) + { + const double compression_correction = + Core::LinAlg::dot(compression_edge.plastic_strain_increment, + elastic_tangent * compression_edge.plastic_strain_increment); + const double tension_correction = Core::LinAlg::dot(tension_edge.plastic_strain_increment, + elastic_tangent * tension_edge.plastic_strain_increment); + return_result = + compression_correction <= tension_correction ? compression_edge : tension_edge; + } + else if (compression_edge.valid) + { + return_result = compression_edge; + } + else if (tension_edge.valid) + { + return_result = tension_edge; + } + else if (compression_edge.numerical_failure || tension_edge.numerical_failure) + { + FOUR_C_THROW("Mohr-Coulomb edge return failed to converge within {} local iterations.", + params_->max_iterations_); + } + } + + // If neither a smooth-face nor an edge return is admissible, return to the apex. + if (!return_result.valid) + { + if (return_result.numerical_failure) + FOUR_C_THROW("Mohr-Coulomb face return failed to converge within {} local iterations.", + params_->max_iterations_); + + const auto trial_deviator = deviator(trial_stress); + return_result.accumulated_plastic_strain = + accumulated_plastic_strain_last_.at(gp) + + std::sqrt(2.0 / 3.0) * Core::LinAlg::norm2(trial_deviator) / (2.0 * shear_modulus); + const double q = 2.0 * std::sqrt(friction_stress_ratio) * + cohesion(*params_, return_result.accumulated_plastic_strain); + const double apex_stress = q / (friction_stress_ratio - 1.0); + return_result.stress = + Core::LinAlg::Tensor{{apex_stress, apex_stress, apex_stress}}; + + const auto stress_correction = trial_stress - return_result.stress; + const auto stress_correction_deviator = deviator(stress_correction); + const double mean_stress_correction = sum_components(stress_correction) / 3.0; + return_result.plastic_strain_increment = + 1.0 / (2.0 * shear_modulus) * stress_correction_deviator + + mean_stress_correction / (3.0 * bulk_modulus) * + Core::LinAlg::Tensor{{1.0, 1.0, 1.0}}; + FOUR_C_ASSERT_ALWAYS(belongs_to_apex_flow_cone(return_result.plastic_strain_increment, + dilatancy_stress_ratio, 100.0 * params_->tolerance_), + "Trial state has no admissible Mohr-Coulomb apex return for the selected dilatation " + "angle."); + + return_result.tangent = Core::LinAlg::Tensor{}; + const auto trial_logarithmic_strain_deviator = deviator(trial_logarithmic_strain); + const double deviator_norm = Core::LinAlg::norm2(trial_logarithmic_strain_deviator); + if (deviator_norm > params_->tolerance_) + { + const double apex_hardening = + 2.0 * std::sqrt(friction_stress_ratio) * + cohesion_derivative(*params_, return_result.accumulated_plastic_strain) / + (friction_stress_ratio - 1.0); + const Core::LinAlg::Tensor accumulated_strain_derivative = + std::sqrt(2.0 / 3.0) / deviator_norm * trial_logarithmic_strain_deviator; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + return_result.tangent(i, j) = apex_hardening * accumulated_strain_derivative(j); + } + return_result.valid = true; + } + } + + Core::LinAlg::SymmetricTensor elastic_lcg{}; + for (int i = 0; i < 3; ++i) + { + const double updated_stretch_square = + stretch_squares[i] * std::exp(-2.0 * return_result.plastic_strain_increment(i)); + elastic_lcg += + updated_stretch_square * Core::LinAlg::self_dyadic(spatial_principal_directions[i]); + stress += return_result.stress(i) * Core::LinAlg::self_dyadic(material_principal_directions[i]); + } + + for (int a = 0; a < 3; ++a) + { + const auto material_projector = Core::LinAlg::self_dyadic(material_principal_directions[a]); + cmat += -2.0 * return_result.stress(a) * + Core::LinAlg::dyadic(material_projector, material_projector); + + for (int b = 0; b < 3; ++b) + { + const auto material_projector_b = Core::LinAlg::self_dyadic(material_principal_directions[b]); + cmat += return_result.tangent(a, b) * + Core::LinAlg::dyadic(material_projector, material_projector_b); + + if (a != b) + { + const double denominator = stretch_squares[a] - stretch_squares[b]; + const double factor = + std::abs(denominator) > 100.0 * std::numeric_limits::epsilon() + ? (return_result.stress(a) * stretch_squares[b] - + return_result.stress(b) * stretch_squares[a]) / + denominator + : 0.5 * (return_result.tangent(b, b) - return_result.tangent(a, b)) - + return_result.stress(b); + const auto mixed_projector = Core::LinAlg::dyadic( + material_principal_directions[a], material_principal_directions[b]); + cmat += factor * Core::LinAlg::assume_symmetry( + Core::LinAlg::dyadic(mixed_projector, mixed_projector)); + cmat += factor * Core::LinAlg::assume_symmetry(Core::LinAlg::dyadic( + mixed_projector, Core::LinAlg::transpose(mixed_projector))); + } + } + } + + inv_plastic_rcg_current_.at(gp) = + Core::LinAlg::assume_symmetry(inverse_deformation_gradient * elastic_lcg * + Core::LinAlg::transpose(inverse_deformation_gradient)); + accumulated_plastic_strain_current_.at(gp) = return_result.accumulated_plastic_strain; + accumulated_plastic_volumetric_strain_current_.at(gp) = + accumulated_plastic_volumetric_strain_last_.at(gp) + + sum_components(return_result.plastic_strain_increment); + dissipated_energy_current_.at(gp) = + dissipated_energy_last_.at(gp) + + Core::LinAlg::dot(return_result.stress, return_result.plastic_strain_increment); +} + +void Mat::PlasticMohrCoulomb::register_output_data_names( + std::unordered_map& names_and_size) const +{ + names_and_size["plastic_strain"] = 6; + names_and_size["accumulated_plastic_strain"] = 1; + names_and_size["accumulated_plastic_volumetric_strain"] = 1; + names_and_size["local_dissipated_energy"] = 1; +} + +bool Mat::PlasticMohrCoulomb::evaluate_output_data( + const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const +{ + if (name == "accumulated_plastic_strain") + { + for (std::size_t gp = 0; gp < accumulated_plastic_strain_current_.size(); ++gp) + data(gp, 0) = accumulated_plastic_strain_current_[gp]; + return true; + } + if (name == "accumulated_plastic_volumetric_strain") + { + for (std::size_t gp = 0; gp < accumulated_plastic_volumetric_strain_current_.size(); ++gp) + data(gp, 0) = accumulated_plastic_volumetric_strain_current_[gp]; + return true; + } + if (name == "local_dissipated_energy") + { + for (std::size_t gp = 0; gp < dissipated_energy_current_.size(); ++gp) + data(gp, 0) = dissipated_energy_current_[gp]; + return true; + } + if (name == "plastic_strain") + { + for (std::size_t gp = 0; gp < inv_plastic_rcg_current_.size(); ++gp) + { + const auto& [eigenvalues, eigenvectors] = Core::LinAlg::eig(inv_plastic_rcg_current_[gp]); + Core::LinAlg::SymmetricTensor logarithmic_plastic_strain{}; + for (int i = 0; i < 3; ++i) + { + FOUR_C_ASSERT_ALWAYS(eigenvalues[i] > 0.0, + "Inverse plastic right Cauchy-Green tensor is not positive definite."); + Core::LinAlg::Tensor direction{}; + for (int j = 0; j < 3; ++j) direction(j) = eigenvectors(j, i); + logarithmic_plastic_strain += + -0.5 * std::log(eigenvalues[i]) * Core::LinAlg::self_dyadic(direction); + } + const Core::LinAlg::Matrix<6, 1> voigt = + Core::LinAlg::make_strain_like_voigt_matrix(logarithmic_plastic_strain); + for (int i = 0; i < 6; ++i) data(gp, i) = voigt(i, 0); + } + return true; + } + return false; +} + +FOUR_C_NAMESPACE_CLOSE diff --git a/src/mat/4C_mat_plasticmohrcoulomb.hpp b/src/mat/4C_mat_plasticmohrcoulomb.hpp new file mode 100644 index 00000000000..81b2db8d830 --- /dev/null +++ b/src/mat/4C_mat_plasticmohrcoulomb.hpp @@ -0,0 +1,118 @@ +// This file is part of 4C multiphysics licensed under the +// GNU Lesser General Public License v3.0 or later. +// +// See the LICENSE.md file in the top-level for license information. +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#ifndef FOUR_C_MAT_PLASTICMOHRCOULOMB_HPP +#define FOUR_C_MAT_PLASTICMOHRCOULOMB_HPP + +#include "4C_config.hpp" + +#include "4C_comm_parobjectfactory.hpp" +#include "4C_mat_so3_material.hpp" +#include "4C_material_parameter_base.hpp" + +FOUR_C_NAMESPACE_OPEN + +namespace Mat +{ + namespace PAR + { + class PlasticMohrCoulomb : public Core::Mat::PAR::Parameter + { + public: + explicit PlasticMohrCoulomb(const Core::Mat::PAR::Parameter::Data& matdata); + + const double youngs_; + const double poisson_ratio_; + const double density_; + const double cohesion_; + const double friction_angle_; + const double dilatancy_angle_; + const double linear_hardening_; + const double saturation_hardening_; + const double hardening_exponent_; + const double tolerance_; + const int max_iterations_; + + std::shared_ptr create_material() override; + }; + } // namespace PAR + + class PlasticMohrCoulombType : public Core::Communication::ParObjectType + { + public: + std::string name() const override { return "PlasticMohrCoulombType"; } + static PlasticMohrCoulombType& instance() { return instance_; } + Core::Communication::ParObject* create(Core::Communication::UnpackBuffer& buffer) override; + + private: + static PlasticMohrCoulombType instance_; + }; + + class PlasticMohrCoulomb : public So3Material + { + public: + PlasticMohrCoulomb(); + explicit PlasticMohrCoulomb(Mat::PAR::PlasticMohrCoulomb* params); + + int unique_par_object_id() const override + { + return PlasticMohrCoulombType::instance().unique_par_object_id(); + } + + void pack(Core::Communication::PackBuffer& data) const override; + void unpack(Core::Communication::UnpackBuffer& buffer) override; + + Core::Materials::MaterialType material_type() const override + { + return Core::Materials::m_plmohrcoulomb; + } + + void valid_kinematics(Solid::KinemType kinem) override + { + if (kinem != Solid::KinemType::nonlinearTotLag) + FOUR_C_THROW("The Mohr-Coulomb material requires nonlinear total-Lagrangian kinematics."); + } + + std::shared_ptr clone() const override + { + return std::make_shared(*this); + } + + Core::Mat::PAR::Parameter* parameter() const override { return params_; } + double density() const override { return params_->density_; } + + void setup(int numgp, const Discret::Elements::Fibers& fibers, + const std::optional& coord_system) override; + void update() override; + + void evaluate(const Core::LinAlg::Tensor* defgrad, + const Core::LinAlg::SymmetricTensor& glstrain, + const Teuchos::ParameterList& params, const EvaluationContext<3>& context, + Core::LinAlg::SymmetricTensor& stress, + Core::LinAlg::SymmetricTensor& cmat, int gp, int eleGID) override; + + void register_output_data_names( + std::unordered_map& names_and_size) const override; + bool evaluate_output_data( + const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const override; + + private: + Mat::PAR::PlasticMohrCoulomb* params_; + std::vector> inv_plastic_rcg_last_; + std::vector> inv_plastic_rcg_current_; + std::vector accumulated_plastic_strain_last_; + std::vector accumulated_plastic_strain_current_; + std::vector accumulated_plastic_volumetric_strain_last_; + std::vector accumulated_plastic_volumetric_strain_current_; + std::vector dissipated_energy_last_; + std::vector dissipated_energy_current_; + }; +} // namespace Mat + +FOUR_C_NAMESPACE_CLOSE + +#endif diff --git a/tests/input_files/mat_mohrcoulomb_return_to_apex.4C.yaml b/tests/input_files/mat_mohrcoulomb_return_to_apex.4C.yaml new file mode 100644 index 00000000000..46c7fcf04f3 --- /dev/null +++ b/tests/input_files/mat_mohrcoulomb_return_to_apex.4C.yaml @@ -0,0 +1,126 @@ +TITLE: + - "One-element structure test for the Mohr-Coulomb apex return" + - "The prescribed homogeneous deformation matches the constitutive material test" +PROBLEM TYPE: + PROBLEMTYPE: "Structure" +IO: + STRUCT_STRESS: "Cauchy" + STRUCT_STRAIN: "LOG" + VERBOSITY: "Standard" +IO/RUNTIME VTK OUTPUT: + INTERVAL_STEPS: 1 +IO/RUNTIME VTK OUTPUT/STRUCTURE: + OUTPUT_STRUCTURE: true + STRESS_STRAIN: true +SOLVER 1: + SOLVER: "UMFPACK" + NAME: "Structure_Solver" +STRUCTURAL DYNAMIC: + DYNAMICTYPE: "Statics" + TIMESTEP: 1 + NUMSTEP: 1 + MAXTIME: 1 + TOLDISP: 1e-12 + TOLRES: 1e-12 + MAXITER: 20 + LOADLIN: true + LINEAR_SOLVER: 1 +MATERIALS: + - MAT: 1 + MAT_Struct_MohrCoulomb: + YOUNG_MODULUS: 1000 + POISSON_RATIO: 0.25 + DENSITY: 0 + COHESION: 10 + FRICTION_ANGLE: 0.5235987755982988 + DILATANCY_ANGLE: 0.3490658503988659 + TOLERANCE: 1e-12 + MAX_ITERATIONS: 100 +FUNCT1: + - SYMBOLIC_FUNCTION_OF_SPACE_TIME: "t" +RESULT DESCRIPTION: + - STRUCTURE: + DIS: "structure" + NODE: 8 + QUANTITY: "stress_xx" + VALUE: 16.527344154661588 + TOLERANCE: 1e-10 + - STRUCTURE: + DIS: "structure" + NODE: 8 + QUANTITY: "stress_yy" + VALUE: 16.527344154661584 + TOLERANCE: 1e-10 + - STRUCTURE: + DIS: "structure" + NODE: 8 + QUANTITY: "stress_zz" + VALUE: 16.527344154661584 + TOLERANCE: 1e-10 +DESIGN SURF DIRICH CONDITIONS: + - E: 1 + NUMDOF: 3 + ONOFF: [1, 0, 0] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 2 + NUMDOF: 3 + ONOFF: [1, 0, 0] + VAL: [0.02531512052442884, 0, 0] + FUNCT: [1, 0, 0] + - E: 3 + NUMDOF: 3 + ONOFF: [0, 1, 0] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 4 + NUMDOF: 3 + ONOFF: [0, 1, 0] + VAL: [0, 0.01099753646217638, 0] + FUNCT: [0, 1, 0] + - E: 5 + NUMDOF: 3 + ONOFF: [0, 0, 1] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 6 + NUMDOF: 3 + ONOFF: [0, 0, 1] + VAL: [0, 0, 0.01099753646217638] + FUNCT: [0, 0, 1] +DSURF-NODE TOPOLOGY: + - "NODE 1 DSURFACE 1" + - "NODE 2 DSURFACE 1" + - "NODE 3 DSURFACE 1" + - "NODE 4 DSURFACE 1" + - "NODE 5 DSURFACE 2" + - "NODE 6 DSURFACE 2" + - "NODE 7 DSURFACE 2" + - "NODE 8 DSURFACE 2" + - "NODE 1 DSURFACE 3" + - "NODE 2 DSURFACE 3" + - "NODE 5 DSURFACE 3" + - "NODE 6 DSURFACE 3" + - "NODE 3 DSURFACE 4" + - "NODE 4 DSURFACE 4" + - "NODE 7 DSURFACE 4" + - "NODE 8 DSURFACE 4" + - "NODE 2 DSURFACE 5" + - "NODE 3 DSURFACE 5" + - "NODE 6 DSURFACE 5" + - "NODE 7 DSURFACE 5" + - "NODE 1 DSURFACE 6" + - "NODE 4 DSURFACE 6" + - "NODE 5 DSURFACE 6" + - "NODE 8 DSURFACE 6" +NODE COORDS: + - "NODE 1 COORD 0 0 1" + - "NODE 2 COORD 0 0 0" + - "NODE 3 COORD 0 1 0" + - "NODE 4 COORD 0 1 1" + - "NODE 5 COORD 1 0 1" + - "NODE 6 COORD 1 0 0" + - "NODE 7 COORD 1 1 0" + - "NODE 8 COORD 1 1 1" +STRUCTURE ELEMENTS: + - "1 SOLID HEX8 1 2 3 4 5 6 7 8 MAT 1 KINEM nonlinear" diff --git a/tests/input_files/mat_mohrcoulomb_return_to_compression_edge.4C.yaml b/tests/input_files/mat_mohrcoulomb_return_to_compression_edge.4C.yaml new file mode 100644 index 00000000000..f3d8f8e0a8c --- /dev/null +++ b/tests/input_files/mat_mohrcoulomb_return_to_compression_edge.4C.yaml @@ -0,0 +1,126 @@ +TITLE: + - "One-element structure test for the Mohr-Coulomb compression-edge return" + - "The prescribed homogeneous deformation matches the constitutive material test" +PROBLEM TYPE: + PROBLEMTYPE: "Structure" +IO: + STRUCT_STRESS: "Cauchy" + STRUCT_STRAIN: "LOG" + VERBOSITY: "Standard" +IO/RUNTIME VTK OUTPUT: + INTERVAL_STEPS: 1 +IO/RUNTIME VTK OUTPUT/STRUCTURE: + OUTPUT_STRUCTURE: true + STRESS_STRAIN: true +SOLVER 1: + SOLVER: "UMFPACK" + NAME: "Structure_Solver" +STRUCTURAL DYNAMIC: + DYNAMICTYPE: "Statics" + TIMESTEP: 1 + NUMSTEP: 1 + MAXTIME: 1 + TOLDISP: 1e-12 + TOLRES: 1e-12 + MAXITER: 20 + LOADLIN: true + LINEAR_SOLVER: 1 +MATERIALS: + - MAT: 1 + MAT_Struct_MohrCoulomb: + YOUNG_MODULUS: 1000 + POISSON_RATIO: 0.25 + DENSITY: 0 + COHESION: 10 + FRICTION_ANGLE: 0.5235987755982988 + DILATANCY_ANGLE: 0.3490658503988659 + TOLERANCE: 1e-12 + MAX_ITERATIONS: 100 +FUNCT1: + - SYMBOLIC_FUNCTION_OF_SPACE_TIME: "t" +RESULT DESCRIPTION: + - STRUCTURE: + DIS: "structure" + NODE: 8 + QUANTITY: "stress_xx" + VALUE: -73.63415802088231 + TOLERANCE: 1e-10 + - STRUCTURE: + DIS: "structure" + NODE: 8 + QUANTITY: "stress_yy" + VALUE: -73.63415802088224 + TOLERANCE: 1e-10 + - STRUCTURE: + DIS: "structure" + NODE: 8 + QUANTITY: "stress_zz" + VALUE: -261.71956603533334 + TOLERANCE: 1e-10 +DESIGN SURF DIRICH CONDITIONS: + - E: 1 + NUMDOF: 3 + ONOFF: [1, 0, 0] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 2 + NUMDOF: 3 + ONOFF: [1, 0, 0] + VAL: [0.02531512052442884, 0, 0] + FUNCT: [1, 0, 0] + - E: 3 + NUMDOF: 3 + ONOFF: [0, 1, 0] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 4 + NUMDOF: 3 + ONOFF: [0, 1, 0] + VAL: [0, 0.01099753646217638, 0] + FUNCT: [0, 1, 0] + - E: 5 + NUMDOF: 3 + ONOFF: [0, 0, 1] + VAL: [0, 0, 0] + FUNCT: [0, 0, 0] + - E: 6 + NUMDOF: 3 + ONOFF: [0, 0, 1] + VAL: [0, 0, -0.1812692469220182] + FUNCT: [0, 0, 1] +DSURF-NODE TOPOLOGY: + - "NODE 1 DSURFACE 1" + - "NODE 2 DSURFACE 1" + - "NODE 3 DSURFACE 1" + - "NODE 4 DSURFACE 1" + - "NODE 5 DSURFACE 2" + - "NODE 6 DSURFACE 2" + - "NODE 7 DSURFACE 2" + - "NODE 8 DSURFACE 2" + - "NODE 1 DSURFACE 3" + - "NODE 2 DSURFACE 3" + - "NODE 5 DSURFACE 3" + - "NODE 6 DSURFACE 3" + - "NODE 3 DSURFACE 4" + - "NODE 4 DSURFACE 4" + - "NODE 7 DSURFACE 4" + - "NODE 8 DSURFACE 4" + - "NODE 2 DSURFACE 5" + - "NODE 3 DSURFACE 5" + - "NODE 6 DSURFACE 5" + - "NODE 7 DSURFACE 5" + - "NODE 1 DSURFACE 6" + - "NODE 4 DSURFACE 6" + - "NODE 5 DSURFACE 6" + - "NODE 8 DSURFACE 6" +NODE COORDS: + - "NODE 1 COORD 0 0 1" + - "NODE 2 COORD 0 0 0" + - "NODE 3 COORD 0 1 0" + - "NODE 4 COORD 0 1 1" + - "NODE 5 COORD 1 0 1" + - "NODE 6 COORD 1 0 0" + - "NODE 7 COORD 1 1 0" + - "NODE 8 COORD 1 1 1" +STRUCTURE ELEMENTS: + - "1 SOLID HEX8 1 2 3 4 5 6 7 8 MAT 1 KINEM nonlinear" diff --git a/tests/list_of_tests.cmake b/tests/list_of_tests.cmake index 6246475cd87..bf73edcb20d 100644 --- a/tests/list_of_tests.cmake +++ b/tests/list_of_tests.cmake @@ -910,6 +910,8 @@ four_c_test(TEST_FILE mat_iso_thermoviscoplast_refJC_log_timint_tsi_monolithic.4 __four_c_test_restart(BASED_ON ${current} SAME_FILE RESTART_STEP 90) four_c_test(TEST_FILE mat_iso_thermoviscoplast_refJC_log_timint_tsi_partitioned.4C.yaml RETURN_AS current) __four_c_test_restart(BASED_ON ${current} SAME_FILE RESTART_STEP 90) +four_c_test(TEST_FILE mat_mohrcoulomb_return_to_apex.4C.yaml) +four_c_test(TEST_FILE mat_mohrcoulomb_return_to_compression_edge.4C.yaml) four_c_test(TEST_FILE mat_muscle_combo_hex.4C.yaml NP 2 RETURN_AS current) __four_c_test_restart(BASED_ON ${current} SAME_FILE NP 2 RESTART_STEP 140) four_c_test(TEST_FILE mat_muscle_combo_hex_act_map_every_timestep.4C.yaml NP 2 RETURN_AS current) diff --git a/unittests/mat/4C_mohrcoulomb_test.cpp b/unittests/mat/4C_mohrcoulomb_test.cpp new file mode 100644 index 00000000000..83c8effc067 --- /dev/null +++ b/unittests/mat/4C_mohrcoulomb_test.cpp @@ -0,0 +1,327 @@ +// This file is part of 4C multiphysics licensed under the +// GNU Lesser General Public License v3.0 or later. +// +// See the LICENSE.md file in the top-level for license information. +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include + +#include "4C_comm_pack_buffer.hpp" +#include "4C_global_data.hpp" +#include "4C_io_input_parameter_container.templates.hpp" +#include "4C_linalg_serialdensematrix.hpp" +#include "4C_linalg_tensor_generators.hpp" +#include "4C_mat_material_factory.hpp" +#include "4C_mat_par_bundle.hpp" +#include "4C_mat_plasticmohrcoulomb.hpp" +#include "4C_material_parameter_base.hpp" +#include "4C_unittest_utils_assertions_test.hpp" +#include "4C_utils_singleton_owner.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + using namespace FourC; + + struct EvaluationResult + { + std::array kirchhoff_stress{}; + Core::LinAlg::SymmetricTensor pk2_stress{}; + Core::LinAlg::SymmetricTensor tangent{}; + Core::LinAlg::SymmetricTensor green_lagrange_strain{}; + }; + + class MohrCoulombTest : public ::testing::Test + { + protected: + void SetUp() override { create_material(); } + + void create_material(double dilatancy_angle = 20.0 * std::numbers::pi / 180.0, + double linear_hardening = 0.0, double saturation_hardening = 0.0, + double hardening_exponent = 0.0, int max_iterations = 100) + { + Core::IO::InputParameterContainer container; + container.add("YOUNG_MODULUS", 1000.0); + container.add("POISSON_RATIO", 0.25); + container.add("DENSITY", 0.0); + container.add("COHESION", 10.0); + container.add("FRICTION_ANGLE", 30.0 * std::numbers::pi / 180.0); + container.add("DILATANCY_ANGLE", dilatancy_angle); + container.add("LINEAR_HARDENING", linear_hardening); + container.add("SATURATION_HARDENING", saturation_hardening); + container.add("HARDENING_EXP", hardening_exponent); + container.add("TOLERANCE", 1.0e-12); + container.add("MAX_ITERATIONS", max_iterations); + + const int material_id = next_material_id_++; + parameter_ = std::shared_ptr( + Mat::make_parameter(material_id, Core::Materials::m_plmohrcoulomb, container)); + Global::Problem& problem = *Global::Problem::instance(); + problem.materials()->set_read_from_problem(0); + problem.materials()->insert(material_id, parameter_); + material_ = std::make_shared( + dynamic_cast(parameter_.get())); + material_->setup(1, {}, {}); + } + + EvaluationResult evaluate(const std::array& logarithmic_strain) + { + Core::LinAlg::Tensor deformation_gradient{}; + EvaluationResult result; + for (int i = 0; i < 3; ++i) + { + deformation_gradient(i, i) = std::exp(logarithmic_strain[i]); + result.green_lagrange_strain(i, i) = 0.5 * (std::exp(2.0 * logarithmic_strain[i]) - 1.0); + } + + Teuchos::ParameterList parameters; + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + material_->evaluate(&deformation_gradient, result.green_lagrange_strain, parameters, context, + result.pk2_stress, result.tangent, 0, 0); + for (int i = 0; i < 3; ++i) + result.kirchhoff_stress[i] = + result.pk2_stress(i, i) * std::exp(2.0 * logarithmic_strain[i]); + std::ranges::sort(result.kirchhoff_stress, std::greater<>()); + return result; + } + + double scalar_output(const std::string& name) + { + Core::LinAlg::SerialDenseMatrix data(1, 1); + EXPECT_TRUE(material_->evaluate_output_data(name, data)); + return data(0, 0); + } + + static double yield_function(const std::array& stress, double cohesion_value) + { + constexpr double friction_angle = 30.0 * std::numbers::pi / 180.0; + const double k = (1.0 + std::sin(friction_angle)) / (1.0 - std::sin(friction_angle)); + return k * stress[0] - stress[2] - 2.0 * std::sqrt(k) * cohesion_value; + } + + std::shared_ptr parameter_; + std::shared_ptr material_; + int next_material_id_ = 1; + Core::Utils::SingletonOwnerRegistry::ScopeGuard guard_; + }; + + TEST_F(MohrCoulombTest, ElasticHenckyResponse) + { + const std::array logarithmic_strain{0.001, -0.0005, -0.001}; + const EvaluationResult result = evaluate(logarithmic_strain); + + constexpr double shear_modulus = 400.0; + constexpr double lame_parameter = 400.0; + std::array expected{}; + const double trace_logarithmic_strain = + logarithmic_strain[0] + logarithmic_strain[1] + logarithmic_strain[2]; + for (int i = 0; i < 3; ++i) + expected[i] = + 2.0 * shear_modulus * logarithmic_strain[i] + lame_parameter * trace_logarithmic_strain; + std::ranges::sort(expected, std::greater<>()); + + for (int i = 0; i < 3; ++i) EXPECT_NEAR(result.kirchhoff_stress[i], expected[i], 1.0e-12); + EXPECT_DOUBLE_EQ(scalar_output("accumulated_plastic_strain"), 0.0); + } + + TEST_F(MohrCoulombTest, ReturnsToSmoothFace) + { + const EvaluationResult result = evaluate({0.025, -0.0453125, -0.2}); + EXPECT_NEAR(yield_function(result.kirchhoff_stress, 10.0), 0.0, 1.0e-9); + EXPECT_GT(result.kirchhoff_stress[0] - result.kirchhoff_stress[1], 1.0); + EXPECT_GT(result.kirchhoff_stress[1] - result.kirchhoff_stress[2], 1.0); + } + + TEST_F(MohrCoulombTest, ReturnsToCompressionEdge) + { + const EvaluationResult result = evaluate({0.025, 0.0109375, -0.2}); + EXPECT_NEAR(yield_function(result.kirchhoff_stress, 10.0), 0.0, 1.0e-9); + EXPECT_NEAR(result.kirchhoff_stress[0], result.kirchhoff_stress[1], 1.0e-9); + } + + TEST_F(MohrCoulombTest, ReturnsToTensionEdge) + { + const EvaluationResult result = evaluate({0.025, -0.0453125, -0.0453125}); + EXPECT_NEAR(yield_function(result.kirchhoff_stress, 10.0), 0.0, 1.0e-9); + EXPECT_NEAR(result.kirchhoff_stress[1], result.kirchhoff_stress[2], 1.0e-9); + } + + TEST_F(MohrCoulombTest, ReturnsToApex) + { + const EvaluationResult result = evaluate({0.025, 0.0109375, 0.0109375}); + EXPECT_NEAR(result.kirchhoff_stress[0], result.kirchhoff_stress[1], 1.0e-9); + EXPECT_NEAR(result.kirchhoff_stress[1], result.kirchhoff_stress[2], 1.0e-9); + EXPECT_NEAR(yield_function(result.kirchhoff_stress, 10.0), 0.0, 1.0e-9); + } + + TEST_F(MohrCoulombTest, NonAssociatedFlowAndOutputs) + { + evaluate({0.025, -0.0453125, -0.2}); + + EXPECT_GT(scalar_output("accumulated_plastic_strain"), 0.0); + EXPECT_GT(scalar_output("accumulated_plastic_volumetric_strain"), 0.0); + EXPECT_GT(scalar_output("local_dissipated_energy"), 0.0); + + Core::LinAlg::SerialDenseMatrix plastic_strain(1, 6); + ASSERT_TRUE(material_->evaluate_output_data("plastic_strain", plastic_strain)); + double norm = 0.0; + for (int i = 0; i < 6; ++i) norm += plastic_strain(0, i) * plastic_strain(0, i); + EXPECT_GT(norm, 0.0); + } + + TEST_F(MohrCoulombTest, AssociativeFlow) + { + create_material(30.0 * std::numbers::pi / 180.0); + const EvaluationResult result = evaluate({0.025, -0.0453125, -0.2}); + + EXPECT_NEAR(yield_function(result.kirchhoff_stress, 10.0), 0.0, 1.0e-9); + EXPECT_GT(scalar_output("accumulated_plastic_volumetric_strain"), 0.0); + } + + TEST_F(MohrCoulombTest, RejectsZeroDilatation) { EXPECT_ANY_THROW(create_material(0.0)); } + + TEST_F(MohrCoulombTest, ReportsLocalSolverFailure) + { + create_material(20.0 * std::numbers::pi / 180.0, 20.0, 30.0, 5.0, 1); + EXPECT_ANY_THROW(evaluate({0.025, -0.0453125, -0.2})); + } + + TEST_F(MohrCoulombTest, LinearAndVoceHardening) + { + create_material(20.0 * std::numbers::pi / 180.0, 20.0, 30.0, 5.0); + const EvaluationResult result = evaluate({0.025, -0.0453125, -0.2}); + const double plastic_strain = scalar_output("accumulated_plastic_strain"); + const double current_cohesion = + 10.0 + 20.0 * plastic_strain + 20.0 * (1.0 - std::exp(-5.0 * plastic_strain)); + + EXPECT_GT(plastic_strain, 0.0); + EXPECT_NEAR(yield_function(result.kirchhoff_stress, current_cohesion), 0.0, 1.0e-8); + } + + TEST_F(MohrCoulombTest, AnalyticalYieldSurfaceMeridians) + { + const char* output_path = std::getenv("FOUR_C_MOHR_COULOMB_YIELD_CURVE_OUTPUT"); + std::ofstream output; + if (output_path != nullptr) + { + output.open(output_path); + ASSERT_TRUE(output.good()); + output << "branch,pressure,analytical_von_mises,calculated_von_mises\n"; + } + + constexpr double bulk_modulus = 2000.0 / 3.0; + constexpr double shear_modulus = 400.0; + constexpr double friction_angle = 30.0 * std::numbers::pi / 180.0; + const double k = (1.0 + std::sin(friction_angle)) / (1.0 - std::sin(friction_angle)); + const double q0 = 2.0 * std::sqrt(k) * 10.0; + + const auto evaluate_meridian = [&](double trial_pressure, bool compression_meridian) + { + const double denominator = compression_meridian ? k + 2.0 : 2.0 * k + 1.0; + const double analytical_von_mises = 3.0 * (q0 + (k - 1.0) * trial_pressure) / denominator; + const double trial_von_mises = 1.25 * analytical_von_mises; + const std::array trial_stress = + compression_meridian + ? std::array{-trial_pressure + trial_von_mises / 3.0, + -trial_pressure + trial_von_mises / 3.0, + -trial_pressure - 2.0 * trial_von_mises / 3.0} + : std::array{-trial_pressure + 2.0 * trial_von_mises / 3.0, + -trial_pressure - trial_von_mises / 3.0, + -trial_pressure - trial_von_mises / 3.0}; + + const double mean_trial_stress = (trial_stress[0] + trial_stress[1] + trial_stress[2]) / 3.0; + std::array logarithmic_strain{}; + for (int i = 0; i < 3; ++i) + logarithmic_strain[i] = (trial_stress[i] - mean_trial_stress) / (2.0 * shear_modulus) + + mean_trial_stress / (3.0 * bulk_modulus); + + const EvaluationResult result = evaluate(logarithmic_strain); + const double calculated_pressure = + -(result.kirchhoff_stress[0] + result.kirchhoff_stress[1] + result.kirchhoff_stress[2]) / + 3.0; + const double calculated_von_mises = result.kirchhoff_stress[0] - result.kirchhoff_stress[2]; + const double analytical_at_calculated_pressure = + 3.0 * (q0 + (k - 1.0) * calculated_pressure) / denominator; + EXPECT_NEAR(calculated_von_mises, analytical_at_calculated_pressure, 1.0e-8); + + if (output.is_open()) + output << (compression_meridian ? "triaxial_compression" : "triaxial_tension") << ',' + << calculated_pressure << ',' << analytical_at_calculated_pressure << ',' + << calculated_von_mises << '\n'; + }; + + for (int point = 0; point <= 20; ++point) + { + const double pressure = 5.0 * point; + evaluate_meridian(pressure, true); + evaluate_meridian(pressure, false); + } + } + + TEST_F(MohrCoulombTest, SmoothFaceConsistentTangent) + { + const std::array logarithmic_strain{0.025, -0.0453125, -0.2}; + const EvaluationResult base = evaluate(logarithmic_strain); + constexpr double perturbation = 1.0e-7; + + for (int j = 0; j < 3; ++j) + { + auto perturbed_green_lagrange_strain = base.green_lagrange_strain; + perturbed_green_lagrange_strain(j, j) += perturbation; + Core::LinAlg::Tensor perturbed_deformation_gradient{}; + for (int i = 0; i < 3; ++i) + perturbed_deformation_gradient(i, i) = + std::sqrt(1.0 + 2.0 * perturbed_green_lagrange_strain(i, i)); + + Core::LinAlg::SymmetricTensor perturbed_stress{}; + Core::LinAlg::SymmetricTensor unused_tangent{}; + Teuchos::ParameterList parameters; + double total_time = 0.0; + double time_step_size = 1.0; + Mat::EvaluationContext<3> context{.total_time = &total_time, + .time_step_size = &time_step_size, + .xi = {}, + .ref_coords = nullptr}; + material_->evaluate(&perturbed_deformation_gradient, perturbed_green_lagrange_strain, + parameters, context, perturbed_stress, unused_tangent, 0, 0); + + for (int i = 0; i < 3; ++i) + { + const double finite_difference = + (perturbed_stress(i, i) - base.pk2_stress(i, i)) / perturbation; + EXPECT_NEAR(base.tangent(i, i, j, j), finite_difference, 2.0e-3); + } + } + } + + TEST_F(MohrCoulombTest, PackAndUnpackPreservesState) + { + evaluate({0.025, -0.0453125, -0.2}); + const double expected_plastic_strain = scalar_output("accumulated_plastic_strain"); + + Core::Communication::PackBuffer pack_buffer; + material_->pack(pack_buffer); + std::vector packed_data; + swap(packed_data, pack_buffer()); + Core::Communication::UnpackBuffer unpack_buffer(packed_data); + Mat::PlasticMohrCoulomb unpacked_material; + unpacked_material.unpack(unpack_buffer); + + Core::LinAlg::SerialDenseMatrix data(1, 1); + ASSERT_TRUE(unpacked_material.evaluate_output_data("accumulated_plastic_strain", data)); + EXPECT_DOUBLE_EQ(data(0, 0), expected_plastic_strain); + } +} // namespace