From 4b66481e33164375bbb789325de644d671812135 Mon Sep 17 00:00:00 2001 From: Theodore Date: Mon, 30 Sep 2024 00:15:53 +0200 Subject: [PATCH] 3D version of subloading surface model for metal (#221) --- CHANGELOG.md | 1 + Example/Material/SubloadingMetal.supan | 50 +++++ Example/Material/SubloadingMetal2.supan | 5 + MSVC/suanPan/suanPan/suanPan.vcxproj | 2 + MSVC/suanPan/suanPan/suanPan.vcxproj.filters | 6 + MSVC/suanPan/suanPan/suanPan.vcxproj.user | 2 +- Material/Material1D/vonMises/Subloading1D.cpp | 3 +- Material/Material3D/CMakeLists.txt | 1 + Material/Material3D/Material3D | 1 + .../Material3D/vonMises/SubloadingMetal.cpp | 182 ++++++++++++++++++ .../Material3D/vonMises/SubloadingMetal.h | 84 ++++++++ Material/MaterialParser.cpp | 23 +++ Toolbox/tensor.cpp | 4 + Toolbox/tensor.h | 2 + 14 files changed, 364 insertions(+), 2 deletions(-) create mode 100644 Example/Material/SubloadingMetal.supan create mode 100644 Example/Material/SubloadingMetal2.supan create mode 100644 Material/Material3D/vonMises/SubloadingMetal.cpp create mode 100644 Material/Material3D/vonMises/SubloadingMetal.h diff --git a/CHANGELOG.md b/CHANGELOG.md index 320828b0d..67c3976f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ 1. add `AFCO1D` material with strain memory [#217](https://github.com/TLCFEM/suanPan/pull/217) 2. update `Catch2` to version `3.7.1` 3. add `Subloading1D` material [#219](https://github.com/TLCFEM/suanPan/pull/219) +4. add `SubloadingMetal` material [#221](https://github.com/TLCFEM/suanPan/pull/221) ## version 3.5 diff --git a/Example/Material/SubloadingMetal.supan b/Example/Material/SubloadingMetal.supan new file mode 100644 index 000000000..af3b82e4e --- /dev/null +++ b/Example/Material/SubloadingMetal.supan @@ -0,0 +1,50 @@ +# A TEST MODEL FOR SUBLOADINGMETAL MATERIAL + +node 1 5 -5 0 +node 2 5 5 0 +node 3 -5 5 0 +node 4 -5 -5 0 +node 5 5 -5 10 +node 6 5 5 10 +node 7 -5 5 10 +node 8 -5 -5 10 +node 9 5 -5 20 +node 10 5 5 20 +node 11 -5 5 20 +node 12 -5 -5 20 + +material SubloadingMetal 1 + +element C3D8 1 1 2 3 4 5 6 7 8 1 G +element C3D8 2 5 6 7 8 9 10 11 12 1 G + +fix 1 1 1 2 5 6 9 10 +fix 2 2 1 4 5 8 9 12 +fix 3 3 1 2 3 4 + +displacement 1 0 -.1 3 11 12 +displacement 2 0 -.12 3 9 10 + +step static 1 +set fixed_step_size 1 +set ini_step_size 1E-2 +set symm_mat 0 + +converger RelIncreDisp 1 1E-10 10 1 + +analyze + +peek element 1 + +# Node 9: +# Coordinate: +# 5.0000e+00 -5.0000e+00 2.0000e+01 +# Displacement: +# 0.0000e+00 0.0000e+00 -1.2000e-01 +# Resistance: +# -2.7793e+02 1.2861e+02 -1.7073e+03 +peek node 9 + +reset +clear +exit \ No newline at end of file diff --git a/Example/Material/SubloadingMetal2.supan b/Example/Material/SubloadingMetal2.supan new file mode 100644 index 000000000..598923d6b --- /dev/null +++ b/Example/Material/SubloadingMetal2.supan @@ -0,0 +1,5 @@ +material SubloadingMetal 1 + +materialTestByLoad3D 1 1 0 0 0 0 0 100 50 50 50 + +exit \ No newline at end of file diff --git a/MSVC/suanPan/suanPan/suanPan.vcxproj b/MSVC/suanPan/suanPan/suanPan.vcxproj index 038003777..a6ab82898 100644 --- a/MSVC/suanPan/suanPan/suanPan.vcxproj +++ b/MSVC/suanPan/suanPan/suanPan.vcxproj @@ -418,6 +418,7 @@ + @@ -910,6 +911,7 @@ + diff --git a/MSVC/suanPan/suanPan/suanPan.vcxproj.filters b/MSVC/suanPan/suanPan/suanPan.vcxproj.filters index 3cb2d2024..fe517a2e3 100644 --- a/MSVC/suanPan/suanPan/suanPan.vcxproj.filters +++ b/MSVC/suanPan/suanPan/suanPan.vcxproj.filters @@ -1588,6 +1588,9 @@ Material\Material1D\vonMises + + Material\Material3D\vonMises + @@ -3040,6 +3043,9 @@ Material\Material1D\vonMises + + Material\Material3D\vonMises + diff --git a/MSVC/suanPan/suanPan/suanPan.vcxproj.user b/MSVC/suanPan/suanPan/suanPan.vcxproj.user index 7d3074fee..a9b0d9d91 100644 --- a/MSVC/suanPan/suanPan/suanPan.vcxproj.user +++ b/MSVC/suanPan/suanPan/suanPan.vcxproj.user @@ -4,7 +4,7 @@ false - -f Subloading1D + -f SubloadingMetal ..\..\..\.dirty WindowsLocalDebugger diff --git a/Material/Material1D/vonMises/Subloading1D.cpp b/Material/Material1D/vonMises/Subloading1D.cpp index 3302cfef6..529d83c89 100644 --- a/Material/Material1D/vonMises/Subloading1D.cpp +++ b/Material/Material1D/vonMises/Subloading1D.cpp @@ -70,6 +70,8 @@ int Subloading1D::update_trial_status(const vec& t_strain) { return SUANPAN_FAIL; } + q = current_q + gamma; + const auto exp_iso = saturation_iso * exp(-m_iso * q); auto y = initial_iso + saturation_iso + k_iso * q - exp_iso; auto dy = k_iso + m_iso * exp_iso; @@ -133,7 +135,6 @@ int Subloading1D::update_trial_status(const vec& t_strain) { gamma -= incre(0); z -= incre(1); - q = current_q + gamma; } } diff --git a/Material/Material3D/CMakeLists.txt b/Material/Material3D/CMakeLists.txt index 133bc31cd..2156dc615 100644 --- a/Material/Material3D/CMakeLists.txt +++ b/Material/Material3D/CMakeLists.txt @@ -44,6 +44,7 @@ set(M3D Material3D/vonMises/NonlinearJ2.cpp Material3D/vonMises/NonlinearPeric.cpp Material3D/vonMises/PolyJ2.cpp + Material3D/vonMises/SubloadingMetal.cpp Material3D/vonMises/TableGurson.cpp Material3D/vonMises/VAFCRP.cpp Material3D/Wrapper/Rotation3D.cpp diff --git a/Material/Material3D/Material3D b/Material/Material3D/Material3D index cf6507a2a..74364d859 100644 --- a/Material/Material3D/Material3D +++ b/Material/Material3D/Material3D @@ -43,6 +43,7 @@ #include "vonMises/NonlinearJ2.h" #include "vonMises/NonlinearPeric.h" #include "vonMises/PolyJ2.h" +#include "vonMises/SubloadingMetal.h" #include "vonMises/TableGurson.h" #include "vonMises/VAFCRP.h" #include "Wrapper/Rotation3D.h" diff --git a/Material/Material3D/vonMises/SubloadingMetal.cpp b/Material/Material3D/vonMises/SubloadingMetal.cpp new file mode 100644 index 000000000..1e8f507d1 --- /dev/null +++ b/Material/Material3D/vonMises/SubloadingMetal.cpp @@ -0,0 +1,182 @@ +/******************************************************************************* + * Copyright (C) 2017-2024 Theodore Chang + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + ******************************************************************************/ + +#include "SubloadingMetal.h" +#include + +const double SubloadingMetal::root_two_third = sqrt(two_third); +const double SubloadingMetal::rate_bound = -log(z_bound); +const mat SubloadingMetal::unit_dev_tensor = tensor::unit_deviatoric_tensor4(); + +vec2 SubloadingMetal::yield_ratio(const double z) { + if(z < z_bound) return {rate_bound, 0.}; + + return {-log(z), -1. / z}; +} + +SubloadingMetal::SubloadingMetal(const unsigned T, DataSubloadingMetal&& D, const double R) + : DataSubloadingMetal{std::move(D)} + , Material3D(T, R) {} + +int SubloadingMetal::initialize(const shared_ptr&) { + trial_stiffness = current_stiffness = initial_stiffness = tensor::isotropic_stiffness(elastic, poissons_ratio); + + initialize_history(14); + + return SUANPAN_SUCCESS; +} + +unique_ptr SubloadingMetal::get_copy() { return make_unique(*this); } + +int SubloadingMetal::update_trial_status(const vec& t_strain) { + incre_strain = (trial_strain = t_strain) - current_strain; + + if(norm(incre_strain) <= datum::eps) return SUANPAN_SUCCESS; + + trial_stress = current_stress + (trial_stiffness = initial_stiffness) * incre_strain; + + trial_history = current_history; + const auto& current_q = current_history(0); + const auto& current_z = current_history(1); + const vec current_alpha(¤t_history(2), 6, false, true); + const vec current_d(¤t_history(8), 6, false, true); + auto& q = trial_history(0); + auto& z = trial_history(1); + vec alpha(&trial_history(2), 6, false, true); + vec d(&trial_history(8), 6, false, true); + + const vec trial_s = tensor::dev(trial_stress); + + auto gamma = 0., ref_error = 0.; + + vec2 residual, incre; + mat22 jacobian; + + const auto current_ratio = yield_ratio(current_z); + + auto counter = 0u; + while(true) { + if(max_iteration == ++counter) { + suanpan_error("Cannot converge within {} iterations.\n", max_iteration); + return SUANPAN_FAIL; + } + + q = current_q + root_two_third * gamma; + + const auto exp_iso = saturation_iso * exp(-m_iso * q); + auto y = initial_iso + saturation_iso + k_iso * q - exp_iso; + auto dy = root_two_third * (k_iso + m_iso * exp_iso); + if(y < 0.) y = dy = 0.; + + const auto exp_kin = saturation_kin * exp(-m_kin * q); + auto a = initial_kin + saturation_kin + k_kin * q - exp_kin; + auto da = root_two_third * (k_kin + m_kin * exp_kin); + if(a < 0.) a = da = 0.; + + const auto bot_alpha = 1. + root_two_third * be * gamma; + const auto bot_d = 1. + root_two_third * ce * gamma; + + const vec pzetapz = current_d / bot_d; + const vec pzetapgamma = root_two_third * (be / bot_alpha / bot_alpha * current_alpha + ce * (1. - z) / bot_d * pzetapz); + + const vec zeta = trial_s - current_alpha / bot_alpha + (z - 1.) * pzetapz; + const auto norm_zeta = tensor::stress::norm(zeta); + const vec n = zeta / norm_zeta; + + vec top_alpha = gamma * be * a * n + current_alpha; + alpha = top_alpha / bot_alpha; + + vec top_d = gamma * ce * ze * y * n + current_d; + d = top_d / bot_d; + + const vec eta = trial_s - gamma * double_shear * n - alpha + (z - 1.) * d; + + residual(0) = tensor::stress::norm(eta) - root_two_third * z * y; + + if(1u == counter && residual(0) < 0.) { + const auto aa = tensor::stress::double_contraction(d) - two_third * y * y; + const auto bb = -tensor::stress::double_contraction(d, eta); + const auto cc = tensor::stress::double_contraction(eta); + const auto sqrt_term = sqrt(bb * bb - aa * cc); + + z = (bb - sqrt_term) / aa; + if(z < 0. || z > 1.) z = (bb + sqrt_term) / aa; + + return SUANPAN_SUCCESS; + } + + const auto trial_ratio = yield_ratio(z); + const auto avg_rate = u * .5 * (current_ratio(0) + trial_ratio(0)); + + residual(1) = z - current_z - root_two_third * gamma * avg_rate; + + jacobian(0, 0) = tensor::stress::double_contraction(n, pzetapgamma) - double_shear - be / bot_alpha * (a + gamma * da - gamma * a * root_two_third * be / bot_alpha) + (z - 1.) * ce * ze / bot_d * (y + gamma * dy - gamma * y * root_two_third * ce / bot_d) - root_two_third * z * dy; + jacobian(0, 1) = tensor::stress::double_contraction(n, pzetapz) + gamma * ce * ze * y / bot_d - root_two_third * y; + + jacobian(1, 0) = -root_two_third * avg_rate; + jacobian(1, 1) = 1. - root_two_third * gamma * u * .5 * trial_ratio(1); + + if(!solve(incre, jacobian, residual)) return SUANPAN_FAIL; + + const auto error = inf_norm(incre); + if(1u == counter) ref_error = error; + suanpan_debug("Local iteration error: {:.5E}.\n", error); + if(error < tolerance * ref_error || ((error < tolerance || inf_norm(residual) < tolerance) && counter > 5u)) { + trial_stress -= gamma * double_shear * n; + + mat right(2, 6, fill::zeros); + right.row(0) = -n.t(); + + const mat left = solve(jacobian, right); + + trial_stiffness -= double_shear * double_shear * n * left.row(0) + double_shear * double_shear * gamma / norm_zeta * (unit_dev_tensor - n * n.t()); + + return SUANPAN_SUCCESS; + } + + gamma -= incre(0); + z -= incre(1); + } +} + +int SubloadingMetal::clear_status() { + current_strain.zeros(); + current_stress.zeros(); + current_history = initial_history; + current_stiffness = initial_stiffness; + return reset_status(); +} + +int SubloadingMetal::commit_status() { + current_strain = trial_strain; + current_stress = trial_stress; + current_history = trial_history; + current_stiffness = trial_stiffness; + return SUANPAN_SUCCESS; +} + +int SubloadingMetal::reset_status() { + trial_strain = current_strain; + trial_stress = current_stress; + trial_history = current_history; + trial_stiffness = current_stiffness; + return SUANPAN_SUCCESS; +} + +void SubloadingMetal::print() { + suanpan_info("A 3D combined hardening material using subloading surface model.\n"); +} diff --git a/Material/Material3D/vonMises/SubloadingMetal.h b/Material/Material3D/vonMises/SubloadingMetal.h new file mode 100644 index 000000000..801029b94 --- /dev/null +++ b/Material/Material3D/vonMises/SubloadingMetal.h @@ -0,0 +1,84 @@ +/******************************************************************************* + * Copyright (C) 2017-2024 Theodore Chang + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + ******************************************************************************/ +/** + * @class SubloadingMetal + * @brief A SubloadingMetal material class. + * @author tlc + * @date 29/09/2024 + * @version 0.1.0 + * @file SubloadingMetal.h + * @addtogroup Material-1D + * @{ + */ + +#ifndef SUBLOADINGMETAL_H +#define SUBLOADINGMETAL_H + +#include + +struct DataSubloadingMetal { + const double elastic; // elastic modulus + const double poissons_ratio; + const double initial_iso; + const double k_iso; + const double saturation_iso; + const double m_iso; + const double initial_kin; + const double k_kin; + const double saturation_kin; + const double m_kin; + const double u; + const double be; + const double ce; + const double ze; +}; + +class SubloadingMetal final : protected DataSubloadingMetal, public Material3D { + static constexpr unsigned max_iteration = 20u; + static constexpr double two_third = 2. / 3.; + static const double root_two_third; + static constexpr double z_bound = 1E-7; + static const double rate_bound; + static const mat unit_dev_tensor; + + static vec2 yield_ratio(double); + + const double double_shear = elastic / (1. + poissons_ratio); // double shear modulus + +public: + SubloadingMetal( + unsigned, // tag + DataSubloadingMetal&&, // data + double = 0. // density + ); + + int initialize(const shared_ptr&) override; + + unique_ptr get_copy() override; + + int update_trial_status(const vec&) override; + + int clear_status() override; + int commit_status() override; + int reset_status() override; + + void print() override; +}; + +#endif + +//! @} diff --git a/Material/MaterialParser.cpp b/Material/MaterialParser.cpp index c43238a49..71fe414c8 100644 --- a/Material/MaterialParser.cpp +++ b/Material/MaterialParser.cpp @@ -2951,6 +2951,28 @@ void new_subloading1d(unique_ptr& return_obj, istringstream& command) return_obj = make_unique(tag, DataSubloading1D{p(0), p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8), p(9), p(10), p(11), p(12)}, density); } +void new_subloadingmetal(unique_ptr& return_obj, istringstream& command) { + unsigned tag; + if(!get_input(command, tag)) { + suanpan_error("A valid tag is required.\n"); + return; + } + + vec p{2E5, .3, 2E2, 0., 2E2, 1E1, 2E2, 0., 2E2, 1E1, 1E1, 1E1, 1E1, .7}; + if(!get_optional_input(command, p)) { + suanpan_error("Valid inputs are required.\n"); + return; + } + + auto density = 0.; + if(!command.eof() && !get_input(command, density)) { + suanpan_error("A valid density is required.\n"); + return; + } + + return_obj = make_unique(tag, DataSubloadingMetal{p(0), p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8), p(9), p(10), p(11), p(12), p(13)}, density); +} + void new_substepping(unique_ptr& return_obj, istringstream& command) { unsigned tag; if(!get_input(command, tag)) { @@ -3523,6 +3545,7 @@ int create_new_material(const shared_ptr& domain, istringstream& com else if(is_equal(material_id, "Stacked")) new_stacked(new_material, command); else if(is_equal(material_id, "SteelBRB")) new_steelbrb(new_material, command); else if(is_equal(material_id, "Subloading1D")) new_subloading1d(new_material, command); + else if(is_equal(material_id, "SubloadingMetal")) new_subloadingmetal(new_material, command); else if(is_equal(material_id, "Substepping")) new_substepping(new_material, command); else if(is_equal(material_id, "TableCDP")) new_tablecdp(new_material, command); else if(is_equal(material_id, "TableGurson")) new_tablegurson(new_material, command); diff --git a/Toolbox/tensor.cpp b/Toolbox/tensor.cpp index 2cb68d4f7..0416710da 100644 --- a/Toolbox/tensor.cpp +++ b/Toolbox/tensor.cpp @@ -379,6 +379,8 @@ double tensor::strain::norm(vec&& in) { throw invalid_argument("need a valid strain vector"); } +double tensor::strain::double_contraction(const vec& a) { return double_contraction(a, a); } + double tensor::strain::double_contraction(const vec& a, const vec& b) { return dot(a % b, norm_weight); } double tensor::strain::double_contraction(vec&& a, vec&& b) { return dot(a % b, norm_weight); } @@ -449,6 +451,8 @@ double tensor::stress::norm(vec&& in) { throw invalid_argument("need a valid stress vector"); } +double tensor::stress::double_contraction(const vec& a) { return double_contraction(a, a); } + double tensor::stress::double_contraction(const vec& a, const vec& b) { return dot(a % b, norm_weight); } double tensor::stress::double_contraction(vec&& a, vec&& b) { return dot(a % b, norm_weight); } diff --git a/Toolbox/tensor.h b/Toolbox/tensor.h index 65873b41f..0913c06f0 100644 --- a/Toolbox/tensor.h +++ b/Toolbox/tensor.h @@ -87,6 +87,7 @@ namespace tensor { vec to_voigt(const mat&); double norm(const vec&); double norm(vec&&); + double double_contraction(const vec&); double double_contraction(const vec&, const vec&); double double_contraction(vec&&, vec&&); } // namespace strain @@ -95,6 +96,7 @@ namespace tensor { vec to_voigt(const mat&); double norm(const vec&); double norm(vec&&); + double double_contraction(const vec&); double double_contraction(const vec&, const vec&); double double_contraction(vec&&, vec&&); } // namespace stress