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