diff --git a/Example/Material/CDPM2.supan b/Example/Material/CDPM2.supan index dad377a8a..47142c666 100644 --- 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 index 1a73218a3..43c557754 100644 --- 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 index a611f69aa..8630c64d6 100644 --- 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 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 index fd29d026e..172e6b7a0 100644 --- 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;