结构动力学仿真-主题015-结构动力响应分析方法

kkchenkx

118人浏览 · 2026-04-03 21:04:01

kkchenkx · 2026-04-03 21:04:01 发布

第十五篇：结构动力响应分析方法

摘要

结构动力响应分析是结构动力学的核心内容，旨在计算结构在各种动力荷载作用下的响应。本教程系统地介绍结构动力响应分析的各种方法，包括时域分析、频域分析、模态分析和随机响应分析等。首先阐述动力响应分析的基本概念和数学模型；然后详细介绍各种分析方法的原理和应用范围；接着讲解数值计算方法，如直接积分法、模态叠加法等；最后通过Python实现多个典型案例，包括时程分析、频率响应分析、模态分析和随机振动分析等。通过本教程的学习，读者将掌握结构动力响应分析的核心方法和实用技术。

关键词

动力响应分析，时域分析，频域分析，模态分析，随机振动，直接积分法，模态叠加法，Python仿真

在这里插入图片描述

1. 引言

1.1 动力响应分析的重要性

结构在使用过程中会受到各种动力荷载的作用，如风荷载、地震作用、车辆荷载、设备振动等。这些动力荷载会引起结构的振动响应，可能导致结构损坏、影响使用功能或降低舒适性。因此，结构动力响应分析是结构设计、安全评估和性能预测的重要手段。

动力响应分析的主要目的：

安全性评估：验证结构在动力荷载作用下是否满足强度和稳定性要求
适用性评估：确保结构在正常使用条件下的振动水平满足舒适度要求
优化设计：通过分析不同设计方案的动力响应，选择最优设计
损伤识别：通过响应分析识别结构的损伤位置和程度
振动控制：为振动控制策略的设计提供依据

1.2 动力响应分析方法分类

结构动力响应分析方法可以从不同角度进行分类：

按分析域分类：

时域分析：直接在时间域内求解运动方程，适用于任意荷载和非线性系统
频域分析：通过傅里叶变换将问题转换到频率域，适用于线性系统和周期性荷载

按求解方法分类：

直接积分法：直接数值积分求解运动方程，如Newmark-β法、Wilson-θ法等
模态叠加法：利用结构的模态特性，将多自由度问题转化为单自由度问题
频率响应法：通过频率域分析求解系统的频率响应函数

按荷载特性分类：

确定性分析：针对确定性荷载（如正弦荷载、脉冲荷载）的响应分析
随机分析：针对随机荷载（如随机风荷载、地震动）的响应分析

按系统特性分类：

线性分析：假设结构为线性系统，适用于小变形情况
非线性分析：考虑结构的非线性特性，适用于大变形或材料非线性情况

1.3 动力响应分析的数学模型

结构动力响应分析的基础是建立结构的运动方程。对于多自由度线性系统，运动方程为：

$Mx¨+Cx˙+Kx=F(t)\mathbf{M} \ddot{\mathbf{x}} + \mathbf{C} \dot{\mathbf{x}} + \mathbf{K} \mathbf{x} = \mathbf{F}(t)$

其中：

$M\mathbf{M}$ 为质量矩阵
$C\mathbf{C}$ 为阻尼矩阵
$K\mathbf{K}$ 为刚度矩阵
$x\mathbf{x}$ 为位移向量
$F(t)\mathbf{F}(t)$ 为荷载向量

对于非线性系统，运动方程为：

$Mx¨+C(x,x˙)x˙+K(x)x=F(t)\mathbf{M} \ddot{\mathbf{x}} + \mathbf{C}(\mathbf{x}, \dot{\mathbf{x}}) \dot{\mathbf{x}} + \mathbf{K}(\mathbf{x}) \mathbf{x} = \mathbf{F}(t)$

2. 时域分析方法

2.1 直接积分法

直接积分法是通过数值积分直接求解运动方程的方法，适用于任意荷载和非线性系统。常用的直接积分法包括：

中心差分法：

基于泰勒展开，使用前一时间步的位移和速度计算当前时间步的加速度
显式方法，不需要解线性方程组
时间步长受稳定性限制

Newmark-β法：

隐式方法，使用当前时间步的位移和速度计算加速度
无条件稳定（当 $β≥1/4\beta \geq 1/4$ 且 $γ≥1/2\gamma \geq 1/2$ ）
广泛应用于工程分析

Wilson-θ法：

基于Newmark法，引入θ参数
无条件稳定（当 $θ≥1.37\theta \geq 1.37$ ）
具有良好的数值稳定性

HHT-α法：

引入α参数控制算法的数值耗散
适用于高频内容丰富的问题

2.2 Newmark-β法的数学推导

Newmark-β法的基本假设为：

$x˙t+Δt=x˙t+Δt[(1−γ)x¨t+γx¨t+Δt]\dot{\mathbf{x}}_{t+\Delta t} = \dot{\mathbf{x}}_t + \Delta t \left[ (1-\gamma) \ddot{\mathbf{x}}_t + \gamma \ddot{\mathbf{x}}_{t+\Delta t} \right]$

$xt+Δt=xt+Δtx˙t+Δt2[(1/2−β)x¨t+βx¨t+Δt]\mathbf{x}_{t+\Delta t} = \mathbf{x}_t + \Delta t \dot{\mathbf{x}}_t + \Delta t^2 \left[ (1/2-\beta) \ddot{\mathbf{x}}_t + \beta \ddot{\mathbf{x}}_{t+\Delta t} \right]$

其中， $γ\gamma$ 和 $β\beta$ 为参数，通常取 $γ=1/2\gamma = 1/2$ ， $β=1/4\beta = 1/4$ （平均加速度法）。

将上述表达式代入运动方程，得到：

$(M+γβΔt2K+γΔtC)x¨t+Δt=Ft+Δt−K(xt+Δtx˙t+Δt2(1/2−β)x¨t)−C(x˙t+Δt(1−γ)x¨t)−Mx¨t\left( \mathbf{M} + \gamma \beta \Delta t^2 \mathbf{K} + \gamma \Delta t \mathbf{C} \right) \ddot{\mathbf{x}}_{t+\Delta t} = \mathbf{F}_{t+\Delta t} - \mathbf{K} \left( \mathbf{x}_t + \Delta t \dot{\mathbf{x}}_t + \Delta t^2 (1/2-\beta) \ddot{\mathbf{x}}_t \right) - \mathbf{C} \left( \dot{\mathbf{x}}_t + \Delta t (1-\gamma) \ddot{\mathbf{x}}_t \right) - \mathbf{M} \ddot{\mathbf{x}}_t$

