From c98741012658c7097ba12e8944759d3a4be58fbd Mon Sep 17 00:00:00 2001 From: Theodore Date: Tue, 29 Oct 2024 04:32:02 +0100 Subject: [PATCH] Subloading surface model with viscosity for metals (#223) --- Example/Material/SubloadingViscous1D.supan | 43 +++++++++++++ Example/Material/SubloadingViscous1D2.supan | 9 +++ MSVC/suanPan/suanPan/suanPan.vcxproj.user | 2 +- Material/Material1D/Material1D | 8 +-- Material/Material1D/vonMises/Subloading1D.cpp | 60 +++++++++++++------ Material/Material1D/vonMises/Subloading1D.h | 20 +++---- Material/MaterialParser.cpp | 45 +++++++++++++- 7 files changed, 152 insertions(+), 35 deletions(-) create mode 100644 Example/Material/SubloadingViscous1D.supan create mode 100644 Example/Material/SubloadingViscous1D2.supan diff --git a/Example/Material/SubloadingViscous1D.supan b/Example/Material/SubloadingViscous1D.supan new file mode 100644 index 000000000..deb822587 --- /dev/null +++ b/Example/Material/SubloadingViscous1D.supan @@ -0,0 +1,43 @@ +# A TEST MODEL FOR SUBLOADINGVISCOUS1D MATERIAL + +node 1 0 0 +node 2 4 0 +node 3 0 -3 + +material SubloadingViscous1D 1 2E5 \ +200 0 0 0 \ +100 10 0 0 \ +4E1 2 3E2 4 \ +2E1 1E1 0 + +element T2D2 1 1 2 1 10 +element T2D2 2 3 2 1 10 + +hdf5recorder 1 Element PEEQ 1 +hdf5recorder 2 Element E 1 + +fix 1 P 1 3 + +displacement 1 0 0.2 2 2 + +step static 1 +set fixed_step_size 1 +set ini_step_size 1E-1 +set symm_mat 0 + +converger RelIncreDisp 1 1E-10 10 1 + +analyze + +# Node 2: +# Coordinate: +# 4.0000e+00 0.0000e+00 +# Displacement: +# -4.6072e-02 2.0000e-01 +# Resistance: +# 4.0400e-09 1.5923e+03 +peek node 2 + +reset +clear +exit \ No newline at end of file diff --git a/Example/Material/SubloadingViscous1D2.supan b/Example/Material/SubloadingViscous1D2.supan new file mode 100644 index 000000000..19b943f8b --- /dev/null +++ b/Example/Material/SubloadingViscous1D2.supan @@ -0,0 +1,9 @@ +material SubloadingViscous1D 1 2E5 \ +200 0 0 0 \ +100 10 0 0 \ +4E1 2 300 4 \ +2E1 1E1 0 + +materialTest1D 1 1E-4 100 200 200 + +exit \ No newline at end of file diff --git a/MSVC/suanPan/suanPan/suanPan.vcxproj.user b/MSVC/suanPan/suanPan/suanPan.vcxproj.user index a9b0d9d91..c23ddea3b 100644 --- a/MSVC/suanPan/suanPan/suanPan.vcxproj.user +++ b/MSVC/suanPan/suanPan/suanPan.vcxproj.user @@ -4,7 +4,7 @@ false - -f SubloadingMetal + -f SubloadingViscous1D ..\..\..\.dirty WindowsLocalDebugger diff --git a/Material/Material1D/Material1D b/Material/Material1D/Material1D index 5100fe137..d4e0c4606 100644 --- a/Material/Material1D/Material1D +++ b/Material/Material1D/Material1D @@ -17,10 +17,10 @@ #include "Elastic/Sinh1D.h" #include "Elastic/Tanh1D.h" #include "Hysteresis/AFC.h" +#include "Hysteresis/BWBN.h" #include "Hysteresis/BilinearOO.h" #include "Hysteresis/BilinearPO.h" #include "Hysteresis/BoucWen.h" -#include "Hysteresis/BWBN.h" #include "Hysteresis/ComplexHysteresis.h" #include "Hysteresis/CoulombFriction.h" #include "Hysteresis/Flag.h" @@ -44,6 +44,9 @@ #include "Viscosity/Nonviscous01.h" #include "Viscosity/Viscosity01.h" #include "Viscosity/Viscosity02.h" +#include "Wrapper/Parallel.h" +#include "Wrapper/Sequential.h" +#include "Wrapper/Uniaxial.h" #include "vonMises/ArmstrongFrederick1D.h" #include "vonMises/Bilinear1D.h" #include "vonMises/BilinearMises1D.h" @@ -56,6 +59,3 @@ #include "vonMises/NonlinearMises1D.h" #include "vonMises/Subloading1D.h" #include "vonMises/VAFCRP1D.h" -#include "Wrapper/Parallel.h" -#include "Wrapper/Sequential.h" -#include "Wrapper/Uniaxial.h" diff --git a/Material/Material1D/vonMises/Subloading1D.cpp b/Material/Material1D/vonMises/Subloading1D.cpp index a09a5c74a..a2042bf06 100644 --- a/Material/Material1D/vonMises/Subloading1D.cpp +++ b/Material/Material1D/vonMises/Subloading1D.cpp @@ -17,6 +17,9 @@ #include "Subloading1D.h" +#include +#include + const double DataSubloading1D::Saturation::root_one_half = sqrt(1.5); const double Subloading1D::rate_bound = -log(z_bound); @@ -31,10 +34,12 @@ Subloading1D::Subloading1D(const unsigned T, DataSubloading1D&& D, const double : DataSubloading1D{std::move(D)} , Material1D(T, R) {} -int Subloading1D::initialize(const shared_ptr&) { +int Subloading1D::initialize(const shared_ptr& D) { + if(nullptr != D) incre_time = &D->get_factory()->modify_incre_time(); + trial_stiffness = current_stiffness = initial_stiffness = elastic; - initialize_history(3u + static_cast(b.size() + c.size())); + initialize_history(4u + static_cast(b.size() + c.size())); return SUANPAN_SUCCESS; } @@ -51,21 +56,25 @@ int Subloading1D::update_trial_status(const vec& t_strain) { trial_history = current_history; const auto& current_q = current_history(1); const auto& current_z = current_history(2); + const auto& current_zv = current_history(3); auto& iteration = trial_history(0); auto& q = trial_history(1); auto& z = trial_history(2); + auto& zv = trial_history(3); + + const vec current_alpha(¤t_history(4), b.size(), false, true); + const vec current_d(¤t_history(4 + b.size()), c.size(), false, true); + vec alpha(&trial_history(4), b.size(), false, true); + vec d(&trial_history(4 + b.size()), c.size(), false, true); - const vec current_alpha(¤t_history(3), b.size(), false, true); - const vec current_d(¤t_history(3 + b.size()), c.size(), false, true); - vec alpha(&trial_history(3), b.size(), false, true); - vec d(&trial_history(3 + b.size()), c.size(), false, true); + const auto norm_mu = mu / (incre_time && *incre_time > 0. ? *incre_time : 1.); iteration = 0.; auto gamma = 0., ref_error = 0.; auto start_z = current_z; - vec2 residual, incre; - mat22 jacobian; + vec3 residual, incre; + mat33 jacobian; auto counter = 0u; while(true) { @@ -90,7 +99,7 @@ int Subloading1D::update_trial_status(const vec& t_strain) { for(auto I = 0llu; I < b.size(); ++I) bottom_alpha(I) = 1. + b[I].r() * gamma; for(auto I = 0llu; I < c.size(); ++I) bottom_d(I) = 1. + c[I].r() * gamma; - const auto n = trial_stress(0) - a * sum(current_alpha / bottom_alpha) + (z - 1.) * y * sum(current_d / bottom_d) > 0. ? 1. : -1.; + const auto n = trial_stress(0) - a * sum(current_alpha / bottom_alpha) + (zv - 1.) * y * sum(current_d / bottom_d) > 0. ? 1. : -1.; for(auto I = 0llu; I < b.size(); ++I) alpha(I) = (b[I].rb() * gamma * n + current_alpha(I)) / bottom_alpha(I); for(auto I = 0llu; I < c.size(); ++I) d(I) = (c[I].rb() * gamma * n + current_d(I)) / bottom_d(I); @@ -101,7 +110,8 @@ int Subloading1D::update_trial_status(const vec& t_strain) { const auto s = (y * sum_d + a * sum_alpha - current_stress(0)) / (trial_stress(0) - current_stress(0)); if(s >= 1.) { // elastic unloading - z = ((trial_stress(0) - a * sum_alpha) / y - sum_d) / (n - sum_d); + zv = ((trial_stress(0) - a * sum_alpha) / y - sum_d) / (n - sum_d); + z = zv / current_zv * current_z; return SUANPAN_SUCCESS; } if(s > 0.) start_z = 0.; @@ -113,15 +123,24 @@ int Subloading1D::update_trial_status(const vec& t_strain) { const auto trial_ratio = yield_ratio(z); const auto avg_rate = u * trial_ratio(0); + const auto fraction_term = (cv * z - zv) * norm_mu * gamma + 1.; + const auto power_term = pow(fraction_term, nv - 1.); - residual(0) = fabs(trial_stress(0) - elastic * gamma * n - a * sum_alpha + (z - 1.) * y * sum_d) - z * y; + residual(0) = fabs(trial_stress(0) - elastic * gamma * n - a * sum_alpha + (zv - 1.) * y * sum_d) - zv * y; residual(1) = z - start_z - gamma * avg_rate; + residual(2) = zv - fraction_term * power_term * z; - jacobian(0, 0) = n * ((z - 1.) * (y * dd + sum_d * dy) - (a * dalpha + sum_alpha * da)) - elastic - z * dy; + jacobian(0, 0) = n * ((zv - 1.) * (y * dd + sum_d * dy) - (a * dalpha + sum_alpha * da)) - elastic - zv * dy; jacobian(0, 1) = n * y * sum_d - y; + jacobian(0, 2) = 0.; jacobian(1, 0) = -avg_rate; - jacobian(1, 1) = 1. - u * gamma * trial_ratio(1); + jacobian(1, 1) = 0.; + jacobian(1, 2) = 1. - u * gamma * trial_ratio(1); + + jacobian(2, 0) = -z * nv * power_term * (cv * z - zv) * norm_mu; + jacobian(2, 1) = 1. + z * nv * power_term * norm_mu * gamma; + jacobian(2, 2) = -power_term * (fraction_term + z * nv * cv * norm_mu * gamma); if(!solve(incre, jacobian, residual)) return SUANPAN_FAIL; @@ -135,14 +154,19 @@ int Subloading1D::update_trial_status(const vec& t_strain) { } iteration = counter; trial_stress -= elastic * gamma * n; - trial_stiffness += elastic / det(jacobian) * elastic * jacobian(1, 1); + trial_stiffness += elastic / det(jacobian) * elastic * det(jacobian.submat(1, 1, 2, 2)); return SUANPAN_SUCCESS; } gamma -= incre(0); - z -= incre(1); - if(z > 1.) z = 1. - datum::eps; - else if(z < 0.) z = 0.; + zv -= incre(1); + z -= incre(2); + if(z < 0.) z = 0.; + else if(z > 1.) z = 1. - datum::eps; + if(is_viscous) { + if(zv < z) zv = z; + else if(zv > cv * z) zv = cv * z - datum::eps; + } } } @@ -171,6 +195,6 @@ int Subloading1D::reset_status() { } void Subloading1D::print() { - suanpan_info("A uniaxial combined hardening material using subloading surface model.\n"); + suanpan_info("A uniaxial combined hardening material using subloading surface model with optional viscosity.\n"); Material1D::print(); } diff --git a/Material/Material1D/vonMises/Subloading1D.h b/Material/Material1D/vonMises/Subloading1D.h index 6495992fd..c7cc87165 100644 --- a/Material/Material1D/vonMises/Subloading1D.h +++ b/Material/Material1D/vonMises/Subloading1D.h @@ -49,18 +49,14 @@ struct DataSubloading1D { }; const double elastic; // elastic modulus - 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 initial_iso, k_iso, saturation_iso, m_iso; + const double initial_kin, k_kin, saturation_kin, m_kin; const double u; + const double cv; + const double mu; + const double nv; - const std::vector b; - const std::vector c; + const std::vector b, c; }; class Subloading1D final : protected DataSubloading1D, public Material1D { @@ -70,6 +66,10 @@ class Subloading1D final : protected DataSubloading1D, public Material1D { static vec2 yield_ratio(double); + const double* incre_time = nullptr; + + const bool is_viscous = mu > 0. && nv > 0.; + public: Subloading1D( unsigned, // tag diff --git a/Material/MaterialParser.cpp b/Material/MaterialParser.cpp index 66bba6af9..b36b2a6d0 100644 --- a/Material/MaterialParser.cpp +++ b/Material/MaterialParser.cpp @@ -2948,11 +2948,47 @@ void new_subloading1d(unique_ptr& return_obj, istringstream& command) return; } - DataSubloading1D para{p(0), p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8), p(9), {{p(10), 1.}}, {{p(11), p(12)}}}; + DataSubloading1D para{p(0), p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8), p(9), 1., 0., 0., {{p(10), 1.}}, {{p(11), p(12)}}}; if(para.m_iso < 0. || para.m_kin < 0.) { suanpan_error("The evolution rate must be positive.\n"); return; } + if(para.cv < 1.) { + suanpan_error("The viscous limit c_v must be greater than unity.\n"); + return; + } + + return_obj = make_unique(tag, std::move(para), density); +} + +void new_subloadingviscous1d(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, 2E2, 0., 2E2, 1E1, 2E2, 0., 2E2, 1E1, 1E1, 1., 0., 0., 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; + } + + DataSubloading1D para{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), 1.}}, {{p(14), p(15)}}}; + if(para.m_iso < 0. || para.m_kin < 0.) { + suanpan_error("The evolution rate must be positive.\n"); + return; + } + if(para.cv < 1.) { + suanpan_error("The viscous limit c_v must be greater than unity.\n"); + return; + } return_obj = make_unique(tag, std::move(para), density); } @@ -2995,11 +3031,15 @@ void new_multisubloading1d(unique_ptr& return_obj, istringstream& comm } } - DataSubloading1D para{p(0), p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8), p(9), std::move(back), std::move(core)}; + DataSubloading1D para{p(0), p(1), p(2), p(3), p(4), p(5), p(6), p(7), p(8), p(9), 1., 0., 0., std::move(back), std::move(core)}; if(para.m_iso < 0. || para.m_kin < 0.) { suanpan_error("The evolution rate must be positive.\n"); return; } + if(para.cv < 1.) { + suanpan_error("The viscous limit c_v must be greater than unity.\n"); + return; + } return_obj = make_unique(tag, std::move(para), p(10)); } @@ -3604,6 +3644,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, "SubloadingViscous1D")) new_subloadingviscous1d(new_material, command); else if(is_equal(material_id, "MultiSubloading1D")) new_multisubloading1d(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);