本文是土木工程设计智能化尝试，基于Python实现D值法自动计算。介绍了D值法，阐述程序流程，包括结构定义、计算αc值、确定反弯点位置等步骤，还绘制了弯矩图。搭建了基于Flask框架和HTML的网页应用，方便用户使用，最后给出演示网页和源码下载地址。
土木工程设计系列-框架结构设计：D值法Python自动计算
框架结构设计-D值法Python自动计算
前言
本文为土木工程设计智能化自动化的尝试，希望存在错误的多多指出。此程序主要讨论以编程的思想实现简单结构下D值法的自动计算，作为互联网时代框架结构设计课程的拓展部分。本文的代码基于python编写，并在最后开发了一个基于网页的应用页面作为demo。
D值法介绍
D值法是由日本抗震结构学家武藤清（1903-1989）在1920年代提出的一种结构分析方法。这种方法专注于通过分布系数简便地计算地震等引起的水平荷载作用下，长方形框架的柱和耐震墙等结构元素上的剪力。它被广泛认为是反弯点法的一个改进，特别是在考虑柱的侧向刚度和反弯点高度的修正上。
武藤清的这项工作在1933年被采纳进日本建筑学会（当时称为建筑学会）的钢筋混凝土结构计算标准中。武藤清不仅在结构工程领域有重要贡献，他的职业生涯还涵盖了广泛的研究领域，包括抗爆研究、强震反应分析、以及高层建筑研究等。他的工作为后来的结构设计和抗震技术发展奠定了重要的基础。D值法的具体理论和计算方法主要基于对结构在受到水平荷载时的响应分析，通过修正反弯点法中的某些关键参数来实现更为精确的结构性能预测。这包括对柱子侧向刚度的特定表示（以D值表示），使得该方法在工程实践中既简便又有效。D值法的计算流程如下所示。
程序流程
程序的整体流程如下图所示：
- 结构定义：这一步骤包括初始化节点（Node）和杆件（Element），以及定义结构的边界条件和荷载情况。
- 结构类初始化：利用节点、杆件、边界条件和荷载信息初始化结构类（Structure），为D值法计算做准备。
- 计算alpha_c值：根据结构的几何形状和材料属性计算每根柱子的alpha_c值，这个值对后续的弯矩和剪力计算很重要。
- 确定反弯点位置和反弯点插值计算：通过插值方法基于结构和荷载条件确定反弯点的位置。
- 计算D值和剪力：利用反弯点位置和alpha_c值计算结构在不同楼层的D值，进而确定剪力分布。
- 弯矩计算：基于D值和剪力，通过结构力学原理计算各节点和杆件的弯矩值。
- 绘制弯矩图：最后，使用matplotlib库根据计算出的弯矩值在结构模型上绘制弯矩分布图。
结构定义
在D值法计算的流程中，首先是结构定义这一部分，这是分析的基础。结构定义主要涉及到初始化节点（Node）、杆件（Element）、边界条件和荷载。这里，我们会根据代码的相关片段，详细介绍每个组件的定义和作用。
节点类（Node）
节点是结构分析中的关键组成部分，表示结构中的关键点位，如梁、柱的连接点。节点类定义了节点的基本属性，包括名称、坐标位置以及与其它节点的连接关系。
class Node:
def __init__(self, name, x, y):
self.x = x
self.y = y
self.name = name
self.link = {'up': None, 'down': None, 'left': None, 'right': None}
# 用于记录各杆的内力弯矩
self.moment = {'up': None, 'down': None, 'left': None, 'right': None}
杆件类（Element）
杆件代表结构中的梁、柱等承载元素。杆件类通过定义两个节点间的连接关系、材料属性（如弹性模量E、截面积A、惯性矩I）以及长度等，来模拟实际结构中的杆件。
class Element:
def __init__(self, p1, p2, E, A, I):
self.p1 = p1 # 点1
self.p2 = p2 # 点2
self.E = E # 弹模
self.A = A # 截面积
self.I = I # 惯性矩
self.l = np.linalg.norm(np.array((p1.x, p1.y)) - np.array((p2.x, p2.y))) # 长度
self.i = self.E*self.I/self.l
self.theta = math.atan2(p2.y - p1.y, p2.x - p1.x)
self.node_link()
结构类（Structure）
结构类将节点和杆件整合在一起，定义了整个待分析的结构。它包含了结构的边界条件和荷载情况，为接下来的分析过程提供了基础数据。
# 定义结构的几何位置
class Structure:
def __init__(self, nodes, elems, bc, load):
# 传入节点字典
self.nodes = nodes
# 传入杆件字典
self.elems = elems
# 边界条件
self.bc = bc
self.fixed()
# 荷载情况
sorted_keys = sorted(load.keys(), reverse=True)
cumulative_sum = 0
transformed_dict = {}
for key in sorted_keys:
cumulative_sum += load[key]
transformed_dict[key] = cumulative_sum
self.load = transformed_dict
# 各竖向轴的杆件
self.axis_ele = defaultdict(list)
self.extract_axis()
# 提取各轴的杆件
def extract_axis(self):
for ele in self.elems:
if self.elems[ele].theta != 0:
self.axis_ele[self.elems[ele].p1.x].append(ele)
# 标记是否为支座
def fixed(self):
for f in self.bc:
self.nodes[f].link['down'] = 'fixed'
这些类和初始化过程为D值法分析提供了必要的数据结构和输入。结构定义是后续分析中不可或缺的一步，确保了分析的准确性和有效性。在完成这一步之后，我们可以继续进行D值法的计算，包括确定反弯点的位置、计算D值、剪力和弯矩等，这将在流程的后续部分进一步介绍。
D值法计算
这一部分包含了核心的结构分析过程。在代码中，这个过程主要由D_value_method类实现，涉及到计算
α
c
\alpha_c
αc值、确定反弯点位置、计算D值、剪力和弯矩等步骤。首先先建立一个类，类里面的变量如下，其中v1.json等文件为规范中对应的表格，可见于附件中：
# D值法的类
class D_value_method:
def __init__(self, structure):
self.structure = structure
self.D_value_info = defaultdict(lambda: defaultdict(dict))
self.n = len(structure.axis_ele[0])
self.D_sum = defaultdict(lambda: 0)
# 记录弯矩值以列表形式记录每根杆的弯矩，从上到下与从左到右的顺序
self.M = defaultdict(list)
# 加载json文件，反弯点表
with open('v0.json', 'r') as file:
self.v0_t = json.load(file)
for n in self.v0_t:
for j in self.v0_t[n]:
x_known = list(map(lambda x:float(x), self.v0_t[n][j].keys()))
y_known = list(map(lambda x:float(x), self.v0_t[n][j].values()))
# 创建插值函数
interp_function = interp1d(x_known, y_known)
self.v0_t[n][j] = interp_function
with open('v1.json', 'r') as file:
data = json.load(file)
# 提取外部键作为一个维度
x = np.array(list(map(float, data.keys())))
# 假设所有外部键具有相同的内部键集合，提取这些内部键作为另一个维度
inner_keys = list(map(float, list(data[next(iter(data))].keys())))
y = np.array(sorted(inner_keys))
# 创建值矩阵
values = np.array([[data[str(x_val)].get(str(y_val), np.nan) for y_val in y] for x_val in x]).T
# 创建 2D 插值函数
self.v1_t = interp2d(x, y, values, kind='linear')
with open('v2.json', 'r') as file:
data = json.load(file)
# 提取外部键作为一个维度
x = np.array(list(map(float, data.keys())))
# 假设所有外部键具有相同的内部键集合，提取这些内部键作为另一个维度
inner_keys = list(map(float, list(data[next(iter(data))].keys())))
y = np.array(sorted(inner_keys))
# 创建值矩阵
values = np.array([[data[str(x_val)].get(str(y_val), np.nan) for y_val in y] for x_val in x]).T
# 创建 2D 插值函数
self.v2_t = interp2d(x, y, values, kind='linear')
with open('v3.json', 'r') as file:
data = json.load(file)
# 提取外部键作为一个维度
x = np.array(list(map(float, data.keys())))
# 假设所有外部键具有相同的内部键集合，提取这些内部键作为另一个维度
inner_keys = list(map(float, list(data[next(iter(data))].keys())))
y = np.array(sorted(inner_keys))
# 创建值矩阵
values = np.array([[data[str(x_val)].get(str(y_val), np.nan) for y_val in y] for x_val in x]).T
# 创建 2D 插值函数
self.v3_t = interp2d(x, y, values, kind='linear')
# 计算alpha_c值
self.alpha_c_cal()
self.v_cal()
self.D_cal()
self.Vim_cal()
self.M_cal()
计算 α c \alpha_c αc值
计算每根柱子的 α c \alpha_c αc值是分析过程的重要一步，它反映了柱子的侧向刚度对整体结构稳定性的影响。 α c \alpha_c αc值的计算需要考虑柱子的线刚度、高度以及与其它杆件的连接情况。
# 计算alphac
def alpha_c_cal(self):
for a in self.structure.axis_ele:
for index, f in enumerate(self.structure.axis_ele[a]):
# 这里的index+1是实际层数
# 该层上部点对应的点位名称
self.D_value_info[a][index + 1]['up_point'] = self.structure.nodes[f[-1]]
self.D_value_info[a][index + 1]['h'] = self.structure.elems[f].l
self.D_value_info[a][index + 1]['i'], \
self.D_value_info[a][index + 1]['alpha_c'], \
self.D_value_info[a][index + 1]['i_1'], \
self.D_value_info[a][index + 1]['i_2'], \
self.D_value_info[a][index + 1]['i_3'], \
self.D_value_info[a][index + 1]['i_4'] = self.get_i_u(self.structure.elems[f].p1,
self.structure.elems[f].p2)
self.D_value_info[a][index + 1]['E'] = self.structure.elems[f].E
self.D_value_info[a][index + 1]['I'] = self.structure.elems[f].I
# 获取i值
# 除底层外其余各层柱的线刚度应乘以 0.9 的修正系数 现浇楼板的中框架的线刚度要乘以2
# 注意此函数只面向一般的情况，若要通用则应该进行更改
def get_i_u(self, p1, p2):
eles = self.structure.elems
if p1.link['down'] != 'fixed':
i_1 = eles[''.join(sorted(p2.name+p2.link['left'].name))].i*2 if p2.link['left']!=None else 0
i_2 = eles[''.join(sorted(p2.name+p2.link['right'].name))].i*2 if p2.link['right']!=None else 0
# 除底层外各层柱要乘以0.9
i_c = eles[p1.name+p2.name].i*0.9
i_3 = eles[''.join(sorted(p1.name+p1.link['left'].name))].i*2 if p1.link['left']!=None else 0
i_4 = eles[''.join(sorted(p1.name+p1.link['right'].name))].i*2 if p1.link['right']!=None else 0
i_u = (i_1+i_2+i_3+i_4)/2/i_c
return i_u, i_u/(2+i_u), i_1, i_2, i_3, i_4
else:
i_1 = eles[''.join(sorted(p2.name+p2.link['left'].name))].i*2 if p2.link['left']!=None else 0
i_2 = eles[''.join(sorted(p2.link['right'].name+p2.name))].i*2 if p2.link['right']!=None else 0
i_c = eles[p1.name+p2.name].i
i_3 = 0
i_4 = 0
i_u = (i_1 + i_2) / i_c
return i_u, (0.5+i_u)/(2+i_u), i_1, i_2, i_3, i_4
确定反弯点的位置
反弯点是弯矩为零的点，确定反弯点的位置对于简化计算和分析结构的内力分布非常关键。这一步骤需要依据结构的几何形状和荷载条件，通过插值方法确定。
# 确定反弯点的位置
def v_cal(self):
for axial_i, axial in self.D_value_info.items():
for f_i, f in axial.items():
self.D_value_info[axial_i][f_i]['v0'] = \
float(self.v0_t[str(self.n)][f"{f_i}.0"](f['i']))
alpha_1 = min((f['i_3']+f['i_4']), (f['i_1']+f['i_2']))\
/max((f['i_3']+f['i_4']), (f['i_1']+f['i_2']))
self.D_value_info[axial_i][f_i]['v1'] = self.v1_t(alpha_1, f['i'])[0]
if f_i == self.n:
self.D_value_info[axial_i][f_i]['v2'] = 0
else:
alpha_2 = axial[f_i+1]['h']/f['h']
self.D_value_info[axial_i][f_i]['v2'] = self.v2_t(alpha_2, f['i'])[0]
if f_i == 1:
self.D_value_info[axial_i][f_i]['v3'] = 0
else:
alpha_3 = axial[f_i - 1]['h'] / f['h']
self.D_value_info[axial_i][f_i]['v3'] = self.v3_t(alpha_3, f['i'])[0]
self.D_value_info[axial_i][f_i]['v'] = self.D_value_info[axial_i][f_i]['v0'] + \
self.D_value_info[axial_i][f_i]['v1'] + \
self.D_value_info[axial_i][f_i]['v2'] + \
self.D_value_info[axial_i][f_i]['v3']
self.D_value_info[axial_i][f_i]['vh'] = self.D_value_info[axial_i][f_i]['v'] * \
self.D_value_info[axial_i][f_i]['h']
计算D值和剪力
D值的计算是D值法的核心，它涉及到利用上一步确定的反弯点位置和 α c \alpha_c αc值，计算结构在不同楼层的D值，进而确定剪力分布。
def D_cal(self):
for axial_i, axial in self.D_value_info.items():
for f_i, f in axial.items():
# 单位变为kn/m
self.D_value_info[axial_i][f_i]['D1'] = 12*f['E']*f['I']/f['h']**3*10e-4
self.D_value_info[axial_i][f_i]['D'] = self.D_value_info[axial_i][f_i]['D1']*f['alpha_c']
# 计算Dsum
self.D_sum[f_i] += self.D_value_info[axial_i][f_i]['D']
计算弯矩
最后一步是计算各节点和杆件的弯矩值。这一步骤是基于前面计算的D值和剪力，通过结构力学原理得到的。
# 计算反弯点处的弯矩
def Vim_cal(self):
for axial_i, axial in self.D_value_info.items():
for f_i, f in axial.items():
# 单位变为kn
self.D_value_info[axial_i][f_i]['Vim'] = f['D']/self.D_sum[f_i]*self.structure.load[f_i]*10e-4
# 计算弯矩
def M_cal(self):
for axial_i, axial in self.D_value_info.items():
for f_i in sorted(axial.keys(), reverse=True):
f = axial[f_i]
self.D_value_info[axial_i][f_i]['M_c_up'] = (f['h'] - f['vh'])*f['Vim']
self.D_value_info[axial_i][f_i]['M_c_down'] = f['vh'] * f['Vim']
# 左右梁弯矩与相对线刚度有关
if f_i == self.n:
# 如果是顶层
self.D_value_info[axial_i][f_i]['M_b_left'] = f['i_1']/(f['i_1']+f['i_2'])*\
self.D_value_info[axial_i][f_i]['M_c_up']
self.D_value_info[axial_i][f_i]['M_b_right'] = f['i_2']/(f['i_1']+f['i_2'])*\
self.D_value_info[axial_i][f_i]['M_c_up']
else:
self.D_value_info[axial_i][f_i]['M_b_left'] = f['i_1'] / (f['i_1'] + f['i_2']) * \
(self.D_value_info[axial_i][f_i]['M_c_up']+
self.D_value_info[axial_i][f_i+1]['M_c_down'])
self.D_value_info[axial_i][f_i]['M_b_right'] = f['i_2'] / (f['i_1'] + f['i_2']) * \
(self.D_value_info[axial_i][f_i]['M_c_up'] +
self.D_value_info[axial_i][f_i + 1]['M_c_down'])
# 配置节点的内力信息
self.D_value_info[axial_i][f_i]['up_point'].moment['down'] \
= self.D_value_info[axial_i][f_i]['M_c_up']
self.D_value_info[axial_i][f_i]['up_point'].link['down'].moment['up'] \
= self.D_value_info[axial_i][f_i]['M_c_down']
self.D_value_info[axial_i][f_i]['up_point'].moment['left'] \
= self.D_value_info[axial_i][f_i]['M_b_left']
self.D_value_info[axial_i][f_i]['up_point'].moment['right'] \
= self.D_value_info[axial_i][f_i]['M_b_right']
绘制弯矩图
这一步骤是D值法分析的可视化部分，通过绘图直观地展示结构在特定荷载作用下的内力分布，尤其是弯矩。这对于理解结构的受力情况和进行设计修改具有重要意义。绘制弯矩图的过程涉及到使用matplotlib库根据计算出的弯矩值在结构模型上绘图。以下是实现这一过程的关键步骤：
# 绘制弯矩图
def draw_moment_diagram(nodes, amp=20):
fig, ax = plt.subplots(figsize=(12, 15), dpi=200) # Higher resolution
# 计算杆件的角度
def calculate_angle(p1, p2):
return np.arctan2(p2['y'] - p1['y'], p2['x'] - p1['x'])
# 值标签偏移以防重叠
def offset_label_position(x, y, theta, offset=0.5):
return x + offset * np.cos(theta), y + offset * np.sin(theta)
# 绘制节点和弯矩
for node in nodes:
ax.plot(node['x'], node['y'], 'ko')
ax.text(node['x'], node['y'], node['name'], fontsize=15, verticalalignment='bottom', horizontalalignment='right')
# 根据方向进行绘制
for direction, linked_node_name in node['link'].items():
if linked_node_name and linked_node_name != 'fixed':
linked_node = next((n for n in nodes if n['name'] == linked_node_name), None)
if linked_node:
theta = calculate_angle(node, linked_node)
# 判断杆件是竖直还是水平的
is_horizontal = np.isclose(theta, 0) or np.isclose(theta, np.pi) or np.isclose(theta, -np.pi)
direction_multiplier = -1 if is_horizontal else 1
moment_value = direction_multiplier * node['moment'].get(direction, 0)
m1_x = node['x'] + (moment_value * np.cos(theta + np.pi / 2) / 1000 * amp)
m1_y = node['y'] + (moment_value * np.sin(theta + np.pi / 2) / 1000 * amp)
next_moment_value = direction_multiplier * linked_node['moment'].get('down' if direction == 'up' else 'up' if direction == 'down' else 'right' if direction == 'left' else 'left', 0)
m2_x = linked_node['x'] - (next_moment_value * np.cos(theta + np.pi / 2) / 1000 * amp)
m2_y = linked_node['y'] - (next_moment_value * np.sin(theta + np.pi / 2) / 1000 * amp)
# 绘图
ax.plot(np.linspace(node['x'], linked_node['x'], 10), np.linspace(node['y'], linked_node['y'], 10), 'k-')
ax.plot(np.linspace(node['x'], m1_x, 10), np.linspace(node['y'], m1_y, 10), 'r-')
ax.plot(np.linspace(m1_x, m2_x, 10), np.linspace(m1_y, m2_y, 10), 'r-')
ax.plot(np.linspace(m2_x, linked_node['x'], 10), np.linspace(m2_y, linked_node['y'], 10), 'r-')
label_x, label_y = offset_label_position(node['x'], node['y'], theta + np.pi / 2)
ax.text(label_x, label_y, f'{abs(moment_value):.2f}', fontsize=12, color='blue', ha='center', va='center')
label_x, label_y = offset_label_position(linked_node['x'], linked_node['y'], theta - np.pi / 2)
ax.text(label_x, label_y, f'{abs(next_moment_value):.2f}', fontsize=12, color='blue', ha='center', va='center')
ax.set_title("弯矩图-D值法", fontsize=20)
plt.show()
这个函数遍历所有的节点，绘制节点和它们之间的连接关系。对于每个节点，根据其连接的其他节点的位置，计算杆件的角度和弯矩方向，然后在图上绘制相应的弯矩线。通过调整amp参数，可以控制弯矩线的放大比例，以便在图上清晰显示。在绘制完节点和杆件之后，该函数还会根据计算出的弯矩值在相应的杆件上绘制弯矩分布。
可视化页面搭建
为了更好地让D值法计算程序使用起来，搭建了便捷的网页程序，其网页的流程图如下所示。网页后端使用了python的flask框架，前端则是用html配合tailwind.css样式库进行。
- 用户访问网页：用户通过浏览器访问D值法自动计算的网页小程序。
- 网页界面加载：HTML前端页面加载，展示用户界面，包括文件上传区域、弯矩图展示区以及结果表格区。
- 用户上传CSV文件：用户选择并上传一个CSV格式的文件，该文件包含结构分析所需的数据。
- 服务器接收文件：后端Flask服务器接收来自前端的文件上传请求。
- 处理文件内容：服务器读取CSV文件内容，将其转化为所需的数据格式以进行D值法的计算。
- 计算D值法结果：后端利用处理后的数据进行D值法的计算，生成所需的弯矩图和结果数据。
- 返回数据到前端：计算结果以JSON格式返回给前端。
- 在网页上展示弯矩图和表格结果：前端接收到后端返回的数据后，使用D3.js等库在网页上动态生成弯矩图和展示计算结果的表格。
- 用户交互（放大/缩小弯矩图）：用户可以通过页面提供的按钮控制弯矩图的放大和缩小，以便更清晰地查看结构的受力情况。
效果展示
演示网页位于：https://asionm.github.io/structure-design/
源码下载
源码可访问 https://github.com/Asionm/structure-design 获取，感兴趣的朋友可以一起继续开发完善这系列的土木工程设计程序。
扫一扫
1645
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
