From 518fb3ba1fadabb5c041a43e7aee4ae8f39e05fc Mon Sep 17 00:00:00 2001 From: jduriez Date: Thu, 7 Dec 2023 12:35:37 +0100 Subject: [PATCH 1/6] 6 new state variables for Mohr Coulomb: 4 flags for failure (tensile or shear; now and/or past) + elastic energy + plastic dissipation --- include/materials/mohr_coulomb.tcc | 43 +++++++++++++++++++++++++++--- 1 file changed, 39 insertions(+), 4 deletions(-) diff --git a/include/materials/mohr_coulomb.tcc b/include/materials/mohr_coulomb.tcc index 3a8a352b7..f108fa7c1 100644 --- a/include/materials/mohr_coulomb.tcc +++ b/include/materials/mohr_coulomb.tcc @@ -70,7 +70,19 @@ mpm::dense_map mpm::MohrCoulomb::initialise_state_variables() { // Theta {"theta", 0.}, // Plastic deviatoric strain - {"pdstrain", 0.}}; + {"pdstrain", 0.}, + // Elastic energy + {"E_el", 0.}, + // Plastic work + {"W_pl", 0.}, + // Flag for current failure because of tensile stress + {"tensile_fail_curr", 0}, + // Flag for current failure because of shear stress + {"shear_fail_curr", 0}, + // Flag for past or current failure because of tensile stress + {"tensile_fail", 0}, + // Flag for past or current failure because of shear stress + {"shear_fail", 0}}; return state_vars; } @@ -78,7 +90,7 @@ mpm::dense_map mpm::MohrCoulomb::initialise_state_variables() { template std::vector mpm::MohrCoulomb::state_variables() const { const std::vector state_vars = { - "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain"}; + "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain", "E_el", "W_pl", "tensile_fail_curr", "shear_fail_curr", "tensile_fail", "shear_fail"}; return state_vars; } @@ -392,10 +404,21 @@ Eigen::Matrix mpm::MohrCoulomb::compute_stress( auto yield_type = this->compute_yield_state(&yield_function, (*state_vars)); // Initialise value of yield function based on stress double yield{std::numeric_limits::max()}; - if (yield_type == mpm::mohrcoulomb::FailureState::Tensile) + if (yield_type == mpm::mohrcoulomb::FailureState::Tensile){ yield = yield_function(0); - if (yield_type == mpm::mohrcoulomb::FailureState::Shear) + (*state_vars).at("tensile_fail_curr") = 1; + } + else (*state_vars).at("tensile_fail_curr") = 0; + + if (yield_type == mpm::mohrcoulomb::FailureState::Shear){ yield = yield_function(1); + (*state_vars).at("shear_fail_curr") = 1; + } + else (*state_vars).at("shear_fail_curr") = 0; + + (*state_vars).at("tensile_fail") = (*state_vars).at("tensile_fail") || (*state_vars).at("tensile_fail_curr"); + (*state_vars).at("shear_fail") = (*state_vars).at("shear_fail") || (*state_vars).at("shear_fail_curr"); + // Compute plastic multiplier based on stress input (Lambda) double softening = 0.; double dp_dq = 0.; @@ -461,5 +484,17 @@ Eigen::Matrix mpm::MohrCoulomb::compute_stress( // Update plastic deviatoric strain (*state_vars).at("pdstrain") += dpdstrain; + // Update elastic energy + Vector6d Id({1, 1, 1, 0, 0, 0}); // Identity in Voigt notation + Vector6d dstress = updated_stress - stress; // Stress increment + double tr_dsig = dstress(0) + dstress(1) + dstress(2); // Trace of the stress increment tensor + Vector6d deps_el = ((1+poisson_ratio_)*dstress - poisson_ratio_*tr_dsig*Id) / youngs_modulus_; // Elastic deformation increment (only true strain coefficients) + Vector6d strain_coefs({1, 1, 1, 2, 2, 2}); // coefficients necessary below + (*state_vars).at("E_el") += updated_stress.cwiseProduct(deps_el.cwiseProduct(strain_coefs)).sum()*ptr->volume(); // with 2*sigma_xy*eps_xy + .. on the out-of-diagonal terms, this is the correct expression sigma_ij:depsilon^el_ij (times volume) + + // Update plastic work + Vector6d deps_pl = dstrain - deps_el.cwiseProduct(strain_coefs); // Plastic deformation increment in Voigt notation (with engineering strain coefficents) + (*state_vars).at("W_pl") += updated_stress.cwiseProduct(deps_pl).sum()*ptr->volume(); + return updated_stress; } From f875e91fcef9c718232587011dd0aca6729554f5 Mon Sep 17 00:00:00 2001 From: jduriez Date: Thu, 7 Dec 2023 14:49:14 +0100 Subject: [PATCH 2/6] A less-demanding way of initializing Vector6d --- include/materials/mohr_coulomb.tcc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/materials/mohr_coulomb.tcc b/include/materials/mohr_coulomb.tcc index f108fa7c1..2b189b535 100644 --- a/include/materials/mohr_coulomb.tcc +++ b/include/materials/mohr_coulomb.tcc @@ -485,11 +485,13 @@ Eigen::Matrix mpm::MohrCoulomb::compute_stress( (*state_vars).at("pdstrain") += dpdstrain; // Update elastic energy - Vector6d Id({1, 1, 1, 0, 0, 0}); // Identity in Voigt notation + Vector6d Id; // Identity in Voigt notation + Id << 1, 1, 1, 0, 0, 0; // the direct brace initialization would require a recent (e.g., Ubuntu 22.04) environment Vector6d dstress = updated_stress - stress; // Stress increment double tr_dsig = dstress(0) + dstress(1) + dstress(2); // Trace of the stress increment tensor Vector6d deps_el = ((1+poisson_ratio_)*dstress - poisson_ratio_*tr_dsig*Id) / youngs_modulus_; // Elastic deformation increment (only true strain coefficients) - Vector6d strain_coefs({1, 1, 1, 2, 2, 2}); // coefficients necessary below + Vector6d strain_coefs; // coefficients necessary below + strain_coefs << 1, 1, 1, 2, 2, 2; (*state_vars).at("E_el") += updated_stress.cwiseProduct(deps_el.cwiseProduct(strain_coefs)).sum()*ptr->volume(); // with 2*sigma_xy*eps_xy + .. on the out-of-diagonal terms, this is the correct expression sigma_ij:depsilon^el_ij (times volume) // Update plastic work From f631f0b0c3416251d06b8086a92d833eda85d7e0 Mon Sep 17 00:00:00 2001 From: jduriez Date: Thu, 7 Dec 2023 15:15:03 +0100 Subject: [PATCH 3/6] Considering new MC state variables when testing MC --- tests/materials/mohr_coulomb_test.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/materials/mohr_coulomb_test.cc b/tests/materials/mohr_coulomb_test.cc index 3c69e8bb6..fc488ed6e 100644 --- a/tests/materials/mohr_coulomb_test.cc +++ b/tests/materials/mohr_coulomb_test.cc @@ -118,7 +118,8 @@ TEST_CASE("MohrCoulomb is checked in 2D (cohesion only, without softening)", REQUIRE(state_variables.at("pdstrain") == Approx(0.).epsilon(Tolerance)); const std::vector state_vars = { - "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain"}; + "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain", "E_el", + "W_pl", "tensile_fail_curr", "shear_fail_curr", "tensile_fail", "shear_fail"}; auto state_vars_test = material->state_variables(); REQUIRE(state_vars == state_vars_test); } From 2cf71b35128dbd297d58d4178a066cefe89a0ca7 Mon Sep 17 00:00:00 2001 From: jduriez Date: Thu, 7 Dec 2023 15:23:06 +0100 Subject: [PATCH 4/6] Applying clang format after updating mohr_coulomb_test.cc --- tests/materials/mohr_coulomb_test.cc | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tests/materials/mohr_coulomb_test.cc b/tests/materials/mohr_coulomb_test.cc index fc488ed6e..efa7c46fa 100644 --- a/tests/materials/mohr_coulomb_test.cc +++ b/tests/materials/mohr_coulomb_test.cc @@ -117,9 +117,19 @@ TEST_CASE("MohrCoulomb is checked in 2D (cohesion only, without softening)", REQUIRE(state_variables.at("theta") == Approx(0.).epsilon(Tolerance)); REQUIRE(state_variables.at("pdstrain") == Approx(0.).epsilon(Tolerance)); - const std::vector state_vars = { - "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain", "E_el", - "W_pl", "tensile_fail_curr", "shear_fail_curr", "tensile_fail", "shear_fail"}; + const std::vector state_vars = {"phi", + "psi", + "cohesion", + "epsilon", + "rho", + "theta", + "pdstrain", + "E_el", + "W_pl", + "tensile_fail_curr", + "shear_fail_curr", + "tensile_fail", + "shear_fail"}; auto state_vars_test = material->state_variables(); REQUIRE(state_vars == state_vars_test); } From f4d183b3c2aebb9d16c7cb2bb1ee6c032a781491 Mon Sep 17 00:00:00 2001 From: jduriez Date: Thu, 14 Dec 2023 14:33:02 +0100 Subject: [PATCH 5/6] New MC state variables e_el and w_pl renamed (as such) and included in regression tests --- include/materials/mohr_coulomb.tcc | 14 +++++++------- tests/materials/mohr_coulomb_test.cc | 27 +++++++++++++++++++++++++-- 2 files changed, 32 insertions(+), 9 deletions(-) diff --git a/include/materials/mohr_coulomb.tcc b/include/materials/mohr_coulomb.tcc index 2b189b535..46bc83353 100644 --- a/include/materials/mohr_coulomb.tcc +++ b/include/materials/mohr_coulomb.tcc @@ -71,10 +71,10 @@ mpm::dense_map mpm::MohrCoulomb::initialise_state_variables() { {"theta", 0.}, // Plastic deviatoric strain {"pdstrain", 0.}, - // Elastic energy - {"E_el", 0.}, - // Plastic work - {"W_pl", 0.}, + // Elastic energy (volume-specific) + {"e_el", 0.}, + // Plastic work (volume-specific) + {"w_pl", 0.}, // Flag for current failure because of tensile stress {"tensile_fail_curr", 0}, // Flag for current failure because of shear stress @@ -90,7 +90,7 @@ mpm::dense_map mpm::MohrCoulomb::initialise_state_variables() { template std::vector mpm::MohrCoulomb::state_variables() const { const std::vector state_vars = { - "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain", "E_el", "W_pl", "tensile_fail_curr", "shear_fail_curr", "tensile_fail", "shear_fail"}; + "phi", "psi", "cohesion", "epsilon", "rho", "theta", "pdstrain", "e_el", "w_pl", "tensile_fail_curr", "shear_fail_curr", "tensile_fail", "shear_fail"}; return state_vars; } @@ -492,11 +492,11 @@ Eigen::Matrix mpm::MohrCoulomb::compute_stress( Vector6d deps_el = ((1+poisson_ratio_)*dstress - poisson_ratio_*tr_dsig*Id) / youngs_modulus_; // Elastic deformation increment (only true strain coefficients) Vector6d strain_coefs; // coefficients necessary below strain_coefs << 1, 1, 1, 2, 2, 2; - (*state_vars).at("E_el") += updated_stress.cwiseProduct(deps_el.cwiseProduct(strain_coefs)).sum()*ptr->volume(); // with 2*sigma_xy*eps_xy + .. on the out-of-diagonal terms, this is the correct expression sigma_ij:depsilon^el_ij (times volume) + (*state_vars).at("e_el") += updated_stress.cwiseProduct(deps_el.cwiseProduct(strain_coefs)).sum(); // with 2*sigma_xy*eps_xy + .. on the out-of-diagonal terms, this is the correct expression sigma_ij:depsilon^el_ij // Update plastic work Vector6d deps_pl = dstrain - deps_el.cwiseProduct(strain_coefs); // Plastic deformation increment in Voigt notation (with engineering strain coefficents) - (*state_vars).at("W_pl") += updated_stress.cwiseProduct(deps_pl).sum()*ptr->volume(); + (*state_vars).at("w_pl") += updated_stress.cwiseProduct(deps_pl).sum(); return updated_stress; } diff --git a/tests/materials/mohr_coulomb_test.cc b/tests/materials/mohr_coulomb_test.cc index efa7c46fa..b9c5eb646 100644 --- a/tests/materials/mohr_coulomb_test.cc +++ b/tests/materials/mohr_coulomb_test.cc @@ -116,6 +116,12 @@ TEST_CASE("MohrCoulomb is checked in 2D (cohesion only, without softening)", REQUIRE(state_variables.at("rho") == Approx(0.).epsilon(Tolerance)); REQUIRE(state_variables.at("theta") == Approx(0.).epsilon(Tolerance)); REQUIRE(state_variables.at("pdstrain") == Approx(0.).epsilon(Tolerance)); + REQUIRE(state_variables.at("e_el") == 0.); + REQUIRE(state_variables.at("w_pl") == 0.); + REQUIRE(state_variables.at("tensile_fail_curr") == 0); + REQUIRE(state_variables.at("tensile_fail") == 0); + REQUIRE(state_variables.at("shear_fail_curr") == 0); + REQUIRE(state_variables.at("shear_fail") == 0); const std::vector state_vars = {"phi", "psi", @@ -124,8 +130,8 @@ TEST_CASE("MohrCoulomb is checked in 2D (cohesion only, without softening)", "rho", "theta", "pdstrain", - "E_el", - "W_pl", + "e_el", + "w_pl", "tensile_fail_curr", "shear_fail_curr", "tensile_fail", @@ -566,6 +572,11 @@ TEST_CASE("MohrCoulomb is checked in 2D (cohesion only, without softening)", REQUIRE(updated_stress(3) == Approx(0.).epsilon(Tolerance)); REQUIRE(updated_stress(4) == Approx(0.).epsilon(Tolerance)); REQUIRE(updated_stress(5) == Approx(0.).epsilon(Tolerance)); + // Check whether failure flag state variables (after being computed in compute_stress) are correct + REQUIRE(state_variables.at("shear_fail") == 1); + REQUIRE(state_variables.at("shear_fail_curr") == 1); + REQUIRE(state_variables.at("tensile_fail") == 0); + REQUIRE(state_variables.at("tensile_fail_curr") == 0); } } } @@ -2866,6 +2877,10 @@ TEST_CASE("MohrCoulomb is checked in 3D (c & phi & psi, with softening)", // Check plastic strain REQUIRE(state_variables.at("pdstrain") == Approx(0.0007751567).epsilon(Tolerance)); + // In the present workflow, all failure flags remain at 0 ?... + // Check plastic work and elastic energy + REQUIRE(state_variables.at("e_el") == Approx(15.33614017380).epsilon(Tolerance)); + REQUIRE(state_variables.at("w_pl") == Approx(4.5135130001).epsilon(Tolerance)); } //! Check for shear failure (pdstrain > pdstrain_residual) @@ -3373,6 +3388,14 @@ TEST_CASE("MohrCoulomb is checked in 3D (c & phi & psi, with softening)", // Check plastic strain REQUIRE(state_variables.at("pdstrain") == Approx(0.0002112522).epsilon(Tolerance)); + // Check failure flags + REQUIRE(state_variables.at("shear_fail") == 1); + REQUIRE(state_variables.at("shear_fail_curr") == 1); + REQUIRE(state_variables.at("tensile_fail") == 0); + REQUIRE(state_variables.at("tensile_fail_curr") == 0); + // Check plastic work and elastic energy + REQUIRE(state_variables.at("e_el") == Approx(1.153960429).epsilon(Tolerance)); + REQUIRE(state_variables.at("w_pl") == Approx(1.0964114497).epsilon(Tolerance)); } //! Check for shear failure (pdstrain > pdstrain_residual) From bc7cfa98d7ce0b84ac873ecd8c0cc2ca447bc540 Mon Sep 17 00:00:00 2001 From: jduriez Date: Thu, 14 Dec 2023 14:50:17 +0100 Subject: [PATCH 6/6] Applying clang format after -1.. --- tests/materials/mohr_coulomb_test.cc | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/tests/materials/mohr_coulomb_test.cc b/tests/materials/mohr_coulomb_test.cc index b9c5eb646..8973f7d07 100644 --- a/tests/materials/mohr_coulomb_test.cc +++ b/tests/materials/mohr_coulomb_test.cc @@ -572,7 +572,8 @@ TEST_CASE("MohrCoulomb is checked in 2D (cohesion only, without softening)", REQUIRE(updated_stress(3) == Approx(0.).epsilon(Tolerance)); REQUIRE(updated_stress(4) == Approx(0.).epsilon(Tolerance)); REQUIRE(updated_stress(5) == Approx(0.).epsilon(Tolerance)); - // Check whether failure flag state variables (after being computed in compute_stress) are correct + // Check whether failure flag state variables (after being computed in + // compute_stress) are correct REQUIRE(state_variables.at("shear_fail") == 1); REQUIRE(state_variables.at("shear_fail_curr") == 1); REQUIRE(state_variables.at("tensile_fail") == 0); @@ -2879,8 +2880,10 @@ TEST_CASE("MohrCoulomb is checked in 3D (c & phi & psi, with softening)", Approx(0.0007751567).epsilon(Tolerance)); // In the present workflow, all failure flags remain at 0 ?... // Check plastic work and elastic energy - REQUIRE(state_variables.at("e_el") == Approx(15.33614017380).epsilon(Tolerance)); - REQUIRE(state_variables.at("w_pl") == Approx(4.5135130001).epsilon(Tolerance)); + REQUIRE(state_variables.at("e_el") == + Approx(15.33614017380).epsilon(Tolerance)); + REQUIRE(state_variables.at("w_pl") == + Approx(4.5135130001).epsilon(Tolerance)); } //! Check for shear failure (pdstrain > pdstrain_residual) @@ -3394,8 +3397,10 @@ TEST_CASE("MohrCoulomb is checked in 3D (c & phi & psi, with softening)", REQUIRE(state_variables.at("tensile_fail") == 0); REQUIRE(state_variables.at("tensile_fail_curr") == 0); // Check plastic work and elastic energy - REQUIRE(state_variables.at("e_el") == Approx(1.153960429).epsilon(Tolerance)); - REQUIRE(state_variables.at("w_pl") == Approx(1.0964114497).epsilon(Tolerance)); + REQUIRE(state_variables.at("e_el") == + Approx(1.153960429).epsilon(Tolerance)); + REQUIRE(state_variables.at("w_pl") == + Approx(1.0964114497).epsilon(Tolerance)); } //! Check for shear failure (pdstrain > pdstrain_residual)