diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h index 7797fd8..67aed65 100644 --- a/include/numsim-materials/materials/drucker_prager_plasticity.h +++ b/include/numsim-materials/materials/drucker_prager_plasticity.h @@ -41,11 +41,9 @@ class drucker_prager_plasticity final m_kappa(base::template add_history_output("equivalent_plastic_strain")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), + m_K_bulk(base::template get_parameter("K_bulk")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), - m_C_e(base::template add_input( - base::template get_parameter("elastic_source"), - "tangent", EdgeKind::Global)), m_strain(base::template add_input( base::template get_parameter("strain_source"), "strain", EdgeKind::Global)), @@ -58,12 +56,15 @@ class drucker_prager_plasticity final { m_yf = yield_fn(base::template get_parameter("eta"), base::template get_parameter("beta"), - base::template get_parameter("K_bulk")); + m_K_bulk); + const auto I = tmech::eye(); + const tensor4 IIvol{tmech::otimes(I, I) / value_type{Dim}}; + m_C_e = value_type{3} * m_K_bulk * IIvol + + value_type{2} * m_G * plasticity_detail::make_IIdev(); } static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source").template add(); para.template insert("hardening_source").template add(); para.template insert("strain_source").template add(); para.template insert("solver_source").template add(); @@ -80,7 +81,7 @@ class drucker_prager_plasticity final } void compute() { - const auto& C_e = m_C_e.get(); + const auto& C_e = m_C_e; const auto kappa_n = m_kappa.old_value(); m_kappa.new_value() = kappa_n; @@ -151,7 +152,16 @@ class drucker_prager_plasticity final const tensor4& C_e, value_type kappa_n, value_type dlambda) { m_eps_p.new_value() = m_eps_p.old_value() + dlambda * ts.N; m_kappa.new_value() = kappa_n + dlambda; - m_stress = tmech::dcontract(C_e, m_strain.get() - m_eps_p.new_value()); + + // sigma = C_e : (eps - eps_p_new) = sig_trial - dl (C_e : N). + // For isotropic C_e that is 2G dev(N) + K tr(N) I -- no rank-4 : rank-2 + // contraction, which measures 33 ns against ~720 for the whole step. The + // flow is non-associative, so tr(N) = beta is generally nonzero and the + // volumetric term does not drop out as it does for J2. + const auto I = tmech::eye(); + const tensor2 Ce_N{value_type{2} * m_G * tmech::dev(ts.N) + + m_K_bulk * tmech::trace(ts.N) * I}; + m_stress = ts.sig - dlambda * Ce_N; // Tangent at trial state (return mapping uses N_trial). // For J2, trial = converged. For DP, they differ. @@ -187,13 +197,20 @@ class drucker_prager_plasticity final const value_type& m_G; const value_type& m_sigma_0; + const value_type& m_K_bulk; material_ref& m_solver; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; yield_fn m_yf{}; + + /// The elastic stiffness, built here rather than read from another material. + /// The tangent collapse requires an isotropic C_e -- C_e : X = 2G X for + /// deviatoric X -- so accepting an arbitrary rank-4 tangent advertised a + /// generality this material cannot honour. K_bulk and G are already + /// parameters, so nothing new is asked of the caller. + tensor4 m_C_e{}; }; diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h index 7322fd4..dee151c 100644 --- a/include/numsim-materials/materials/j2_plasticity.h +++ b/include/numsim-materials/materials/j2_plasticity.h @@ -45,9 +45,10 @@ namespace numsim::materials { /// may be nonlinear (see exponential_isotropic_hardening). /// /// Parameters: -/// "name", "elastic_source", "hardening_source", "strain_source", -/// "solver_source", "G", "sigma_0" -/// -- the same set small_strain_plasticity takes, so this is a drop-in. +/// "name", "hardening_source", "strain_source", "solver_source", +/// "K", "G", "sigma_0" +/// +/// No elastic_source: the stiffness is built from K and G here. template class j2_plasticity final : public material_base, Traits> { @@ -71,11 +72,9 @@ class j2_plasticity final "equivalent_plastic_strain")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), + m_K(base::template get_parameter("K")), m_solver(base::template add_material_ref( base::template get_parameter("solver_source"))), - m_C_e(base::template add_input( - base::template get_parameter("elastic_source"), - "tangent", EdgeKind::Global)), m_strain(base::template add_input( base::template get_parameter("strain_source"), "strain", EdgeKind::Global)), @@ -85,26 +84,41 @@ class j2_plasticity final m_dH(base::template add_input( base::template get_parameter("hardening_source"), "hardening_modulus", EdgeKind::Local)), - m_IIdev(plasticity_detail::make_IIdev()) + m_IIdev(plasticity_detail::make_IIdev()), + m_C_e(build_elastic_tangent(m_K, m_G, m_IIdev)) {} + /// The elastic stiffness, built here rather than read from another material. + /// + /// The closed forms below REQUIRE an isotropic C_e -- that is what makes + /// C_e : N = 2G N and N : C_e : N = 3G true. Accepting an arbitrary rank-4 + /// tangent from an elastic_source advertised a generality this material + /// cannot honour, and dragged a linear_elasticity into every plasticity graph + /// whose own "stress" output (C : eps, ignoring eps_p) is meaningless once + /// yielding starts. + static tensor4 build_elastic_tangent(value_type K, value_type G, + const tensor4& IIdev) { + const auto I = tmech::eye(); + return value_type{3} * K * (tmech::otimes(I, I) / value_type{Dim}) + + value_type{2} * G * IIdev; + } + static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source") - .template add(); para.template insert("hardening_source") .template add(); para.template insert("strain_source") .template add(); para.template insert("solver_source") .template add(); + para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); return para; } void compute() { - const auto& C_e = m_C_e.get(); + const auto& C_e = m_C_e; const auto kappa_n = m_kappa.old_value(); m_kappa.new_value() = kappa_n; @@ -143,8 +157,13 @@ class j2_plasticity final m_eps_p.new_value() = m_eps_p.old_value() + dlambda * N; m_kappa.new_value() = kappa_n + dlambda; - m_stress = tmech::dcontract( - C_e, tensor2(m_strain.get() - m_eps_p.new_value())); + + // sigma = C_e : (eps - eps_p_new) + // = C_e : (eps - eps_p_old) - dl (C_e : N) + // = sig_trial - 2G dl N, since C_e : N = 2G N. + // The same identity the tangent uses. Avoids a second rank-4 : rank-2 + // contraction, which measures 33 ns against ~245 for the whole step. + m_stress = sig_trial - value_type{2} * m_G * dlambda * N; m_H.update_source(); const auto GG = m_G * m_G; @@ -162,14 +181,15 @@ class j2_plasticity final const value_type& m_G; const value_type& m_sigma_0; + const value_type& m_K; material_ref& m_solver; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; const tensor4 m_IIdev; + const tensor4 m_C_e; }; } // namespace numsim::materials diff --git a/include/numsim-materials/materials/rk_plasticity.h b/include/numsim-materials/materials/rk_plasticity.h index cf5b0b1..05d40bd 100644 --- a/include/numsim-materials/materials/rk_plasticity.h +++ b/include/numsim-materials/materials/rk_plasticity.h @@ -37,16 +37,14 @@ class rk_plasticity final m_tangent(base::template add_output("tangent")), m_eps_p(base::template add_history_output("plastic_strain")), m_kappa(base::template add_history_output("equivalent_plastic_strain")), + m_K(base::template get_parameter("K")), m_G(base::template get_parameter("G")), m_sigma_0(base::template get_parameter("sigma_0")), m_tol(base::template get_parameter("tolerance")), m_max_iter(base::template get_parameter("max_iter")), m_tableau(base::template get_parameter("tableau")), - m_elastic_source(base::template get_parameter("elastic_source")), m_hardening_source(base::template get_parameter("hardening_source")), m_strain_source(base::template get_parameter("strain_source")), - m_C_e(base::template add_input( - m_elastic_source, "tangent", EdgeKind::Global)), m_strain(base::template add_input( m_strain_source, "strain", EdgeKind::Global)), m_H(base::template add_input( @@ -54,6 +52,11 @@ class rk_plasticity final m_dH(base::template add_input( m_hardening_source, "hardening_modulus", EdgeKind::Local)) { + const auto I = tmech::eye(); + const tensor4 IIvol{tmech::otimes(I, I) / value_type{Dim}}; + m_C_e = value_type{3} * m_K * IIvol + + value_type{2} * m_G * plasticity_detail::make_IIdev(); + const int s = m_tableau->stages(); m_dlambda.resize(s, value_type{0}); m_N_stage.resize(s); @@ -67,9 +70,9 @@ class rk_plasticity final static input_parameter_controller parameters() { input_parameter_controller para{base::parameters()}; - para.template insert("elastic_source").template add(); para.template insert("hardening_source").template add(); para.template insert("strain_source").template add(); + para.template insert("K").template add(); para.template insert("G").template add(); para.template insert("sigma_0").template add(); para.template insert("tolerance") @@ -80,7 +83,7 @@ class rk_plasticity final } void compute() { - const auto& C_e = m_C_e.get(); + const auto& C_e = m_C_e; const auto& eps = m_strain.get(); const auto kappa_n = m_kappa.old_value(); const auto eps_p_n = m_eps_p.old_value(); @@ -181,20 +184,26 @@ class rk_plasticity final history_property& m_eps_p; history_property& m_kappa; + const value_type& m_K; const value_type& m_G; const value_type& m_sigma_0; const value_type& m_tol; const int& m_max_iter; const butcher_tableau* m_tableau; - const std::string& m_elastic_source; const std::string& m_hardening_source; const std::string& m_strain_source; - const input_property& m_C_e; const input_property& m_strain; const input_property& m_H; const input_property& m_dH; + /// Built here, not read from an elastic_source. compute_tangent's collapse + /// (C_e : X = 2G X) and effective_modulus(G) = 3G both require an isotropic + /// C_e, so taking an arbitrary rank-4 tangent promised more than this class + /// can deliver -- and pulled a linear_elasticity into the graph whose own + /// "stress" output is meaningless once eps_p is nonzero. + tensor4 m_C_e{}; + yield_fn m_yf{}; std::vector m_dlambda; std::vector m_N_stage; diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp index b32ed3f..b32db4b 100644 --- a/tests/plot_data.cpp +++ b/tests/plot_data.cpp @@ -78,10 +78,10 @@ run_result run_j2(T increment, int steps, p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", K_val); p.insert("G", G_val); p.insert("sigma_0", sigma_0); ctx.create(p); @@ -144,7 +144,6 @@ run_result run_dp(T increment, int steps, p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp index aa8f4fc..c495c32 100644 --- a/tests/test_drucker_prager.cpp +++ b/tests/test_drucker_prager.cpp @@ -116,7 +116,6 @@ class DruckerPragerTest : public ::testing::Test { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -221,7 +220,6 @@ class DPTangentTest : public ::testing::Test { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -295,7 +293,6 @@ T run_dp_max_tangent_error(T increment, int steps) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -383,7 +380,6 @@ T max_tangent_error(std::vector direction, T increment, int steps) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -481,7 +477,6 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) { p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); @@ -540,7 +535,6 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) { ctx.create>(p); p.clear(); p.insert("name", "dp"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp index 03ff086..a7b3429 100644 --- a/tests/test_j2_plasticity.cpp +++ b/tests/test_j2_plasticity.cpp @@ -56,10 +56,10 @@ class J2PlasticityTest : public ::testing::Test { // J2 plasticity — solver passed as pointer p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); - p.insert("hardening_source", "hardening"); + p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", K); p.insert("G", G); p.insert("sigma_0", sigma_0); ctx.create>(p); @@ -168,10 +168,10 @@ class J2TangentTest : public ::testing::Test { p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); - p.insert("hardening_source", "hardening"); + p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); p.insert("solver_source", "solver"); + p.insert("K", T{166.67}); p.insert("G", T{76.92}); p.insert("sigma_0", T{50.0}); ctx.create>(p); @@ -258,9 +258,9 @@ class RKPlasticityTest : public ::testing::Test { p.clear(); p.insert("name", "j2"); - p.insert("elastic_source", "elastic"); p.insert("hardening_source", "hardening"); p.insert("strain_source", "stepper"); + p.insert("K", T{166.67}); p.insert("G", T{76.92}); p.insert("sigma_0", T{50.0}); p.insert("tableau", &m_tab);