双线性弹塑性模型(五)

本节内容为多杆结构的弹塑性有限元计算。

对于弹塑性材料， ,其中

含多个杆单元的结构，需要分别判断每个单元的弹塑性状态，确定是 或者 参与计算。

[算例]

如图所示，两个并联的杆，一段固定，另一端另一端施加轴向力P。荷载逐渐增加到 ，然后逐渐卸载至0。杆一。

杆二。

import matplotlib.pyplot as plt
import numpy as np

import numpy as np
import harden

E1 = 10000; H1 = 1111.11; sYield1 = 5 
E2 = 5000; H2 = 555.55;   sYield2 = 7.5

Et1 = 1000; Et2 = 500

mp1 = [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 = 0

RES = 0
nincr = 0
max_iter = 20

Load =[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 = 0

X = 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 math

def Sgn(x):
    if x > 0 :
        return 1
    elif x < 0 :
        return -1
    else :
        return 0

def 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

迭代路径