本节内容为多杆结构的弹塑性有限元计算。对于弹塑性材料， ,其中 含多个杆单元的结构，需要分别判断每个单元的弹塑性状态，确定是 或者 参与计算。[算例]如图所示，两个并联的杆，一段固定，另一端另一端施加轴向力P。荷载逐渐增加到 ，然后逐渐卸载至0。杆一。杆二。import matplotlib.pyplot as pltimport numpy as npimport numpy as npimport hardenE1 = 10000; H1 = 1111.11; sYield1 = 5 E2 = 5000; H2 = 555.55; sYield2 = 7.5Et1 = 1000; Et2 = 500mp1 = [E1, H1, sYield1]; mp2 = [E2, H2, sYield2]; nStress1 = 0; nAlpha1=0; neps1 = 0; eps_new1=0; nStress2 = 0; nAlpha2=0; neps2 = 0; eps_new2=0; A1 = 0.75; L1 = 100; A2 = 1.25; L2 = 100; tol = 1.0E-5; u = 0RES = 0nincr = 0max_iter = 20Load =[0,1,2,3,4,5,6,7,8,9,10,11,12,13, 14,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0]N = len(Load)flag1 = 0; flag2 = 0X = np.zeros( (N) )Y1 = np.zeros( (N) )Y2 = np.zeros( (N) )for i in range(N): nincr += 1 print('第{}增量步：' .format(nincr) ) P = Load[i] RES = P - nStress1 * A1 - nStress2 * A2 du = 0 niter = 0 conv = 2e12 print('迭代步 位移 不平衡力 收敛参数') while ( conv > tol and niter < max_iter ): niter += 1 Eep1 = E1; Eep2 = E2 if flag1 == 1: Eep1 = Et1 if flag2 == 1: Eep2 = Et2 K = Eep1 * A1 / L1 + Eep2 * A2 / L2 delta_u = RES / K du = du + delta_u delta_eps = delta_u / L1 Stress_new1, Alpha_new1, eps_new1, flag1 = harden.KinematicHard1D(mp1, delta_eps,nStress1,nAlpha1,neps1) Stress_new2, Alpha_new2, eps_new2, flag2 = harden.KinematicHard1D(mp2, delta_eps,nStress2,nAlpha2,neps2) RES = P - Stress_new1*A1 - Stress_new2*A2 conv = RES**2 / (1 + P**2) nStress1 = Stress_new1; nStress2 = Stress_new2 neps1 = eps_new1; neps2 = eps_new2 nAlpha1 = Alpha_new1; nAlpha2 = Alpha_new2 print(format(niter, '>3d'), format(u, '>20.12f'), format(nStress1, '>22.8f'), format(nStress2, '>22.8f'), format(conv, '>22.14e')) u = u + du X[i] = u; Y1[i] = nStress1; Y2[i] = nStress2 print(X)print(Y1)print(Y2)fig, axs = plt.subplots(1, 1, figsize=(7,12) )axs.plot(X, Y1, label="M1", linewidth = 2, color = "deeppink",marker='o', markersize=6) axs.plot(X, Y2, label="M1", linewidth = 2, color = "green",marker='o', markersize=6) axs.set_xlabel('$ Displacement(mm) $', fontsize = 18)axs.set_ylabel('$ Stress(MPa) $', fontsize = 18)plt.ylim([-5,9])plt.legend(['Bar1','Bar2']) fig.savefig('./f428.png', dpi = 300) #保存图片 plt.show()##### harden模块import mathdef Sgn(x): if x > 0 : return 1 elif x < 0 : return -1 else : return 0def KinematicHard1D(MP,deps, stressN,alphaN, epsN): E = MP[0] H = MP[1] Y0 = MP[2] stress_tr = stressN + E*deps eta_tr = stress_tr - alphaN f_tr = math.fabs(eta_tr ) - Y0 if f_tr < 0: stress = stress_tr alpha = alphaN ep = epsN flag = 0 # 处于弹性状态 else: dep = f_tr / ( E + H ) stress = stress_tr - Sgn(eta_tr) * E * dep alpha = alphaN + Sgn(eta_tr) * H * dep ep = epsN + dep flag = 1 # 处于塑性状态 return stress, alpha, ep, flag 迭代路径来源：数值分析与有限元编程
