ABAQUS混凝土塑性损伤模型参数生成

def cal_con_c(s_cr, ac, n, fck, e0, ec, epc, scp): # 混凝土受压时应力应变曲线定义/求解各非弹性应变对应的应力、应变、损伤  sin_cr = s_cr - fck / ec # 计算抗压强度峰值位置的非弹性应变，用于后续迭代求解各非弹性应变值位置时，判别迭代区间。  pc = fck / (e0 * s_cr)
   s0 = Cal_con_c0(s_cr, ac, n, fck, e0, epc).s0 if scp < 10**(-8):
      dc = 0.0  return fck*epc, s0, dc elif scp <= sin_cr: # 求取上升段的应力、应变及损伤参数。  res = 1  m = 0  s_min = s0
      s_max = s_cr
      s = s_cr / 2  while (abs(res) > 10 ** (-8)) and m < 100:
         res = s - pc * n * e0 * s / (n - 1  (s / s_cr) ** n)/ec - scp if res > 0:
            s_max = s
            s = (s   s_min) / 2  else:
            s_min = s
            s = (s   s_max) / 2  m  = 1  stress_c = pc * n * e0 * s / (n - 1  (s / s_cr) ** n)
      dc = 1 - (stress_c/s/ec)**0.5  return stress_c, s, dc else: # 求取下降段的应力、应变及损伤参数。  res = 1  m = 0  s_min = s_cr
      s_max = 5  s = s_cr*2  while (abs(res) > 10 ** (-8)) and m < 100:
         res = s - pc * e0 * s / (ac*(s/s_cr - 1)**2  s/s_cr)/ec - scp if res > 0:
            s_max = s
            s = (s   s_min) / 2  else:
            s_min = s
            s = (s   s_max) / 2  m  = 1  stress_c = pc * e0 * s / (ac*(s/s_cr - 1)**2  s/s_cr)
      dc = 1 - (stress_c/s/ec)**0.5  return stress_c, s, dc def cal_con_t(ftk, s_tr, at, e0, ec, stp): # 混凝土受拉时应力应变曲线定义/求解各非弹性应变对应的应力、应变、损伤  if stp < 10**(-10):
      dt = 0.0  return ftk, s_tr, dt else:
      res = 1  m = 0  s_min = s_tr
      s_max = 10  s = s_tr * 2  pt = ftk/(e0*s_tr) while (abs(res) > 10 ** (-8)) and m < 100:
         res = s - pt * e0 * s / (at * (s / s_tr - 1) ** 1.7  s / s_tr) / ec - stp if res > 0:
            s_max = s
            s = (s   s_min) / 2  else:
            s_min = s
            s = (s   s_max) / 2  m  = 1  stress_t = pt * e0 * s / (at * (s / s_tr - 1) ** 1.7  s / s_tr)
      dt = 1 - (stress_t / s / ec)**0.5  return stress_t, s, dt