# Python桁架分析3步实现力学求解与可视化 ## 1. 问题解构 在Python中实现一个完整的桁架结构力学分析，通常需要完成三个核心步骤： 1. **结构建模与参数输入**：定义节点坐标、杆件连接关系、材料属性（如弹性模量EA）和荷载条件。 2. **力学方程建立与求解**：组装整体刚度矩阵，施加边界条件，求解节点位移，进而计算杆件内力。 3. **结果可视化**：绘制原始结构及变形后的结构，直观展示分析结果。 ## 2. 方案推演与代码实现 下面，我们结合具体的Python代码示例，分三步详细说明如何实现这一过程。 ### 第一步：结构建模与参数输入 我们需要定义节点、杆件单元，并初始化材料属性和荷载。这里，我们采用面向对象的方法来组织代码。 ```python import numpy as np import matplotlib.pyplot as plt # 定义节点类，存储节点ID和坐标 class Node: def __init__(self, node_id, x, y): self.id = node_id self.x = x self.y = y # 定义杆件单元类，存储单元ID、连接的节点ID、材料属性，并计算单元刚度矩阵 class TrussElement: def __init__(self, elem_id, node_i, node_j, EA): self.id = elem_id self.node_i = node_i # 节点i的ID self.node_j = node_j # 节点j的ID self.EA = EA # 轴向刚度 def compute_length(self, nodes_dict): """计算杆件长度 [ref_3]""" ni = nodes_dict[self.node_i] nj = nodes_dict[self.node_j] return np.sqrt((nj.x - ni.x)**2 + (nj.y - ni.y)**2) def compute_stiffness_matrix(self, nodes_dict): """计算局部坐标系下的单元刚度矩阵，并转换到全局坐标系 [ref_2][ref_6]""" ni = nodes_dict[self.node_i] nj = nodes_dict[self.node_j] L = self.compute_length(nodes_dict) c = (nj.x - ni.x) / L # 方向余弦 cosθ s = (nj.y - ni.y) / L # 方向正弦 sinθ # 局部刚度矩阵 (2x2 对于杆件轴向) k_local = (self.EA / L) * np.array([[1, -1], [-1, 1]]) # 转换矩阵 T T = np.array([[c, s, 0, 0], [0, 0, c, s]]) # 全局刚度矩阵 (4x4) k_global = T.T @ k_local @ T return k_global ``` ### 第二步：力学方程建立与求解 这一步是核心，需要组装总刚度矩阵，施加边界条件和荷载，并求解线性方程组。 ```python class TrussSolver: def __init__(self): self.nodes = {} # 存储所有节点对象，键为节点ID self.elements = [] # 存储所有杆件单元对象 self.EA = 1000 # 默认轴向刚度，单位：力 [ref_1] self.loads = {} # 存储节点荷载，键为节点ID，值为[Fx, Fy] self.constraints = {} # 存储约束条件，键为节点ID，值为[约束x, 约束y] (True/False) def add_node(self, node_id, x, y): self.nodes[node_id] = Node(node_id, x, y) def add_element(self, elem_id, node_i_id, node_j_id): elem = TrussElement(elem_id, node_i_id, node_j_id, self.EA) self.elements.append(elem) def apply_load(self, node_id, fx, fy): """在指定节点施加荷载 [ref_5]""" self.loads[node_id] = [fx, fy] def apply_constraint(self, node_id, fix_x, fix_y): """施加约束，True表示固定 [ref_2]""" self.constraints[node_id] = [fix_x, fix_y] def solve(self): """组装总刚度矩阵，施加边界条件，求解位移 [ref_2][ref_6]""" n_nodes = len(self.nodes) n_dof = 2 * n_nodes # 总自由度 (每个节点x,y方向) K_global = np.zeros((n_dof, n_dof)) F_global = np.zeros(n_dof) # 1. 组装总刚度矩阵 for elem in self.elements: k_elem = elem.compute_stiffness_matrix(self.nodes) # 获取单元在总刚度矩阵中的位置索引 i = int(elem.node_i) j = int(elem.node_j) dof_indices = [2*i, 2*i+1, 2*j, 2*j+1] for row_local, row_global in enumerate(dof_indices): for col_local, col_global in enumerate(dof_indices): K_global[row_global, col_global] += k_elem[row_local, col_local] # 2. 组装荷载向量 for node_id, load in self.loads.items(): idx = int(node_id) F_global[2*idx] = load[0] F_global[2*idx+1] = load[1] # 3. 处理边界条件（置1法） dof_to_keep = [] # 需要保留的自由度索引 for node_id, constraint in self.constraints.items(): idx = int(node_id) if constraint[0]: # x方向固定 K_global[2*idx, :] = 0 K_global[:, 2*idx] = 0 K_global[2*idx, 2*idx] = 1 F_global[2*idx] = 0 else: dof_to_keep.append(2*idx) if constraint[1]: # y方向固定 K_global[2*idx+1, :] = 0 K_global[:, 2*idx+1] = 0 K_global[2*idx+1, 2*idx+1] = 1 F_global[2*idx+1] = 0 else: dof_to_keep.append(2*idx+1) # 对于未约束的自由度，也加入保留列表 for i in range(n_dof): if i not in dof_to_keep and (i % 2 == 0 or i % 2 == 1): # 检查该自由度是否未被任何约束处理（即既不在constraints中明确固定，也未因约束而被排除） node_id_for_dof = i // 2 if str(node_id_for_dof) not in self.constraints: dof_to_keep.append(i) else: # 如果在约束字典中，但该特定方向未固定，则应该已经在上面添加了 pass # 去除重复并排序 dof_to_keep = sorted(list(set(dof_to_keep))) # 4. 求解缩减后的方程组 K_reduced = K_global[np.ix_(dof_to_keep, dof_to_keep)] F_reduced = F_global[dof_to_keep] U_reduced = np.linalg.solve(K_reduced, F_reduced) # 5. 将解得的位移赋值回完整的位移向量 U_full = np.zeros(n_dof) for i, dof in enumerate(dof_to_keep): U_full[dof] = U_reduced[i] # 6. 计算杆件内力 internal_forces = {} for elem in self.elements: ni = self.nodes[elem.node_i] nj = self.nodes[elem.node_j] L = elem.compute_length(self.nodes) c = (nj.x - ni.x) / L s = (nj.y - ni.y) / L # 节点位移 ui = U_full[2*int(elem.node_i)] vi = U_full[2*int(elem.node_i)+1] uj = U_full[2*int(elem.node_j)] vj = U_full[2*int(elem.node_j)+1] # 局部位移差 delta_L = c*(uj-ui) + s*(vj-vi) # 轴力 (拉力为正) N = (elem.EA / L) * delta_L internal_forces[elem.id] = N return U_full, internal_forces ``` ### 第三步：结果可视化 使用Matplotlib绘制结构变形前后的形态，并用颜色或箭头表示内力。 ```python def plot_results(self, U, scale_factor=50): """ 绘制桁架结构变形图。 scale_factor: 位移放大系数，以便观察变形。 """ plt.figure(figsize=(10, 6)) # 绘制原始结构 for elem in self.elements: ni = self.nodes[elem.node_i] nj = self.nodes[elem.node_j] plt.plot([ni.x, nj.x], [ni.y, nj.y], 'k-', linewidth=2, alpha=0.5, label='Original' if elem.id=='0' else "") # 绘制变形后的结构 for elem in self.elements: ni = self.nodes[elem.node_i] nj = self.nodes[elem.node_j] # 获取位移 ui, vi = U[2*int(elem.node_i)], U[2*int(elem.node_i)+1] uj, vj = U[2*int(elem.node_j)], U[2*int(elem.node_j)+1] # 绘制变形后的杆件 plt.plot([ni.x + scale_factor*ui, nj.x + scale_factor*uj], [ni.y + scale_factor*vi, nj.y + scale_factor*vj], 'r--', linewidth=1.5, label='Deformed' if elem.id=='0' else "") # 绘制节点 for node_id, node in self.nodes.items(): idx = int(node_id) plt.plot(node.x, node.y, 'bo', markersize=8) # 如果有位移，绘制位移向量 if abs(U[2*idx]) > 1e-10 or abs(U[2*idx+1]) > 1e-10: plt.arrow(node.x, node.y, scale_factor*U[2*idx], scale_factor*U[2*idx+1], head_width=0.05, head_length=0.1, fc='g', ec='g') plt.xlabel('X') plt.ylabel('Y') plt.title('Truss Structure Analysis (Original vs. Deformed)') plt.axis('equal') plt.grid(True, alpha=0.3) plt.legend() plt.show() ``` ## 3. 完整应用案例 下面，我们创建一个简单的两杆桁架算例，演示上述三步流程。 ```python # 主程序示例 [ref_1] if __name__ == "__main__": # 初始化求解器 solver = TrussSolver() # 第一步：建立模型 # 定义节点 (ID, x, y) solver.add_node('0', 0, 0) solver.add_node('1', 1, 0) solver.add_node('2', 0.5, 1) # 定义杆件单元 (单元ID, 起始节点ID, 终止节点ID) solver.add_element('0', '0', '2') solver.add_element('1', '1', '2') # 施加荷载 (在节点2施加竖向力) solver.apply_load('2', fx=0, fy=-10) # 向下10个单位力 [ref_1] # 施加约束 (节点0和节点1为铰接支座) solver.apply_constraint('0', fix_x=True, fix_y=True) solver.apply_constraint('1', fix_x=True, fix_y=True) # 第二步：求解 displacements, internal_forces = solver.solve() # 打印结果 print("节点位移:") for node_id, node in solver.nodes.items(): idx = int(node_id) print(f" 节点 {node_id}: ux = {displacements[2*idx]:.6f}, uy = {displacements[2*idx+1]:.6f}") print("\n杆件内力:") for elem_id, force in internal_forces.items(): print(f" 杆件 {elem_id}: 轴力 N = {force:.6f}") # 第三步：可视化 solver.plot_results(displacements, scale_factor=100) ``` ## 4. 关键技术与要点总结 下表对比了实现桁架分析时不同环节的核心技术与注意事项： | 分析环节 | 核心技术/公式 | 注意事项 | 参考来源 | | :--- | :--- | :--- | :--- | | **单元刚度矩阵** | `k_global = T.T @ k_local @ T`，其中`T`为坐标转换矩阵 | 确保方向余弦计算正确，长度不为零 | [ref_2], [ref_6] | | **总刚组装** | 将单元刚度矩阵按节点自由度编号“对号入座”叠加到总刚中 | 注意节点编号与自由度索引的映射关系 | [ref_2] | | **边界条件处理** | 采用“置1法”或“划行划列法”引入支座约束 | 处理约束后，方程组必须是正定的 | [ref_2] | | **方程求解** | 使用`np.linalg.solve`求解线性方程组`K * U = F` | 对于大型问题，需考虑稀疏矩阵求解器 | [ref_1] | | **内力计算** | `N = (EA / L) * (c*Δu + s*Δv)` | 结果为正值表示拉力，负值表示压力 | [ref_3], [ref_6] | | **结果可视化** | 使用Matplotlib绘制线条和箭头 | 变形通常需要放大显示，原始与变形图需清晰区分 | [ref_1] | 通过以上三步，我们实现了一个从建模、求解到可视化的完整Python桁架分析流程。此框架具有良好的扩展性，可以方便地添加更多节点、杆件、不同类型的荷载和约束，以分析更复杂的桁架结构。对于更复杂的静不定结构或动力分析，可以在总刚组装和求解步骤上进行扩展[ref_4][ref_5]。