0. 背景

最近在梳理混凝土塔筒的计算指南，在求解弯矩曲率关系以及MN相关曲线时，需要混凝土的本构关系作为输入条件。

1. 代码

这段代码还是比较简单的。不过需要注意的是，我把受拉和受压两种状态统一了起来，即规定受压为正，受拉为负。

def constitutive_GB(epsilon, Ec, f_cr, f_tr):'''按照 GB50010-2010混凝土结构设计规范 附录C的本构关系编制，由应变给出应力:param epsilon: 应变，压为正，拉为负:param Ec: 混凝土弹性模量，Mpa:param f_cr: 混凝土单轴抗压强度代表值， Mpa:param f_tr: 混凝土单轴抗拉强度代表值， Mpa:return: 应力， Mpa'''alpha_t = 0.312*f_tr*f_tr   #epsilon_tr = f_tr**0.54*65/10**6rho_t = f_tr / Ec / epsilon_tralpha_c = 0.157 * f_cr**0.785 - 0.905epsilon_cr = (700 + 172 * m.sqrt(f_cr))/10**6# epsilon_cu是极限压应变，不知道后面会不会用上，先保留# epsilon_cu = epsilon_cr / 2 / alpha_c * (1 + 2*alpha_c + m.sqrt(1 + 4*alpha_c))n = Ec * epsilon_cr / (Ec * epsilon_cr - f_cr)rho_c = f_cr / Ec / epsilon_crx = epsilon/epsilon_cr if epsilon >= 0 else epsilon/epsilon_trif x > 1:    # 受压状态d = 1 - rho_c / (alpha_c*(x-1)**2 + x)elif 0 <= x < 1:    # 受压状态d = 1 - rho_c*n / (n - 1 + x**n)elif -1 <= x < 0:   # 受拉状态d = 1 - rho_t * (1.2 + 0.2*x**5)     # 受拉时，应变为负值，x也为负值else:d = 1 - rho_t / (alpha_t * (-x-1)**1.7 -x)sigma = (1 - d)*Ec*epsilonreturn sigma

2. 结果测试

如下，以C70混凝土的相关参数，取应变值为0.0001670进行计算，结果与本文所附的Excel结果一致。说明这段代码应该是没有问题的。

# 受压时
epsilon = 0.0001670
Ec = 37000
f_cr = 31.8
f_tr = 2.14
sigma = constitutive_GB(epsilon, Ec, f_cr, f_tr)
# (sigma, d, n, rho_c)
# (6.12871040573435, 0.008138791756862052, 2.060440315282528, 0.5146668444686883)
# C70, 应变0.000167时结果与已有文件一致，

# 受拉时
epsilon = -9.8025*10**(-5) * 0.5
Ec = 37000
f_cr = 31.8
f_tr = 2.14sigma = constitutive_GB(epsilon, Ec, f_cr, f_tr)# (sigma, d, rho_t) 
# (-1.2773126624557158, 0.295649806678817, 0.5900315775743109)
# C70, 应变-9.8025*10**(-5) * 0.5时结果与已有文件一致(注意！受拉时应变应该为负值！！！)