第十七篇：结构地震响应分析与设计反应谱
摘要
地震是结构工程中最具破坏性的自然灾害之一，对结构进行地震响应分析是确保结构抗震安全性的关键环节。本教程系统地介绍结构地震响应分析的基本理论、计算方法和工程应用。首先阐述地震动的基本特性和地震波的数学描述；然后详细介绍反应谱理论及其计算方法；接着讲解设计反应谱的构建原理和规范要求；最后通过Python实现地震波生成、反应谱计算、设计反应谱构建以及多自由度体系的地震响应分析。通过本教程的学习，读者将掌握结构抗震设计的核心理论和实用技术。
关键词
地震响应，反应谱，设计反应谱，地震波，抗震设计，Python仿真
1. 引言
1.1 地震工程的重要性
地震是一种突发性强、破坏性大的自然灾害，对人类社会造成巨大的人员伤亡和财产损失。结构抗震设计是土木工程领域最重要的研究方向之一，其目标是使结构在遭遇地震时能够保持足够的承载能力和变形能力，避免倒塌，保护生命安全。
地震工程的主要研究内容：
- 地震动特性研究：地震波的传播规律和地面运动特性
- 结构地震响应分析：结构在地震作用下的动力响应计算
- 抗震设计方法：基于性能的结构抗震设计理论和方法
- 减震隔震技术：通过特殊装置减小结构地震响应
1.2 地震动的基本特性
地震动是指地震引起的地面运动，通常用地面加速度、速度和位移时程来描述。
地震动的主要参数：
- 峰值加速度（PGA）：地面加速度的最大值，是衡量地震动强度的基本指标
- 峰值速度（PGV）：地面速度的最大值，与结构的破坏程度密切相关
- 峰值位移（PGD）：地面位移的最大值，对长周期结构影响较大
- 持时：强震动的持续时间，影响结构的累积损伤
- 频谱特性：地震动的频率成分，决定对不同周期结构的破坏作用
地震动的随机性：
地震动具有显著的随机性，即使在同一地点、相同震级的地震，其地面运动特性也可能差异很大。这种随机性来源于：
- 震源机制的不确定性
- 传播路径的复杂性
- 局部场地条件的影响
1.3 反应谱理论的意义
反应谱是地震工程中最核心的概念之一，它描述了单自由度体系在不同周期下的最大响应。
反应谱的物理意义：
- 反应谱反映了地震动对不同周期结构的破坏潜力
- 通过反应谱可以快速估算结构的地震响应
- 设计反应谱是抗震规范的基础
反应谱的优势：
- 将复杂的时程分析简化为谱分析
- 便于工程设计和规范制定
- 可以考虑地震动的统计特性
2. 地震波的数学描述
2.1 地震波的类型
根据传播方式，地震波可分为体波和面波：
体波：
- P波（纵波）：质点振动方向与波的传播方向一致，传播速度最快
- S波（横波）：质点振动方向与波的传播方向垂直，传播速度次之
面波：
- 瑞利波：沿地表传播，质点做椭圆运动
- 洛夫波：沿地表传播，质点做水平横向振动
2.2 地震动的数学模型
确定性模型：
用明确的数学函数描述地震动时程，常用的有：
- 简谐波模型
- 脉冲模型
- 合成地震波模型
随机模型：
将地震动视为随机过程，常用的有：
- 白噪声模型
- 过滤白噪声模型（Kanai-Tajimi模型）
- 非平稳随机模型
2.3 Kanai-Tajimi模型
Kanai-Tajimi模型是描述地震动最常用的随机模型，它将地震动视为白噪声通过滤波器后的输出。
功率谱密度函数：
Su¨g(ω)=1+4ζg2(ω/ωg)2[1−(ω/ωg)2]2+4ζg2(ω/ωg)2S0S_{\ddot{u}_g}(\omega) = \frac{1 + 4\zeta_g^2(\omega/\omega_g)^2}{[1-(\omega/\omega_g)^2]^2 + 4\zeta_g^2(\omega/\omega_g)^2} S_0Su¨g(ω)=[1−(ω/ωg)2]2+4ζg2(ω/ωg)21+4ζg2(ω/ωg)2S0
其中：
- ωg\omega_gωg：场地特征频率
- ζg\zeta_gζg：场地阻尼比
- S0S_0S0：白噪声强度
不同场地条件：
- 硬土：ωg=8π\omega_g = 8\piωg=8π rad/s, ζg=0.6\zeta_g = 0.6ζg=0.6
- 中硬土：ωg=5π\omega_g = 5\piωg=5π rad/s, ζg=0.6\zeta_g = 0.6ζg=0.6
- 软土：ωg=2.5π\omega_g = 2.5\piωg=2.5π rad/s, ζg=0.85\zeta_g = 0.85ζg=0.85
3. 反应谱理论
3.1 反应谱的定义
反应谱是指在给定地震动作用下，单自由度体系的最大响应（位移、速度或加速度）与体系自振周期之间的关系曲线。
位移反应谱：
Sd(T)=maxt∣x(t)∣S_d(T) = \max_t |x(t)|Sd(T)=tmax∣x(t)∣
速度反应谱：
Sv(T)=maxt∣x˙(t)∣S_v(T) = \max_t |\dot{x}(t)|Sv(T)=tmax∣x˙(t)∣
加速度反应谱：
Sa(T)=maxt∣x¨(t)+u¨g(t)∣S_a(T) = \max_t |\ddot{x}(t) + \ddot{u}_g(t)|Sa(T)=tmax∣x¨(t)+u¨g(t)∣
3.2 反应谱的计算方法
对于单自由度体系，运动方程为：
mx¨+cx˙+kx=−mu¨gm\ddot{x} + c\dot{x} + kx = -m\ddot{u}_gmx¨+cx˙+kx=−mu¨g
或：
x¨+2ζωnx˙+ωn2x=−u¨g\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2x = -\ddot{u}_gx¨+2ζωnx˙+ωn2x=−u¨g
计算步骤：
- 选择一系列周期值 TiT_iTi
- 对每个周期，计算相应的 omegan=2π/Ti\\ omega_n = 2\pi/T_iomegan=2π/Ti
- 数值求解运动方程，得到响应时程
- 提取最大响应值
- 绘制反应谱曲线
3.3 反应谱的特性
反应谱的形状特征：
- 短周期段（T < 0.1s）：加速度反应谱接近常数，等于地面最大加速度
- 中周期段（0.1s < T < T_g）：速度反应谱接近常数
- 长周期段（T > T_g）：位移反应谱接近常数
特征周期 TgT_gTg：
特征周期是反应谱形状发生显著变化的周期点，与场地条件密切相关：
- 硬土：Tg≈0.2T_g \approx 0.2Tg≈0.2 s
- 中硬土：Tg≈0.4T_g \approx 0.4Tg≈0.4 s
- 软土：Tg≈0.6T_g \approx 0.6Tg≈0.6 s
3.4 伪反应谱
伪反应谱是工程实践中常用的简化概念：
伪速度反应谱：
PSv=ωnSdPS_v = \omega_n S_dPSv=ωnSd
伪加速度反应谱：
PSa=ωn2SdPS_a = \omega_n^2 S_dPSa=ωn2Sd
对于小阻尼情况，伪反应谱与真实反应谱非常接近。
4. 设计反应谱
4.1 设计反应谱的概念
设计反应谱是基于大量地震记录统计分析得到的平滑化反应谱，用于工程抗震设计。
设计反应谱的特点：
- 基于统计平均，具有代表性
- 形状平滑，便于工程应用
- 考虑场地条件的影响
- 具有规定的超越概率
4.2 中国规范设计反应谱
根据《建筑抗震设计规范》（GB 50011），设计反应谱的表达式为：
上升段（0<T<0.10 < T < 0.10<T<0.1 s）：
α=η2αmax\alpha = \eta_2 \alpha_{max}α=η2αmax
平台段（0.1≤T≤Tg0.1 \leq T \leq T_g0.1≤T≤Tg）：
α=αmax\alpha = \alpha_{max}α=αmax
下降段（Tg<T≤5TgT_g < T \leq 5T_gTg<T≤5Tg）：
α=(TgT)γαmax\alpha = \left(\frac{T_g}{T}\right)^\gamma \alpha_{max}α=(TTg)γαmax
长周期段（5Tg<T≤6.05T_g < T \leq 6.05Tg<T≤6.0 s）：
α=[η20.2γ−η1(T−5Tg)]αmax\alpha = [\eta_2 0.2^\gamma - \eta_1(T - 5T_g)] \alpha_{max}α=[η20.2γ−η1(T−5Tg)]αmax
其中：
- αmax\alpha_{max}αmax：地震影响系数最大值
- TgT_gTg：特征周期
- γ\gammaγ：下降段衰减指数
- η1\eta_1η1：直线下降段斜率调整系数
- η2\eta_2η2：阻尼调整系数
4.3 地震影响系数最大值
地震影响系数最大值与设防烈度和设计地震分组有关：
| 设防烈度 | 6度 | 7度 | 8度 | 9度 |
|---|---|---|---|---|
| αmax\alpha_{max}αmax | 0.04 | 0.08(0.12) | 0.16(0.24) | 0.32 |
4.4 特征周期
特征周期与场地类别和设计地震分组有关：
| 场地类别 | 第一组 | 第二组 | 第三组 |
|---|---|---|---|
| I0 | 0.20 | 0.25 | 0.30 |
| I1 | 0.25 | 0.30 | 0.35 |
| II | 0.35 | 0.40 | 0.45 |
| III | 0.45 | 0.55 | 0.65 |
| IV | 0.65 | 0.75 | 0.90 |
5. 多自由度体系地震响应分析
5.1 振型分解法
对于多自由度体系，地震响应可以通过振型分解法计算：
运动方程：
Mx¨+Cx˙+Kx=−M1u¨g\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = -\mathbf{M}\mathbf{1}\ddot{u}_gMx¨+Cx˙+Kx=−M1u¨g
振型分解：
x(t)=∑i=1nϕiqi(t)\mathbf{x}(t) = \sum_{i=1}^{n} \phi_i q_i(t)x(t)=i=1∑nϕiqi(t)
广义坐标方程：
q¨i+2ζiωiq˙i+ωi2qi=−γiu¨g\ddot{q}_i + 2\zeta_i\omega_i\dot{q}_i + \omega_i^2 q_i = -\gamma_i \ddot{u}_gq¨i+2ζiωiq˙i+ωi2qi=−γiu¨g
其中，γi=ϕiTM1ϕiTMϕi\gamma_i = \frac{\phi_i^T \mathbf{M} \mathbf{1}}{\phi_i^T \mathbf{M} \phi_i}γi=ϕiTMϕiϕiTM1 为振型参与系数。
5.2 振型组合方法
由于各振型的最大响应不会同时发生，需要采用适当的组合方法：
SRSS法（平方和开方法）：
R=∑i=1nRi2R = \sqrt{\sum_{i=1}^{n} R_i^2}R=i=1∑nRi2
适用于频率分离较好的情况。
CQC法（完全二次组合法）：
R=∑i=1n∑j=1nρijRiRjR = \sqrt{\sum_{i=1}^{n}\sum_{j=1}^{n} \rho_{ij} R_i R_j}R=i=1∑nj=1∑nρijRiRj
适用于频率密集的情况，其中 ρij\rho_{ij}ρij 为互相关系数。
5.3 底部剪力法
对于规则结构，可以采用简化的底部剪力法：
总水平地震作用：
FEk=α1GeqF_{Ek} = \alpha_1 G_{eq}FEk=α1Geq
其中：
- α1\alpha_1α1：相应于结构基本周期的水平地震影响系数
- GeqG_{eq}Geq：结构等效总重力荷载
各层水平地震作用：
Fi=GiHi∑j=1nGjHjFEk(1−δn)F_i = \frac{G_i H_i}{\sum_{j=1}^{n} G_j H_j} F_{Ek}(1 - \delta_n)Fi=∑j=1nGjHjGiHiFEk(1−δn)
顶部附加地震作用：
ΔFn=δnFEk\Delta F_n = \delta_n F_{Ek}ΔFn=δnFEk
6. Python实现案例
6.1 案例1：地震波生成与处理
本案例实现基于Kanai-Tajimi模型的地震波生成和频谱分析。
# -*- coding: utf-8 -*-
"""
主题017：结构地震响应分析与设计反应谱 - 案例1
地震波生成与处理
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题017'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("=" * 70)
print("主题017：结构地震响应分析与设计反应谱")
print("案例1：地震波生成与处理")
print("=" * 70)
# ============================================================
# Kanai-Tajimi模型参数
# ============================================================
print("\nKanai-Tajimi模型参数：")
print("-" * 70)
# 不同场地条件
site_conditions = {
'硬土': {
'omega_g': 8*np.pi, 'zeta_g': 0.6, 'S0': 0.01},
'中硬土': {
'omega_g': 5*np.pi, 'zeta_g': 0.6, 'S0': 0.02},
'软土': {
'omega_g': 2.5*np.pi, 'zeta_g': 0.85, 'S0': 0.05}
}
for site, params in site_conditions.items():
print(f" {
site}: ω_g = {
params['omega_g']:.4f} rad/s, ζ_g = {
params['zeta_g']}, S0 = {
params['S0']}")
# ============================================================
# 生成地震波
# ============================================================
def generate_earthquake(dt, t_end, site_type='中硬土', seed=None):
"""
基于Kanai-Tajimi模型生成地震波
"""
if seed is not None:
np.random.seed(seed)
params = site_conditions[site_type]
omega_g = params['omega_g']
zeta_g = params['zeta_g']
S0 = params['S0']
t = np.arange(0, t_end, dt)
n = len(t)
# 生成白噪声
white_noise = np.random.randn(n)
# 构建Kanai-Tajimi滤波器
# H(ω) = (ω_g^2 + 2iζ_gω_gω) / (ω_g^2 - ω^2 + 2iζ_gω_gω)
# 使用状态空间实现
# 滤波器系数
a = [1, 2*zeta_g*omega_g, omega_g**2]
b = [omega_g**2, 2*zeta_g*omega_g, 0]
# 应用滤波器
filtered = signal.lfilter(b, a, white_noise)
# 应用包络函数（模拟非平稳性）
t_mid = t_end / 2
envelope = np.zeros(n)
for i, ti in enumerate(t):
if ti < t_mid:
envelope[i] = (ti / t_mid)**2
else:
envelope[i] = np.exp(-2 * (ti - t_mid) / t_mid)
# 合成地震波
acc = filtered * envelope * np.sqrt(S0 / dt)
# 标准化（调整PGA到目标值）
target_pga = 0.2 * 9.8 # 目标PGA = 0.2g
current_pga = np.max(np.abs(acc))
acc = acc * target_pga / current_pga
return t, acc
# ============================================================
# 计算速度和位移
# ============================================================
def integrate_acceleration(t, acc):
"""通过积分计算速度和位移"""
dt = t[1] - t[0]
# 高通滤波去除趋势
sos = signal.butter(4, 0.1, 'high', fs=1/dt, output='sos')
acc_filtered = signal.sosfilt(sos, acc)
# 积分得到速度
vel = np.cumsum(acc_filtered) * dt
vel = signal.sosfilt(sos, vel) # 再次滤波
# 积分得到位移
disp = np.cumsum(vel) * dt
disp = signal.sosfilt(sos, disp) # 再次滤波
return vel, disp
# ============================================================
# 生成不同场地的地震波
# ============================================================
print("\n生成地震波...")
print("-" * 70)
dt = 0.01
t_end = 30.0
earthquake_data = {
}
for site in site_conditions.keys():
print(f" 生成{
site}地震波...")
t, acc = generate_earthquake(dt, t_end, site, seed=42)
vel, disp = integrate_acceleration(t, acc)
earthquake_data[site] = {
't': t,
'acc': acc,
'vel': vel,
'disp': disp,
'pga': np.max(np.abs(acc)),
'pgv': np.max(np.abs(vel)),
'pgd': np.max(np.abs(disp))
}
# ============================================================
# 打印地震波参数
# ============================================================
print("\n地震波参数：")
print("-" * 70)
print(f"{
'场地':<10} {
'PGA (m/s²)':<15} {
'PGV (m/s)':<15} {
'PGD (m)':<15}")
print("-" * 70)
for site, data in earthquake_data.items():
print(f"{
site:<10} {
data['pga']:<15.4f} {
data['pgv']:<15.4f} {
data['pgd']:<15.4f}")
# ============================================================
# 绘制地震波时程
# ============================================================
print("\n正在生成可视化图表...")
fig1, axes = plt.subplots(3, 1, figsize=(14, 12))
colors = ['steelblue', 'forestgreen', 'coral']
for i, (site, data) in enumerate(earthquake_data.items()):
# 加速度
axes[0].plot(data['t'], data['acc'], color=colors[i], linewidth=1, label=f'{
site}')
# 速度
axes[1].plot(data['t'], data['vel'], color=colors[i], linewidth=1, label=f'{
site}')
# 位移
axes[2].plot(data['t'], data['disp'], color=colors[i], linewidth=1, label=f'{
site}')
axes[0].set_ylabel('Acceleration (m/s²)', fontsize=12)
axes[0].set_title('Ground Acceleration Time History', fontsize=14, fontweight='bold')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[1].set_ylabel('Velocity (m/s)', fontsize=12)
axes[1].set_title('Ground Velocity Time History', fontsize=14, fontweight='bold')
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[2].set_xlabel('Time (s)', fontsize=12)
axes[2].set_ylabel('Displacement (m)', fontsize=12)
axes[2].set_title('Ground Displacement Time History', fontsize=14, fontweight='bold')
axes[2].legend(fontsize=11)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{
output_dir}/earthquake_time_history.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 地震波时程图已保存")
# ============================================================
# 傅里叶谱分析
# ============================================================
print("\n进行傅里叶谱分析...")
fig2, axes = plt.subplots(1, 2, figsize=(14, 5))
for i, (site, data) in enumerate(earthquake_data.items()):
# FFT
n = len(data['t'])
freqs = np.fft.fftfreq(n, dt)[:n//2]
fft_values = np.fft.fft(data['acc'])
fft_magnitude = np.abs(fft_values)[:n//2] * 2 / n
# 傅里叶幅值谱
axes[0].semilogy(freqs, fft_magnitude, color=colors[i], linewidth=1.5, label=f'{
site}')
# 功率谱密度
f_psd, psd = signal.welch(data['acc'], fs=1/dt, nperseg=1024)
axes[1].semilogy(f_psd, psd, color=colors[i], linewidth=1.5, label=f'{
site}')
axes[0].set_xlabel('Frequency (Hz)', fontsize=12)
axes[0].set_ylabel('Fourier Amplitude', fontsize=12)
axes[0].set_title('Fourier Amplitude Spectrum', fontsize=14, fontweight='bold')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 25)
axes[1].set_xlabel('Frequency (Hz)', fontsize=12)
axes[1].set_ylabel('PSD (m²/s⁴/Hz)', fontsize=12)
axes[1].set_title('Power Spectral Density', fontsize=14, fontweight='bold')
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 25)
plt.tight_layout()
plt.savefig(f'{
output_dir}/earthquake_spectrum.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 地震波频谱图已保存")
# ============================================================
# 生成地震波动画
# ============================================================
print("\n正在生成地震波动画...")
fig3, axes = plt.subplots(3, 1, figsize=(12, 10))
site = '中硬土'
data = earthquake_data[site]
# 初始化线条
line_acc, = axes[0].plot([], [], 'b-', linewidth=1.5)
line_vel, = axes[1].plot([], [], 'g-', linewidth=1.5)
line_disp, = axes[2].plot([], [], 'r-', linewidth=1.5)
# 设置坐标轴
axes[0].set_xlim(0, t_end)
axes[0].set_ylim(np.min(data['acc'])*1.1, np.max(data['acc'])*1.1)
axes[0].set_ylabel('Acceleration (m/s²)', fontsize=11)
axes[0].set_title(f'Earthquake Ground Motion ({
site})', fontsize=13, fontweight='bold')
axes[0].grid(True, alpha=0.3)
axes[1].set_xlim(0, t_end)
axes[1].set_ylim(np.min(data['vel'])*1.1, np.max(data['vel'])*1.1)
axes[1].set_ylabel('Velocity (m/s)', fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[2].set_xlim(0, t_end)
axes[2].set_ylim(np.min(data['disp'])*1.1, np.max(data['disp'])*1.1)
axes[2].set_xlabel('Time (s)', fontsize=11)
axes[2].set_ylabel('Displacement (m)', fontsize=11)
axes[2].grid(True, alpha=0.3)
from matplotlib.animation import FuncAnimation
def update(frame):
idx = frame * 10
if idx > len(data['t']):
idx = len(data['t'])
line_acc.set_data(data['t'][:idx], data['acc'][:idx])
line_vel.set_data(data['t'][:idx], data['vel'][:idx])
line_disp.set_data(data['t'][:idx], data['disp'][:idx])
return line_acc, line_vel, line_disp
n_frames = len(data['t']) // 10
anim = FuncAnimation(fig3, update, frames=n_frames, interval=50, blit=True)
anim.save(f'{
output_dir}/earthquake_animation.gif', writer='pillow', fps=20, dpi=100)
plt.close()
print(" 地震波动画已保存")
print("\n" + "=" * 70)
print("地震波生成与处理完成！")
print("=" * 70)
6.2 案例2：反应谱计算
本案例实现单自由度体系地震反应谱的计算。
# -*- coding: utf-8 -*-
"""
主题017：结构地震响应分析与设计反应谱 - 案例2
反应谱计算
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题017'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("=" * 70)
print("主题017：结构地震响应分析与设计反应谱")
print("案例2：反应谱计算")
print("=" * 70)
# ============================================================
# 生成地震波
# ============================================================
def generate_earthquake(dt, t_end, omega_g=5*np.pi, zeta_g=0.6, S0=0.02, seed=None):
"""生成地震波"""
if seed is not None:
np.random.seed(seed)
t = np.arange(0, t_end, dt)
n = len(t)
white_noise = np.random.randn(n)
freqs = np.fft.fftfreq(n, dt)
omega = 2 * np.pi * freqs
omega_sq = omega**2
numerator = omega_g**4 + 4*zeta_g**2*omega_g**2*omega_sq
denominator = (omega_g**2 - omega_sq)**2 + 4*zeta_g**2*omega_g**2*omega_sq
psd = np.zeros_like(omega)
mask = denominator > 1e-10
psd[mask] = S0 * numerator[mask] / denominator[mask]
white_fft = np.fft.fft(white_noise)
filtered_fft = white_fft * np.sqrt(psd * n * dt / 2)
filtered = np.real(np.fft.ifft(filtered_fft))
t_mid = t_end / 2
envelope = np.zeros(n)
for i, ti in enumerate(t):
if ti < t_mid:
envelope[i] = (ti / t_mid)**2
else:
envelope[i] = np.exp(-2 * (ti - t_mid) / t_mid)
acc = filtered * envelope
target_pga = 0.2 * 9.8
current_pga = np.max(np.abs(acc))
if current_pga > 0:
acc = acc * target_pga / current_pga
return t, acc
# ============================================================
# 反应谱计算函数
# ============================================================
def compute_response_spectrum(acc, dt, T_range, zeta=0.05):
"""
计算反应谱
参数:
acc: 地面加速度时程
dt: 时间步长
T_range: 周期范围
zeta: 阻尼比
返回:
Sd: 位移反应谱
Sv: 速度反应谱
Sa: 加速度反应谱
PSv: 伪速度反应谱
PSa: 伪加速度反应谱
"""
n_T = len(T_range)
Sd = np.zeros(n_T)
Sv = np.zeros(n_T)
Sa = np.zeros(n_T)
PSv = np.zeros(n_T)
PSa = np.zeros(n_T)
for i, T in enumerate(T_range):
if T == 0:
Sd[i] = 0
Sv[i] = 0
Sa[i] = np.max(np.abs(acc))
PSv[i] = 0
PSa[i] = np.max(np.abs(acc))
else:
omega_n = 2 * np.pi / T
# Newmark-beta方法参数（平均加速度法）
gamma = 0.5
beta = 0.25
# 初始化
n_steps = len(acc)
x = np.zeros(n_steps)
v = np.zeros(n_steps)
a = np.zeros(n_steps)
# 等效刚度
k_eff = omega_n**2 + gamma / (beta * dt) * 2 * zeta * omega_n + 1 / (beta * dt**2)
for j in range(n_steps - 1):
# 等效荷载
p_eff = -acc[j+1] + (1/(beta*dt**2))*x[j] + (1/(beta*dt))*v[j] + (1/(2*beta)-1)*a[j]
# 计算位移
x[j+1] = p_eff / k_eff
# 计算速度和加速度
v[j+1] = gamma/(beta*dt)*(x[j+1]-x[j]) + (1-gamma/beta)*v[j] + dt*(1-gamma/(2*beta))*a[j]
a[j+1] = 1/(beta*dt**2)*(x[j+1]-x[j]) - 1/(beta*dt)*v[j] - (1/(2*beta)-1)*a[j]
# 计算绝对加速度
abs_acc = a + acc
# 提取最大值
Sd[i] = np.max(np.abs(x))
Sv[i] = np.max(np.abs(v))
Sa[i] = np.max(np.abs(abs_acc))
PSv[i] = omega_n * Sd[i]
PSa[i] = omega_n**2 * Sd[i]
return Sd, Sv, Sa, PSv, PSa
# ============================================================
# 生成地震波并计算反应谱
# ============================================================
print("\n生成地震波...")
print("-" * 70)
dt = 0.01
t_end = 30.0
t, acc = generate_earthquake(dt, t_end, seed=42)
print(f" 时间步长: {
dt} s")
print(f" 持续时间: {
t_end} s")
print(f" PGA: {
np.max(np.abs(acc)):.4f} m/s²")
# 计算反应谱
print("\n计算反应谱...")
print("-" * 70)
T_range = np.linspace(0.01, 6.0, 200)
zeta_values = [0.02, 0.05, 0.10]
response_spectra = {
}
for zeta in zeta_values:
print(f" 计算阻尼比 ζ = {
zeta} 的反应谱...")
Sd, Sv, Sa, PSv, PSa = compute_response_spectrum(acc, dt, T_range, zeta)
response_spectra[zeta] = {
'Sd': Sd, 'Sv': Sv, 'Sa': Sa, 'PSv': PSv, 'PSa': PSa
}
# ============================================================
# 绘制反应谱
# ============================================================
print("\n正在生成可视化图表...")
fig1, axes = plt.subplots(2, 2, figsize=(14, 10))
colors = ['steelblue', 'forestgreen', 'coral']
for i, zeta in enumerate(zeta_values):
data = response_spectra[zeta]
# 位移反应谱
axes[0, 0].plot(T_range, data['Sd']*1000, color=colors[i], linewidth=2, label=f'ζ={
zeta}')
# 速度反应谱
axes[0, 1].plot(T_range, data['Sv'], color=colors[i], linewidth=2, label=f'ζ={
zeta}')
# 加速度反应谱
axes[1, 0].plot(T_range, data['Sa']/9.8, color=colors[i], linewidth=2, label=f'ζ={
zeta}')
# 伪加速度反应谱
axes[1, 1].plot(T_range, data['PSa']/9.8, color=colors[i], linewidth=2, label=f'ζ={
zeta}')
axes[0, 0].set_xlabel('Period (s)', fontsize=12)
axes[0, 0].set_ylabel('Displacement (mm)', fontsize=12)
axes[0, 0].set_title('Displacement Response Spectrum', fontsize=14, fontweight='bold')
axes[0, 0].legend(fontsize=11)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].set_xlabel('Period (s)', fontsize=12)
axes[0, 1].set_ylabel('Velocity (m/s)', fontsize=12)
axes[0, 1].set_title('Velocity Response Spectrum', fontsize=14, fontweight='bold')
axes[0, 1].legend(fontsize=11)
axes[0, 1].grid(True, alpha=0.3)
axes[1, 0].set_xlabel('Period (s)', fontsize=12)
axes[1, 0].set_ylabel('Acceleration (g)', fontsize=12)
axes[1, 0].set_title('Absolute Acceleration Response Spectrum', fontsize=14, fontweight='bold')
axes[1, 0].legend(fontsize=11)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].set_xlabel('Period (s)', fontsize=12)
axes[1, 1].set_ylabel('Pseudo Acceleration (g)', fontsize=12)
axes[1, 1].set_title('Pseudo Acceleration Response Spectrum', fontsize=14, fontweight='bold')
axes[1, 1].legend(fontsize=11)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{
output_dir}/response_spectrum.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 反应谱图已保存")
# ============================================================
# 绘制真实反应谱与伪反应谱对比
# ============================================================
fig2, axes = plt.subplots(1, 2, figsize=(14, 5))
zeta = 0.05
data = response_spectra[zeta]
# 速度反应谱对比
axes[0].plot(T_range, data['Sv'], 'b-', linewidth=2, label='True Velocity')
axes[0].plot(T_range, data['PSv'], 'r--', linewidth=2, label='Pseudo Velocity')
axes[0].set_xlabel('Period (s)', fontsize=12)
axes[0].set_ylabel('Velocity (m/s)', fontsize=12)
axes[0].set_title('True vs Pseudo Velocity Response Spectrum (ζ=0.05)', fontsize=13, fontweight='bold')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# 加速度反应谱对比
axes[1].plot(T_range, data['Sa']/9.8, 'b-', linewidth=2, label='True Acceleration')
axes[1].plot(T_range, data['PSa']/9.8, 'r--', linewidth=2, label='Pseudo Acceleration')
axes[1].set_xlabel('Period (s)', fontsize=12)
axes[1].set_ylabel('Acceleration (g)', fontsize=12)
axes[1].set_title('True vs Pseudo Acceleration Response Spectrum (ζ=0.05)', fontsize=13, fontweight='bold')
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{
output_dir}/true_vs_pseudo_spectrum.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 真实与伪反应谱对比图已保存")
# ============================================================
# 生成反应谱动画
# ============================================================
print("\n正在生成反应谱动画...")
fig3, axes = plt.subplots(2, 2, figsize=(14, 10))
zeta = 0.05
data = response_spectra[zeta]
# 初始化线条
line_sd, = axes[0, 0].plot([], [], 'b-', linewidth=2)
line_sv, = axes[0, 1].plot([], [], 'g-', linewidth=2)
line_sa, = axes[1, 0].plot([], [], 'r-', linewidth=2)
line_psa, = axes[1, 1].plot([], [], 'm-', linewidth=2)
# 设置坐标轴
axes[0, 0].set_xlim(0, 6)
axes[0, 0].set_ylim(0, np.max(data['Sd'])*1000*1.1)
axes[0, 0].set_xlabel('Period (s)', fontsize=11)
axes[0, 0].set_ylabel('Displacement (mm)', fontsize=11)
axes[0, 0].set_title('Displacement Response Spectrum', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].set_xlim(0, 6)
axes[0, 1].set_ylim(0, np.max(data['Sv'])*1.1)
axes[0, 1].set_xlabel('Period (s)', fontsize=11)
axes[0, 1].set_ylabel('Velocity (m/s)', fontsize=11)
axes[0, 1].set_title('Velocity Response Spectrum', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)
axes[1, 0].set_xlim(0, 6)
axes[1, 0].set_ylim(0, np.max(data['Sa']/9.8)*1.1)
axes[1, 0].set_xlabel('Period (s)', fontsize=11)
axes[1, 0].set_ylabel('Acceleration (g)', fontsize=11)
axes[1, 0].set_title('Acceleration Response Spectrum', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].set_xlim(0, 6)
axes[1, 1].set_ylim(0, np.max(data['PSa']/9.8)*1.1)
axes[1, 1].set_xlabel('Period (s)', fontsize=11)
axes[1, 1].set_ylabel('Pseudo Acceleration (g)', fontsize=11)
axes[1, 1].set_title('Pseudo Acceleration Response Spectrum', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)
from matplotlib.animation import FuncAnimation
def update(frame):
idx = frame + 1
if idx > len(T_range):
idx = len(T_range)
line_sd.set_data(T_range[:idx], data['Sd'][:idx]*1000)
line_sv.set_data(T_range[:idx], data['Sv'][:idx])
line_sa.set_data(T_range[:idx], data['Sa'][:idx]/9.8)
line_psa.set_data(T_range[:idx], data['PSa'][:idx]/9.8)
return line_sd, line_sv, line_sa, line_psa
n_frames = len(T_range)
anim = FuncAnimation(fig3, update, frames=n_frames, interval=30, blit=True)
anim.save(f'{
output_dir}/response_spectrum_animation.gif', writer='pillow', fps=30, dpi=100)
plt.close()
print(" 反应谱动画已保存")
# ============================================================
# 打印关键结果
# ============================================================
print("\n反应谱关键结果：")
print("-" * 70)
for zeta in zeta_values:
data = response_spectra[zeta]
max_sd_idx = np.argmax(data['Sd'])
max_sa_idx = np.argmax(data['Sa'])
print(f"\n阻尼比 ζ = {
zeta}:")
print(f" 最大位移: {
data['Sd'][max_sd_idx]*1000:.4f} mm (T={
T_range[max_sd_idx]:.2f}s)")
print(f" 最大加速度: {
data['Sa'][max_sa_idx]/9.8:.4f} g (T={
T_range[max_sa_idx]:.2f}s)")
print("\n" + "=" * 70)
print("反应谱计算完成！")
print("=" * 70)
6.3 案例3：设计反应谱（中国规范GB 50011）
本案例实现中国《建筑抗震设计规范》（GB 50011）的设计反应谱。
# -*- coding: utf-8 -*-
"""
主题017：结构地震响应分析与设计反应谱 - 案例3
设计反应谱（中国规范GB 50011）
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题017'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("=" * 70)
print("主题017：结构地震响应分析与设计反应谱")
print("案例3：设计反应谱（中国规范GB 50011）")
print("=" * 70)
# ============================================================
# 中国规范设计反应谱参数
# ============================================================
print("\n中国规范设计反应谱参数：")
print("-" * 70)
# 地震影响系数最大值（多遇地震）
alpha_max_table = {
'6度': 0.04,
'7度': 0.08,
'7度(0.15g)': 0.12,
'8度': 0.16,
'8度(0.30g)': 0.24,
'9度': 0.32
}
# 特征周期（第一组）
Tg_table = {
'I0': 0.20,
'I1': 0.25,
'II': 0.35,
'III': 0.45,
'IV': 0.65
}
print(" 地震影响系数最大值（多遇地震）：")
for intensity, alpha_max in alpha_max_table.items():
print(f" {
intensity}: αmax = {
alpha_max}")
print("\n 特征周期（第一组）：")
for site, Tg in Tg_table.items():
print(f" {
site}类场地: Tg = {
Tg} s")
# ============================================================
# 设计反应谱计算函数
# ============================================================
def compute_design_spectrum(T, alpha_max, Tg, zeta=0.05):
"""
计算中国规范设计反应谱
参数:
T: 周期数组
alpha_max: 地震影响系数最大值
Tg: 特征周期
zeta: 阻尼比
返回:
alpha: 地震影响系数
"""
# 阻尼调整系数
gamma = 0.9 + (0.05 - zeta) / (0.3 + 6 * zeta)
eta2 = 1 + (0.05 - zeta) / (0.08 + 1.6 * zeta)
eta1 = 0.02 + (0.05 - zeta) / (4 + 32 * zeta)
# 确保系数在合理范围
gamma = max(0.5, min(gamma, 1.0))
eta2 = max(0.5, min(eta2, 1.0))
eta1 = max(0.0, eta1)
alpha = np.zeros_like(T)
for i, t in enumerate(T):
if t < 0.1:
# 上升段
alpha[i] = (0.45 + 5.5 * t) * alpha_max
elif t <= Tg:
# 平台段
alpha[i] = alpha_max
elif t <= 5 * Tg:
# 下降段
alpha[i] = (Tg / t) ** gamma * alpha_max
elif t <= 6.0:
# 长周期段
alpha[i] = (eta2 * 0.2 ** gamma - eta1 * (t - 5 * Tg)) * alpha_max
else:
alpha[i] = 0
return alpha
# ============================================================
# 生成设计反应谱
# ============================================================
print("\n生成设计反应谱...")
print("-" * 70)
T = np.linspace(0.01, 6.0, 200)
# 不同设防烈度的设计反应谱
fig1, axes = plt.subplots(1, 2, figsize=(14, 6))
colors_intensity = ['steelblue', 'forestgreen', 'coral', 'purple', 'brown', 'crimson']
Tg = 0.35 # II类场地
for i, (intensity, alpha_max) in enumerate(alpha_max_table.items()):
alpha = compute_design_spectrum(T, alpha_max, Tg)
axes[0].plot(T, alpha, color=colors_intensity[i], linewidth=2, label=f'{
intensity}')
axes[0].set_xlabel('Period (s)', fontsize=12)
axes[0].set_ylabel('Seismic Influence Coefficient α', fontsize=12)
axes[0].set_title('Design Response Spectrum - Different Intensities (Site II)', fontsize=13, fontweight='bold')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 6)
# 不同场地类别的设计反应谱
alpha_max = 0.16 # 8度设防
colors_site = ['steelblue', 'forestgreen', 'coral', 'purple', 'brown']
for i, (site, Tg_val) in enumerate(Tg_table.items()):
alpha = compute_design_spectrum(T, alpha_max, Tg_val)
axes[1].plot(T, alpha, color=colors_site[i], linewidth=2, label=f'Site {
site}')
axes[1].set_xlabel('Period (s)', fontsize=12)
axes[1].set_ylabel('Seismic Influence Coefficient α', fontsize=12)
axes[1].set_title('Design Response Spectrum - Different Site Classes (Intensity 8)', fontsize=13, fontweight='bold')
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 6)
plt.tight_layout()
plt.savefig(f'{
output_dir}/design_spectrum_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 不同设防烈度设计反应谱已保存")
# 不同阻尼比的设计反应谱
fig2, ax = plt.subplots(figsize=(10, 6))
alpha_max = 0.16
Tg = 0.35
zeta_values = [0.02, 0.05, 0.10, 0.20]
colors_zeta = ['steelblue', 'forestgreen', 'coral', 'purple']
for i, zeta in enumerate(zeta_values):
alpha = compute_design_spectrum(T, alpha_max, Tg, zeta)
ax.plot(T, alpha, color=colors_zeta[i], linewidth=2, label=f'ζ={
zeta}')
ax.set_xlabel('Period (s)', fontsize=12)
ax.set_ylabel('Seismic Influence Coefficient α', fontsize=12)
ax.set_title('Design Response Spectrum - Different Damping Ratios', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 6)
plt.tight_layout()
plt.savefig(f'{
output_dir}/design_spectrum_damping.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 不同阻尼比设计反应谱已保存")
# 综合对比图
fig3, axes = plt.subplots(2, 2, figsize=(14, 10))
# 子图1：不同设防烈度
for i, (intensity, alpha_max) in enumerate(alpha_max_table.items()):
alpha = compute_design_spectrum(T, alpha_max, Tg)
axes[0, 0].plot(T, alpha, color=colors_intensity[i], linewidth=2, label=f'{
intensity}')
axes[0, 0].set_title('Different Seismic Intensities', fontsize=12, fontweight='bold')
axes[0, 0].set_xlabel('Period (s)')
axes[0, 0].set_ylabel('α')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
# 子图2：不同场地类别
for i, (site, Tg_val) in enumerate(Tg_table.items()):
alpha = compute_design_spectrum(T, alpha_max, Tg_val)
axes[0, 1].plot(T, alpha, color=colors_site[i], linewidth=2, label=f'Site {
site}')
axes[0, 1].set_title('Different Site Classes', fontsize=12, fontweight='bold')
axes[0, 1].set_xlabel('Period (s)')
axes[0, 1].set_ylabel('α')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)
# 子图3：不同阻尼比
for i, zeta in enumerate(zeta_values):
alpha = compute_design_spectrum(T, alpha_max, Tg, zeta)
axes[1, 0].plot(T, alpha, color=colors_zeta[i], linewidth=2, label=f'ζ={
zeta}')
axes[1, 0].set_title('Different Damping Ratios', fontsize=12, fontweight='bold')
axes[1, 0].set_xlabel('Period (s)')
axes[1, 0].set_ylabel('α')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)
# 子图4：规范反应谱四段特征
alpha = compute_design_spectrum(T, alpha_max, Tg)
axes[1, 1].plot(T, alpha, 'b-', linewidth=2.5)
axes[1, 1].axvline(x=0.1, color='r', linestyle='--', alpha=0.7, label='Rising end')
axes[1, 1].axvline(x=Tg, color='g', linestyle='--', alpha=0.7, label=f'Tg={
Tg}s')
axes[1, 1].axvline(x=5*Tg, color='m', linestyle='--', alpha=0.7, label=f'5Tg={
5*Tg}s')
axes[1, 1].set_title('Four Segments of Design Spectrum', fontsize=12, fontweight='bold')
axes[1, 1].set_xlabel('Period (s)')
axes[1, 1].set_ylabel('α')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{
output_dir}/design_spectrum_comprehensive.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 综合对比图已保存")
# ============================================================
# 生成设计反应谱动画
# ============================================================
print("\n正在生成设计反应谱动画...")
fig4, ax = plt.subplots(figsize=(10, 6))
alpha = compute_design_spectrum(T, alpha_max, Tg)
line, = ax.plot([], [], 'b-', linewidth=2.5)
ax.set_xlim(0, 6)
ax.set_ylim(0, alpha_max * 1.1)
ax.set_xlabel('Period (s)', fontsize=12)
ax.set_ylabel('Seismic Influence Coefficient α', fontsize=12)
ax.set_title(f'Design Response Spectrum Animation (Intensity 8, Site II)', fontsize=13, fontweight='bold')
ax.grid(True, alpha=0.3)
# 添加特征周期标记
ax.axvline(x=Tg, color='r', linestyle='--', alpha=0.7, label=f'Characteristic Period Tg={
Tg}s')
ax.legend(fontsize=10)
from matplotlib.animation import FuncAnimation
def update(frame):
idx = frame + 1
if idx > len(T):
idx = len(T)
line.set_data(T[:idx], alpha[:idx])
return line,
n_frames = len(T)
anim = FuncAnimation(fig4, update, frames=n_frames, interval=30, blit=True)
anim.save(f'{
output_dir}/design_spectrum_animation.gif', writer='pillow', fps=30, dpi=100)
plt.close()
print(" 设计反应谱动画已保存")
# ============================================================
# 打印关键结果
# ============================================================
print("\n设计反应谱关键结果：")
print("-" * 70)
print(" 设防烈度与地震影响系数最大值：")
for intensity, alpha_max in alpha_max_table.items():
print(f" {
intensity}: αmax = {
alpha_max}")
print(f"\n 场地类别与特征周期（第一组）：")
for site, Tg in Tg_table.items():
print(f" {
site}类: Tg = {
Tg} s")
print("\n" + "=" * 70)
print("设计反应谱计算完成！")
print("=" * 70)
6.4 案例4：多自由度体系地震响应分析
本案例实现多自由度体系的地震响应分析，包括振型分解法和底部剪力法。
# -*- coding: utf-8 -*-
"""
主题017：结构地震响应分析与设计反应谱 - 案例4
多自由度体系地震响应分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.linalg import eigh
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题017'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("=" * 70)
print("主题017：结构地震响应分析与设计反应谱")
print("案例4：多自由度体系地震响应分析")
print("=" * 70)
# ============================================================
# 生成地震波
# ============================================================
def generate_earthquake(dt, t_end, omega_g=5*np.pi, zeta_g=0.6, S0=0.02, seed=None):
"""生成地震波"""
if seed is not None:
np.random.seed(seed)
t = np.arange(0, t_end, dt)
n = len(t)
white_noise = np.random.randn(n)
freqs = np.fft.fftfreq(n, dt)
omega = 2 * np.pi * freqs
omega_sq = omega**2
numerator = omega_g**4 + 4*zeta_g**2*omega_g**2*omega_sq
denominator = (omega_g**2 - omega_sq)**2 + 4*zeta_g**2*omega_g**2*omega_sq
psd = np.zeros_like(omega)
mask = denominator > 1e-10
psd[mask] = S0 * numerator[mask] / denominator[mask]
white_fft = np.fft.fft(white_noise)
filtered_fft = white_fft * np.sqrt(psd * n * dt / 2)
filtered = np.real(np.fft.ifft(filtered_fft))
t_mid = t_end / 2
envelope = np.zeros(n)
for i, ti in enumerate(t):
if ti < t_mid:
envelope[i] = (ti / t_mid)**2
else:
envelope[i] = np.exp(-2 * (ti - t_mid) / t_mid)
acc = filtered * envelope
target_pga = 0.2 * 9.8
current_pga = np.max(np.abs(acc))
if current_pga > 0:
acc = acc * target_pga / current_pga
return t, acc
# ============================================================
# 定义多自由度体系
# ============================================================
print("\n定义多自由度体系参数：")
print("-" * 70)
# 5层剪切型框架
n_dof = 5
m = 1000.0 # 每层质量 (kg)
k = 50000.0 # 每层刚度 (N/m)
zeta = 0.05 # 阻尼比
# 质量矩阵
M = np.eye(n_dof) * m
# 刚度矩阵（三对角矩阵）
K = np.zeros((n_dof, n_dof))
for i in range(n_dof):
if i == 0:
K[i, i] = 2 * k
K[i, i+1] = -k
elif i == n_dof - 1:
K[i, i] = k
K[i, i-1] = -k
else:
K[i, i] = 2 * k
K[i, i-1] = -k
K[i, i+1] = -k
print(f" 层数: {
n_dof}")
print(f" 每层质量: {
m} kg")
print(f" 每层刚度: {
k} N/m")
print(f" 阻尼比: {
zeta}")
# ============================================================
# 模态分析
# ============================================================
print("\n进行模态分析...")
print("-" * 70)
# 求解特征值问题
eigenvalues, eigenvectors = eigh(K, M)
# 固有频率
omega_n = np.sqrt(eigenvalues)
f_n = omega_n / (2 * np.pi)
T_n = 1 / f_n
print(" 模态频率和周期：")
for i in range(n_dof):
print(f" 模态{
i+1}: f = {
f_n[i]:.4f} Hz, T = {
T_n[i]:.4f} s")
# 振型归一化（质量归一化）
phi = np.zeros_like(eigenvectors)
for i in range(n_dof):
phi[:, i] = eigenvectors[:, i] / np.sqrt(eigenvectors[:, i].T @ M @ eigenvectors[:, i])
# 振型参与系数
gamma = np.zeros(n_dof)
for i in range(n_dof):
gamma[i] = phi[:, i].T @ M @ np.ones(n_dof)
print("\n 振型参与系数：")
for i in range(n_dof):
print(f" 模态{
i+1}: γ = {
gamma[i]:.4f}")
# 参与质量比
participation_mass = gamma**2 / np.sum(gamma**2) * 100
print("\n 参与质量比：")
for i in range(n_dof):
print(f" 模态{
i+1}: {
participation_mass[i]:.2f}%")
# ============================================================
# 生成地震波
# ============================================================
print("\n生成地震波...")
print("-" * 70)
dt = 0.01
t_end = 30.0
t, acc = generate_earthquake(dt, t_end, seed=42)
print(f" 时间步长: {
dt} s")
print(f" 持续时间: {
t_end} s")
print(f" PGA: {
np.max(np.abs(acc)):.4f} m/s²")
# ============================================================
# 振型分解法计算地震响应
# ============================================================
print("\n使用振型分解法计算地震响应...")
print("-" * 70)
# Newmark-beta参数
gamma_nb = 0.5
beta_nb = 0.25
# 计算各模态响应
n_steps = len(t)
modal_disp = np.zeros((n_dof, n_steps))
modal_vel = np.zeros((n_dof, n_steps))
modal_acc = np.zeros((n_dof, n_steps))
for mode in range(n_dof):
omega = omega_n[mode]
zeta_mode = zeta
gamma_mode = gamma[mode]
# 广义坐标方程
q = np.zeros(n_steps)
q_dot = np.zeros(n_steps)
q_ddot = np.zeros(n_steps)
# 等效刚度
k_eff = omega**2 + gamma_nb / (beta_nb * dt) * 2 * zeta_mode * omega + 1 / (beta_nb * dt**2)
for i in range(n_steps - 1):
# 等效荷载
p_eff = -gamma_mode * acc[i+1] + (1/(beta_nb*dt**2))*q[i] + (1/(beta_nb*dt))*q_dot[i] + (1/(2*beta_nb)-1)*q_ddot[i]
# 计算广义位移
q[i+1] = p_eff / k_eff
# 计算广义速度和加速度
q_dot[i+1] = gamma_nb/(beta_nb*dt)*(q[i+1]-q[i]) + (1-gamma_nb/beta_nb)*q_dot[i] + dt*(1-gamma_nb/(2*beta_nb))*q_ddot[i]
q_ddot[i+1] = 1/(beta_nb*dt**2)*(q[i+1]-q[i]) - 1/(beta_nb*dt)*q_dot[i] - (1/(2*beta_nb)-1)*q_ddot[i]
modal_disp[mode, :] = q
modal_vel[mode, :] = q_dot
modal_acc[mode, :] = q_ddot
# SRSS法组合各模态响应
n_modes_to_use = 3 # 使用前3阶模态
disp = np.zeros((n_dof, n_steps))
vel = np.zeros((n_dof, n_steps))
acc_total = np.zeros((n_dof, n_steps))
for i in range(n_dof):
for j in range(n_steps):
# SRSS组合
disp_sq = sum((phi[i, mode] * modal_disp[mode, j])**2 for mode in range(n_modes_to_use))
vel_sq = sum((phi[i, mode] * modal_vel[mode, j])**2 for mode in range(n_modes_to_use))
acc_sq = sum((phi[i, mode] * modal_acc[mode, j])**2 for mode in range(n_modes_to_use))
disp[i, j] = np.sqrt(disp_sq)
vel[i, j] = np.sqrt(vel_sq)
acc_total[i, j] = np.sqrt(acc_sq)
print(f" 使用{
n_modes_to_use}阶模态进行SRSS组合")
# ============================================================
# 绘制振型
# ============================================================
print("\n正在生成可视化图表...")
fig1, axes = plt.subplots(1, 2, figsize=(14, 6))
# 振型图
for i in range(n_dof):
mode_shape = phi[:, i]
# 归一化到最大值为1
mode_shape = mode_shape / np.max(np.abs(mode_shape))
axes[0].plot(mode_shape, np.arange(n_dof)+1, 'o-', linewidth=2, markersize=8, label=f'Mode {
i+1}')
axes[0].set_xlabel('Normalized Mode Shape', fontsize=12)
axes[0].set_ylabel('Story', fontsize=12)
axes[0].set_title('Mode Shapes', fontsize=14, fontweight='bold')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0.5, n_dof+0.5)
axes[0].invert_yaxis()
# 参与质量比
axes[1].bar(range(1, n_dof+1), participation_mass, color='steelblue', alpha=0.7)
axes[1].set_xlabel('Mode Number', fontsize=12)
axes[1].set_ylabel('Participation Mass Ratio (%)', fontsize=12)
axes[1].set_title('Modal Participation Mass Ratio', fontsize=14, fontweight='bold')
axes[1].grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig(f'{
output_dir}/modal_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 模态分析图已保存")
# ============================================================
# 绘制楼层响应时程
# ============================================================
fig2, axes = plt.subplots(3, 1, figsize=(14, 12))
# 选择几个楼层显示
display_floors = [0, 2, 4] # 第1、3、5层
colors = ['steelblue', 'forestgreen', 'coral']
for i, floor in enumerate(display_floors):
axes[0].plot(t, disp[floor, :]*1000, color=colors[i], linewidth=1, label=f'Floor {
floor+1}')
axes[1].plot(t, vel[floor, :], color=colors[i], linewidth=1, label=f'Floor {
floor+1}')
axes[2].plot(t, acc_total[floor, :], color=colors[i], linewidth=1, label=f'Floor {
floor+1}')
axes[0].set_ylabel('Displacement (mm)', fontsize=12)
axes[0].set_title('Floor Displacement Time History', fontsize=14, fontweight='bold')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[1].set_ylabel('Velocity (m/s)', fontsize=12)
axes[1].set_title('Floor Velocity Time History', fontsize=14, fontweight='bold')
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[2].set_xlabel('Time (s)', fontsize=12)
axes[2].set_ylabel('Acceleration (m/s²)', fontsize=12)
axes[2].set_title('Floor Acceleration Time History', fontsize=14, fontweight='bold')
axes[2].legend(fontsize=11)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{
output_dir}/floor_response.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 楼层响应时程图已保存")
# ============================================================
# 绘制楼层最大响应包络
# ============================================================
fig3, axes = plt.subplots(1, 3, figsize=(15, 5))
floors = np.arange(1, n_dof+1)
max_disp = np.max(np.abs(disp), axis=1) * 1000 # mm
max_vel = np.max(np.abs(vel), axis=1) # m/s
max_acc = np.max(np.abs(acc_total), axis=1) # m/s²
axes[0].plot(max_disp, floors, 'o-', linewidth=2, markersize=8, color='steelblue')
axes[0].set_xlabel('Max Displacement (mm)', fontsize=11)
axes[0].set_ylabel('Floor', fontsize=11)
axes[0].set_title('Maximum Floor Displacement', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0.5, n_dof+0.5)
axes[0].invert_yaxis()
axes[1].plot(max_vel, floors, 'o-', linewidth=2, markersize=8, color='forestgreen')
axes[1].set_xlabel('Max Velocity (m/s)', fontsize=11)
axes[1].set_ylabel('Floor', fontsize=11)
axes[1].set_title('Maximum Floor Velocity', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(0.5, n_dof+0.5)
axes[1].invert_yaxis()
axes[2].plot(max_acc, floors, 'o-', linewidth=2, markersize=8, color='coral')
axes[2].set_xlabel('Max Acceleration (m/s²)', fontsize=11)
axes[2].set_ylabel('Floor', fontsize=11)
axes[2].set_title('Maximum Floor Acceleration', fontsize=12, fontweight='bold')
axes[2].grid(True, alpha=0.3)
axes[2].set_ylim(0.5, n_dof+0.5)
axes[2].invert_yaxis()
plt.tight_layout()
plt.savefig(f'{
output_dir}/max_response_envelope.png', dpi=150, bbox_inches='tight')
plt.close()
print(" 最大响应包络图已保存")
# ============================================================
# 生成结构地震响应动画
# ============================================================
print("\n正在生成结构地震响应动画...")
fig4, ax = plt.subplots(figsize=(10, 10))
# 绘制结构轮廓
def draw_building(ax, displacements, floor_height=3.0):
"""绘制建筑结构"""
ax.clear()
# 楼层高度
heights = np.arange(n_dof) * floor_height
# 建筑宽度
width = 6.0
# 绘制变形后的结构
for i in range(n_dof):
x = displacements[i] * 1000 # 转换为mm
y = heights[i]
# 绘制楼层
rect = plt.Rectangle((x - width/2, y), width, floor_height*0.8,
linewidth=1, edgecolor='black', facecolor='lightblue', alpha=0.7)
ax.add_patch(rect)
# 标注楼层号
ax.text(x + width/2 + 0.5, y + floor_height*0.4, f'{
i+1}F', fontsize=10)
# 绘制柱子
for i in range(n_dof):
if i == 0:
# 底层柱
x_left = displacements[i] * 1000 - width/2
x_right = displacements[i] * 1000 + width/2
y_bottom = 0
y_top = heights[i]
ax.plot([x_left, x_left], [y_bottom, y_top], 'k-', linewidth=2)
ax.plot([x_right, x_right], [y_bottom, y_top], 'k-', linewidth=2)
else:
# 上层柱
x_left_bottom = displacements[i-1] * 1000 - width/2
x_left_top = displacements[i] * 1000 - width/2
x_right_bottom = displacements[i-1] * 1000 + width/2
x_right_top = displacements[i] * 1000 + width/2
y_bottom = heights[i-1] + floor_height*0.8
y_top = heights[i]
ax.plot([x_left_bottom, x_left_top], [y_bottom, y_top], 'k-', linewidth=2)
ax.plot([x_right_bottom, x_right_top], [y_bottom, y_top], 'k-', linewidth=2)
ax.set_xlim(-20, 20)
ax.set_ylim(-1, n_dof * floor_height + 1)
ax.set_xlabel('Displacement (mm)', fontsize=12)
ax.set_ylabel('Height (m)', fontsize=12)
ax.set_title('Building Seismic Response Animation', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
from matplotlib.animation import FuncAnimation
def update(frame):
idx = frame * 5
if idx >= n_steps:
idx = n_steps - 1
# 获取当前时刻的位移
current_disp = disp[:, idx]
draw_building(ax, current_disp)
return []
n_frames = n_steps // 5
anim = FuncAnimation(fig4, update, frames=n_frames, interval=50, blit=True)
anim.save(f'{
output_dir}/building_response_animation.gif', writer='pillow', fps=20, dpi=100)
plt.close()
print(" 结构地震响应动画已保存")
# ============================================================
# 底部剪力法对比
# ============================================================
print("\n底部剪力法计算...")
print("-" * 70)
# 计算底部剪力
G = np.sum(np.diag(M)) * 9.8 # 总重力荷载
alpha_1 = 0.16 # 假设8度设防，II类场地
F_ek = alpha_1 * G * 0.85 # 等效总重力荷载系数0.85
print(f" 结构总质量: {
np.sum(np.diag(M)):.2f} kg")
print(f" 总重力荷载: {
G:.2f} N")
print(f" 水平地震影响系数: α1 = {
alpha_1}")
print(f" 底部总剪力: F_Ek = {
F_ek:.2f} N")
# 各层地震作用
heights = np.arange(1, n_dof+1) * 3.0 # 假设层高3m
F_i = np.zeros(n_dof)
for i in range(n_dof):
F_i[i] = np.diag(M)[i] * 9.8 * heights[i] / np.sum(np.diag(M) * 9.8 * heights) * F_ek
print(f"\n 各层水平地震作用：")
for i in range(n_dof):
print(f" 第{
i+1}层: F_{
i+1} = {
F_i[i]:.2f} N")
# ============================================================
# 打印关键结果
# ============================================================
print("\n地震响应关键结果：")
print("-" * 70)
print(" 楼层最大位移：")
for i in range(n_dof):
print(f" 第{
i+1}层: {
max_disp[i]:.4f} mm")
print("\n 楼层最大加速度：")
for i in range(n_dof):
print(f" 第{
i+1}层: {
max_acc[i]:.4f} m/s²")
print("\n 顶层最大位移: {:.4f} mm".format(max_disp[-1]))
print(" 顶层最大加速度: {:.4f} m/s²".format(max_acc[-1]))
print("\n" + "=" * 70)
print("多自由度体系地震响应分析完成！")
print("=" * 70)
7. 结果分析与讨论
7.1 地震波特性分析
通过案例1的地震波生成和频谱分析，可以观察到：
-
场地条件对地震波的影响：
- 硬土场地的高频成分更丰富，PGA较大
- 软土场地的低频成分更显著，持时更长
- 中硬土场地的频谱特性介于两者之间
-
包络函数的作用：
- 模拟地震动的非平稳特性
- 强震段出现在地震波的中部
- 符合实际地震动的能量分布规律
7.2 反应谱特性分析
通过案例2的反应谱计算，可以得出以下结论：
-
阻尼比的影响：
- 阻尼比越大，反应谱值越小
- 高阻尼对长周期结构的响应减小更明显
- 工程设计中通常取ζ=0.05作为标准值
-
反应谱的形状特征：
- 短周期段：加速度反应谱接近常数
- 中周期段：速度反应谱接近常数
- 长周期段：位移反应谱接近常数
-
伪反应谱的适用性：
- 对于小阻尼体系，伪反应谱与真实反应谱非常接近
- 伪加速度反应谱便于计算地震作用
- 伪速度反应谱与能量输入相关
7.3 设计反应谱的工程意义
通过案例3的设计反应谱分析，可以理解：
-
设防烈度的影响：
- 设防烈度越高，地震影响系数越大
- 9度设防的地震作用约为6度设防的8倍
- 设计反应谱的形状与设防烈度无关
-
场地类别的影响：
- 软土场地的特征周期较长
- 长周期结构在软土场地上响应更大
- 场地条件是抗震设计的重要考虑因素
-
阻尼比的影响：
- 阻尼比越大，设计反应谱值越小
- 阻尼调整系数反映了这一影响
- 增加结构阻尼是有效的减震措施
7.4 多自由度体系地震响应
通过案例4的多自由度体系分析，可以得出：
-
振型参与特性：
- 第一阶振型的参与质量比最大（约88%）
- 前3阶振型累计参与质量超过99%
- 高阶振型对总响应的贡献较小
-
楼层响应分布：
- 楼层位移随高度增加而增大
- 顶层位移最大，底层剪力最大
- 加速度响应在各楼层分布相对均匀
-
设计方法对比：
- 振型分解法结果精确，但计算复杂
- 底部剪力法简便，适用于规则结构
- 两种方法的结果具有可比性
8. 工程应用与展望
8.1 抗震设计方法
基于反应谱理论的抗震设计方法包括：
- 底部剪力法：适用于高度不超过40m、质量和刚度沿高度分布均匀的规则结构
- 振型分解反应谱法：适用于大多数建筑结构，考虑多振型的贡献
- 时程分析法：用于特别不规则结构、超高层结构或重要建筑
8.2 减震隔震技术
现代抗震工程的发展趋势：
- 基础隔震：通过隔震支座延长结构周期，避开地震能量集中的频段
- 消能减震：利用阻尼器耗散地震能量，减小结构响应
- 调谐质量阻尼器（TMD）：通过调谐质量系统吸收主结构的振动能量
8.3 基于性能的抗震设计
新一代抗震设计理念的核心理念：
-
多级设防目标：
- 小震不坏：结构保持弹性
- 中震可修：结构可修复
- 大震不倒：保证生命安全
-
性能化设计：根据建筑的重要性确定性能目标
-
位移控制设计：以层间位移角作为控制指标
8.4 未来发展方向
地震工程领域的研究前沿：
- 地震动预测：基于地震学原理的强地面运动预测
- 结构健康监测：实时监测结构状态，评估震后损伤
- 韧性抗震设计：不仅保证生命安全，还要实现震后快速恢复
- 人工智能应用：机器学习在地震响应预测和损伤识别中的应用
9. 结论
本教程系统地介绍了结构地震响应分析与设计反应谱的理论和方法，通过Python实现了地震波生成、反应谱计算、设计反应谱构建以及多自由度体系地震响应分析。主要结论如下：
- 地震动特性：地震动具有显著的随机性和非平稳性，Kanai-Tajimi模型可以有效描述其统计特性
- 反应谱理论：反应谱是连接地震动特性和结构响应的桥梁，是抗震设计的核心工具
- 设计反应谱：基于统计分析和规范要求的设计反应谱，为工程抗震设计提供了科学依据
- 多自由度体系：振型分解法是分析多自由度体系地震响应的有效方法，SRSS组合适用于大多数工程结构
通过本教程的学习，读者应该能够：
- 理解地震动的基本特性和数学描述
- 掌握反应谱的计算方法和物理意义
- 应用设计反应谱进行结构抗震设计
- 使用Python进行地震响应分析和可视化