Add prototype

TLCFEM · Nov 2, 2024 · 3cd4cfd · 3cd4cfd
1 parent 36fc9cf
commit 3cd4cfd
Show file tree

Hide file tree

Showing 6 changed files with 153 additions and 25 deletions.
diff --git a/Example/Material/SubloadingViscous1D.supan 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
diff --git a/Example/Material/SubloadingViscous1D2.supan 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
diff --git a/MSVC/suanPan/suanPan/suanPan.vcxproj.user b/MSVC/suanPan/suanPan/suanPan.vcxproj.user
@@ -4,7 +4,7 @@
     <ShowAllFiles>false</ShowAllFiles>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
-    <LocalDebuggerCommandArguments>-f SubloadingMetal</LocalDebuggerCommandArguments>
+    <LocalDebuggerCommandArguments>-f SubloadingViscous1D</LocalDebuggerCommandArguments>
     <LocalDebuggerWorkingDirectory>..\..\..\.dirty</LocalDebuggerWorkingDirectory>
     <DebuggerFlavor>WindowsLocalDebugger</DebuggerFlavor>
   </PropertyGroup>

diff --git a/Material/Material1D/vonMises/Subloading1D.cpp b/Material/Material1D/vonMises/Subloading1D.cpp
@@ -17,6 +17,9 @@
 
 #include "Subloading1D.h"
 
+#include <Domain/DomainBase.h>
+#include <Domain/Factory.hpp>
+
 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<DomainBase>&) {
+int Subloading1D::initialize(const shared_ptr<DomainBase>& D) {
+    if(nullptr != D) incre_time = &D->get_factory()->modify_incre_time();
+
     trial_stiffness = current_stiffness = initial_stiffness = elastic;
 
-    initialize_history(3u + static_cast<unsigned>(b.size() + c.size()));
+    initialize_history(4u + static_cast<unsigned>(b.size() + c.size()));
 
     return SUANPAN_SUCCESS;
 }
@@ -44,7 +49,9 @@ unique_ptr<Material> Subloading1D::get_copy() { return make_unique<Subloading1D>
 int Subloading1D::update_trial_status(const vec& t_strain) {
     incre_strain = (trial_strain = t_strain) - current_strain;
 
-    if(fabs(incre_strain(0)) <= datum::eps) return SUANPAN_SUCCESS;
+    const auto zero_increment = fabs(incre_strain(0)) <= datum::eps;
+
+    if(!is_viscous && zero_increment) return SUANPAN_SUCCESS;
 
     trial_stress = current_stress + (trial_stiffness = initial_stiffness) * incre_strain;
 
@@ -54,18 +61,21 @@ int Subloading1D::update_trial_status(const vec& t_strain) {
     auto& iteration = trial_history(0);
     auto& q = trial_history(1);
     auto& z = trial_history(2);
+    auto& zv = trial_history(3);
 
-    const vec current_alpha(&current_history(3), b.size(), false, true);
-    const vec current_d(&current_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 vec current_alpha(&current_history(4), b.size(), false, true);
+    const vec current_d(&current_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 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) {
@@ -97,12 +107,16 @@ int Subloading1D::update_trial_status(const vec& t_strain) {
 
         const auto sum_alpha = sum(alpha), sum_d = sum(d);
 
-        if(1u == counter) {
+        if(1u == counter && !zero_increment) {
             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);
-                return SUANPAN_SUCCESS;
+                zv = ((trial_stress(0) - a * sum_alpha) / y - sum_d) / (n - sum_d);
+                if(zv < z) {
+                    // elastic unloading
+                    z = zv;
+                    return SUANPAN_SUCCESS;
+                }
+                // return SUANPAN_SUCCESS;
             }
             if(s > 0.) start_z = 0.;
         }
@@ -113,15 +127,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 + (z - 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, 1) = n * y * sum_d - y;
+        jacobian(0, 0) = n * ((z - 1.) * (y * dd + sum_d * dy) - (a * dalpha + sum_alpha * da)) - elastic - zv * dy;
+        jacobian(0, 1) = -y;
+        jacobian(0, 2) = n * y * sum_d;
 
         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 +158,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 +199,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
@@ -52,6 +52,9 @@ struct DataSubloading1D {
     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<Saturation> b, c;
 };
@@ -63,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
@@ -2948,11 +2948,47 @@ void new_subloading1d(unique_ptr<Material>& 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<Subloading1D>(tag, std::move(para), density);
+}
+
+void new_subloadingviscous1d(unique_ptr<Material>& 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<Subloading1D>(tag, std::move(para), density);
 }
@@ -2995,11 +3031,15 @@ void new_multisubloading1d(unique_ptr<Material>& 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<Subloading1D>(tag, std::move(para), p(10));
 }
@@ -3604,6 +3644,7 @@ int create_new_material(const shared_ptr<DomainBase>& 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);