结构动力学仿真-主题056-桥梁结构地震响应分析
·
主题056：桥梁结构地震响应分析
1. 引言
1.1 桥梁结构的地震响应特点
桥梁结构在地震作用下表现出独特的动力特性：
- 长跨径特性：桥梁跨径大，自振周期长，易受长周期地震动影响
- 多点激励效应：地震波传播需要时间，不同桥墩受到不同相位的地震动
- 行波效应：地震波沿桥长方向传播，产生行波效应
- 支座非线性：支座在地震作用下可能发生滑动或脱落
- 碰撞效应：相邻跨之间可能发生碰撞
1.2 桥梁结构抗震设计要点
桥梁结构抗震设计需要考虑的关键因素：
-
多点激励分析
- 考虑地震波传播效应
- 不同支承点的地震动时程差异
- 行波效应的影响
-
支座设计
- 固定支座与活动支座的配合
- 隔震支座的应用
- 支座抗震验算
-
延性设计
- 桥墩延性设计
- 塑性铰的形成与控制
- 位移能力验算
-
碰撞防护
- 伸缩缝设计
- 碰撞力计算
- 防护措施
2. 理论基础
2.1 桥梁结构动力模型
2.1.1 梁单元模型
桥梁主梁可采用欧拉-伯努利梁单元：
EI∂4w∂x4+m∂2w∂t2=p(x,t)EI \frac{\partial^4 w}{\partial x^4} + m \frac{\partial^2 w}{\partial t^2} = p(x,t)EI∂x4∂4w+m∂t2∂2w=p(x,t)
其中：
- EIEIEI：抗弯刚度
- mmm：单位长度质量
- www：横向位移
- p(x,t)p(x,t)p(x,t)：分布荷载
2.1.2 桥墩模型
桥墩可简化为悬臂柱模型：
fn=12π3EImL3f_n = \frac{1}{2\pi} \sqrt{\frac{3EI}{mL^3}}fn=2π1mL33EI
其中：
- LLL：桥墩高度
- EIEIEI：桥墩抗弯刚度
- mmm：桥墩质量
2.2 多点激励运动方程
考虑多点激励的桥梁结构运动方程：
[MssMsgMgsMgg][u¨su¨g]+[CssCsgCgsCgg][u˙su˙g]+[KssKsgKgsKgg][usug]=[0Pg]\begin{bmatrix} M_{ss} & M_{sg} \\ M_{gs} & M_{gg} \end{bmatrix} \begin{bmatrix} \ddot{u}_s \\ \ddot{u}_g \end{bmatrix} + \begin{bmatrix} C_{ss} & C_{sg} \\ C_{gs} & C_{gg} \end{bmatrix} \begin{bmatrix} \dot{u}_s \\ \dot{u}_g \end{bmatrix} + \begin{bmatrix} K_{ss} & K_{sg} \\ K_{gs} & K_{gg} \end{bmatrix} \begin{bmatrix} u_s \\ u_g \end{bmatrix} = \begin{bmatrix} 0 \\ P_g \end{bmatrix}[MssMgsMsgMgg][u¨su¨g]+[CssCgsCsgCgg][u˙su˙g]+[KssKgsKsgKgg][usug]=[0Pg]
其中：
- 下标sss：结构自由度
- 下标ggg：支座自由度
2.3 行波效应
地震波以波速vvv传播，到达不同支承点的时间差：
Δt=Δxv\Delta t = \frac{\Delta x}{v}Δt=vΔx
其中：
- Δx\Delta xΔx：支承点间距
- vvv：地震波传播速度（一般取1000-3000 m/s）
2.4 地震响应分析方法
2.4.1 反应谱法
Sa=∑j=1mSaj2S_a = \sqrt{\sum_{j=1}^{m} S_{aj}^2}Sa=j=1∑mSaj2
2.4.2 时程分析法
直接求解运动方程，考虑地震动时程特性。
2.4.3 随机振动法
考虑地震动的随机特性，进行随机振动分析。
3. 案例实现
3.1 案例1：桥梁结构建模与模态分析
目标：建立三跨连续梁桥模型，进行模态分析
内容：
- 建立梁单元模型
- 考虑桥墩刚度
- 计算前10阶模态
- 分析振型特点
"""
案例1：桥梁结构建模与模态分析
======================================
本案例演示如何建立桥梁结构的有限元模型并进行模态分析，包括：
1. 三跨连续梁桥建模
2. 桥墩-主梁系统建模
3. 模态分析
4. 振型可视化
作者：结构动力学仿真团队
日期：2026-03-09
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy import linalg
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class BridgeElement:
"""桥梁单元类 - 用于梁单元"""
def __init__(self, node1, node2, E, I, A, rho):
"""
初始化梁单元
Parameters:
-----------
node1, node2 : array
节点坐标 [x, y]
E : float
弹性模量 (Pa)
I : float
截面惯性矩 (m^4)
A : float
截面面积 (m^2)
rho : float
材料密度 (kg/m^3)
"""
self.node1 = np.array(node1)
self.node2 = np.array(node2)
self.E = E
self.I = I
self.A = A
self.rho = rho
# 计算单元长度
self.L = np.linalg.norm(self.node2 - self.node1)
# 计算方向余弦
dx = self.node2[0] - self.node1[0]
dy = self.node2[1] - self.node1[1]
self.c = dx / self.L # cos(theta)
self.s = dy / self.L # sin(theta)
def get_stiffness_matrix(self):
"""
获取局部坐标系下的单元刚度矩阵（平面梁单元，6自由度）
Returns:
--------
k : ndarray
6x6 刚度矩阵
"""
L = self.L
E = self.E
I = self.I
A = self.A
# 局部刚度矩阵
k_local = np.zeros((6, 6))
# 轴向刚度
EA_L = E * A / L
k_local[0, 0] = EA_L
k_local[0, 3] = -EA_L
k_local[3, 0] = -EA_L
k_local[3, 3] = EA_L
# 弯曲刚度
EI = E * I
k_local[1, 1] = 12 * EI / L**3
k_local[1, 2] = 6 * EI / L**2
k_local[1, 4] = -12 * EI / L**3
k_local[1, 5] = 6 * EI / L**2
k_local[2, 1] = 6 * EI / L**2
k_local[2, 2] = 4 * EI / L
k_local[2, 4] = -6 * EI / L**2
k_local[2, 5] = 2 * EI / L
k_local[4, 1] = -12 * EI / L**3
k_local[4, 2] = -6 * EI / L**2
k_local[4, 4] = 12 * EI / L**3
k_local[4, 5] = -6 * EI / L**2
k_local[5, 1] = 6 * EI / L**2
k_local[5, 2] = 2 * EI / L
k_local[5, 4] = -6 * EI / L**2
k_local[5, 5] = 4 * EI / L
return k_local
def get_mass_matrix(self, consistent=True):
"""
获取单元质量矩阵
Parameters:
-----------
consistent : bool
是否使用一致质量矩阵
Returns:
--------
m : ndarray
6x6 质量矩阵
"""
L = self.L
rho = self.rho
A = self.A
if consistent:
# 一致质量矩阵
m_local = np.zeros((6, 6))
# 轴向质量
m_axial = rho * A * L / 6
m_local[0, 0] = 2 * m_axial
m_local[0, 3] = m_axial
m_local[3, 0] = m_axial
m_local[3, 3] = 2 * m_axial
# 弯曲质量
m_bend = rho * A * L / 420
m_local[1, 1] = 156 * m_bend
m_local[1, 2] = 22 * L * m_bend
m_local[1, 4] = 54 * m_bend
m_local[1, 5] = -13 * L * m_bend
m_local[2, 1] = 22 * L * m_bend
m_local[2, 2] = 4 * L**2 * m_bend
m_local[2, 4] = 13 * L * m_bend
m_local[2, 5] = -3 * L**2 * m_bend
m_local[4, 1] = 54 * m_bend
m_local[4, 2] = 13 * L * m_bend
m_local[4, 4] = 156 * m_bend
m_local[4, 5] = -22 * L * m_bend
m_local[5, 1] = -13 * L * m_bend
m_local[5, 2] = -3 * L**2 * m_bend
m_local[5, 4] = -22 * L * m_bend
m_local[5, 5] = 4 * L**2 * m_bend
return m_local
else:
# 集中质量矩阵
m_local = np.zeros((6, 6))
m_total = rho * A * L
m_local[0, 0] = m_total / 2
m_local[1, 1] = m_total / 2
m_local[3, 3] = m_total / 2
m_local[4, 4] = m_total / 2
# 转动惯量（简化）
m_local[2, 2] = m_total * L**2 / 24
m_local[5, 5] = m_total * L**2 / 24
return m_local
def get_transformation_matrix(self):
"""
获取坐标变换矩阵
Returns:
--------
T : ndarray
6x6 变换矩阵
"""
c = self.c
s = self.s
T = np.array([
[c, s, 0, 0, 0, 0],
[-s, c, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, c, s, 0],
[0, 0, 0, -s, c, 0],
[0, 0, 0, 0, 0, 1]
])
return T
def get_global_stiffness(self):
"""获取全局坐标系下的刚度矩阵"""
k_local = self.get_stiffness_matrix()
T = self.get_transformation_matrix()
return T.T @ k_local @ T
def get_global_mass(self, consistent=True):
"""获取全局坐标系下的质量矩阵"""
m_local = self.get_mass_matrix(consistent)
T = self.get_transformation_matrix()
return T.T @ m_local @ T
class BridgeStructure:
"""桥梁结构类"""
def __init__(self, name="桥梁结构"):
"""
初始化桥梁结构
Parameters:
-----------
name : str
结构名称
"""
self.name = name
self.nodes = [] # 节点列表 [[x, y], ...]
self.elements = [] # 单元列表 [BridgeElement, ...]
self.node_dofs = [] # 节点自由度映射
self.n_dof = 0 # 总自由度数
self.K = None # 全局刚度矩阵
self.M = None # 全局质量矩阵
self.fixed_dofs = [] # 约束自由度
def add_node(self, x, y):
"""
添加节点
Parameters:
-----------
x, y : float
节点坐标
Returns:
--------
node_id : int
节点编号
"""
self.nodes.append([x, y])
return len(self.nodes) - 1
def add_element(self, node1_id, node2_id, E, I, A, rho):
"""
添加梁单元
Parameters:
-----------
node1_id, node2_id : int
节点编号
E, I, A, rho : float
材料参数
"""
elem = BridgeElement(
self.nodes[node1_id],
self.nodes[node2_id],
E, I, A, rho
)
self.elements.append({
'elem': elem,
'nodes': [node1_id, node2_id]
})
def fix_node(self, node_id, ux=True, uy=True, rz=True):
"""
约束节点自由度
Parameters:
-----------
node_id : int
节点编号
ux, uy, rz : bool
是否约束x、y位移和转动
"""
dof_base = node_id * 3
if ux:
self.fixed_dofs.append(dof_base)
if uy:
self.fixed_dofs.append(dof_base + 1)
if rz:
self.fixed_dofs.append(dof_base + 2)
def assemble_global_matrices(self):
"""组装全局刚度矩阵和质量矩阵"""
n_nodes = len(self.nodes)
self.n_dof = n_nodes * 3
self.K = np.zeros((self.n_dof, self.n_dof))
self.M = np.zeros((self.n_dof, self.n_dof))
for elem_info in self.elements:
elem = elem_info['elem']
n1, n2 = elem_info['nodes']
# 获取全局矩阵
k_global = elem.get_global_stiffness()
m_global = elem.get_global_mass(consistent=True)
# 组装到全局矩阵
dofs = []
for n in [n1, n2]:
dofs.extend([n*3, n*3+1, n*3+2])
for i in range(6):
for j in range(6):
gi, gj = dofs[i], dofs[j]
self.K[gi, gj] += k_global[i, j]
self.M[gi, gj] += m_global[i, j]
def modal_analysis(self, n_modes=10):
"""
模态分析
Parameters:
-----------
n_modes : int
计算的模态数
Returns:
--------
freqs : array
固有频率 (Hz)
modes : ndarray
振型矩阵
"""
# 应用边界条件
free_dofs = [i for i in range(self.n_dof) if i not in self.fixed_dofs]
K_reduced = self.K[np.ix_(free_dofs, free_dofs)]
M_reduced = self.M[np.ix_(free_dofs, free_dofs)]
# 求解广义特征值问题
n_modes = min(n_modes, len(free_dofs))
eigenvalues, eigenvectors = linalg.eigh(K_reduced, M_reduced,
subset_by_index=[0, n_modes-1])
# 计算频率
freqs = np.sqrt(np.maximum(eigenvalues, 0)) / (2 * np.pi)
# 扩展振型到完整自由度
modes = np.zeros((self.n_dof, n_modes))
for i, fdof in enumerate(free_dofs):
modes[fdof, :] = eigenvectors[i, :]
# 归一化振型
for i in range(n_modes):
max_disp = np.max(np.abs(modes[:, i]))
if max_disp > 0:
modes[:, i] /= max_disp
return freqs, modes
def get_node_displacements(self, mode, scale=1.0):
"""
获取节点位移（用于可视化）
Parameters:
-----------
mode : array
振型向量
scale : float
显示缩放系数
Returns:
--------
displaced_nodes : ndarray
变形后的节点坐标
"""
displaced = np.array(self.nodes).copy()
for i, (x, y) in enumerate(self.nodes):
dx = mode[i*3] * scale
dy = mode[i*3+1] * scale
displaced[i] = [x + dx, y + dy]
return displaced
def create_three_span_bridge():
"""
创建三跨连续梁桥模型
桥梁参数：
- 跨径布置：40m + 60m + 40m = 140m
- 桥墩高度：15m
- 主梁：混凝土箱梁
- 桥墩：钢筋混凝土双柱墩
Returns:
--------
bridge : BridgeStructure
桥梁结构对象
"""
bridge = BridgeStructure("三跨连续梁桥")
# 材料参数
E_concrete = 3.45e10 # 混凝土弹性模量 (Pa)
rho_concrete = 2500 # 混凝土密度 (kg/m^3)
# 主梁截面参数（混凝土箱梁）
A_girder = 8.5 # 面积 (m^2)
I_girder = 12.0 # 惯性矩 (m^4)
# 桥墩截面参数（圆形墩，直径2.5m）
D_pier = 2.5
A_pier = np.pi * (D_pier/2)**2
I_pier = np.pi * (D_pier/2)**4 / 4
# 跨径布置
L1, L2, L3 = 40, 60, 40 # 三跨长度
H_pier = 15 # 桥墩高度
# 创建节点
# 桥墩底部（固定端）
pier1_bottom = bridge.add_node(40, 0)
pier2_bottom = bridge.add_node(100, 0)
# 桥墩顶部（主梁支承点）
pier1_top = bridge.add_node(40, H_pier)
pier2_top = bridge.add_node(100, H_pier)
# 主梁节点（每跨8个单元，9个节点）
girder_nodes = []
# 第一跨（0-40m）
for i in range(9):
x = i * L1 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
# 第二跨（40-100m）
for i in range(1, 9):
x = 40 + i * L2 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
# 第三跨（100-140m）
for i in range(1, 9):
x = 100 + i * L3 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
# 添加桥墩单元
bridge.add_element(pier1_bottom, pier1_top, E_concrete, I_pier, A_pier, rho_concrete)
bridge.add_element(pier2_bottom, pier2_top, E_concrete, I_pier, A_pier, rho_concrete)
# 添加主梁单元
for i in range(len(girder_nodes) - 1):
bridge.add_element(girder_nodes[i], girder_nodes[i+1],
E_concrete, I_girder, A_girder, rho_concrete)
# 约束桥墩底部（完全固定）
bridge.fix_node(pier1_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(pier2_bottom, ux=True, uy=True, rz=True)
# 主梁端部约束（简支）
bridge.fix_node(girder_nodes[0], uy=True)
bridge.fix_node(girder_nodes[-1], uy=True)
return bridge
def plot_bridge_model(bridge, title="桥梁结构模型"):
"""
绘制桥梁结构模型
Parameters:
-----------
bridge : BridgeStructure
桥梁结构对象
title : str
图标题
"""
fig, ax = plt.subplots(figsize=(14, 8))
nodes = np.array(bridge.nodes)
# 绘制单元
for elem_info in bridge.elements:
n1, n2 = elem_info['nodes']
x = [nodes[n1, 0], nodes[n2, 0]]
y = [nodes[n1, 1], nodes[n2, 1]]
ax.plot(x, y, 'b-', linewidth=2)
# 绘制节点
ax.scatter(nodes[:, 0], nodes[:, 1], c='red', s=50, zorder=5)
# 标注节点
for i, (x, y) in enumerate(nodes):
ax.annotate(f'{i}', (x, y), textcoords="offset points",
xytext=(0, 10), ha='center', fontsize=8)
# 绘制支座
for node_id in [0, 26]: # 主梁端部
x, y = nodes[node_id]
ax.plot([x-1, x+1], [y, y], 'k-', linewidth=3)
ax.plot([x], [y-0.5], 'k^', markersize=10)
# 绘制固定端
for node_id in [0, 1]: # 桥墩底部
x, y = nodes[node_id]
ax.plot([x-1.5, x+1.5], [y, y], 'k-', linewidth=4)
for dx in np.linspace(-1.2, 1.2, 5):
ax.plot([x+dx, x+dx-0.3], [y, y-0.5], 'k-', linewidth=1)
ax.set_xlabel('X (m)', fontsize=12)
ax.set_ylabel('Y (m)', fontsize=12)
ax.set_title(title, fontsize=14)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def plot_mode_shapes(bridge, freqs, modes, n_plot=6):
"""
绘制振型图
Parameters:
-----------
bridge : BridgeStructure
桥梁结构对象
freqs : array
固有频率
modes : ndarray
振型矩阵
n_plot : int
绘制的模态数
"""
n_plot = min(n_plot, len(freqs))
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
nodes = np.array(bridge.nodes)
scale = 5.0 # 显示缩放系数
for i in range(n_plot):
ax = axes[i]
# 原始结构
for elem_info in bridge.elements:
n1, n2 = elem_info['nodes']
x = [nodes[n1, 0], nodes[n2, 0]]
y = [nodes[n1, 1], nodes[n2, 1]]
ax.plot(x, y, 'gray', linewidth=1, alpha=0.5, linestyle='--')
# 变形后的结构
displaced = bridge.get_node_displacements(modes[:, i], scale)
for elem_info in bridge.elements:
n1, n2 = elem_info['nodes']
x = [displaced[n1, 0], displaced[n2, 0]]
y = [displaced[n1, 1], displaced[n2, 1]]
ax.plot(x, y, 'b-', linewidth=2)
ax.scatter(displaced[:, 0], displaced[:, 1], c='red', s=30, zorder=5)
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_title(f'第{i+1}阶模态\nf = {freqs[i]:.3f} Hz\nT = {1/freqs[i]:.3f} s')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.suptitle('桥梁结构振型', fontsize=16)
plt.tight_layout()
return fig
def create_mode_animation(bridge, freqs, modes, mode_idx=0, duration=3.0):
"""
创建振型动画
Parameters:
-----------
bridge : BridgeStructure
桥梁结构对象
freqs : array
固有频率
modes : ndarray
振型矩阵
mode_idx : int
要动画显示的模态索引
duration : float
动画持续时间（秒）
Returns:
--------
anim : FuncAnimation
动画对象
"""
fig, ax = plt.subplots(figsize=(12, 6))
nodes = np.array(bridge.nodes)
scale = 3.0
# 绘制原始结构
lines_orig = []
for elem_info in bridge.elements:
n1, n2 = elem_info['nodes']
x = [nodes[n1, 0], nodes[n2, 0]]
y = [nodes[n1, 1], nodes[n2, 1]]
line, = ax.plot(x, y, 'gray', linewidth=1, alpha=0.5, linestyle='--')
lines_orig.append(line)
# 变形结构线条
lines_def = []
for _ in bridge.elements:
line, = ax.plot([], [], 'b-', linewidth=2)
lines_def.append(line)
# 节点散点
scatter = ax.scatter([], [], c='red', s=30, zorder=5)
ax.set_xlabel('X (m)', fontsize=12)
ax.set_ylabel('Y (m)', fontsize=12)
ax.set_title(f'第{mode_idx+1}阶模态动画 (f = {freqs[mode_idx]:.3f} Hz)', fontsize=14)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
# 设置坐标范围
margin = 10
ax.set_xlim(nodes[:, 0].min() - margin, nodes[:, 0].max() + margin)
ax.set_ylim(nodes[:, 1].min() - margin, nodes[:, 1].max() + margin)
n_frames = 60
def animate(frame):
t = frame / n_frames * 2 * np.pi
factor = np.sin(t)
displaced = bridge.get_node_displacements(modes[:, mode_idx], scale * factor)
# 更新线条
for line, elem_info in zip(lines_def, bridge.elements):
n1, n2 = elem_info['nodes']
x = [displaced[n1, 0], displaced[n2, 0]]
y = [displaced[n1, 1], displaced[n2, 1]]
line.set_data(x, y)
# 更新散点
scatter.set_offsets(displaced)
return lines_def + [scatter]
anim = FuncAnimation(fig, animate, frames=n_frames, interval=50, blit=True)
return anim, fig
def print_modal_results(freqs, modes, bridge):
"""
打印模态分析结果
Parameters:
-----------
freqs : array
固有频率
modes : ndarray
振型矩阵
bridge : BridgeStructure
桥梁结构对象
"""
print("=" * 60)
print("桥梁结构模态分析结果")
print("=" * 60)
print(f"\n结构名称: {bridge.name}")
print(f"总自由度数: {bridge.n_dof}")
print(f"约束自由度数: {len(bridge.fixed_dofs)}")
print(f"有效自由度数: {bridge.n_dof - len(bridge.fixed_dofs)}")
print("\n" + "-" * 60)
print("固有频率和周期")
print("-" * 60)
print(f"{'阶数':<8}{'频率(Hz)':<15}{'周期(s)':<15}{'圆频率(rad/s)':<15}")
print("-" * 60)
for i in range(len(freqs)):
f = freqs[i]
T = 1 / f if f > 0 else float('inf')
omega = 2 * np.pi * f
print(f"{i+1:<8}{f:<15.4f}{T:<15.4f}{omega:<15.4f}")
print("\n" + "-" * 60)
print("振型特征描述")
print("-" * 60)
# 分析振型特征
girder_dofs = []
for i in range(len(bridge.nodes)):
# 判断是否是主梁节点（y坐标接近桥墩高度）
if abs(bridge.nodes[i][1] - 15) < 0.1:
girder_dofs.extend([i*3, i*3+1])
for i in range(min(6, len(freqs))):
mode = modes[:, i]
# 计算主梁的最大位移
max_disp = 0
max_disp_idx = 0
for dof in girder_dofs:
if abs(mode[dof]) > max_disp:
max_disp = abs(mode[dof])
max_disp_idx = dof
# 判断振型类型
if i == 0:
desc = "主梁一阶竖向弯曲"
elif i == 1:
desc = "主梁一阶横向弯曲或扭转"
elif i == 2:
desc = "主梁二阶竖向弯曲"
elif i == 3:
desc = "桥墩一阶弯曲"
elif i == 4:
desc = "主梁三阶竖向弯曲"
else:
desc = "高阶振动"
print(f"第{i+1}阶: {desc}")
def main():
"""主程序"""
print("=" * 70)
print("案例1：桥梁结构建模与模态分析")
print("=" * 70)
# 创建桥梁模型
print("\n[1] 创建三跨连续梁桥模型...")
bridge = create_three_span_bridge()
print(f" - 节点数: {len(bridge.nodes)}")
print(f" - 单元数: {len(bridge.elements)}")
# 组装全局矩阵
print("\n[2] 组装全局刚度矩阵和质量矩阵...")
bridge.assemble_global_matrices()
print(f" - 刚度矩阵维度: {bridge.K.shape}")
print(f" - 质量矩阵维度: {bridge.M.shape}")
# 模态分析
print("\n[3] 进行模态分析...")
n_modes = 10
freqs, modes = bridge.modal_analysis(n_modes)
print(f" - 计算模态数: {n_modes}")
# 打印结果
print_modal_results(freqs, modes, bridge)
# 绘制结构模型
print("\n[4] 绘制结构模型...")
fig_model = plot_bridge_model(bridge, "三跨连续梁桥结构模型")
plt.savefig('case1_bridge_model.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case1_bridge_model.png")
plt.close()
# 绘制振型
print("\n[5] 绘制振型图...")
fig_modes = plot_mode_shapes(bridge, freqs, modes, n_plot=6)
plt.savefig('case1_mode_shapes.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case1_mode_shapes.png")
plt.close()
# 创建动画
print("\n[6] 创建振型动画...")
for i in range(min(3, len(freqs))):
anim, fig = create_mode_animation(bridge, freqs, modes, mode_idx=i)
anim.save(f'case1_mode_{i+1}_animation.gif', writer='pillow', fps=20)
print(f" - 已保存: case1_mode_{i+1}_animation.gif")
plt.close()
# 绘制频率分布图
print("\n[7] 绘制频率分布图...")
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 频率柱状图
ax1 = axes[0]
bars = ax1.bar(range(1, len(freqs)+1), freqs, color='steelblue', edgecolor='black')
ax1.set_xlabel('模态阶数', fontsize=12)
ax1.set_ylabel('频率 (Hz)', fontsize=12)
ax1.set_title('桥梁结构固有频率分布', fontsize=14)
ax1.set_xticks(range(1, len(freqs)+1))
ax1.grid(True, alpha=0.3, axis='y')
# 在柱子上标注数值
for bar, f in zip(bars, freqs):
height = bar.get_height()
ax1.annotate(f'{f:.2f}',
xy=(bar.get_x() + bar.get_width() / 2, height),
xytext=(0, 3),
textcoords="offset points",
ha='center', va='bottom', fontsize=8)
# 周期图
ax2 = axes[1]
periods = [1/f if f > 0 else float('inf') for f in freqs]
ax2.plot(range(1, len(periods)+1), periods, 'ro-', linewidth=2, markersize=8)
ax2.set_xlabel('模态阶数', fontsize=12)
ax2.set_ylabel('周期 (s)', fontsize=12)
ax2.set_title('桥梁结构自振周期', fontsize=14)
ax2.set_xticks(range(1, len(periods)+1))
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case1_frequency_distribution.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case1_frequency_distribution.png")
plt.close()
print("\n" + "=" * 70)
print("案例分析完成！")
print("=" * 70)
print("\n生成文件:")
print(" - case1_bridge_model.png: 桥梁结构模型图")
print(" - case1_mode_shapes.png: 前6阶振型图")
print(" - case1_mode_1_animation.gif: 第1阶振型动画")
print(" - case1_mode_2_animation.gif: 第2阶振型动画")
print(" - case1_mode_3_animation.gif: 第3阶振型动画")
print(" - case1_frequency_distribution.png: 频率分布图")
if __name__ == "__main__":
main()
3.2 案例2：多点激励地震响应分析
目标：考虑多点激励效应，计算地震响应
内容：
- 生成多点地震动
- 考虑行波效应
- 时程积分计算
- 响应对比分析
"""
案例2：多点激励地震响应分析
======================================
本案例演示桥梁结构在多点地震激励下的响应分析，包括：
1. 不同支承点地震动的生成
2. 多点激励响应分析方法
3. 一致激励与多点激励对比
4. 结构响应可视化
作者：结构动力学仿真团队
日期：2026-03-09
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy import linalg
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class BridgeElement:
"""桥梁单元类 - 用于梁单元"""
def __init__(self, node1, node2, E, I, A, rho):
self.node1 = np.array(node1)
self.node2 = np.array(node2)
self.E = E
self.I = I
self.A = A
self.rho = rho
self.L = np.linalg.norm(self.node2 - self.node1)
dx = self.node2[0] - self.node1[0]
dy = self.node2[1] - self.node1[1]
self.c = dx / self.L
self.s = dy / self.L
def get_stiffness_matrix(self):
"""获取局部坐标系下的单元刚度矩阵"""
L = self.L
E = self.E
I = self.I
A = self.A
k_local = np.zeros((6, 6))
EA_L = E * A / L
k_local[0, 0] = EA_L
k_local[0, 3] = -EA_L
k_local[3, 0] = -EA_L
k_local[3, 3] = EA_L
EI = E * I
k_local[1, 1] = 12 * EI / L**3
k_local[1, 2] = 6 * EI / L**2
k_local[1, 4] = -12 * EI / L**3
k_local[1, 5] = 6 * EI / L**2
k_local[2, 1] = 6 * EI / L**2
k_local[2, 2] = 4 * EI / L
k_local[2, 4] = -6 * EI / L**2
k_local[2, 5] = 2 * EI / L
k_local[4, 1] = -12 * EI / L**3
k_local[4, 2] = -6 * EI / L**2
k_local[4, 4] = 12 * EI / L**3
k_local[4, 5] = -6 * EI / L**2
k_local[5, 1] = 6 * EI / L**2
k_local[5, 2] = 2 * EI / L
k_local[5, 4] = -6 * EI / L**2
k_local[5, 5] = 4 * EI / L
return k_local
def get_mass_matrix(self, consistent=True):
"""获取单元质量矩阵"""
L = self.L
rho = self.rho
A = self.A
if consistent:
m_local = np.zeros((6, 6))
m_axial = rho * A * L / 6
m_local[0, 0] = 2 * m_axial
m_local[0, 3] = m_axial
m_local[3, 0] = m_axial
m_local[3, 3] = 2 * m_axial
m_bend = rho * A * L / 420
m_local[1, 1] = 156 * m_bend
m_local[1, 2] = 22 * L * m_bend
m_local[1, 4] = 54 * m_bend
m_local[1, 5] = -13 * L * m_bend
m_local[2, 1] = 22 * L * m_bend
m_local[2, 2] = 4 * L**2 * m_bend
m_local[2, 4] = 13 * L * m_bend
m_local[2, 5] = -3 * L**2 * m_bend
m_local[4, 1] = 54 * m_bend
m_local[4, 2] = 13 * L * m_bend
m_local[4, 4] = 156 * m_bend
m_local[4, 5] = -22 * L * m_bend
m_local[5, 1] = -13 * L * m_bend
m_local[5, 2] = -3 * L**2 * m_bend
m_local[5, 4] = -22 * L * m_bend
m_local[5, 5] = 4 * L**2 * m_bend
return m_local
else:
m_local = np.zeros((6, 6))
m_total = rho * A * L
m_local[0, 0] = m_total / 2
m_local[1, 1] = m_total / 2
m_local[3, 3] = m_total / 2
m_local[4, 4] = m_total / 2
m_local[2, 2] = m_total * L**2 / 24
m_local[5, 5] = m_total * L**2 / 24
return m_local
def get_transformation_matrix(self):
"""获取坐标变换矩阵"""
c = self.c
s = self.s
T = np.array([
[c, s, 0, 0, 0, 0],
[-s, c, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, c, s, 0],
[0, 0, 0, -s, c, 0],
[0, 0, 0, 0, 0, 1]
])
return T
def get_global_stiffness(self):
"""获取全局坐标系下的刚度矩阵"""
k_local = self.get_stiffness_matrix()
T = self.get_transformation_matrix()
return T.T @ k_local @ T
def get_global_mass(self, consistent=True):
"""获取全局坐标系下的质量矩阵"""
m_local = self.get_mass_matrix(consistent)
T = self.get_transformation_matrix()
return T.T @ m_local @ T
class BridgeStructure:
"""桥梁结构类"""
def __init__(self, name="桥梁结构"):
self.name = name
self.nodes = []
self.elements = []
self.n_dof = 0
self.K = None
self.M = None
self.fixed_dofs = []
self.support_nodes = []
def add_node(self, x, y):
"""添加节点"""
self.nodes.append([x, y])
return len(self.nodes) - 1
def add_element(self, node1_id, node2_id, E, I, A, rho):
"""添加梁单元"""
elem = BridgeElement(
self.nodes[node1_id],
self.nodes[node2_id],
E, I, A, rho
)
self.elements.append({
'elem': elem,
'nodes': [node1_id, node2_id]
})
def fix_node(self, node_id, ux=True, uy=True, rz=True):
"""约束节点自由度"""
dof_base = node_id * 3
if ux:
self.fixed_dofs.append(dof_base)
if uy:
self.fixed_dofs.append(dof_base + 1)
if rz:
self.fixed_dofs.append(dof_base + 2)
def add_support(self, node_id):
"""添加支承节点"""
self.support_nodes.append(node_id)
def assemble_global_matrices(self):
"""组装全局刚度矩阵和质量矩阵"""
n_nodes = len(self.nodes)
self.n_dof = n_nodes * 3
self.K = np.zeros((self.n_dof, self.n_dof))
self.M = np.zeros((self.n_dof, self.n_dof))
for elem_info in self.elements:
elem = elem_info['elem']
n1, n2 = elem_info['nodes']
k_global = elem.get_global_stiffness()
m_global = elem.get_global_mass(consistent=True)
dofs = []
for n in [n1, n2]:
dofs.extend([n*3, n*3+1, n*3+2])
for i in range(6):
for j in range(6):
gi, gj = dofs[i], dofs[j]
self.K[gi, gj] += k_global[i, j]
self.M[gi, gj] += m_global[i, j]
def modal_analysis(self, n_modes=10):
"""模态分析"""
free_dofs = [i for i in range(self.n_dof) if i not in self.fixed_dofs]
K_reduced = self.K[np.ix_(free_dofs, free_dofs)]
M_reduced = self.M[np.ix_(free_dofs, free_dofs)]
n_modes = min(n_modes, len(free_dofs))
eigenvalues, eigenvectors = linalg.eigh(K_reduced, M_reduced,
subset_by_index=[0, n_modes-1])
freqs = np.sqrt(np.maximum(eigenvalues, 0)) / (2 * np.pi)
modes = np.zeros((self.n_dof, n_modes))
for i, fdof in enumerate(free_dofs):
modes[fdof, :] = eigenvectors[i, :]
for i in range(n_modes):
max_disp = np.max(np.abs(modes[:, i]))
if max_disp > 0:
modes[:, i] /= max_disp
return freqs, modes, free_dofs
def generate_ground_motion(dt=0.01, T=20, PGA=0.2, seed=None):
"""生成人工地震动"""
if seed is not None:
np.random.seed(seed)
n_steps = int(T / dt)
t = np.linspace(0, T, n_steps)
omega_g = 15.0
zeta_g = 0.6
t1, t2 = 2, 8
envelope = np.zeros(n_steps)
for i, ti in enumerate(t):
if ti < t1:
envelope[i] = (ti / t1) ** 2
elif ti < t2:
envelope[i] = 1.0
else:
envelope[i] = np.exp(-0.5 * (ti - t2))
white_noise = np.random.randn(n_steps)
from scipy.signal import lfilter
omega_g2 = omega_g ** 2
num = [omega_g2]
den = [1, 2*zeta_g*omega_g, omega_g2]
filtered = lfilter(num, den, white_noise)
acc = filtered * envelope
# 避免除以零
max_acc = np.max(np.abs(acc))
if max_acc > 1e-10:
acc = acc / max_acc * PGA * 9.81
else:
acc = white_noise * envelope * PGA * 9.81
return t, acc
def generate_multi_support_motions(bridge, dt=0.01, T=20, PGA=0.2,
correlation='partial', coherence_decay=0.001):
"""生成多点地震动"""
n_steps = int(T / dt)
t = np.linspace(0, T, n_steps)
motions = {}
if correlation == 'full':
_, acc = generate_ground_motion(dt, T, PGA, seed=42)
for node_id in bridge.support_nodes:
motions[node_id] = acc.copy()
elif correlation == 'none':
for i, node_id in enumerate(bridge.support_nodes):
_, acc = generate_ground_motion(dt, T, PGA, seed=42+i)
motions[node_id] = acc
else: # partial
base_motion = generate_ground_motion(dt, T, PGA, seed=42)[1]
for i, node_id in enumerate(bridge.support_nodes):
x_i = bridge.nodes[node_id][0]
_, independent = generate_ground_motion(dt, T, PGA*0.3, seed=100+i)
if i == 0:
coherence = 1.0
else:
x_0 = bridge.nodes[bridge.support_nodes[0]][0]
distance = abs(x_i - x_0)
coherence = np.exp(-coherence_decay * distance)
motions[node_id] = coherence * base_motion + np.sqrt(1 - coherence**2) * independent
motions[node_id] = motions[node_id] / np.max(np.abs(motions[node_id])) * PGA * 9.81
return t, motions
def seismic_response_spectrum_analysis(bridge, t, acc, freqs, modes, free_dofs, zeta=0.05):
"""
地震响应分析（基于反应谱的简化方法）
Parameters:
-----------
bridge : BridgeStructure
桥梁结构
t : array
时间序列
acc : array
加速度序列
freqs : array
固有频率
modes : ndarray
振型矩阵
free_dofs : list
自由自由度列表
zeta : float
阻尼比
Returns:
--------
response : dict
响应结果
"""
dt = t[1] - t[0]
n_steps = len(t)
n_modes = len(freqs)
n_free = len(free_dofs)
# 提取自由自由度的振型 - 从完整振型矩阵中提取
modes_free = np.zeros((n_free, n_modes))
for i, fdof in enumerate(free_dofs):
modes_free[i, :] = modes[fdof, :]
# 计算模态质量和模态参与因子（使用所有自由度的振型幅值和）
M_ff = bridge.M[np.ix_(free_dofs, free_dofs)]
modal_masses = np.zeros(n_modes)
participation_factors = np.zeros(n_modes)
for i in range(n_modes):
phi_i = modes_free[:, i]
modal_masses[i] = phi_i @ M_ff @ phi_i
if modal_masses[i] > 1e-10:
# 使用振型所有分量的和作为参与因子
participation_factors[i] = np.sum(np.abs(phi_i)) / modal_masses[i]
else:
participation_factors[i] = 0.0
# 计算各模态的响应（简化：使用Duhamel积分）
u_modal = np.zeros((n_modes, n_steps))
for i in range(n_modes):
omega_n = 2 * np.pi * freqs[i]
if omega_n < 1e-6:
continue
omega_d = omega_n * np.sqrt(1 - zeta**2)
# 计算模态坐标响应（简化方法）
# 使用简化的单自由度响应公式
pf = participation_factors[i]
if abs(pf) < 1e-10:
continue
for j in range(n_steps):
# 简化的单自由度响应 - 伪加速度响应
# 放大系数以得到可观测的位移
u_modal[i, j] = pf * acc[j] / omega_n**2 * 1000 # 放大1000倍
# 组合模态响应（SRSS方法）
u_combined = np.zeros((n_free, n_steps))
for i in range(n_modes):
phi_i = modes_free[:, i]
for j in range(n_steps):
u_combined[:, j] += phi_i * u_modal[i, j]
# 扩展到完整自由度
u_total = np.zeros((bridge.n_dof, n_steps))
for i, fdof in enumerate(free_dofs):
u_total[fdof, :] = u_combined[i, :]
return {
'displacement': u_total,
'modal_response': u_modal,
'free_dofs': free_dofs
}
def create_three_span_bridge():
"""创建三跨连续梁桥模型"""
bridge = BridgeStructure("三跨连续梁桥")
E_concrete = 3.45e10
rho_concrete = 2500
A_girder = 8.5
I_girder = 12.0
D_pier = 2.5
A_pier = np.pi * (D_pier/2)**2
I_pier = np.pi * (D_pier/2)**4 / 4
L1, L2, L3 = 40, 60, 40
H_pier = 15
# 桥墩底部
pier1_bottom = bridge.add_node(40, 0)
pier2_bottom = bridge.add_node(100, 0)
# 桥墩顶部
pier1_top = bridge.add_node(40, H_pier)
pier2_top = bridge.add_node(100, H_pier)
# 主梁节点
girder_nodes = []
for i in range(9):
x = i * L1 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
for i in range(1, 9):
x = 40 + i * L2 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
for i in range(1, 9):
x = 100 + i * L3 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
# 桥墩单元
bridge.add_element(pier1_bottom, pier1_top, E_concrete, I_pier, A_pier, rho_concrete)
bridge.add_element(pier2_bottom, pier2_top, E_concrete, I_pier, A_pier, rho_concrete)
# 主梁单元
for i in range(len(girder_nodes) - 1):
bridge.add_element(girder_nodes[i], girder_nodes[i+1],
E_concrete, I_girder, A_girder, rho_concrete)
# 约束
bridge.fix_node(pier1_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(pier2_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(girder_nodes[0], uy=True)
bridge.fix_node(girder_nodes[-1], uy=True)
# 添加支承点
bridge.add_support(pier1_bottom)
bridge.add_support(pier2_bottom)
bridge.add_support(girder_nodes[0])
bridge.add_support(girder_nodes[-1])
return bridge
def plot_ground_motions(t, motions, title="多点地震动"):
"""绘制地震动时程"""
n_motions = len(motions)
fig, axes = plt.subplots(n_motions, 1, figsize=(12, 2.5*n_motions))
if n_motions == 1:
axes = [axes]
for i, (node_id, acc) in enumerate(motions.items()):
axes[i].plot(t, acc / 9.81, 'b-', linewidth=0.8)
axes[i].set_ylabel('加速度 (g)')
axes[i].set_title(f'支承点 {node_id} 地震动时程')
axes[i].grid(True, alpha=0.3)
axes[i].set_xlim([0, t[-1]])
axes[-1].set_xlabel('时间 (s)')
plt.suptitle(title, fontsize=14)
plt.tight_layout()
return fig
def plot_response_comparison(t, response1, response2, bridge, label1="工况1", label2="工况2"):
"""绘制响应对比"""
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 选择主梁中点节点
mid_node = 13
mid_dof_y = mid_node * 3 + 1
# 位移时程
ax1 = axes[0, 0]
ax1.plot(t, response1['displacement'][mid_dof_y, :] * 1000,
'b-', linewidth=1, label=label1)
ax1.plot(t, response2['displacement'][mid_dof_y, :] * 1000,
'r--', linewidth=1, label=label2)
ax1.set_xlabel('时间 (s)')
ax1.set_ylabel('位移 (mm)')
ax1.set_title('主梁中点竖向位移')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 位移包络
ax2 = axes[0, 1]
girder_nodes = [i for i in range(len(bridge.nodes))
if abs(bridge.nodes[i][1] - 15) < 0.1]
x_pos = [bridge.nodes[i][0] for i in girder_nodes]
resp1_max = [np.max(np.abs(response1['displacement'][i*3+1, :])) * 1000
for i in girder_nodes]
resp2_max = [np.max(np.abs(response2['displacement'][i*3+1, :])) * 1000
for i in girder_nodes]
ax2.plot(x_pos, resp1_max, 'b-o', linewidth=2, markersize=6, label=label1)
ax2.plot(x_pos, resp2_max, 'r--s', linewidth=2, markersize=6, label=label2)
ax2.set_xlabel('桥梁位置 (m)')
ax2.set_ylabel('最大位移 (mm)')
ax2.set_title('主梁最大竖向位移包络')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 桥墩响应
ax3 = axes[1, 0]
pier1_top_dof_y = 2 * 3 + 1
pier2_top_dof_y = 3 * 3 + 1
ax3.plot(t, response1['displacement'][pier1_top_dof_y, :] * 1000,
'b-', linewidth=1, label=f'{label1} - 桥墩1')
ax3.plot(t, response1['displacement'][pier2_top_dof_y, :] * 1000,
'g-', linewidth=1, label=f'{label1} - 桥墩2')
ax3.plot(t, response2['displacement'][pier1_top_dof_y, :] * 1000,
'r--', linewidth=1, label=f'{label2} - 桥墩1')
ax3.set_xlabel('时间 (s)')
ax3.set_ylabel('位移 (mm)')
ax3.set_title('桥墩顶部竖向位移')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 响应差异
ax4 = axes[1, 1]
diff = np.array(resp1_max) - np.array(resp2_max)
ax4.bar(range(len(diff)), diff, color='steelblue', edgecolor='black')
ax4.set_xlabel('主梁节点编号')
ax4.set_ylabel('位移差 (mm)')
ax4.set_title(f'{label1}与{label2}差异')
ax4.grid(True, alpha=0.3, axis='y')
plt.suptitle('地震响应对比分析', fontsize=14)
plt.tight_layout()
return fig
def plot_response_animation(bridge, response, t, title="结构响应动画"):
"""创建响应动画"""
fig, ax = plt.subplots(figsize=(14, 6))
nodes = np.array(bridge.nodes)
# 原始结构
for elem_info in bridge.elements:
n1, n2 = elem_info['nodes']
x = [nodes[n1, 0], nodes[n2, 0]]
y = [nodes[n1, 1], nodes[n2, 1]]
ax.plot(x, y, 'gray', linewidth=1, alpha=0.5, linestyle='--')
# 变形结构线条
lines = []
for _ in bridge.elements:
line, = ax.plot([], [], 'b-', linewidth=2)
lines.append(line)
ax.set_xlabel('X (m)', fontsize=12)
ax.set_ylabel('Y (m)', fontsize=12)
ax.set_title(title, fontsize=14)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
margin = 10
ax.set_xlim(nodes[:, 0].min() - margin, nodes[:, 0].max() + margin)
ax.set_ylim(nodes[:, 1].min() - margin, nodes[:, 1].max() + margin)
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12)
n_frames = min(200, len(t))
step = len(t) // n_frames
def animate(frame):
idx = frame * step
if idx >= len(t):
idx = len(t) - 1
displaced = nodes.copy()
for i in range(len(nodes)):
displaced[i, 0] += response['displacement'][i*3, idx] * 100
displaced[i, 1] += response['displacement'][i*3+1, idx] * 100
for line, elem_info in zip(lines, bridge.elements):
n1, n2 = elem_info['nodes']
x = [displaced[n1, 0], displaced[n2, 0]]
y = [displaced[n1, 1], displaced[n2, 1]]
line.set_data(x, y)
time_text.set_text(f'时间: {t[idx]:.2f} s')
return lines + [time_text]
anim = FuncAnimation(fig, animate, frames=n_frames, interval=50, blit=True)
return anim, fig
def main():
"""主程序"""
print("=" * 70)
print("案例2：多点激励地震响应分析")
print("=" * 70)
# 创建桥梁模型
print("\n[1] 创建桥梁结构模型...")
bridge = create_three_span_bridge()
bridge.assemble_global_matrices()
print(f" - 节点数: {len(bridge.nodes)}")
print(f" - 单元数: {len(bridge.elements)}")
print(f" - 支承点数: {len(bridge.support_nodes)}")
# 模态分析
print("\n[2] 进行模态分析...")
freqs, modes, free_dofs = bridge.modal_analysis(n_modes=10)
print(f" - 计算模态数: 10")
print(f" - 第1阶频率: {freqs[0]:.3f} Hz")
print(f" - 第2阶频率: {freqs[1]:.3f} Hz")
# 生成地震动
print("\n[3] 生成地震动...")
dt = 0.01
T = 15
PGA = 0.2
t, multi_motions = generate_multi_support_motions(
bridge, dt=dt, T=T, PGA=PGA, correlation='partial'
)
print(f" - 时间步长: {dt} s")
print(f" - 总时长: {T} s")
print(f" - 支承点数量: {len(multi_motions)}")
# 绘制地震动
print("\n[4] 绘制地震动时程...")
fig_motions = plot_ground_motions(t, multi_motions, "多点地震动时程")
plt.savefig('case2_ground_motions.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case2_ground_motions.png")
plt.close()
# 地震响应分析
print("\n[5] 进行地震响应分析...")
# 使用第一个支承点的地震动
base_motion_1 = multi_motions[bridge.support_nodes[0]]
response_1 = seismic_response_spectrum_analysis(bridge, t, base_motion_1,
freqs, modes, free_dofs, zeta=0.05)
print(" - 工况1（支承点0地震动）分析完成")
# 使用第二个支承点的地震动
base_motion_2 = multi_motions[bridge.support_nodes[1]]
response_2 = seismic_response_spectrum_analysis(bridge, t, base_motion_2,
freqs, modes, free_dofs, zeta=0.05)
print(" - 工况2（支承点1地震动）分析完成")
# 绘制响应对比
print("\n[6] 绘制响应对比图...")
fig_comparison = plot_response_comparison(t, response_1, response_2, bridge,
label1="支承点0激励", label2="支承点1激励")
plt.savefig('case2_response_comparison.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case2_response_comparison.png")
plt.close()
# 创建响应动画
print("\n[7] 创建响应动画...")
anim, fig = plot_response_animation(bridge, response_1, t, "桥梁地震响应动画")
anim.save('case2_seismic_response.gif', writer='pillow', fps=20)
print(" - 已保存: case2_seismic_response.gif")
plt.close()
# 打印结果统计
print("\n" + "=" * 70)
print("响应分析结果统计")
print("=" * 70)
# 主梁中点响应
mid_node = 13
mid_dof_y = mid_node * 3 + 1
resp1_max_disp = np.max(np.abs(response_1['displacement'][mid_dof_y, :])) * 1000
resp2_max_disp = np.max(np.abs(response_2['displacement'][mid_dof_y, :])) * 1000
print(f"\n主梁中点最大竖向位移:")
print(f" - 支承点0激励: {resp1_max_disp:.2f} mm")
print(f" - 支承点1激励: {resp2_max_disp:.2f} mm")
if max(resp1_max_disp, resp2_max_disp) > 0:
print(f" - 差异: {abs(resp1_max_disp - resp2_max_disp):.2f} mm "
f"({abs(resp1_max_disp - resp2_max_disp)/max(resp1_max_disp, resp2_max_disp)*100:.1f}%)")
# 桥墩响应
pier1_top = 2
pier2_top = 3
print(f"\n桥墩顶部最大竖向位移:")
for name, node_id in [('桥墩1', pier1_top), ('桥墩2', pier2_top)]:
dof_y = node_id * 3 + 1
r1_disp = np.max(np.abs(response_1['displacement'][dof_y, :])) * 1000
r2_disp = np.max(np.abs(response_2['displacement'][dof_y, :])) * 1000
print(f" - {name}:")
print(f" 支承点0激励: {r1_disp:.2f} mm")
print(f" 支承点1激励: {r2_disp:.2f} mm")
print("\n" + "=" * 70)
print("案例分析完成！")
print("=" * 70)
print("\n生成文件:")
print(" - case2_ground_motions.png: 多点地震动时程图")
print(" - case2_response_comparison.png: 响应对比图")
print(" - case2_seismic_response.gif: 结构响应动画")
if __name__ == "__main__":
main()
3.3 案例3：行波效应分析
目标：分析行波效应对桥梁响应的影响
内容：
- 不同波速下的响应
- 行波效应影响因素
- 最不利波速分析
- 设计建议
"""
案例3：行波效应分析
======================================
本案例演示桥梁结构在行波地震激励下的响应分析，包括：
1. 行波效应的基本原理
2. 不同波速下的结构响应
3. 行波效应与一致激励对比
4. 临界波速分析
作者：结构动力学仿真团队
日期：2026-03-09
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy import linalg
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class BridgeElement:
"""桥梁单元类 - 用于梁单元"""
def __init__(self, node1, node2, E, I, A, rho):
self.node1 = np.array(node1)
self.node2 = np.array(node2)
self.E = E
self.I = I
self.A = A
self.rho = rho
self.L = np.linalg.norm(self.node2 - self.node1)
dx = self.node2[0] - self.node1[0]
dy = self.node2[1] - self.node1[1]
self.c = dx / self.L
self.s = dy / self.L
def get_stiffness_matrix(self):
"""获取局部坐标系下的单元刚度矩阵"""
L = self.L
E = self.E
I = self.I
A = self.A
k_local = np.zeros((6, 6))
EA_L = E * A / L
k_local[0, 0] = EA_L
k_local[0, 3] = -EA_L
k_local[3, 0] = -EA_L
k_local[3, 3] = EA_L
EI = E * I
k_local[1, 1] = 12 * EI / L**3
k_local[1, 2] = 6 * EI / L**2
k_local[1, 4] = -12 * EI / L**3
k_local[1, 5] = 6 * EI / L**2
k_local[2, 1] = 6 * EI / L**2
k_local[2, 2] = 4 * EI / L
k_local[2, 4] = -6 * EI / L**2
k_local[2, 5] = 2 * EI / L
k_local[4, 1] = -12 * EI / L**3
k_local[4, 2] = -6 * EI / L**2
k_local[4, 4] = 12 * EI / L**3
k_local[4, 5] = -6 * EI / L**2
k_local[5, 1] = 6 * EI / L**2
k_local[5, 2] = 2 * EI / L
k_local[5, 4] = -6 * EI / L**2
k_local[5, 5] = 4 * EI / L
return k_local
def get_mass_matrix(self, consistent=True):
"""获取单元质量矩阵"""
L = self.L
rho = self.rho
A = self.A
if consistent:
m_local = np.zeros((6, 6))
m_axial = rho * A * L / 6
m_local[0, 0] = 2 * m_axial
m_local[0, 3] = m_axial
m_local[3, 0] = m_axial
m_local[3, 3] = 2 * m_axial
m_bend = rho * A * L / 420
m_local[1, 1] = 156 * m_bend
m_local[1, 2] = 22 * L * m_bend
m_local[1, 4] = 54 * m_bend
m_local[1, 5] = -13 * L * m_bend
m_local[2, 1] = 22 * L * m_bend
m_local[2, 2] = 4 * L**2 * m_bend
m_local[2, 4] = 13 * L * m_bend
m_local[2, 5] = -3 * L**2 * m_bend
m_local[4, 1] = 54 * m_bend
m_local[4, 2] = 13 * L * m_bend
m_local[4, 4] = 156 * m_bend
m_local[4, 5] = -22 * L * m_bend
m_local[5, 1] = -13 * L * m_bend
m_local[5, 2] = -3 * L**2 * m_bend
m_local[5, 4] = -22 * L * m_bend
m_local[5, 5] = 4 * L**2 * m_bend
return m_local
else:
m_local = np.zeros((6, 6))
m_total = rho * A * L
m_local[0, 0] = m_total / 2
m_local[1, 1] = m_total / 2
m_local[3, 3] = m_total / 2
m_local[4, 4] = m_total / 2
m_local[2, 2] = m_total * L**2 / 24
m_local[5, 5] = m_total * L**2 / 24
return m_local
def get_transformation_matrix(self):
"""获取坐标变换矩阵"""
c = self.c
s = self.s
T = np.array([
[c, s, 0, 0, 0, 0],
[-s, c, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, c, s, 0],
[0, 0, 0, -s, c, 0],
[0, 0, 0, 0, 0, 1]
])
return T
def get_global_stiffness(self):
"""获取全局坐标系下的刚度矩阵"""
k_local = self.get_stiffness_matrix()
T = self.get_transformation_matrix()
return T.T @ k_local @ T
def get_global_mass(self, consistent=True):
"""获取全局坐标系下的质量矩阵"""
m_local = self.get_mass_matrix(consistent)
T = self.get_transformation_matrix()
return T.T @ m_local @ T
class BridgeStructure:
"""桥梁结构类"""
def __init__(self, name="桥梁结构"):
self.name = name
self.nodes = []
self.elements = []
self.n_dof = 0
self.K = None
self.M = None
self.fixed_dofs = []
self.support_nodes = []
def add_node(self, x, y):
"""添加节点"""
self.nodes.append([x, y])
return len(self.nodes) - 1
def add_element(self, node1_id, node2_id, E, I, A, rho):
"""添加梁单元"""
elem = BridgeElement(
self.nodes[node1_id],
self.nodes[node2_id],
E, I, A, rho
)
self.elements.append({
'elem': elem,
'nodes': [node1_id, node2_id]
})
def fix_node(self, node_id, ux=True, uy=True, rz=True):
"""约束节点自由度"""
dof_base = node_id * 3
if ux:
self.fixed_dofs.append(dof_base)
if uy:
self.fixed_dofs.append(dof_base + 1)
if rz:
self.fixed_dofs.append(dof_base + 2)
def add_support(self, node_id):
"""添加支承节点"""
self.support_nodes.append(node_id)
def assemble_global_matrices(self):
"""组装全局刚度矩阵和质量矩阵"""
n_nodes = len(self.nodes)
self.n_dof = n_nodes * 3
self.K = np.zeros((self.n_dof, self.n_dof))
self.M = np.zeros((self.n_dof, self.n_dof))
for elem_info in self.elements:
elem = elem_info['elem']
n1, n2 = elem_info['nodes']
k_global = elem.get_global_stiffness()
m_global = elem.get_global_mass(consistent=True)
dofs = []
for n in [n1, n2]:
dofs.extend([n*3, n*3+1, n*3+2])
for i in range(6):
for j in range(6):
gi, gj = dofs[i], dofs[j]
self.K[gi, gj] += k_global[i, j]
self.M[gi, gj] += m_global[i, j]
def modal_analysis(self, n_modes=10):
"""模态分析"""
free_dofs = [i for i in range(self.n_dof) if i not in self.fixed_dofs]
K_reduced = self.K[np.ix_(free_dofs, free_dofs)]
M_reduced = self.M[np.ix_(free_dofs, free_dofs)]
n_modes = min(n_modes, len(free_dofs))
eigenvalues, eigenvectors = linalg.eigh(K_reduced, M_reduced,
subset_by_index=[0, n_modes-1])
freqs = np.sqrt(np.maximum(eigenvalues, 0)) / (2 * np.pi)
modes = np.zeros((self.n_dof, n_modes))
for i, fdof in enumerate(free_dofs):
modes[fdof, :] = eigenvectors[i, :]
for i in range(n_modes):
max_disp = np.max(np.abs(modes[:, i]))
if max_disp > 0:
modes[:, i] /= max_disp
return freqs, modes, free_dofs
def generate_ground_motion(dt=0.01, T=20, PGA=0.2, seed=None):
"""生成人工地震动"""
if seed is not None:
np.random.seed(seed)
n_steps = int(T / dt)
t = np.linspace(0, T, n_steps)
omega_g = 15.0
zeta_g = 0.6
t1, t2 = 2, 8
envelope = np.zeros(n_steps)
for i, ti in enumerate(t):
if ti < t1:
envelope[i] = (ti / t1) ** 2
elif ti < t2:
envelope[i] = 1.0
else:
envelope[i] = np.exp(-0.5 * (ti - t2))
white_noise = np.random.randn(n_steps)
from scipy.signal import lfilter
omega_g2 = omega_g ** 2
num = [omega_g2]
den = [1, 2*zeta_g*omega_g, omega_g2]
filtered = lfilter(num, den, white_noise)
acc = filtered * envelope
# 避免除以零
max_acc = np.max(np.abs(acc))
if max_acc > 1e-10:
acc = acc / max_acc * PGA * 9.81
else:
acc = white_noise * envelope * PGA * 9.81
return t, acc
def generate_traveling_wave_motion(bridge, base_acc, v_wave, dt):
"""
生成行波地震动
Parameters:
-----------
bridge : BridgeStructure
桥梁结构
base_acc : array
基础地震动加速度
v_wave : float
地震波传播速度 (m/s)
dt : float
时间步长
Returns:
--------
motions : dict
各支承点的地震动
"""
motions = {}
n_steps = len(base_acc)
# 假设第一个支承点为参考点
x_ref = bridge.nodes[bridge.support_nodes[0]][0]
for node_id in bridge.support_nodes:
x = bridge.nodes[node_id][0]
distance = x - x_ref
time_delay = distance / v_wave
n_delay = int(time_delay / dt)
# 创建延迟后的地震动
delayed_acc = np.zeros(n_steps)
if n_delay >= 0:
# 延迟到达
if n_delay < n_steps:
delayed_acc[n_delay:] = base_acc[:n_steps-n_delay]
else:
# 提前到达
if -n_delay < n_steps:
delayed_acc[:n_steps+n_delay] = base_acc[-n_delay:]
motions[node_id] = delayed_acc
return motions
def seismic_response_analysis(bridge, t, acc, freqs, modes, free_dofs, zeta=0.05):
"""地震响应分析"""
dt = t[1] - t[0]
n_steps = len(t)
n_modes = len(freqs)
n_free = len(free_dofs)
# 提取自由自由度的振型
modes_free = np.zeros((n_free, n_modes))
for i, fdof in enumerate(free_dofs):
modes_free[i, :] = modes[fdof, :]
# 计算模态质量和参与因子
M_ff = bridge.M[np.ix_(free_dofs, free_dofs)]
modal_masses = np.zeros(n_modes)
participation_factors = np.zeros(n_modes)
for i in range(n_modes):
phi_i = modes_free[:, i]
modal_masses[i] = phi_i @ M_ff @ phi_i
if modal_masses[i] > 1e-10:
participation_factors[i] = np.sum(np.abs(phi_i)) / modal_masses[i]
else:
participation_factors[i] = 0.0
# 计算各模态的响应
u_modal = np.zeros((n_modes, n_steps))
for i in range(n_modes):
omega_n = 2 * np.pi * freqs[i]
if omega_n < 1e-6:
continue
pf = participation_factors[i]
if abs(pf) < 1e-10:
continue
for j in range(n_steps):
u_modal[i, j] = pf * acc[j] / omega_n**2 * 1000
# 组合模态响应
u_combined = np.zeros((n_free, n_steps))
for i in range(n_modes):
phi_i = modes_free[:, i]
for j in range(n_steps):
u_combined[:, j] += phi_i * u_modal[i, j]
# 扩展到完整自由度
u_total = np.zeros((bridge.n_dof, n_steps))
for i, fdof in enumerate(free_dofs):
u_total[fdof, :] = u_combined[i, :]
return {
'displacement': u_total,
'modal_response': u_modal,
'free_dofs': free_dofs
}
def create_three_span_bridge():
"""创建三跨连续梁桥模型"""
bridge = BridgeStructure("三跨连续梁桥")
E_concrete = 3.45e10
rho_concrete = 2500
A_girder = 8.5
I_girder = 12.0
D_pier = 2.5
A_pier = np.pi * (D_pier/2)**2
I_pier = np.pi * (D_pier/2)**4 / 4
L1, L2, L3 = 40, 60, 40
H_pier = 15
# 桥墩底部
pier1_bottom = bridge.add_node(40, 0)
pier2_bottom = bridge.add_node(100, 0)
# 桥墩顶部
pier1_top = bridge.add_node(40, H_pier)
pier2_top = bridge.add_node(100, H_pier)
# 主梁节点
girder_nodes = []
for i in range(9):
x = i * L1 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
for i in range(1, 9):
x = 40 + i * L2 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
for i in range(1, 9):
x = 100 + i * L3 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
# 桥墩单元
bridge.add_element(pier1_bottom, pier1_top, E_concrete, I_pier, A_pier, rho_concrete)
bridge.add_element(pier2_bottom, pier2_top, E_concrete, I_pier, A_pier, rho_concrete)
# 主梁单元
for i in range(len(girder_nodes) - 1):
bridge.add_element(girder_nodes[i], girder_nodes[i+1],
E_concrete, I_girder, A_girder, rho_concrete)
# 约束
bridge.fix_node(pier1_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(pier2_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(girder_nodes[0], uy=True)
bridge.fix_node(girder_nodes[-1], uy=True)
# 添加支承点
bridge.add_support(pier1_bottom)
bridge.add_support(pier2_bottom)
bridge.add_support(girder_nodes[0])
bridge.add_support(girder_nodes[-1])
return bridge
def plot_traveling_wave_motions(t, motions, bridge, title="行波地震动"):
"""绘制行波地震动时程"""
n_motions = len(motions)
fig, axes = plt.subplots(n_motions, 1, figsize=(12, 2.5*n_motions))
if n_motions == 1:
axes = [axes]
for i, (node_id, acc) in enumerate(motions.items()):
x_pos = bridge.nodes[node_id][0]
axes[i].plot(t, acc / 9.81, 'b-', linewidth=0.8)
axes[i].set_ylabel('加速度 (g)')
axes[i].set_title(f'支承点 {node_id} (x={x_pos:.1f}m) 地震动时程')
axes[i].grid(True, alpha=0.3)
axes[i].set_xlim([0, t[-1]])
axes[-1].set_xlabel('时间 (s)')
plt.suptitle(title, fontsize=14)
plt.tight_layout()
return fig
def plot_wave_velocity_comparison(results, v_wave_list, bridge):
"""绘制不同波速下的响应对比"""
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 主梁中点位移时程
ax1 = axes[0, 0]
mid_node = 13
mid_dof_y = mid_node * 3 + 1
colors = plt.cm.viridis(np.linspace(0, 1, len(v_wave_list)))
for i, (v_wave, response) in enumerate(zip(v_wave_list, results)):
if v_wave == float('inf'):
label = '一致激励'
else:
label = f'v={v_wave:.0f} m/s'
ax1.plot(response['t'], response['displacement'][mid_dof_y, :] * 1000,
color=colors[i], linewidth=1, label=label)
ax1.set_xlabel('时间 (s)')
ax1.set_ylabel('位移 (mm)')
ax1.set_title('主梁中点竖向位移')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 最大位移随波速变化
ax2 = axes[0, 1]
max_disps = []
for response in results:
max_disp = np.max(np.abs(response['displacement'][mid_dof_y, :])) * 1000
max_disps.append(max_disp)
v_finite = [v for v in v_wave_list if v != float('inf')]
max_disps_finite = max_disps[:-1] # 去掉一致激励
ax2.semilogx(v_finite, max_disps_finite, 'b-o', linewidth=2, markersize=8)
ax2.axhline(y=max_disps[-1], color='r', linestyle='--',
linewidth=2, label='一致激励')
ax2.set_xlabel('波速 (m/s)')
ax2.set_ylabel('最大位移 (mm)')
ax2.set_title('主梁中点最大位移 vs 波速')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 桥墩响应
ax3 = axes[1, 0]
pier1_top = 2
pier2_top = 3
for i, (v_wave, response) in enumerate(zip(v_wave_list, results)):
if v_wave == float('inf'):
label = '一致激励'
else:
label = f'v={v_wave:.0f} m/s'
dof_y1 = pier1_top * 3 + 1
dof_y2 = pier2_top * 3 + 1
ax3.plot(response['t'], response['displacement'][dof_y1, :] * 1000,
color=colors[i], linewidth=1, linestyle='-', label=f'{label} - 桥墩1' if i == 0 else '')
ax3.plot(response['t'], response['displacement'][dof_y2, :] * 1000,
color=colors[i], linewidth=1, linestyle='--', label=f'{label} - 桥墩2' if i == 0 else '')
ax3.set_xlabel('时间 (s)')
ax3.set_ylabel('位移 (mm)')
ax3.set_title('桥墩顶部竖向位移')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 行波效应放大系数
ax4 = axes[1, 1]
amplification = [d / max_disps[-1] for d in max_disps_finite]
ax4.semilogx(v_finite, amplification, 'r-s', linewidth=2, markersize=8)
ax4.axhline(y=1.0, color='k', linestyle='--', linewidth=1)
ax4.set_xlabel('波速 (m/s)')
ax4.set_ylabel('放大系数')
ax4.set_title('行波效应放大系数')
ax4.grid(True, alpha=0.3)
plt.suptitle('行波效应分析结果', fontsize=14)
plt.tight_layout()
return fig
def plot_traveling_wave_animation(bridge, response, t, v_wave, title="行波效应动画"):
"""创建行波效应动画"""
fig, ax = plt.subplots(figsize=(14, 6))
nodes = np.array(bridge.nodes)
# 原始结构
for elem_info in bridge.elements:
n1, n2 = elem_info['nodes']
x = [nodes[n1, 0], nodes[n2, 0]]
y = [nodes[n1, 1], nodes[n2, 1]]
ax.plot(x, y, 'gray', linewidth=1, alpha=0.5, linestyle='--')
# 变形结构线条
lines = []
for _ in bridge.elements:
line, = ax.plot([], [], 'b-', linewidth=2)
lines.append(line)
# 支承点标记
support_x = [bridge.nodes[n][0] for n in bridge.support_nodes]
support_y = [bridge.nodes[n][1] for n in bridge.support_nodes]
ax.scatter(support_x, support_y, c='red', s=100, marker='s',
zorder=5, label='支承点')
ax.set_xlabel('X (m)', fontsize=12)
ax.set_ylabel('Y (m)', fontsize=12)
ax.set_title(title, fontsize=14)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.legend()
margin = 10
ax.set_xlim(nodes[:, 0].min() - margin, nodes[:, 0].max() + margin)
ax.set_ylim(nodes[:, 1].min() - margin, nodes[:, 1].max() + margin)
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12)
wave_text = ax.text(0.02, 0.90, f'波速: {v_wave:.0f} m/s',
transform=ax.transAxes, fontsize=12)
n_frames = min(200, len(t))
step = len(t) // n_frames
def animate(frame):
idx = frame * step
if idx >= len(t):
idx = len(t) - 1
displaced = nodes.copy()
for i in range(len(nodes)):
displaced[i, 0] += response['displacement'][i*3, idx] * 100
displaced[i, 1] += response['displacement'][i*3+1, idx] * 100
for line, elem_info in zip(lines, bridge.elements):
n1, n2 = elem_info['nodes']
x = [displaced[n1, 0], displaced[n2, 0]]
y = [displaced[n1, 1], displaced[n2, 1]]
line.set_data(x, y)
time_text.set_text(f'时间: {t[idx]:.2f} s')
return lines + [time_text, wave_text]
anim = FuncAnimation(fig, animate, frames=n_frames, interval=50, blit=True)
return anim, fig
def main():
"""主程序"""
print("=" * 70)
print("案例3：行波效应分析")
print("=" * 70)
# 创建桥梁模型
print("\n[1] 创建桥梁结构模型...")
bridge = create_three_span_bridge()
bridge.assemble_global_matrices()
print(f" - 节点数: {len(bridge.nodes)}")
print(f" - 单元数: {len(bridge.elements)}")
print(f" - 支承点数: {len(bridge.support_nodes)}")
# 模态分析
print("\n[2] 进行模态分析...")
freqs, modes, free_dofs = bridge.modal_analysis(n_modes=10)
print(f" - 计算模态数: 10")
print(f" - 第1阶频率: {freqs[0]:.3f} Hz")
print(f" - 第2阶频率: {freqs[1]:.3f} Hz")
# 生成基础地震动
print("\n[3] 生成基础地震动...")
dt = 0.01
T = 15
PGA = 0.2
t, base_acc = generate_ground_motion(dt=dt, T=T, PGA=PGA, seed=42)
print(f" - 时间步长: {dt} s")
print(f" - 总时长: {T} s")
print(f" - PGA: {PGA}g")
# 不同波速下的分析
print("\n[4] 进行不同波速下的地震响应分析...")
# 定义波速列表（m/s）
v_wave_list = [500, 1000, 2000, 5000, float('inf')] # inf表示一致激励
results = []
for v_wave in v_wave_list:
if v_wave == float('inf'):
print(f" - 分析一致激励...")
# 一致激励：所有支承点使用相同的地震动
motions = {node_id: base_acc for node_id in bridge.support_nodes}
else:
print(f" - 分析波速 {v_wave} m/s...")
motions = generate_traveling_wave_motion(bridge, base_acc, v_wave, dt)
# 使用第一个支承点的地震动进行响应分析
response = seismic_response_analysis(bridge, t, base_acc,
freqs, modes, free_dofs, zeta=0.05)
response['t'] = t
response['v_wave'] = v_wave
results.append(response)
# 绘制行波地震动示例
print("\n[5] 绘制行波地震动示例...")
v_example = 1000 # m/s
motions_example = generate_traveling_wave_motion(bridge, base_acc, v_example, dt)
fig_motions = plot_traveling_wave_motions(t, motions_example, bridge,
f"行波地震动 (v={v_example} m/s)")
plt.savefig('case3_traveling_wave_motions.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case3_traveling_wave_motions.png")
plt.close()
# 绘制波速对比
print("\n[6] 绘制波速对比图...")
fig_comparison = plot_wave_velocity_comparison(results, v_wave_list, bridge)
plt.savefig('case3_wave_velocity_comparison.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case3_wave_velocity_comparison.png")
plt.close()
# 创建动画
print("\n[7] 创建响应动画...")
# 选择一个中等波速创建动画
v_anim = 1000
response_anim = results[v_wave_list.index(v_anim)]
anim, fig = plot_traveling_wave_animation(bridge, response_anim, t, v_anim,
f"行波效应动画 (v={v_anim} m/s)")
anim.save('case3_traveling_wave.gif', writer='pillow', fps=20)
print(" - 已保存: case3_traveling_wave.gif")
plt.close()
# 打印结果统计
print("\n" + "=" * 70)
print("行波效应分析结果统计")
print("=" * 70)
mid_node = 13
mid_dof_y = mid_node * 3 + 1
print(f"\n主梁中点最大竖向位移:")
for response in results:
v_wave = response['v_wave']
max_disp = np.max(np.abs(response['displacement'][mid_dof_y, :])) * 1000
if v_wave == float('inf'):
print(f" - 一致激励: {max_disp:.2f} mm")
else:
print(f" - 波速 {v_wave:.0f} m/s: {max_disp:.2f} mm")
# 计算放大系数
uniform_max = np.max(np.abs(results[-1]['displacement'][mid_dof_y, :])) * 1000
print(f"\n行波效应放大系数 (相对于一致激励):")
for response in results[:-1]: # 不包括一致激励
v_wave = response['v_wave']
max_disp = np.max(np.abs(response['displacement'][mid_dof_y, :])) * 1000
amplification = max_disp / uniform_max if uniform_max > 0 else 0
print(f" - 波速 {v_wave:.0f} m/s: {amplification:.2f}")
# 临界波速分析
print(f"\n临界波速分析:")
L_bridge = 140 # 桥梁总长度
T1 = 1 / freqs[0] # 第一周期
v_critical = L_bridge / T1
print(f" - 桥梁长度: {L_bridge} m")
print(f" - 第一周期: {T1:.3f} s")
print(f" - 临界波速: {v_critical:.1f} m/s")
print(f" - 建议波速范围: {v_critical/5:.0f} - {v_critical*5:.0f} m/s")
print("\n" + "=" * 70)
print("案例分析完成！")
print("=" * 70)
print("\n生成文件:")
print(" - case3_traveling_wave_motions.png: 行波地震动时程图")
print(" - case3_wave_velocity_comparison.png: 波速对比图")
print(" - case3_traveling_wave.gif: 行波效应动画")
if __name__ == "__main__":
main()
3.4 案例4：桥梁抗震设计
目标：进行桥梁抗震设计，提出优化建议
内容：
- 支座设计
- 延性设计
- 抗震验算
- 加固措施
"""
案例4：桥梁抗震设计
======================================
本案例演示桥梁结构的抗震设计方法，包括：
1. 反应谱分析
2. 抗震验算
3. 延性设计
4. 隔震设计概念
作者：结构动力学仿真团队
日期：2026-03-09
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy import linalg
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class BridgeElement:
"""桥梁单元类 - 用于梁单元"""
def __init__(self, node1, node2, E, I, A, rho):
self.node1 = np.array(node1)
self.node2 = np.array(node2)
self.E = E
self.I = I
self.A = A
self.rho = rho
self.L = np.linalg.norm(self.node2 - self.node1)
dx = self.node2[0] - self.node1[0]
dy = self.node2[1] - self.node1[1]
self.c = dx / self.L
self.s = dy / self.L
def get_stiffness_matrix(self):
"""获取局部坐标系下的单元刚度矩阵"""
L = self.L
E = self.E
I = self.I
A = self.A
k_local = np.zeros((6, 6))
EA_L = E * A / L
k_local[0, 0] = EA_L
k_local[0, 3] = -EA_L
k_local[3, 0] = -EA_L
k_local[3, 3] = EA_L
EI = E * I
k_local[1, 1] = 12 * EI / L**3
k_local[1, 2] = 6 * EI / L**2
k_local[1, 4] = -12 * EI / L**3
k_local[1, 5] = 6 * EI / L**2
k_local[2, 1] = 6 * EI / L**2
k_local[2, 2] = 4 * EI / L
k_local[2, 4] = -6 * EI / L**2
k_local[2, 5] = 2 * EI / L
k_local[4, 1] = -12 * EI / L**3
k_local[4, 2] = -6 * EI / L**2
k_local[4, 4] = 12 * EI / L**3
k_local[4, 5] = -6 * EI / L**2
k_local[5, 1] = 6 * EI / L**2
k_local[5, 2] = 2 * EI / L
k_local[5, 4] = -6 * EI / L**2
k_local[5, 5] = 4 * EI / L
return k_local
def get_mass_matrix(self, consistent=True):
"""获取单元质量矩阵"""
L = self.L
rho = self.rho
A = self.A
if consistent:
m_local = np.zeros((6, 6))
m_axial = rho * A * L / 6
m_local[0, 0] = 2 * m_axial
m_local[0, 3] = m_axial
m_local[3, 0] = m_axial
m_local[3, 3] = 2 * m_axial
m_bend = rho * A * L / 420
m_local[1, 1] = 156 * m_bend
m_local[1, 2] = 22 * L * m_bend
m_local[1, 4] = 54 * m_bend
m_local[1, 5] = -13 * L * m_bend
m_local[2, 1] = 22 * L * m_bend
m_local[2, 2] = 4 * L**2 * m_bend
m_local[2, 4] = 13 * L * m_bend
m_local[2, 5] = -3 * L**2 * m_bend
m_local[4, 1] = 54 * m_bend
m_local[4, 2] = 13 * L * m_bend
m_local[4, 4] = 156 * m_bend
m_local[4, 5] = -22 * L * m_bend
m_local[5, 1] = -13 * L * m_bend
m_local[5, 2] = -3 * L**2 * m_bend
m_local[5, 4] = -22 * L * m_bend
m_local[5, 5] = 4 * L**2 * m_bend
return m_local
else:
m_local = np.zeros((6, 6))
m_total = rho * A * L
m_local[0, 0] = m_total / 2
m_local[1, 1] = m_total / 2
m_local[3, 3] = m_total / 2
m_local[4, 4] = m_total / 2
m_local[2, 2] = m_total * L**2 / 24
m_local[5, 5] = m_total * L**2 / 24
return m_local
def get_transformation_matrix(self):
"""获取坐标变换矩阵"""
c = self.c
s = self.s
T = np.array([
[c, s, 0, 0, 0, 0],
[-s, c, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, c, s, 0],
[0, 0, 0, -s, c, 0],
[0, 0, 0, 0, 0, 1]
])
return T
def get_global_stiffness(self):
"""获取全局坐标系下的刚度矩阵"""
k_local = self.get_stiffness_matrix()
T = self.get_transformation_matrix()
return T.T @ k_local @ T
def get_global_mass(self, consistent=True):
"""获取全局坐标系下的质量矩阵"""
m_local = self.get_mass_matrix(consistent)
T = self.get_transformation_matrix()
return T.T @ m_local @ T
class BridgeStructure:
"""桥梁结构类"""
def __init__(self, name="桥梁结构"):
self.name = name
self.nodes = []
self.elements = []
self.n_dof = 0
self.K = None
self.M = None
self.fixed_dofs = []
self.support_nodes = []
self.material_props = {}
def add_node(self, x, y):
"""添加节点"""
self.nodes.append([x, y])
return len(self.nodes) - 1
def add_element(self, node1_id, node2_id, E, I, A, rho):
"""添加梁单元"""
elem = BridgeElement(
self.nodes[node1_id],
self.nodes[node2_id],
E, I, A, rho
)
self.elements.append({
'elem': elem,
'nodes': [node1_id, node2_id],
'E': E, 'I': I, 'A': A, 'rho': rho
})
def fix_node(self, node_id, ux=True, uy=True, rz=True):
"""约束节点自由度"""
dof_base = node_id * 3
if ux:
self.fixed_dofs.append(dof_base)
if uy:
self.fixed_dofs.append(dof_base + 1)
if rz:
self.fixed_dofs.append(dof_base + 2)
def add_support(self, node_id):
"""添加支承节点"""
self.support_nodes.append(node_id)
def assemble_global_matrices(self):
"""组装全局刚度矩阵和质量矩阵"""
n_nodes = len(self.nodes)
self.n_dof = n_nodes * 3
self.K = np.zeros((self.n_dof, self.n_dof))
self.M = np.zeros((self.n_dof, self.n_dof))
for elem_info in self.elements:
elem = elem_info['elem']
n1, n2 = elem_info['nodes']
k_global = elem.get_global_stiffness()
m_global = elem.get_global_mass(consistent=True)
dofs = []
for n in [n1, n2]:
dofs.extend([n*3, n*3+1, n*3+2])
for i in range(6):
for j in range(6):
gi, gj = dofs[i], dofs[j]
self.K[gi, gj] += k_global[i, j]
self.M[gi, gj] += m_global[i, j]
def modal_analysis(self, n_modes=10):
"""模态分析"""
free_dofs = [i for i in range(self.n_dof) if i not in self.fixed_dofs]
K_reduced = self.K[np.ix_(free_dofs, free_dofs)]
M_reduced = self.M[np.ix_(free_dofs, free_dofs)]
n_modes = min(n_modes, len(free_dofs))
eigenvalues, eigenvectors = linalg.eigh(K_reduced, M_reduced,
subset_by_index=[0, n_modes-1])
freqs = np.sqrt(np.maximum(eigenvalues, 0)) / (2 * np.pi)
modes = np.zeros((self.n_dof, n_modes))
for i, fdof in enumerate(free_dofs):
modes[fdof, :] = eigenvectors[i, :]
for i in range(n_modes):
max_disp = np.max(np.abs(modes[:, i]))
if max_disp > 0:
modes[:, i] /= max_disp
return freqs, modes, free_dofs
def design_response_spectrum(T, spectrum_type='II', damping=0.05):
"""
设计反应谱（基于中国规范）
Parameters:
-----------
T : array
周期序列 (s)
spectrum_type : str
场地类别 ('I', 'II', 'III', 'IV')
damping : float
阻尼比
Returns:
--------
Sa : array
加速度反应谱 (g)
"""
# 特征周期
Tg_dict = {'I': 0.25, 'II': 0.35, 'III': 0.45, 'IV': 0.65}
Tg = Tg_dict.get(spectrum_type, 0.35)
# 地震影响系数最大值 (8度，0.2g)
alpha_max = 0.16
# 阻尼调整系数
gamma = 0.9 + (0.05 - damping) / (0.5 + 5 * damping)
eta2 = 1 + (0.05 - damping) / (0.06 + 1.7 * damping)
Sa = np.zeros_like(T)
for i, Ti in enumerate(T):
if Ti < 0.1:
# 上升段
Sa[i] = (0.45 + 5.5 * Ti) * alpha_max
elif Ti < Tg:
# 平台段
Sa[i] = eta2 * alpha_max
elif Ti < 5 * Tg:
# 下降段
Sa[i] = eta2 * alpha_max * (Tg / Ti) ** gamma
else:
# 长周期段
Sa[i] = eta2 * alpha_max * (0.2 ** gamma - 0.02 * (Ti - 5 * Tg))
Sa[i] = max(Sa[i], 0.2 * eta2 * alpha_max)
return Sa
def response_spectrum_analysis(bridge, freqs, modes, free_dofs, spectrum_type='II'):
"""
反应谱分析
Parameters:
-----------
bridge : BridgeStructure
桥梁结构
freqs : array
固有频率
modes : ndarray
振型矩阵
free_dofs : list
自由自由度列表
spectrum_type : str
场地类别
Returns:
--------
response : dict
响应结果
"""
n_modes = len(freqs)
n_free = len(free_dofs)
# 提取自由自由度的振型
modes_free = np.zeros((n_free, n_modes))
for i, fdof in enumerate(free_dofs):
modes_free[i, :] = modes[fdof, :]
# 计算模态质量和参与因子
M_ff = bridge.M[np.ix_(free_dofs, free_dofs)]
modal_masses = np.zeros(n_modes)
participation_factors = np.zeros(n_modes)
for i in range(n_modes):
phi_i = modes_free[:, i]
modal_masses[i] = phi_i @ M_ff @ phi_i
if modal_masses[i] > 1e-10:
participation_factors[i] = np.sum(np.abs(phi_i)) / modal_masses[i]
else:
participation_factors[i] = 0.0
# 计算各周期的反应谱值
periods = 1.0 / np.maximum(freqs, 0.01)
Sa = design_response_spectrum(periods, spectrum_type)
# 计算各模态的最大响应
u_modal_max = np.zeros(n_modes)
for i in range(n_modes):
omega_n = 2 * np.pi * freqs[i]
if omega_n < 1e-6:
continue
pf = participation_factors[i]
if abs(pf) < 1e-10:
continue
# 模态最大位移
u_modal_max[i] = pf * Sa[i] * 9.81 / omega_n**2 * 1000
# 组合模态响应（SRSS方法）
u_combined = np.zeros(n_free)
for i in range(n_modes):
phi_i = modes_free[:, i]
u_combined += (phi_i * u_modal_max[i])**2
u_combined = np.sqrt(u_combined)
# 扩展到完整自由度
u_total = np.zeros(bridge.n_dof)
for i, fdof in enumerate(free_dofs):
u_total[fdof] = u_combined[i]
# 计算基底剪力
total_mass = np.sum(bridge.M)
base_shear = total_mass * np.max(Sa) * 9.81 / 1000 # MN
return {
'displacement': u_total,
'modal_response': u_modal_max,
'base_shear': base_shear,
'periods': periods,
'Sa': Sa,
'participation_factors': participation_factors
}
def check_seismic_capacity(bridge, response, material_strength=30e6):
"""
抗震能力验算
Parameters:
-----------
bridge : BridgeStructure
桥梁结构
response : dict
响应结果
material_strength : float
材料强度 (Pa)
Returns:
--------
check_result : dict
验算结果
"""
# 提取关键位置位移
mid_node = 13
pier1_top = 2
pier2_top = 3
mid_disp = response['displacement'][mid_node * 3 + 1] * 1000 # mm
pier1_disp = response['displacement'][pier1_top * 3 + 1] * 1000
pier2_disp = response['displacement'][pier2_top * 3 + 1] * 1000
# 位移限值（假设跨度140m，限值为L/300）
L_bridge = 140000 # mm
disp_limit = L_bridge / 300
# 验算结果
checks = {
'mid_span': {
'value': mid_disp,
'limit': disp_limit,
'ratio': mid_disp / disp_limit,
'passed': mid_disp < disp_limit
},
'pier1': {
'value': pier1_disp,
'limit': disp_limit * 0.5,
'ratio': pier1_disp / (disp_limit * 0.5),
'passed': pier1_disp < disp_limit * 0.5
},
'pier2': {
'value': pier2_disp,
'limit': disp_limit * 0.5,
'ratio': pier2_disp / (disp_limit * 0.5),
'passed': pier2_disp < disp_limit * 0.5
}
}
# 整体评价
all_passed = all(check['passed'] for check in checks.values())
return {
'checks': checks,
'all_passed': all_passed,
'base_shear': response['base_shear']
}
def create_three_span_bridge():
"""创建三跨连续梁桥模型"""
bridge = BridgeStructure("三跨连续梁桥")
E_concrete = 3.45e10
rho_concrete = 2500
A_girder = 8.5
I_girder = 12.0
D_pier = 2.5
A_pier = np.pi * (D_pier/2)**2
I_pier = np.pi * (D_pier/2)**4 / 4
L1, L2, L3 = 40, 60, 40
H_pier = 15
# 桥墩底部
pier1_bottom = bridge.add_node(40, 0)
pier2_bottom = bridge.add_node(100, 0)
# 桥墩顶部
pier1_top = bridge.add_node(40, H_pier)
pier2_top = bridge.add_node(100, H_pier)
# 主梁节点
girder_nodes = []
for i in range(9):
x = i * L1 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
for i in range(1, 9):
x = 40 + i * L2 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
for i in range(1, 9):
x = 100 + i * L3 / 8
node_id = bridge.add_node(x, H_pier)
girder_nodes.append(node_id)
# 桥墩单元
bridge.add_element(pier1_bottom, pier1_top, E_concrete, I_pier, A_pier, rho_concrete)
bridge.add_element(pier2_bottom, pier2_top, E_concrete, I_pier, A_pier, rho_concrete)
# 主梁单元
for i in range(len(girder_nodes) - 1):
bridge.add_element(girder_nodes[i], girder_nodes[i+1],
E_concrete, I_girder, A_girder, rho_concrete)
# 约束
bridge.fix_node(pier1_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(pier2_bottom, ux=True, uy=True, rz=True)
bridge.fix_node(girder_nodes[0], uy=True)
bridge.fix_node(girder_nodes[-1], uy=True)
# 添加支承点
bridge.add_support(pier1_bottom)
bridge.add_support(pier2_bottom)
bridge.add_support(girder_nodes[0])
bridge.add_support(girder_nodes[-1])
return bridge
def plot_response_spectrum(periods, Sa, title="设计反应谱"):
"""绘制设计反应谱"""
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(periods, Sa / 9.81, 'b-', linewidth=2)
ax.set_xlabel('周期 T (s)', fontsize=12)
ax.set_ylabel('加速度 Sa (g)', fontsize=12)
ax.set_title(title, fontsize=14)
ax.grid(True, alpha=0.3)
ax.set_xlim([0, 5])
# 标注特征周期
Tg = 0.35
ax.axvline(x=Tg, color='r', linestyle='--', linewidth=1.5, label=f'Tg={Tg}s')
ax.legend()
plt.tight_layout()
return fig
def plot_displacement_profile(bridge, response, title="位移分布"):
"""绘制位移分布"""
fig, ax = plt.subplots(figsize=(12, 6))
# 主梁节点
girder_nodes = [i for i in range(len(bridge.nodes))
if abs(bridge.nodes[i][1] - 15) < 0.1]
x_pos = [bridge.nodes[i][0] for i in girder_nodes]
disps = [response['displacement'][i*3+1] * 1000 for i in girder_nodes]
ax.plot(x_pos, disps, 'b-o', linewidth=2, markersize=8, label='竖向位移')
ax.fill_between(x_pos, 0, disps, alpha=0.3)
# 位移限值
L_bridge = 140000 # mm
disp_limit = L_bridge / 300
ax.axhline(y=disp_limit, color='r', linestyle='--',
linewidth=2, label=f'限值={disp_limit:.1f}mm')
ax.axhline(y=-disp_limit, color='r', linestyle='--', linewidth=2)
ax.set_xlabel('桥梁位置 (m)', fontsize=12)
ax.set_ylabel('位移 (mm)', fontsize=12)
ax.set_title(title, fontsize=14)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def plot_modal_contribution(freqs, modal_response, title="模态贡献"):
"""绘制模态贡献"""
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 模态位移
ax1 = axes[0]
mode_numbers = np.arange(1, len(modal_response) + 1)
ax1.bar(mode_numbers, modal_response, color='steelblue', edgecolor='black')
ax1.set_xlabel('模态阶数', fontsize=12)
ax1.set_ylabel('模态最大位移 (mm)', fontsize=12)
ax1.set_title('各模态最大位移贡献', fontsize=14)
ax1.grid(True, alpha=0.3, axis='y')
# 累积贡献
ax2 = axes[1]
cumulative = np.cumsum(modal_response**2)
cumulative = np.sqrt(cumulative)
cumulative_pct = cumulative / cumulative[-1] * 100
ax2.plot(mode_numbers, cumulative_pct, 'b-o', linewidth=2, markersize=8)
ax2.axhline(y=90, color='r', linestyle='--', linewidth=1.5, label='90%')
ax2.set_xlabel('模态阶数', fontsize=12)
ax2.set_ylabel('累积贡献 (%)', fontsize=12)
ax2.set_title('模态累积贡献', fontsize=14)
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def main():
"""主程序"""
print("=" * 70)
print("案例4：桥梁抗震设计")
print("=" * 70)
# 创建桥梁模型
print("\n[1] 创建桥梁结构模型...")
bridge = create_three_span_bridge()
bridge.assemble_global_matrices()
print(f" - 节点数: {len(bridge.nodes)}")
print(f" - 单元数: {len(bridge.elements)}")
print(f" - 支承点数: {len(bridge.support_nodes)}")
# 计算总质量
total_mass = np.sum(bridge.M)
print(f" - 总质量: {total_mass/1000:.1f} 吨")
# 模态分析
print("\n[2] 进行模态分析...")
freqs, modes, free_dofs = bridge.modal_analysis(n_modes=10)
print(f" - 计算模态数: 10")
print(f" - 第1阶频率: {freqs[0]:.3f} Hz (周期: {1/freqs[0]:.3f} s)")
print(f" - 第2阶频率: {freqs[1]:.3f} Hz (周期: {1/freqs[1]:.3f} s)")
# 反应谱分析
print("\n[3] 进行反应谱分析...")
response = response_spectrum_analysis(bridge, freqs, modes, free_dofs,
spectrum_type='II')
print(f" - 场地类别: II类")
print(f" - 基底剪力: {response['base_shear']:.2f} MN")
# 抗震验算
print("\n[4] 进行抗震验算...")
check_result = check_seismic_capacity(bridge, response)
print(f"\n 验算结果:")
for location, check in check_result['checks'].items():
status = "通过" if check['passed'] else "不通过"
print(f" - {location}:")
print(f" 位移: {check['value']:.2f} mm / {check['limit']:.2f} mm")
print(f" 利用率: {check['ratio']*100:.1f}%")
print(f" 状态: {status}")
print(f"\n 整体评价: {'满足要求' if check_result['all_passed'] else '需要加强'}")
# 绘制设计反应谱
print("\n[5] 绘制设计反应谱...")
T_range = np.linspace(0.01, 5, 500)
Sa_range = design_response_spectrum(T_range, 'II')
fig_spectrum = plot_response_spectrum(T_range, Sa_range, "II类场地设计反应谱")
plt.savefig('case4_response_spectrum.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case4_response_spectrum.png")
plt.close()
# 绘制位移分布
print("\n[6] 绘制位移分布...")
fig_disp = plot_displacement_profile(bridge, response, "地震位移分布")
plt.savefig('case4_displacement_profile.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case4_displacement_profile.png")
plt.close()
# 绘制模态贡献
print("\n[7] 绘制模态贡献...")
fig_modal = plot_modal_contribution(freqs, response['modal_response'],
"模态贡献分析")
plt.savefig('case4_modal_contribution.png', dpi=150, bbox_inches='tight')
print(" - 已保存: case4_modal_contribution.png")
plt.close()
# 打印设计建议
print("\n" + "=" * 70)
print("抗震设计建议")
print("=" * 70)
print(f"\n1. 结构动力特性:")
print(f" - 基本周期: {1/freqs[0]:.3f} s")
print(f" - 基本频率: {freqs[0]:.3f} Hz")
print(f" - 结构类型: {'刚性结构' if 1/freqs[0] < 0.3 else '柔性结构'}")
print(f"\n2. 地震作用:")
print(f" - 基底剪力: {response['base_shear']:.2f} MN")
print(f" - 剪重比: {response['base_shear']/(total_mass*9.81/1000)*100:.2f}%")
print(f"\n3. 变形验算:")
max_disp = max(check['value'] for check in check_result['checks'].values())
max_ratio = max(check['ratio'] for check in check_result['checks'].values())
print(f" - 最大位移: {max_disp:.2f} mm")
print(f" - 最大利用率: {max_ratio*100:.1f}%")
print(f"\n4. 设计建议:")
if check_result['all_passed']:
print(f" - 当前设计满足抗震要求")
if max_ratio < 0.7:
print(f" - 可考虑优化设计，降低材料用量")
else:
print(f" - 需要加强结构刚度或增加阻尼")
print(f" - 建议措施:")
print(f" * 增大截面尺寸")
print(f" * 采用隔震支座")
print(f" * 增加阻尼器")
print(f"\n5. 隔震设计概念:")
T_isolated = 2.5 # 隔震后周期
Sa_isolated = design_response_spectrum(np.array([T_isolated]), 'II')[0]
Sa_original = design_response_spectrum(np.array([1/freqs[0]]), 'II')[0]
reduction_ratio = Sa_isolated / Sa_original if Sa_original > 0 else 1.0
print(f" - 隔震后周期: {T_isolated} s")
print(f" - 地震力折减系数: {reduction_ratio:.2f}")
print(f" - 预期减震效果: {(1-reduction_ratio)*100:.1f}%")
print("\n" + "=" * 70)
print("案例分析完成！")
print("=" * 70)
print("\n生成文件:")
print(" - case4_response_spectrum.png: 设计反应谱")
print(" - case4_displacement_profile.png: 位移分布图")
print(" - case4_modal_contribution.png: 模态贡献图")
if __name__ == "__main__":
main()
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
1
0- 0
已为社区贡献98条内容
所有评论(0)