Improve stability using second order integration

TLCFEM · Oct 15, 2023 · 2a2b850 · 2a2b850
1 parent c62da07
commit 2a2b850
Show file tree

Hide file tree

Showing 4 changed files with 43 additions and 27 deletions.
diff --git a/Example/Material/CDPM2.supan b/Example/Material/CDPM2.supan
@@ -48,31 +48,31 @@ peek element 1
 # Displacement:
 #   0.0000e+00  0.0000e+00 -1.0000e-01
 # Resistance:
-#  -9.1743e-15  9.7543e-15 -2.0043e+01
+#  -6.8989e-13  6.8141e-13 -2.0141e+01
 # 
 # Node 10:
 # Coordinate:
 #   5.0000e+00  5.0000e+00  2.0000e+01
 # Displacement:
-#   0.0000e+00  3.6694e-02 -1.0000e-01
+#   0.0000e+00  3.6699e-02 -1.0000e-01
 # Resistance:
-#  -9.5924e-15  6.3983e-16 -2.0043e+01
+#  -6.1428e-13 -1.2614e-13 -2.0141e+01
 # 
 # Node 11:
 # Coordinate:
 #  -5.0000e+00  5.0000e+00  2.0000e+01
 # Displacement:
-#  -3.6694e-02  3.6694e-02 -1.0000e-01
+#  -3.6699e-02  3.6699e-02 -1.0000e-01
 # Resistance:
-#   1.5466e-15  3.8572e-15 -2.0043e+01
+#   1.2746e-13 -1.2966e-13 -2.0141e+01
 # 
 # Node 12:
 # Coordinate:
 #  -5.0000e+00 -5.0000e+00  2.0000e+01
 # Displacement:
-#  -3.6694e-02  0.0000e+00 -1.0000e-01
+#  -3.6699e-02  0.0000e+00 -1.0000e-01
 # Resistance:
-#  -1.7442e-15  5.5802e-15 -2.0043e+01
+#   1.2714e-13  6.0664e-13 -2.0141e+01
 peek node 9 10 11 12
 
 # save recorder 1

diff --git a/Example/Material/CDPM2PS.supan b/Example/Material/CDPM2PS.supan
@@ -32,17 +32,17 @@ analyze
 # Coordinate:
 #   1.0000e+00  1.0000e+00
 # Displacement:
-#  -8.0000e-04  9.1915e-04
+#  -8.0000e-04  8.8256e-04
 # Resistance:
-#  -2.4764e-01 -2.7409e-16
+#  -2.7424e-01  1.1102e-16
 # 
 # Node 4:
 # Coordinate:
 #   0.0000e+00  1.0000e+00
 # Displacement:
-#  -8.0000e-04 -2.0510e-04
+#  -8.0000e-04 -2.1286e-04
 # Resistance:
-#  -5.2001e-01 -2.6021e-16
+#  -6.4287e-01 -1.3878e-17
 peek node 3 4
 
 reset

diff --git a/Material/Material3D/Concrete/CDPM2.cpp b/Material/Material3D/Concrete/CDPM2.cpp
@@ -356,6 +356,7 @@ int CDPM2::update_trial_status(const vec& t_strain) {
     const auto trial_p = hydro_stress;
     const vec n = dev_stress / trial_s;
 
+    auto current_gs = 0., current_gp = 0., current_gg = 0.;
     auto gamma = 0., s = trial_s, p = trial_p;
 
     mat jacobian(4, 4, fill::none), left(4, 6, fill::zeros);
@@ -390,31 +391,44 @@ int CDPM2::update_trial_status(const vec& t_strain) {
 
         compute_plasticity(s, p, kp, data);
 
-        if(1u == counter && f < 0.) break;
+        if(1u == counter) {
+            if(f < 0.) break;
+
+            const vec current_effective_stress = initial_stiffness * (current_strain - plastic_strain);
+            const auto current_s = tensor::stress::norm(tensor::dev(current_effective_stress));
+            const auto current_p = tensor::mean3(current_effective_stress);
+
+            podarray<double> current_data(15);
+            compute_plasticity(current_s, current_p, current_kp, current_data);
+
+            current_gs = current_data(4);
+            current_gp = current_data(5);
+            current_gg = current_data(6);
+        }
 
         residual(0) = f;
-        residual(1) = s + double_shear * gamma * gs - trial_s;
-        residual(2) = p + bulk * gamma * gp - trial_p;
-        residual(3) = xh * (current_kp - kp) + gamma * gg;
+        residual(1) = s + shear * gamma * (gs + current_gs) - trial_s;
+        residual(2) = p + half_bulk * gamma * (gp + current_gp) - trial_p;
+        residual(3) = xh * (current_kp - kp) + .5 * gamma * (gg + current_gg);
 
         jacobian(0, 1) = pfps;
         jacobian(0, 2) = pfpp;
         jacobian(0, 3) = pfpkp;
 
-        jacobian(1, 0) = double_shear * gs;
-        jacobian(1, 1) = double_shear * gamma * pgsps + 1.;
-        jacobian(1, 2) = double_shear * gamma * pgspp;
-        jacobian(1, 3) = double_shear * gamma * pgspkp;
+        jacobian(1, 0) = shear * (gs + current_gs);
+        jacobian(1, 1) = shear * gamma * pgsps + 1.;
+        jacobian(1, 2) = shear * gamma * pgspp;
+        jacobian(1, 3) = shear * gamma * pgspkp;
 
-        jacobian(2, 0) = bulk * gp;
-        jacobian(2, 1) = bulk * gamma * pgpps;
-        jacobian(2, 2) = bulk * gamma * pgppp + 1.;
-        jacobian(2, 3) = bulk * gamma * pgppkp;
+        jacobian(2, 0) = half_bulk * (gp + current_gp);
+        jacobian(2, 1) = half_bulk * gamma * pgpps;
+        jacobian(2, 2) = half_bulk * gamma * pgppp + 1.;
+        jacobian(2, 3) = half_bulk * gamma * pgppkp;
 
-        jacobian(3, 0) = gg;
-        jacobian(3, 1) = gamma / gg * (gs * pgsps + gp / 3. * pgpps);
-        jacobian(3, 2) = gamma / gg * (gs * pgspp + gp / 3. * pgppp) + (current_kp - kp) * dxhdp;
-        jacobian(3, 3) = gamma / gg * (gs * pgspkp + gp / 3. * pgppkp) - xh;
+        jacobian(3, 0) = .5 * (gg + current_gg);
+        jacobian(3, 1) = .5 * gamma / gg * (gs * pgsps + gp / 3. * pgpps);
+        jacobian(3, 2) = .5 * gamma / gg * (gs * pgspp + gp / 3. * pgppp) + (current_kp - kp) * dxhdp;
+        jacobian(3, 3) = .5 * gamma / gg * (gs * pgspkp + gp / 3. * pgppkp) - xh;
 
         if(!solve(incre, jacobian, residual, solve_opts::equilibrate + solve_opts::refine)) return SUANPAN_FAIL;
 

diff --git a/Material/Material3D/Concrete/CDPM2.h b/Material/Material3D/Concrete/CDPM2.h
@@ -74,7 +74,9 @@ class CDPM2 final : protected DataCDPM2, public Material3D {
     static const mat unit_dev_tensor;
 
     const double double_shear = elastic_modulus / (1. + poissons_ratio);
+    const double shear = .5 * double_shear;
     const double bulk = elastic_modulus / (3. - 6. * poissons_ratio);
+    const double half_bulk = .5 * bulk;
 
     const DamageType damage_type = DamageType::ANISOTROPIC;