2.3 直接积分法的稳定性分析

数值稳定性：

显式方法（如中心差分法）：有条件稳定，时间步长必须小于临界时间步长
隐式方法（如Newmark-β法）：可无条件稳定，时间步长仅受精度影响

临界时间步长：
对于中心差分法，临界时间步长为：

$Δtcr=2ωmax\Delta t_{cr} = \frac{2}{\omega_{max}}$

其中， $ωmax\omega_{max}$ 为系统的最大固有频率。

3. 频域分析方法

3.1 傅里叶变换

频域分析的基础是傅里叶变换，将时域信号转换为频域表示：

$F(ω)=∫−∞∞f(t)e−iωtdtF(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt$

逆变换为：

$\frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} d\omega$

3.2 频率响应函数

对于线性系统，输入和输出之间的关系可以用频率响应函数（FRF）表示：

$H(ω)=X(ω)F(ω)H(\omega) = \frac{X(\omega)}{F(\omega)}$

其中， $X(ω)X(\omega)$ 为输出的傅里叶变换， $F(ω)F(\omega)$ 为输入的傅里叶变换。

对于单自由度系统，频率响应函数为：

$H(ω)=1−mω2+icω+kH(\omega) = \frac{1}{-m\omega^2 + ic\omega + k}$

对于多自由度系统，频率响应函数矩阵为：

$H(ω)=(−ω2M+iωC+K)−1\mathbf{H}(\omega) = \left( -\omega^2 \mathbf{M} + i\omega \mathbf{C} + \mathbf{K} \right)^{-1}$

3.3 频域分析的应用

频率响应分析：

计算系统在不同频率下的响应幅值和相位
识别系统的共振频率和阻尼特性
评估系统对谐波荷载的响应

随机振动分析：

通过功率谱密度（PSD）分析随机荷载的响应
计算响应的统计特性（均值、方差、概率分布）

4. 模态分析方法

4.1 模态叠加法

模态叠加法利用结构的模态特性，将多自由度问题转化为单自由度问题：

模态分析：求解结构的固有频率和模态振型
坐标变换：将物理坐标转换为模态坐标
单自由度求解：在模态坐标下求解每个模态的响应
叠加：将各模态的响应转换回物理坐标并叠加

4.2 模态分析的数学基础

结构的自由振动方程为：

$Mx¨+Kx=0\mathbf{M} \ddot{\mathbf{x}} + \mathbf{K} \mathbf{x} = 0$

假设解的形式为 $x=ϕeiωt\mathbf{x} = \boldsymbol{\phi} e^{i\omega t}$ ，代入得到特征值问题：

$(K−ω2M)ϕ=0\left( \mathbf{K} - \omega^2 \mathbf{M} \right) \boldsymbol{\phi} = 0$

求解得到固有频率 $ω1,ω2,…,ωn\omega_1, \omega_2, \dots, \omega_n$ 和对应的模态振型 $ϕ1,ϕ2,…,ϕn\boldsymbol{\phi}_1, \boldsymbol{\phi}_2, \dots, \boldsymbol{\phi}_n$ 。

4.3 模态叠加法的应用

线性系统：

适用于线性系统的动力响应分析
计算效率高，特别是对于大型系统
可以只考虑贡献较大的前几阶模态

非线性系统：

模态叠加法不适用于强非线性系统
对于弱非线性系统，可以采用非线性模态叠加法

5. 随机响应分析

5.1 随机过程基础

随机振动分析处理的是随机荷载作用下的结构响应。随机过程的基本特性包括：

均值：随机过程的平均水平
自相关函数：描述不同时刻随机变量之间的相关性
功率谱密度（PSD）：描述随机过程在频率域的能量分布

5.2 随机响应的计算方法

频域方法：

利用功率谱密度和频率响应函数计算响应的功率谱密度
适用于线性系统
计算效率高

时域方法：

生成随机荷载的时间历程样本
对每个样本进行时程分析
统计分析响应的统计特性
适用于非线性系统

5.3 响应统计特性的计算

对于线性系统，响应的统计特性可以通过以下公式计算：

响应功率谱密度： $SX(ω)=∣H(ω)∣2SF(ω)S_X(\omega) = |H(\omega)|^2 S_F(\omega)$
响应方差： $σX2=∫0∞SX(ω)dω\sigma_X^2 = \int_{0}^{\infty} S_X(\omega) d\omega$
响应均方根： $RMSX=σX2\text{RMS}_X = \sqrt{\sigma_X^2}$

6. Python实现案例

6.1 案例1：直接积分法（Newmark-β法）

本案例实现Newmark-β法求解结构的动力响应。

# -*- coding: utf-8 -*-
"""
主题015：结构动力响应分析方法 - 案例1
直接积分法（Newmark-β法）
"""
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仿真领域\工程仿真\结构动力学仿真\主题015'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("主题015：结构动力响应分析方法")
print("案例1：直接积分法（Newmark-β法）")
print("=" * 70)

# ============================================================
# 多自由度结构模型
# ============================================================
print("\n多自由度结构模型：")
print("-" * 70)

# 质量矩阵
M = np.array([[1000, 0, 0],
              [0, 1000, 0],
              [0, 0, 1000]])

# 刚度矩阵
K = np.array([[40000, -20000, 0],
              [-20000, 40000, -20000],
              [0, -20000, 20000]])

# 阻尼矩阵（瑞利阻尼）
alpha = 0.1
beta = 0.005
C = alpha * M + beta * K

print("  质量矩阵 M (kg):")
print(M)
print("\n  刚度矩阵 K (N/m):")
print(K)
print("\n  阻尼矩阵 C (N·s/m):")
print(C)

# 计算固有频率
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(M) @ K)
idx = np.argsort(eigenvalues)
omega_n = np.sqrt(eigenvalues[idx])
f_n = omega_n / (2 * np.pi)

print(f"\n固有频率：")
for i, (wn, fn) in enumerate(zip(omega_n, f_n)):
    print(f"  第{i+1}阶: ω = {wn:.4f} rad/s, f = {fn:.4f} Hz")

