双线性弹塑性模型(三)

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  #alpha不变
        ep = epsN       #塑性应变不增加，为0
        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

E = 200E3; H = 25E3; sYield = 250
Et = E*H / (E+H)
mp = [E, H, sYield]
#上一步的应力，塑性应变以及alpha值
nS = 200
nep = 1E-4
nA = 2.5

delta_eps = -0.003

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

print(Snew, Anew, epnew, flag)