文章目录
- 0. 背景
- 1. 代码
- 2. 结果测试
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时结果与已有文件一致(注意!受拉时应变应该为负值!!!)