# ============================================================
# Newmark-β法
# ============================================================
def newmark_beta(M, C, K, F, x0, v0, a0, dt, t_end, gamma=0.5, beta=0.25):
    """
    Newmark-β法求解运动方程
    
    参数：
    M: 质量矩阵
    C: 阻尼矩阵
    K: 刚度矩阵
    F: 荷载时程函数
    x0: 初始位移
    v0: 初始速度
    a0: 初始加速度
    dt: 时间步长
    t_end: 结束时间
    gamma: Newmark参数
    beta: Newmark参数
    
    返回：
    t: 时间数组
    x: 位移时程
    v: 速度时程
    a: 加速度时程
    """
    n_dof = len(x0)
    n_steps = int(t_end / dt)
    
    t = np.linspace(0, t_end, n_steps + 1)
    x = np.zeros((n_steps + 1, n_dof))
    v = np.zeros((n_steps + 1, n_dof))
    a = np.zeros((n_steps + 1, n_dof))
    
    # 初始条件
    x[0] = x0
    v[0] = v0
    a[0] = a0
    
    # 构建有效刚度矩阵
    K_eff = K + gamma / (beta * dt) * C + 1 / (beta * dt**2) * M
    
    # 预计算常数
    a1 = 1 / (beta * dt**2)
    a2 = gamma / (beta * dt)
    a3 = 1 / (beta * dt)
    a4 = 1 / (2 * beta) - 1
    a5 = gamma / beta - 1
    a6 = dt * (gamma / (2 * beta) - 1)
    
    for i in range(n_steps):
        # 计算有效荷载
        F_t = F(t[i+1])
        F_eff = F_t + M @ (a1 * x[i] + a3 * v[i] + a4 * a[i]) + C @ (a2 * x[i] + a5 * v[i] + a6 * a[i])
        
        # 求解位移
        x[i+1] = np.linalg.solve(K_eff, F_eff)
        
        # 计算速度和加速度
        a[i+1] = a1 * (x[i+1] - x[i]) - a3 * v[i] - a4 * a[i]
        v[i+1] = v[i] + dt * ((1 - gamma) * a[i] + gamma * a[i+1])
    
    return t, x, v, a

# ============================================================
# 荷载函数
# ============================================================
def seismic_load(t):
    """地震荷载时程"""
    # 生成人工地震波
    if t < 0:
        return np.array([0, 0, 0])
    
    # 模拟El Centro地震波
    omega = 2 * np.pi * 1.0  # 主频1Hz
    envelope = np.exp(-0.1 * t) * np.sin(np.pi * t / 10)**2
    if t > 10:
        envelope = 0
    
    acc = 0.1 * np.sin(omega * t) + 0.05 * np.sin(3 * omega * t) + 0.02 * np.sin(5 * omega * t)
    acc *= envelope
    
    # 转换为地震力
    return -M @ np.array([acc, acc, acc])

# ============================================================
# 求解动力响应
# ============================================================
print("\n求解动力响应：")
print("-" * 70)

# 初始条件
x0 = np.zeros(3)
v0 = np.zeros(3)
a0 = np.linalg.solve(M, seismic_load(0) - C @ v0 - K @ x0)

# 时间参数
dt = 0.01
t_end = 20.0

print(f"  时间步长: {dt} s")
print(f"  计算时间: {t_end} s")
print(f"  总步数: {int(t_end/dt)} 步")

# 使用Newmark-β法求解
t, x, v, a = newmark_beta(M, C, K, seismic_load, x0, v0, a0, dt, t_end)

# 计算加速度响应（转换为m/s²）
a_accel = np.zeros_like(a)
for i in range(len(a)):
    a_accel[i] = a[i] + np.array([seismic_load(t[i])[0]/M[0,0], 
                                   seismic_load(t[i])[1]/M[1,1], 
                                   seismic_load(t[i])[2]/M[2,2]])

# ============================================================
# 结果分析
# ============================================================
print("\n响应结果分析：")
print("-" * 70)

for i in range(3):
    peak_disp = np.max(np.abs(x[:, i]))
    peak_vel = np.max(np.abs(v[:, i]))
    peak_acc = np.max(np.abs(a_accel[:, i]))
    
    print(f"\n第{i+1}层：")
    print(f"  峰值位移: {peak_disp*1000:.2f} mm")
    print(f"  峰值速度: {peak_vel:.2f} m/s")
    print(f"  峰值加速度: {peak_acc:.2f} m/s²")

# ============================================================
# 绘制结果
# ============================================================
print("\n正在生成可视化图表...")

fig1, axes = plt.subplots(3, 3, figsize=(16, 12))

