双线性弹塑性模型(四)

本节内容为在牛顿-拉夫逊方法中集成基于随动硬化模型的当前应力计算。

对于非线性的问题，一般将其线性化为

一次迭代得到的是位移增量，如图所示

接下来要将位移增量转化为应变增量，以一维杆结构为例，其应变增量

其中    为杆初始长度。

[算例]

一根各向同性杆，一端固定，另一端施加轴向力做拉伸试验.荷载逐渐增加到    ，然后逐渐卸载至0。，    。

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

#########  以上为harden模块


import matplotlib.pyplot as plt
import numpy as np

import numpy as np
import harden

E = 70000; H = 10000; sYield = 250
Et = E*H / (E+H)
mp = [E, H, sYield]

nS = 0
nep = 0
nA = 0

A = 100
L = 1000

tol = 1.0E-1
u = 0
Res = 0
nincr = 0
max_iter = 20
epnew = 0

Load = [0,5000, 10000, 15000,20000,25000,30000,25000,20000,
        15000,10000,5000,0 ]

N = len(Load)
flag = 0

X = np.zeros( (N) )
Y = np.zeros( (N) )

for i in range(N):
    nincr += 1
    print('第{}增量步：' .format(nincr) )

    P = Load[i]

    Res = P - nS * A
    du = 0
    
    niter = 0
    conv = 2e12
    print('迭代步            位移                   不平衡力                    收敛参数')
    
    while ( conv > tol and niter < max_iter ):
        niter += 1

        Eep = E
        if flag == 1:
            Eep = Et

        K = Eep*A/L
        delta_u = Res / K
        du = du + delta_u
        delta_eps = delta_u / L

        Snew, Anew, epnew, flag = harden.KinematicHard1D(mp,delta_eps,nS,nA,nep)

        Res = P - Snew*A
        conv = Res**2 / (1 + P**2)
        nS = Snew
        nep = epnew
        nA  = Anew
        print(format(niter, '>3d'), format(u, '>20.12f'), format(nS, '>26.14f'), format(conv, '>28.16e'))

    u = u + du
    
    X[i] = u
    Y[i] = nS

print(X)
print(Y)

fig, axs = plt.subplots(1, 1,  figsize=(8,6) )
axs.plot(X, Y, label="M1", linewidth = 3, color = "deeppink",marker='o', markersize=12) 
axs.set_xlabel('$Displacement(mm)$', fontsize = 18)
axs.set_ylabel('$Stress(MPa)$', fontsize = 18)

fig.savefig('./f358.png', dpi = 300) #保存图片 
plt.show()