Subloading surface model with viscosity for metals (#223)

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

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

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>

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"

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;
 }
@@ -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(&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 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 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();
 }

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<Saturation> b;
-    const std::vector<Saturation> c;
+    const std::vector<Saturation> 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

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);