for i in range(3):
    # 位移时程
    axes[i, 0].plot(t, x[:, i]*1000, 'b-', linewidth=1.5)
    axes[i, 0].set_xlabel('Time (s)', fontsize=11)
    axes[i, 0].set_ylabel('Displacement (mm)', fontsize=11)
    axes[i, 0].set_title(f'Floor {i+1} Displacement', fontsize=12, fontweight='bold')
    axes[i, 0].grid(True, alpha=0.3)
    
    # 速度时程
    axes[i, 1].plot(t, v[:, i], 'g-', linewidth=1.5)
    axes[i, 1].set_xlabel('Time (s)', fontsize=11)
    axes[i, 1].set_ylabel('Velocity (m/s)', fontsize=11)
    axes[i, 1].set_title(f'Floor {i+1} Velocity', fontsize=12, fontweight='bold')
    axes[i, 1].grid(True, alpha=0.3)
    
    # 加速度时程
    axes[i, 2].plot(t, a_accel[:, i], 'r-', linewidth=1.5)
    axes[i, 2].set_xlabel('Time (s)', fontsize=11)
    axes[i, 2].set_ylabel('Acceleration (m/s²)', fontsize=11)
    axes[i, 2].set_title(f'Floor {i+1} Acceleration', fontsize=12, fontweight='bold')
    axes[i, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/newmark_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  Newmark分析图已保存")

# 绘制层间位移
interstory_disp = np.zeros((len(t), 2))
for i in range(len(t)):
    interstory_disp[i, 0] = x[i, 1] - x[i, 0]  # 1-2层间
    interstory_disp[i, 1] = x[i, 2] - x[i, 1]  # 2-3层间

fig2, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(t, interstory_disp[:, 0]*1000, 'b-', linewidth=1.5)
axes[0].set_xlabel('Time (s)', fontsize=11)
axes[0].set_ylabel('Interstory Displacement (mm)', fontsize=11)
axes[0].set_title('1-2 Interstory Displacement', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, interstory_disp[:, 1]*1000, 'r-', linewidth=1.5)
axes[1].set_xlabel('Time (s)', fontsize=11)
axes[1].set_ylabel('Interstory Displacement (mm)', fontsize=11)
axes[1].set_title('2-3 Interstory Displacement', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/interstory_disp.png', dpi=150, bbox_inches='tight')
plt.close()
print("  层间位移图已保存")

# ============================================================
# 不同时间步长的影响
# ============================================================
print("\n不同时间步长的影响分析：")
print("-" * 70)

dt_values = [0.005, 0.01, 0.02, 0.05]
peak_disps = []

for dt_test in dt_values:
    t_test, x_test, v_test, a_test = newmark_beta(M, C, K, seismic_load, x0, v0, a0, dt_test, t_end)
    peak_disp = np.max(np.abs(x_test[:, 2]))  # 顶层峰值位移
    peak_disps.append(peak_disp)
    print(f"  dt = {dt_test:.3f} s: 顶层峰值位移 = {peak_disp*1000:.2f} mm")

# 绘制时间步长影响
fig3, ax = plt.subplots(figsize=(10, 6))

ax.plot(dt_values, np.array(peak_disps)*1000, 'b-o', linewidth=2, markersize=8)
ax.set_xlabel('Time Step (s)', fontsize=11)
ax.set_ylabel('Peak Displacement (mm)', fontsize=11)
ax.set_title('Effect of Time Step on Peak Displacement', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/time_step_effect.png', dpi=150, bbox_inches='tight')
plt.close()
print("  时间步长影响图已保存")

print("\n" + "=" * 70)
print("Newmark-β法动力响应分析完成！")
print("=" * 70)

6.2 案例2：频率响应分析

本案例实现频率响应分析，计算结构在谐波荷载作用下的响应。

# -*- coding: utf-8 -*-
"""
主题015：结构动力响应分析方法 - 案例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仿真领域\工程仿真\结构动力学仿真\主题015'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 70)
print("案例2：频率响应分析")
print("=" * 70)

# ============================================================
# 多自由度结构模型
# ============================================================
print("\n多自由度结构模型：")
print("-" * 70)

# 质量矩阵
M = np.array([[1000, 0, 0],
              [0, 1000, 0],
              [0, 0, 1000]])

# 刚度矩阵
K = np.array([[40000, -20000, 0],
              [-20000, 40000, -20000],
              [0, -20000, 20000]])

# 阻尼矩阵（瑞利阻尼）
alpha = 0.1
beta = 0.005
C = alpha * M + beta * K

print("  质量矩阵 M (kg):")
print(M)
print("\n  刚度矩阵 K (N/m):")
print(K)
print("\n  阻尼矩阵 C (N·s/m):")
print(C)

# 计算固有频率
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(M) @ K)
idx = np.argsort(eigenvalues)
omega_n = np.sqrt(eigenvalues[idx])
f_n = omega_n / (2 * np.pi)

print(f"\n固有频率：")
for i, (wn, fn) in enumerate(zip(omega_n, f_n)):
    print(f"  第{i+1}阶: ω = {wn:.4f} rad/s, f = {fn:.4f} Hz")

# ============================================================
# 频率响应函数计算
# ============================================================
def frequency_response_function(M, C, K, omega_range):
    """
    计算频率响应函数
    
    参数：
    M: 质量矩阵
    C: 阻尼矩阵
    K: 刚度矩阵
    omega_range: 频率范围
    
    返回：
    H: 频率响应函数矩阵
    """
    n_dof = M.shape[0]
    n_freq = len(omega_range)
    H = np.zeros((n_dof, n_dof, n_freq), dtype=complex)
    
    for i, omega in enumerate(omega_range):
        # 构建动力刚度矩阵
        K_dyn = -omega**2 * M + 1j * omega * C + K
        # 计算频率响应函数（逆矩阵）
        H[:, :, i] = np.linalg.inv(K_dyn)
    
    return H

# ============================================================
# 频率范围
# ============================================================
print("\n频率响应分析：")
print("-" * 70)

# 频率范围（覆盖前3阶固有频率）
f_min = 0.1
top_freq = f_n[-1] * 2
f_range = np.logspace(np.log10(f_min), np.log10(top_freq), 200)
omega_range = 2 * np.pi * f_range

print(f"  频率范围: {f_min:.1f} ~ {top_freq:.1f} Hz")
print(f"  频率点数: {len(f_range)}")

# 计算频率响应函数
H = frequency_response_function(M, C, K, omega_range)

# 荷载分布（作用于顶层）
F0 = 1000.0  # 荷载幅值
load_distribution = np.array([0, 0, 1])  # 仅顶层受荷

# 计算响应
response_amp = np.zeros((3, len(f_range)))
response_phase = np.zeros((3, len(f_range)))

for i in range(len(f_range)):
    # 计算响应
    X = H[:, :, i] @ (F0 * load_distribution)
    response_amp[:, i] = np.abs(X)
    response_phase[:, i] = np.angle(X) * 180 / np.pi  # 转换为度

# ============================================================
# 结果分析
# ============================================================
print("\n频率响应分析结果：")
print("-" * 70)

for i in range(3):
    # 找到峰值响应频率
    peak_idx = np.argmax(response_amp[i])
    peak_freq = f_range[peak_idx]
    peak_amp = response_amp[i, peak_idx]
    
    print(f"\n第{i+1}层：")
    print(f"  峰值响应频率: {peak_freq:.2f} Hz")
    print(f"  峰值响应幅值: {peak_amp*1000:.2f} mm")
    print(f"  对应固有频率: {f_n[i]:.2f} Hz")

# ============================================================
# 绘制频率响应曲线
# ============================================================
print("\n正在生成可视化图表...")

fig1, axes = plt.subplots(3, 2, figsize=(16, 12))

for i in range(3):
    # 幅值响应
    axes[i, 0].semilogx(f_range, response_amp[i]*1000, 'b-', linewidth=2)
    axes[i, 0].axvline(f_n[i], color='r', linestyle='--', alpha=0.7, label=f'f{i+1} = {f_n[i]:.2f} Hz')
    axes[i, 0].set_xlabel('Frequency (Hz)', fontsize=11)
    axes[i, 0].set_ylabel('Amplitude (mm)', fontsize=11)
    axes[i, 0].set_title(f'Floor {i+1} Amplitude Response', fontsize=12, fontweight='bold')
    axes[i, 0].legend()
    axes[i, 0].grid(True, alpha=0.3)
    
    # 相位响应
    axes[i, 1].semilogx(f_range, response_phase[i], 'g-', linewidth=2)
    axes[i, 1].axvline(f_n[i], color='r', linestyle='--', alpha=0.7, label=f'f{i+1} = {f_n[i]:.2f} Hz')
    axes[i, 1].set_xlabel('Frequency (Hz)', fontsize=11)
    axes[i, 1].set_ylabel('Phase (deg)', fontsize=11)
    axes[i, 1].set_title(f'Floor {i+1} Phase Response', fontsize=12, fontweight='bold')
    axes[i, 1].legend()
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/frequency_response.png', dpi=150, bbox_inches='tight')
plt.close()
print("  频率响应图已保存")

# ============================================================
# 三维频率响应图
# ============================================================
fig2 = plt.figure(figsize=(12, 8))
ax = fig2.add_subplot(111, projection='3d')

# 网格
X, Y = np.meshgrid(f_range, np.arange(1, 4))
Z = response_amp.T * 1000  # 转换为mm

# 绘制表面
surf = ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)

ax.set_xlabel('Frequency (Hz)', fontsize=11)
ax.set_ylabel('Floor', fontsize=11)
ax.set_zlabel('Amplitude (mm)', fontsize=11)
ax.set_title('3D Frequency Response Surface', fontsize=12, fontweight='bold')

# 添加颜色条
fig2.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

plt.tight_layout()
plt.savefig(f'{output_dir}/frequency_response_3d.png', dpi=150, bbox_inches='tight')
plt.close()
print("  3D频率响应图已保存")

# ============================================================
# 不同阻尼比的影响
# ============================================================
print("\n不同阻尼比的影响分析：")
print("-" * 70)

damping_ratios = [0.02, 0.05, 0.10, 0.20]
response_amp_damping = []

for zeta in damping_ratios:
    # 重新计算阻尼矩阵
    C_zeta = 2 * zeta * np.sqrt(M[0,0] * K[0,0]) * np.eye(3)
    # 计算频率响应
    H_zeta = frequency_response_function(M, C_zeta, K, omega_range)
    # 计算顶层响应
    X_zeta = np.abs(H_zeta[2, 2, :] * F0)
    response_amp_damping.append(X_zeta)
    
    print(f"  阻尼比 ζ = {zeta:.2f}: 顶层峰值响应 = {np.max(X_zeta)*1000:.2f} mm")

# 绘制阻尼比影响
fig3, ax = plt.subplots(figsize=(10, 6))

for i, zeta in enumerate(damping_ratios):
    ax.semilogx(f_range, response_amp_damping[i]*1000, linewidth=2, label=f'ζ = {zeta:.2f}')

ax.axvline(f_n[2], color='r', linestyle='--', alpha=0.7, label=f'f3 = {f_n[2]:.2f} Hz')
ax.set_xlabel('Frequency (Hz)', fontsize=11)
ax.set_ylabel('Amplitude (mm)', fontsize=11)
ax.set_title('Effect of Damping Ratio on Frequency Response', fontsize=12, fontweight='bold')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/damping_effect.png', dpi=150, bbox_inches='tight')
plt.close()
print("  阻尼比影响图已保存")

print("\n" + "=" * 70)
print("频率响应分析完成！")
print("=" * 70)

6.3 案例3：模态叠加法

本案例实现模态叠加法求解结构的动力响应。

# -*- coding: utf-8 -*-
"""
主题015：结构动力响应分析方法 - 案例3
模态叠加法
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题015'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 70)
print("案例3：模态叠加法")
print("=" * 70)

# ============================================================
# 多自由度结构模型
# ============================================================
print("\n多自由度结构模型：")
print("-" * 70)

# 质量矩阵
M = np.array([[1000, 0, 0],
              [0, 1000, 0],
              [0, 0, 1000]])

# 刚度矩阵
K = np.array([[40000, -20000, 0],
              [-20000, 40000, -20000],
              [0, -20000, 20000]])

# 阻尼矩阵（瑞利阻尼）
alpha = 0.1
beta = 0.005
C = alpha * M + beta * K

print("  质量矩阵 M (kg):")
print(M)
print("\n  刚度矩阵 K (N/m):")
print(K)
print("\n  阻尼矩阵 C (N·s/m):")
print(C)

# ============================================================
# 模态分析
# ============================================================
print("\n模态分析：")
print("-" * 70)

# 求解特征值问题
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(M) @ K)

# 排序（从小到大）
idx = np.argsort(eigenvalues)
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# 计算固有频率和阻尼比
omega_n = np.sqrt(np.real(eigenvalues))
f_n = omega_n / (2 * np.pi)

# 模态阻尼比（假设各阶模态阻尼比相同）
zeta_n = np.zeros(3)
for i in range(3):
    # 计算模态阻尼
    phi = eigenvectors[:, i].reshape(-1, 1)
    c_modal = phi.T @ C @ phi
    m_modal = phi.T @ M @ phi
    zeta_n[i] = c_modal[0, 0] / (2 * np.sqrt(m_modal[0, 0] * eigenvalues[i]))

print("  模态振型矩阵 Φ:")
print(eigenvectors)

print(f"\n固有频率和阻尼比：")
for i in range(3):
    print(f"  第{i+1}阶: ω = {omega_n[i]:.4f} rad/s, f = {f_n[i]:.4f} Hz, ζ = {zeta_n[i]:.4f}")

# 模态质量矩阵
M_modal = eigenvectors.T @ M @ eigenvectors
print("\n  模态质量矩阵 M_modal:")
print(M_modal)

# 模态刚度矩阵
K_modal = eigenvectors.T @ K @ eigenvectors
print("\n  模态刚度矩阵 K_modal:")
print(K_modal)

# 模态阻尼矩阵
C_modal = eigenvectors.T @ C @ eigenvectors
print("\n  模态阻尼矩阵 C_modal:")
print(C_modal)

# ============================================================
# 模态叠加法
# ============================================================
def modal_superposition(M, C, K, phi, omega_n, zeta_n, F_func, t):
    """
    模态叠加法求解动力响应
    
    参数：
    M: 质量矩阵
    C: 阻尼矩阵
    K: 刚度矩阵
    phi: 模态振型矩阵
    omega_n: 固有频率
    zeta_n: 模态阻尼比
    F_func: 荷载时程函数
    t: 时间数组
    
    返回：
    x: 位移时程
    """
    n_dof = M.shape[0]
    n_modes = phi.shape[1]
    n_steps = len(t)
    
    x = np.zeros((n_steps, n_dof))
    
    # 对每个模态求解
    for i in range(n_modes):
        print(f"  求解第{i+1}阶模态...")
        
        # 模态质量、刚度、阻尼
        m_i = phi[:, i].T @ M @ phi[:, i]
        k_i = phi[:, i].T @ K @ phi[:, i]
        c_i = phi[:, i].T @ C @ phi[:, i]
        
        # 模态荷载
        def modal_force(t_val):
            F = F_func(t_val)
            return phi[:, i].T @ F
        
        # 单自由度运动方程
        def sdof_system(y, t_val):
            q, q_dot = y
            F_modal = modal_force(t_val)
            q_ddot = (F_modal - c_i * q_dot - k_i * q) / m_i
            return [q_dot, q_ddot]
        
        # 初始条件
        y0 = [0, 0]
        
        # 求解单自由度方程
        solution = odeint(sdof_system, y0, t)
        q = solution[:, 0]
        
        # 转换为物理坐标并叠加
        x_i = np.outer(q, phi[:, i])
        x += x_i
    
    return x

# ============================================================
# 荷载函数
# ============================================================
def pulse_load(t):
    """脉冲荷载"""
    if t < 0.5:
        return np.array([0, 0, 5000 * np.sin(np.pi * t / 0.5)])
    else:
        return np.array([0, 0, 0])

def harmonic_load(t):
    """谐波荷载"""
    omega = 2 * np.pi * 2.0  # 2Hz
    return np.array([0, 0, 1000 * np.sin(omega * t)])

# ============================================================
# 求解动力响应
# ============================================================
print("\n求解动力响应：")
print("-" * 70)

# 时间参数
t = np.linspace(0, 10, 500)

# 脉冲荷载响应
print("  求解脉冲荷载响应...")
x_pulse = modal_superposition(M, C, K, eigenvectors, omega_n, zeta_n, pulse_load, t)

# 谐波荷载响应
print("  求解谐波荷载响应...")
x_harmonic = modal_superposition(M, C, K, eigenvectors, omega_n, zeta_n, harmonic_load, t)

# ============================================================
# 结果分析
# ============================================================
print("\n响应结果分析：")
print("-" * 70)

print("脉冲荷载响应：")
for i in range(3):
    peak_disp = np.max(np.abs(x_pulse[:, i]))
    print(f"  第{i+1}层峰值位移: {peak_disp*1000:.2f} mm")

print("\n谐波荷载响应：")
for i in range(3):
    peak_disp = np.max(np.abs(x_harmonic[:, i]))
    print(f"  第{i+1}层峰值位移: {peak_disp*1000:.2f} mm")

# ============================================================
# 绘制结果
# ============================================================
print("\n正在生成可视化图表...")

fig1, axes = plt.subplots(3, 2, figsize=(16, 12))

for i in range(3):
    # 脉冲荷载响应
    axes[i, 0].plot(t, x_pulse[:, i]*1000, 'b-', linewidth=1.5)
    axes[i, 0].set_xlabel('Time (s)', fontsize=11)
    axes[i, 0].set_ylabel('Displacement (mm)', fontsize=11)
    axes[i, 0].set_title(f'Floor {i+1} - Pulse Load Response', fontsize=12, fontweight='bold')
    axes[i, 0].grid(True, alpha=0.3)
    
    # 谐波荷载响应
    axes[i, 1].plot(t, x_harmonic[:, i]*1000, 'r-', linewidth=1.5)
    axes[i, 1].set_xlabel('Time (s)', fontsize=11)
    axes[i, 1].set_ylabel('Displacement (mm)', fontsize=11)
    axes[i, 1].set_title(f'Floor {i+1} - Harmonic Load Response', fontsize=12, fontweight='bold')
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/modal_superposition.png', dpi=150, bbox_inches='tight')
plt.close()
print("  模态叠加法响应图已保存")

# ============================================================
# 模态参与分析
# ============================================================
print("\n模态参与分析：")
print("-" * 70)

# 计算模态参与因子
modal_participation = np.zeros(3)
for i in range(3):
    modal_participation[i] = np.sum(eigenvectors[:, i])
    print(f"  第{i+1}阶模态参与因子: {modal_participation[i]:.4f}")

# 计算有效质量
modal_mass = np.zeros(3)
for i in range(3):
    modal_mass[i] = (np.sum(eigenvectors[:, i]))**2 / M_modal[i, i]
    print(f"  第{i+1}阶有效质量: {modal_mass[i]:.2f} kg")

# 绘制模态参与和有效质量
fig2, axes = plt.subplots(1, 2, figsize=(14, 5))

modes = ['Mode 1', 'Mode 2', 'Mode 3']

axes[0].bar(modes, modal_participation, color='steelblue', alpha=0.8)
axes[0].set_xlabel('Mode', fontsize=11)
axes[0].set_ylabel('Participation Factor', fontsize=11)
axes[0].set_title('Modal Participation Factors', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3, axis='y')

axes[1].bar(modes, modal_mass, color='coral', alpha=0.8)
axes[1].set_xlabel('Mode', fontsize=11)
axes[1].set_ylabel('Effective Mass (kg)', fontsize=11)
axes[1].set_title('Modal Effective Mass', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig(f'{output_dir}/modal_participation.png', dpi=150, bbox_inches='tight')
plt.close()
print("  模态参与分析图已保存")

# ============================================================
# 不同模态数量的影响
# ============================================================
print("\n不同模态数量的影响分析：")
print("-" * 70)

n_modes_list = [1, 2, 3]
peak_disps = []

for n_modes in n_modes_list:
    # 使用前n_modes阶模态
    phi_truncated = eigenvectors[:, :n_modes]
    omega_truncated = omega_n[:n_modes]
    zeta_truncated = zeta_n[:n_modes]
    
    # 求解响应
    x_truncated = modal_superposition(M, C, K, phi_truncated, omega_truncated, zeta_truncated, pulse_load, t)
    
    # 计算顶层峰值位移
    peak_disp = np.max(np.abs(x_truncated[:, 2]))
    peak_disps.append(peak_disp)
    
    print(f"  前{n_modes}阶模态: 顶层峰值位移 = {peak_disp*1000:.2f} mm")

# 绘制模态数量影响
fig3, ax = plt.subplots(figsize=(10, 6))

ax.plot(n_modes_list, np.array(peak_disps)*1000, 'b-o', linewidth=2, markersize=8)
ax.set_xlabel('Number of Modes', fontsize=11)
ax.set_ylabel('Peak Displacement (mm)', fontsize=11)
ax.set_title('Effect of Number of Modes on Response', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/modal_number_effect.png', dpi=150, bbox_inches='tight')
plt.close()
print("  模态数量影响图已保存")

print("\n" + "=" * 70)
print("模态叠加法分析完成！")
print("=" * 70)

6.4 案例4：随机振动分析

本案例实现随机振动分析，计算结构在随机荷载作用下的响应。

# -*- coding: utf-8 -*-
"""
主题015：结构动力响应分析方法 - 案例4
随机振动分析
"""
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仿真领域\工程仿真\结构动力学仿真\主题015'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 70)
print("案例4：随机振动分析")
print("=" * 70)

# ============================================================
# 单自由度结构模型
# ============================================================
print("\n单自由度结构模型：")
print("-" * 70)

m = 1000.0  # 质量 (kg)
k = 40000.0  # 刚度 (N/m)
zeta = 0.05  # 阻尼比

omega_n = np.sqrt(k / m)  # 固有频率
c = 2 * zeta * m * omega_n  # 阻尼系数

print(f"  质量 m = {m} kg")
print(f"  刚度 k = {k} N/m")
print(f"  阻尼比 ζ = {zeta}")
print(f"  固有频率 ω_n = {omega_n:.4f} rad/s")
print(f"  阻尼系数 c = {c:.2f} N·s/m")

# ============================================================
# 随机荷载生成
# ============================================================
def generate_random_load(t, dt, seed=42):
    """生成随机荷载时程"""
    np.random.seed(seed)
    n = len(t)
    
    # 生成白噪声
    white_noise = np.random.randn(n)
    
    # 带通滤波（模拟风荷载或地震动）
    sos = signal.butter(4, [0.1, 10], 'bandpass', fs=1/dt, output='sos')
    filtered = signal.sosfilt(sos, white_noise)
    
    # 调整功率谱密度
    freq = np.fft.fftfreq(n, dt)
    psd = np.zeros(n)
    
    # 模拟风荷载PSD（随频率增加而减小）
    for i in range(n):
        f = abs(freq[i])
        if f > 0:
            psd[i] = 1.0 / (1 + f**2)
    
    # 应用PSD
    filtered_fft = np.fft.fft(filtered)
    filtered_fft *= np.sqrt(psd)
    random_load = np.real(np.fft.ifft(filtered_fft))
    
    # 标准化
    random_load = random_load / np.std(random_load) * 1000  # 标准差为1000N
    
    return random_load

# ============================================================
# 时间参数
# ============================================================
dt = 0.01
t_end = 60.0
t = np.arange(0, t_end, dt)

print(f"\n时间参数：")
print(f"  时间步长: {dt} s")
print(f"  持续时间: {t_end} s")
print(f"  总点数: {len(t)}")

# 生成随机荷载
F = generate_random_load(t, dt)

# ============================================================
# 求解随机响应（时域方法）
# ============================================================
def sdof_random_response(m, c, k, F, dt):
    """求解单自由度系统的随机响应"""
    n = len(F)
    x = np.zeros(n)
    v = np.zeros(n)
    a = np.zeros(n)
    
    # 初始条件
    x[0] = 0
    v[0] = 0
    a[0] = (F[0] - c * v[0] - k * x[0]) / m
    
    # 使用中心差分法
    for i in range(1, n):
        # 预测加速度
        a[i] = (F[i] - c * v[i-1] - k * x[i-1]) / m
        
        # 更新速度和位移
        v[i] = v[i-1] + a[i] * dt
        x[i] = x[i-1] + v[i-1] * dt + 0.5 * a[i-1] * dt**2
    
    return x, v, a

print("\n求解随机响应：")
print("-" * 70)

x, v, a = sdof_random_response(m, c, k, F, dt)

# ============================================================
# 统计分析
# ============================================================
print("\n响应统计分析：")
print("-" * 70)

# 位移统计
x_mean = np.mean(x)
x_std = np.std(x)
x_rms = np.sqrt(np.mean(x**2))
x_max = np.max(np.abs(x))

# 速度统计
v_mean = np.mean(v)
v_std = np.std(v)
v_rms = np.sqrt(np.mean(v**2))
v_max = np.max(np.abs(v))

# 加速度统计
a_mean = np.mean(a)
a_std = np.std(a)
a_rms = np.sqrt(np.mean(a**2))
a_max = np.max(np.abs(a))

print(f"位移统计：")
print(f"  均值: {x_mean*1000:.2f} mm")
print(f"  标准差: {x_std*1000:.2f} mm")
print(f"  RMS: {x_rms*1000:.2f} mm")
print(f"  最大值: {x_max*1000:.2f} mm")

print(f"\n速度统计：")
print(f"  均值: {v_mean:.2f} m/s")
print(f"  标准差: {v_std:.2f} m/s")
print(f"  RMS: {v_rms:.2f} m/s")
print(f"  最大值: {v_max:.2f} m/s")

print(f"\n加速度统计：")
print(f"  均值: {a_mean:.2f} m/s²")
print(f"  标准差: {a_std:.2f} m/s²")
print(f"  RMS: {a_rms:.2f} m/s²")
print(f"  最大值: {a_max:.2f} m/s²")

# ============================================================
# 功率谱密度分析
# ============================================================
print("\n功率谱密度分析：")
print("-" * 70)

# 计算荷载和响应的功率谱密度
freq_F, Pxx_F = signal.welch(F, fs=1/dt, nperseg=1024)
freq_x, Pxx_x = signal.welch(x, fs=1/dt, nperseg=1024)
freq_v, Pxx_v = signal.welch(v, fs=1/dt, nperseg=1024)
freq_a, Pxx_a = signal.welch(a, fs=1/dt, nperseg=1024)

# 计算理论频率响应函数
omega = 2 * np.pi * freq_F
H = 1 / (-m * omega**2 + 1j * c * omega + k)
H_mag = np.abs(H)

# 理论响应PSD
Pxx_x_theory = H_mag**2 * Pxx_F

print(f"  荷载PSD峰值频率: {freq_F[np.argmax(Pxx_F)]:.2f} Hz")
print(f"  位移PSD峰值频率: {freq_x[np.argmax(Pxx_x)]:.2f} Hz")
print(f"  理论峰值频率: {omega_n/(2*np.pi):.2f} Hz")

# ============================================================
# 绘制结果
# ============================================================
print("\n正在生成可视化图表...")

fig1, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 荷载时程
axes[0, 0].plot(t, F, 'k-', linewidth=1)
axes[0, 0].set_xlabel('Time (s)', fontsize=11)
axes[0, 0].set_ylabel('Force (N)', fontsize=11)
axes[0, 0].set_title('Random Load Time History', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3)

# 2. 位移时程
axes[0, 1].plot(t, x*1000, 'b-', linewidth=1)
axes[0, 1].set_xlabel('Time (s)', fontsize=11)
axes[0, 1].set_ylabel('Displacement (mm)', fontsize=11)
axes[0, 1].set_title('Displacement Time History', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)

# 3. 速度时程
axes[1, 0].plot(t, v, 'g-', linewidth=1)
axes[1, 0].set_xlabel('Time (s)', fontsize=11)
axes[1, 0].set_ylabel('Velocity (m/s)', fontsize=11)
axes[1, 0].set_title('Velocity Time History', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)

# 4. 加速度时程
axes[1, 1].plot(t, a, 'r-', linewidth=1)
axes[1, 1].set_xlabel('Time (s)', fontsize=11)
axes[1, 1].set_ylabel('Acceleration (m/s²)', fontsize=11)
axes[1, 1].set_title('Acceleration Time History', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/random_time_history.png', dpi=150, bbox_inches='tight')
plt.close()
print("  随机时程图已保存")

# 功率谱密度图
fig2, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 荷载PSD
axes[0, 0].semilogx(freq_F, Pxx_F, 'k-', linewidth=2)
axes[0, 0].set_xlabel('Frequency (Hz)', fontsize=11)
axes[0, 0].set_ylabel('PSD (N²/Hz)', fontsize=11)
axes[0, 0].set_title('Load Power Spectral Density', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3)

# 2. 位移PSD
axes[0, 1].semilogx(freq_x, Pxx_x, 'b-', linewidth=2, label='Simulation')
axes[0, 1].semilogx(freq_F, Pxx_x_theory, 'r--', linewidth=2, label='Theory')
axes[0, 1].axvline(omega_n/(2*np.pi), color='g', linestyle='--', alpha=0.7, label=f'f_n = {omega_n/(2*np.pi):.2f} Hz')
axes[0, 1].set_xlabel('Frequency (Hz)', fontsize=11)
axes[0, 1].set_ylabel('PSD (m²/Hz)', fontsize=11)
axes[0, 1].set_title('Displacement Power Spectral Density', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 3. 速度PSD
axes[1, 0].semilogx(freq_v, Pxx_v, 'g-', linewidth=2)
axes[1, 0].axvline(omega_n/(2*np.pi), color='r', linestyle='--', alpha=0.7, label=f'f_n = {omega_n/(2*np.pi):.2f} Hz')
axes[1, 0].set_xlabel('Frequency (Hz)', fontsize=11)
axes[1, 0].set_ylabel('PSD (m²/s²/Hz)', fontsize=11)
axes[1, 0].set_title('Velocity Power Spectral Density', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 4. 加速度PSD
axes[1, 1].semilogx(freq_a, Pxx_a, 'r-', linewidth=2)
axes[1, 1].axvline(omega_n/(2*np.pi), color='g', linestyle='--', alpha=0.7, label=f'f_n = {omega_n/(2*np.pi):.2f} Hz')
axes[1, 1].set_xlabel('Frequency (Hz)', fontsize=11)
axes[1, 1].set_ylabel('PSD (m²/s⁴/Hz)', fontsize=11)
axes[1, 1].set_title('Acceleration Power Spectral Density', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/random_psd.png', dpi=150, bbox_inches='tight')
plt.close()
print("  功率谱密度图已保存")

# ============================================================
# 概率分布分析
# ============================================================
print("\n概率分布分析：")
print("-" * 70)

# 计算位移的概率分布
hist_x, bins_x = np.histogram(x, bins=50, density=True)
mid_bins_x = (bins_x[:-1] + bins_x[1:]) / 2

# 理论正态分布
x_normal = np.linspace(np.min(x), np.max(x), 100)
y_normal = (1 / (x_std * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((x_normal - x_mean) / x_std)**2)

print(f"  位移分布均值: {x_mean*1000:.2f} mm")
print(f"  位移分布标准差: {x_std*1000:.2f} mm")

# 绘制概率分布图
fig3, ax = plt.subplots(figsize=(10, 6))

ax.hist(x*1000, bins=50, density=True, alpha=0.7, color='steelblue', label='Histogram')
ax.plot(x_normal*1000, y_normal*1000, 'r-', linewidth=2, label='Normal Distribution')
ax.set_xlabel('Displacement (mm)', fontsize=11)
ax.set_ylabel('Probability Density', fontsize=11)
ax.set_title('Displacement Probability Distribution', fontsize=12, fontweight='bold')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/displacement_distribution.png', dpi=150, bbox_inches='tight')
plt.close()
print("  概率分布图已保存")

print("\n" + "=" * 70)
print("随机振动分析完成！")
print("=" * 70)

7. 总结与展望

7.1 本教程总结

本教程系统地介绍了结构动力响应分析的各种方法，包括：

时域分析方法：详细讲解了直接积分法（如Newmark-β法）的原理和实现，适用于任意荷载和非线性系统。
频域分析方法：介绍了傅里叶变换、频率响应函数等概念，适用于线性系统和周期性荷载。
模态分析方法：讲解了模态叠加法的原理和应用，利用结构的模态特性提高计算效率。
随机振动分析：介绍了随机过程基础、功率谱密度分析和响应统计特性计算。

通过四个Python案例实现，读者可以掌握：

Newmark-β法求解动力响应
频率响应分析识别结构共振特性
模态叠加法高效计算多自由度系统响应
随机振动分析评估结构在随机荷载下的性能

7.2 应用前景

结构动力响应分析在工程实践中具有广泛的应用前景：

结构设计：优化结构参数，提高结构的动力性能
抗震设计：评估结构在地震作用下的响应
风振分析：预测结构在风荷载作用下的振动水平
设备振动：分析设备振动对结构的影响
振动控制：为振动控制策略的设计提供依据

7.3 未来发展方向

结构动力响应分析的未来发展方向包括：

非线性动力分析：考虑材料非线性、几何非线性等因素
多物理场耦合分析：如结构-流体耦合、结构-土壤耦合等
不确定性分析：考虑参数不确定性对响应的影响
智能算法应用：如机器学习在动力响应预测中的应用
实时监测与分析：结合传感器技术实现结构响应的实时监测和分析

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

ArkUI 项目结构设计：小项目 vs 大项目

AtomGit开源社区

当工业AI原生企业撞上制造业——全面提质增效时代来临

大制造贯穿从订单到计划、生产、交付的全流程，通过速赢+AI、卓越+AI打造覆盖计划、生产、质量、仓储、能源、设备六大核心域的制造智能体，链接研发与制造工程、构建端到端的制造闭环。而广域铭岛在更复杂的工业场景中，同样用可量化的收益回答了工业AI的现实价值。工业AI的核心，不是单点算法，而是覆盖产品需求、设计、工艺验证的大研发，以及从订单到计划、生产、交付的大制造。传统工厂引入AI往往只是局部点状提效