主题034：阻尼模型参数识别与应用

目录

1. 概述
2. 阻尼的基本概念
3. 瑞利阻尼模型
4. 模态阻尼模型
5. 非比例阻尼系统
6. 阻尼参数识别方法
7. 工程应用案例
8. Python仿真实现

---

概述

阻尼是结构动力学中描述能量耗散的重要参数，它反映了结构在振动过程中由于材料内摩擦、节点摩擦、空气阻力等因素导致的能量损失。准确识别和应用阻尼模型对于结构地震响应分析、振动控制和疲劳寿命预测具有重要意义。

本主题将系统介绍：

• 阻尼的物理机制和数学模型
• 瑞利阻尼的参数识别方法
• 模态阻尼的应用
• 非比例阻尼系统的处理
• 工程实际中的阻尼参数取值

---

阻尼的基本概念

2.1 阻尼的物理来源

结构中的能量耗散主要来自以下几个方面：

材料阻尼

• 材料内部晶格缺陷的运动
• 分子链的摩擦
• 应力-应变滞回曲线包围的面积

结构阻尼

• 节点和连接处的摩擦
• 构件之间的相对滑移
• 非结构构件的耗能

流体阻尼

• 空气阻力
• 液体粘滞效应
• 涡激振动

辐射阻尼

• 地基的波辐射
• 能量向无限域传播

2.2 阻尼的数学描述

对于单自由度系统，阻尼力可以表示为：

$$F_d=c\cdot \dot{u}$$

其中：

• $c$ 为阻尼系数
• $\dot{u}$ 为速度

运动方程为：

$$m\ddot{u}+c\dot{u}+ku=F(t)$$

2.3 阻尼比与临界阻尼

临界阻尼系数：

$$c_{cr}=2\sqrt{km}=2m{\omega }_n$$

阻尼比：

$$\xi =\frac{c}{c_{cr}}=\frac{c}{2m{\omega }_n}$$

不同阻尼比下系统的响应特性：

阻尼比范围	系统特性	应用场景
$\xi <1$	欠阻尼，振荡衰减	大多数工程结构
$\xi =1$	临界阻尼，最快回到平衡	减震装置设计
$\xi >1$	过阻尼，缓慢回到平衡	特殊隔震系统

2.4 工程结构的典型阻尼比

根据大量实测数据，不同类型结构的阻尼比典型值如下：

钢结构

• 焊接钢框架：1.0% - 2.0%
• 螺栓连接钢框架：2.0% - 3.0%
• 钢桥：0.5% - 1.5%

混凝土结构

• 钢筋混凝土框架：2.0% - 5.0%
• 预应力混凝土：1.0% - 2.0%
• 混凝土坝：5.0% - 10.0%

其他结构

• 木框架：3.0% - 5.0%
• 砌体结构：3.0% - 6.0%
• 索结构：0.3% - 1.0%

---

瑞利阻尼模型

3.1 瑞利阻尼的基本形式

瑞利阻尼（Rayleigh Damping）是最常用的阻尼模型之一，其阻尼矩阵与质量矩阵和刚度矩阵成线性关系：

$$\mathbf{C}=\alpha \mathbf{M}+\beta \mathbf{K}$$

其中：

• $\alpha$ 为质量比例系数
• $\beta$ 为刚度比例系数

3.2 瑞利阻尼的模态特性

对于第 $i$ 阶模态，瑞利阻尼对应的模态阻尼比为：

$${\xi }_i=\frac{\alpha }{2{\omega }_i}+\frac{\beta {\omega }_i}{2}$$

这个公式表明：

• 质量比例项 $\alpha$ 在高频模态贡献较小
• 刚度比例项 $\beta$ 在低频模态贡献较小
• 阻尼比随频率变化呈"V"字形曲线

3.3 瑞利阻尼参数的确定

方法一：指定两阶模态的阻尼比

假设已知第 $i$ 阶和第 $j$ 阶模态的阻尼比 ${\xi }_i$ 和 ${\xi }_j$，则：

$$\left[\begin{array}{c}{\xi }_i\\ {\xi }_j\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1/{\omega }_i& {\omega }_i\\ 1/{\omega }_j& {\omega }_j\end{array}\right]\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]$$

解得：

$$\alpha =\frac{2{\omega }_i{\omega }_j({\xi }_i{\omega }_j-{\xi }_j{\omega }_i)}{{\omega }_j^2-{\omega }_i^2}$$

$$\beta =\frac{2({\xi }_j{\omega }_j-{\xi }_i{\omega }_i)}{{\omega }_j^2-{\omega }_i^2}$$

方法二：最小二乘法拟合

当需要拟合多阶模态的阻尼比时，可以采用最小二乘法：

$$\underset{\alpha ,\beta }{\min }\sum _{i=1}^n{\left({\xi }_i-\frac{\alpha }{2{\omega }_i}-\frac{\beta {\omega }_i}{2}\right)}^2$$

3.4 瑞利阻尼的优缺点

优点：

• 数学形式简单，易于实现
• 保持模态正交性，便于模态叠加分析
• 计算效率高

缺点：

• 难以准确拟合宽频带内的阻尼特性
• 高频模态阻尼比可能过大或过小
• 不能反映复杂的阻尼机制

---

模态阻尼模型

4.1 模态阻尼的概念

模态阻尼（Modal Damping）假设各阶模态之间相互独立，每阶模态具有独立的阻尼比：

$${\xi }_i^{(modal)}={\xi }_i$$

在模态坐标下，运动方程解耦为：

$${\ddot{q}}_i+2{\xi }_i{\omega }_i{\dot{q}}_i+{\omega }_i^2{q}_i=\frac{{\varphi }_i^{T}F(t)}{M_i^{\ast }}$$

4.2 模态阻尼矩阵的构建

模态阻尼矩阵可以通过以下方式构建：

$$\mathbf{C}=\mathbf{M}\left(\sum _{i=1}^n\frac{2{\xi }_i{\omega }_i}{M_i^{\ast }}{\varphi }_i{\varphi }_i^T\right)\mathbf{M}$$

其中：

• $M_i^{\ast }={\varphi }_i^T\mathbf{M}{\varphi }_i$ 为模态质量
• ${\varphi }_i$ 为模态振型

4.3 模态阻尼的应用

适用条件：

• 阻尼较小（通常 $\xi <0.2$）
• 模态之间耦合较弱
• 响应主要由前几阶模态控制

实施步骤：

1. 进行模态分析，获取频率和振型
2. 为每阶模态指定合适的阻尼比
3. 在模态坐标下求解运动方程
4. 通过模态叠加获得物理坐标响应

---

非比例阻尼系统

5.1 非比例阻尼的定义

当阻尼矩阵不能满足以下正交条件时，称为非比例阻尼：

$${\varphi }_i^T\mathbf{C}{\varphi }_j\ne 0,i\ne j$$

非比例阻尼的常见情况：

• 结构-土相互作用
• 隔震/减震装置
• 流体-结构相互作用
• 局部附加阻尼

5.2 非比例阻尼的处理方法

方法1：复模态分析

将运动方程转化为状态空间形式：

$$\left[\begin{array}{cc}\mathbf{M}& 0\\ 0& \mathbf{M}\end{array}\right]\left[\begin{array}{c}\ddot{u}\\ \dot{u}\end{array}\right]+\left[\begin{array}{cc}\mathbf{C}& \mathbf{K}\\ -\mathbf{M}& 0\end{array}\right]\left[\begin{array}{c}\dot{u}\\ u\end{array}\right]=\left[\begin{array}{c}F(t)\\ 0\end{array}\right]$$

求解广义特征值问题，得到复频率和复振型。

方法2：直接积分法

直接在物理坐标下进行时程积分，无需解耦：

$$\mathbf{M}\ddot{u}+\mathbf{C}\dot{u}+\mathbf{K}u=F(t)$$

常用的积分方法包括：

• Newmark-β法
• Wilson-θ法
• 中心差分法

5.3 非比例阻尼的识别

对于已知的非比例阻尼系统，可以通过实验模态分析识别阻尼矩阵：

$$\mathbf{C}=\mathbf{M}\left(\sum _{r=1}^{n}2{\xi }_r{\omega }_r\frac{{\varphi }_r{\varphi }_r^T}{M_r^{\ast }}\right)\mathbf{M}+{\mathbf{C}}_{non}$$

其中 ${\mathbf{C}}_{non}$ 为非比例部分。

---

阻尼参数识别方法

6.1 自由振动衰减法

原理：利用自由振动的振幅衰减计算阻尼比。

对数衰减率：

$$\delta =\ln \frac{u_n}{u_{n+m}}=\frac{2m\pi \xi }{\sqrt{1-{\xi }^2}}$$

当 $\xi \ll 1$ 时：

$$\xi \approx \frac{\delta }{2m\pi }$$

实施步骤：

1. 测量结构自由振动的位移时程
2. 确定相邻（或间隔m个）振幅峰值
3. 计算对数衰减率
4. 求解阻尼比

6.2 半功率带宽法

原理：利用频响函数的半功率带宽计算阻尼比。

$$\xi =\frac{\Delta \omega }{2{\omega }_n}$$

其中 $\Delta \omega$ 为半功率点之间的频率差。

实施步骤：

1. 测量或计算频响函数
2. 确定共振峰值
3. 找到峰值两侧半功率点（$1/\sqrt{2}$ 倍峰值）
4. 计算带宽并求解阻尼比

6.3 模态参数识别法

特征系统实现算法（ERA）

利用脉冲响应数据构建Hankel矩阵：

$$\mathbf{H}(k)=\left[\begin{array}{cccc}h(k)& h(k+1)& \cdots & h(k+s-1)\\ h(k+1)& h(k+2)& \cdots & h(k+s)\\ ⋮& ⋮& \ddots & ⋮\\ h(k+r-1)& h(k+r)& \cdots & h(k+r+s-2)\end{array}\right]$$

通过奇异值分解提取模态参数。

随机子空间识别（SSI）

适用于环境激励下的模态参数识别，包括：

• 数据驱动随机子空间（DATA-SSI）
• 协方差驱动随机子空间（COV-SSI）

6.4 优化识别方法

最小二乘优化

目标函数：

$$J=\sum _{i=1}^N{\left(y_i^{exp}-y_i^{sim}(\theta )\right)}^2$$

其中：

• $y_i^{exp}$ 为实验测量值
• $y_i^{sim}$ 为仿真计算值
• $\theta$ 为待识别的阻尼参数

贝叶斯方法

利用贝叶斯定理更新参数的后验分布：

$$p(\theta |D)=\frac{p(D|\theta )p(\theta )}{p(D)}$$

其中：

• $p(\theta |D)$ 为后验分布
• $p(D|\theta )$ 为似然函数
• $p(\theta )$ 为先验分布

---

工程应用案例

7.1 高层建筑阻尼比识别

项目背景：某50层钢筋混凝土框架-核心筒结构

识别方法：环境振动测试 + 随机子空间识别

识别结果：

模态阶数	频率(Hz)	阻尼比(%)	振型描述
1	0.25	2.1	X向平动
2	0.28	2.3	Y向平动
3	0.45	1.8	扭转
4	0.72	1.5	X向二阶
5	0.78	1.6	Y向二阶

工程应用：

• 基于识别的阻尼比进行时程分析
• 评估结构的风振响应
• 验证设计假设的合理性

7.2 桥梁结构阻尼模型建立

项目背景：某大跨度斜拉桥

阻尼特性：

• 钢主梁：低阻尼（0.5%-1.0%）
• 拉索：极低阻尼（0.1%-0.3%）
• 桥塔：中等阻尼（1.5%-2.5%）

建模策略：

1. 对不同构件采用不同的瑞利阻尼参数
2. 拉索采用附加粘滞阻尼器模型
3. 整体组装非比例阻尼矩阵

分析结果：

• 考虑真实阻尼的响应比常阻尼假设小15-20%
• 拉索振动得到有效控制

7.3 隔震结构阻尼设计

项目背景：某医院建筑采用铅芯橡胶支座隔震

阻尼组成：

• 上部结构：钢筋混凝土框架，2.5%
• 隔震层：铅芯橡胶支座，15-20%
• 整体等效阻尼比：约12%

设计要点：

1. 隔震层提供大部分阻尼
2. 上部结构保持弹性
3. 阻尼器参数优化设计

效果评估：

• 基底剪力降低60%
• 上部结构加速度降低70%
• 满足罕遇地震下的性能目标

---

Python仿真实现

本主题将通过以下四个案例演示阻尼模型的参数识别与应用：

案例1：瑞利阻尼参数识别

• 根据目标模态阻尼比识别瑞利阻尼系数
• 分析阻尼比随频率的变化规律
• 评估瑞利阻尼的拟合精度

案例2：模态阻尼与时程分析

• 实现模态阻尼的时程分析
• 对比不同阻尼模型的响应差异
• 分析模态阻尼的适用范围

案例3：非比例阻尼系统分析

• 建立非比例阻尼矩阵
• 复模态分析实现
• 对比比例阻尼与非比例阻尼的差异

案例4：阻尼对结构响应的影响

• 参数化研究阻尼比的影响
• 分析阻尼对地震响应的敏感性
• 工程阻尼比取值建议

---

---

完整Python代码实现

以下是本主题的完整Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例1：瑞利阻尼参数识别
=====================================
演示瑞利阻尼系数的识别方法和阻尼比随频率的变化规律

本案例演示：
1. 根据目标模态阻尼比识别瑞利阻尼系数
2. 分析阻尼比随频率的变化规律
3. 评估瑞利阻尼的拟合精度
4. 最小二乘法拟合多阶模态阻尼比
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("案例1：瑞利阻尼参数识别")
print("=" * 70)

# ============================================================================
# 1. 创建示例结构模型
# ============================================================================
print("\n\n【1. 创建示例结构模型】")
print("=" * 70)

# 创建10层剪切型结构
n_dof = 10
m_floor = 10000  # kg
k_story = 5e7   # N/m

M = np.eye(n_dof) * m_floor

# 组装刚度矩阵
K = np.zeros((n_dof, n_dof))
for i in range(n_dof):
    if i == 0:
        K[i, i] = k_story
        if n_dof > 1:
            K[i, i+1] = -k_story
    elif i == n_dof - 1:
        K[i, i] = k_story
        K[i, i-1] = -k_story
    else:
        K[i, i] = 2 * k_story
        K[i, i-1] = -k_story
        K[i, i+1] = -k_story

# 模态分析 - 使用广义特征值求解
try:
    from scipy.linalg import eig
    eigvals, eigvecs = eig(K, M)
    idx = np.argsort(np.real(eigvals))
    eigvals_sorted = np.real(eigvals[idx])
except ImportError:
    # 手动求解广义特征值问题
    L = np.linalg.cholesky(M)
    L_inv = np.linalg.inv(L)
    K_tilde = L_inv @ K @ L_inv.T
    eigvals, eigvecs_tilde = np.linalg.eigh(K_tilde)
    eigvecs = L_inv.T @ eigvecs_tilde
    eigvals_sorted = eigvals

# 过滤掉负值或接近零的特征值
eigvals_sorted = np.where(eigvals_sorted > 1e-6, eigvals_sorted, 1e-6)
omega = np.sqrt(eigvals_sorted)  # 圆频率
f = omega / (2 * np.pi)  # 频率
T = 2 * np.pi / omega  # 周期
phi = eigvecs[:, idx] if 'idx' in dir() else eigvecs

print(f"结构模态参数:")
print(f"{'模态':<6}{'频率(Hz)':<12}{'周期(s)':<12}{'圆频率(rad/s)':<15}")
print("-" * 50)
for i in range(min(5, n_dof)):
    print(f"{i+1:<6}{f[i]:<12.3f}{T[i]:<12.3f}{omega[i]:<15.3f}")

# ============================================================================
# 2. 瑞利阻尼参数识别 - 两阶模态法
# ============================================================================
print("\n\n【2. 瑞利阻尼参数识别 - 两阶模态法】")
print("=" * 70)

def identify_rayleigh_damping(omega_i, omega_j, xi_i, xi_j):
    """
    根据两阶模态的阻尼比识别瑞利阻尼系数
    
    参数:
        omega_i, omega_j: 两阶模态的圆频率
        xi_i, xi_j: 两阶模态的阻尼比
    
    返回:
        alpha, beta: 瑞利阻尼系数
    """
    # 构建方程组
    # xi_i = alpha/(2*omega_i) + beta*omega_i/2
    # xi_j = alpha/(2*omega_j) + beta*omega_j/2
    
    A = np.array([[1/(2*omega_i), omega_i/2],
                  [1/(2*omega_j), omega_j/2]])
    b = np.array([xi_i, xi_j])
    
    solution = np.linalg.solve(A, b)
    alpha = solution[0]
    beta = solution[1]
    
    return alpha, beta

def rayleigh_damping_ratio(omega, alpha, beta):
    """
    计算瑞利阻尼在各频率下的阻尼比
    """
    return alpha / (2 * omega) + beta * omega / 2

# 情况1：基于第1阶和第3阶模态（控制低阶模态）
mode_i, mode_j = 0, 2  # 第1阶和第3阶（0-indexed）
xi_target = 0.05  # 目标阻尼比 5%

alpha1, beta1 = identify_rayleigh_damping(omega[mode_i], omega[mode_j], xi_target, xi_target)

print(f"\n情况1：基于第{mode_i+1}阶和第{mode_j+1}阶模态")
print(f"  目标阻尼比: {xi_target*100:.1f}%")
print(f"  识别得到的瑞利阻尼系数:")
print(f"    alpha = {alpha1:.6f}")
print(f"    beta  = {beta1:.8f}")

# 计算各阶模态的阻尼比
xi_rayleigh1 = rayleigh_damping_ratio(omega, alpha1, beta1)

print(f"\n  各阶模态阻尼比:")
print(f"{'模态':<6}{'目标阻尼比':<15}{'瑞利阻尼':<15}{'误差(%)':<10}")
print("-" * 50)
for i in range(min(5, n_dof)):
    target = xi_target if i in [mode_i, mode_j] else 0
    error = (xi_rayleigh1[i] - target) / xi_target * 100 if target > 0 else 0
    print(f"{i+1:<6}{target*100:<15.2f}{xi_rayleigh1[i]*100:<15.2f}{error:<10.2f}")

# 情况2：基于第1阶和第5阶模态（覆盖更宽频带）
mode_i, mode_j = 0, 4
alpha2, beta2 = identify_rayleigh_damping(omega[mode_i], omega[mode_j], xi_target, xi_target)
xi_rayleigh2 = rayleigh_damping_ratio(omega, alpha2, beta2)

print(f"\n情况2：基于第{mode_i+1}阶和第{mode_j+1}阶模态")
print(f"  alpha = {alpha2:.6f}")
print(f"  beta  = {beta2:.8f}")

# 情况3：基于第2阶和第4阶模态
mode_i, mode_j = 1, 3
alpha3, beta3 = identify_rayleigh_damping(omega[mode_i], omega[mode_j], xi_target, xi_target)
xi_rayleigh3 = rayleigh_damping_ratio(omega, alpha3, beta3)

print(f"\n情况3：基于第{mode_i+1}阶和第{mode_j+1}阶模态")
print(f"  alpha = {alpha3:.6f}")
print(f"  beta  = {beta3:.8f}")

# ============================================================================
# 3. 最小二乘法拟合多阶模态
# ============================================================================
print("\n\n【3. 最小二乘法拟合多阶模态】")
print("=" * 70)

def identify_rayleigh_least_squares(omega_modes, xi_targets):
    """
    使用最小二乘法识别瑞利阻尼系数
    
    参数:
        omega_modes: 模态频率数组
        xi_targets: 目标阻尼比数组
    
    返回:
        alpha, beta: 瑞利阻尼系数
    """
    # 确保输入是有效的数值
    omega_modes = np.array(omega_modes)
    xi_targets = np.array(xi_targets)
    
    # 过滤掉无效值
    valid_mask = (omega_modes > 1e-6) & np.isfinite(omega_modes) & np.isfinite(xi_targets)
    omega_valid = omega_modes[valid_mask]
    xi_valid = xi_targets[valid_mask]
    
    if len(omega_valid) < 2:
        # 如果数据不足，返回默认值
        return 0.0, 0.0, np.array([0.0])
    
    # 构建设计矩阵
    # xi = alpha/(2*omega) + beta*omega/2
    A = np.column_stack([1/(2*omega_valid), omega_valid/2])
    
    # 使用正规方程求解 (更稳定)
    AtA = A.T @ A
    Atb = A.T @ xi_valid
    
    try:
        solution = np.linalg.solve(AtA, Atb)
        residuals = np.sum((A @ solution - xi_valid)**2)
    except np.linalg.LinAlgError:
        # 如果矩阵奇异，使用伪逆
        solution = np.linalg.pinv(AtA) @ Atb
        residuals = np.sum((A @ solution - xi_valid)**2)
    
    alpha = solution[0]
    beta = solution[1]
    
    return alpha, beta, np.array([residuals])

# 使用前5阶模态进行拟合
n_modes_fit = 5
omega_fit = omega[:n_modes_fit]
xi_targets_fit = np.ones(n_modes_fit) * xi_target

alpha_ls, beta_ls, residuals = identify_rayleigh_least_squares(omega_fit, xi_targets_fit)
xi_rayleigh_ls = rayleigh_damping_ratio(omega, alpha_ls, beta_ls)

print(f"使用最小二乘法拟合前{n_modes_fit}阶模态")
print(f"  alpha = {alpha_ls:.6f}")
print(f"  beta  = {beta_ls:.8f}")
print(f"  残差平方和: {residuals[0]:.10f}" if len(residuals) > 0 else "  残差: 极小")

print(f"\n  拟合结果对比:")
print(f"{'模态':<6}{'目标':<12}{'拟合值':<12}{'误差(%)':<10}")
print("-" * 45)
for i in range(n_modes_fit):
    error = (xi_rayleigh_ls[i] - xi_targets_fit[i]) / xi_targets_fit[i] * 100
    print(f"{i+1:<6}{xi_targets_fit[i]*100:<12.2f}{xi_rayleigh_ls[i]*100:<12.2f}{error:<10.2f}")

# ============================================================================
# 4. 阻尼比随频率变化分析
# ============================================================================
print("\n\n【4. 阻尼比随频率变化分析】")
print("=" * 70)

# 生成连续频率范围
omega_continuous = np.linspace(omega[0]*0.5, omega[-1]*1.5, 500)

xi_continuous1 = rayleigh_damping_ratio(omega_continuous, alpha1, beta1)
xi_continuous2 = rayleigh_damping_ratio(omega_continuous, alpha2, beta2)
xi_continuous3 = rayleigh_damping_ratio(omega_continuous, alpha3, beta3)
xi_continuous_ls = rayleigh_damping_ratio(omega_continuous, alpha_ls, beta_ls)

# 分析高频和低频行为
print(f"\n瑞利阻尼的频率特性:")
print(f"{'方法':<20}{'低频阻尼比(%)':<15}{'高频阻尼比(%)':<15}{'最小阻尼比(%)':<15}")
print("-" * 70)

methods = [
    ("情况1(1,3阶)", alpha1, beta1),
    ("情况2(1,5阶)", alpha2, beta2),
    ("情况3(2,4阶)", alpha3, beta3),
    ("最小二乘", alpha_ls, beta_ls)
]

for name, alpha, beta in methods:
    xi_low = rayleigh_damping_ratio(omega[0]*0.5, alpha, beta) * 100
    xi_high = rayleigh_damping_ratio(omega[-1]*1.5, alpha, beta) * 100
    
    # 找到最小阻尼比
    xi_cont = rayleigh_damping_ratio(omega_continuous, alpha, beta)
    xi_min = np.min(xi_cont) * 100
    
    print(f"{name:<20}{xi_low:<15.2f}{xi_high:<15.2f}{xi_min:<15.2f}")

# ============================================================================
# 5. 可视化
# ============================================================================
print("\n\n【5. 生成可视化图形】")
print("=" * 70)

fig = plt.figure(figsize=(16, 12))

# 子图1: 阻尼比-频率曲线（情况1）
ax1 = fig.add_subplot(3, 3, 1)
ax1.plot(omega_continuous/(2*np.pi), xi_continuous1*100, 'b-', linewidth=2, label='Rayleigh Damping')
ax1.scatter(f[:5], xi_rayleigh1[:5]*100, c='red', s=100, zorder=5, label='Modal Damping')
ax1.axhline(y=xi_target*100, color='gray', linestyle='--', alpha=0.5, label='Target 5%')
ax1.set_xlabel('Frequency (Hz)', fontsize=10)
ax1.set_ylabel('Damping Ratio (%)', fontsize=10)
ax1.set_title('Case 1: Based on Mode 1 & 3', fontsize=11, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_ylim([0, 15])

# 子图2: 阻尼比-频率曲线（情况2）
ax2 = fig.add_subplot(3, 3, 2)
ax2.plot(omega_continuous/(2*np.pi), xi_continuous2*100, 'b-', linewidth=2, label='Rayleigh Damping')
ax2.scatter(f[:5], xi_rayleigh2[:5]*100, c='red', s=100, zorder=5, label='Modal Damping')
ax2.axhline(y=xi_target*100, color='gray', linestyle='--', alpha=0.5, label='Target 5%')
ax2.set_xlabel('Frequency (Hz)', fontsize=10)
ax2.set_ylabel('Damping Ratio (%)', fontsize=10)
ax2.set_title('Case 2: Based on Mode 1 & 5', fontsize=11, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_ylim([0, 15])

# 子图3: 阻尼比-频率曲线（情况3）
ax3 = fig.add_subplot(3, 3, 3)
ax3.plot(omega_continuous/(2*np.pi), xi_continuous3*100, 'b-', linewidth=2, label='Rayleigh Damping')
ax3.scatter(f[:5], xi_rayleigh3[:5]*100, c='red', s=100, zorder=5, label='Modal Damping')
ax3.axhline(y=xi_target*100, color='gray', linestyle='--', alpha=0.5, label='Target 5%')
ax3.set_xlabel('Frequency (Hz)', fontsize=10)
ax3.set_ylabel('Damping Ratio (%)', fontsize=10)
ax3.set_title('Case 3: Based on Mode 2 & 4', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_ylim([0, 15])

# 子图4: 最小二乘拟合结果
ax4 = fig.add_subplot(3, 3, 4)
ax4.plot(omega_continuous/(2*np.pi), xi_continuous_ls*100, 'b-', linewidth=2, label='Rayleigh Damping')
ax4.scatter(f[:n_modes_fit], xi_rayleigh_ls[:n_modes_fit]*100, c='red', s=100, zorder=5, label='Fitted Modes')
ax4.axhline(y=xi_target*100, color='gray', linestyle='--', alpha=0.5, label='Target 5%')
ax4.set_xlabel('Frequency (Hz)', fontsize=10)
ax4.set_ylabel('Damping Ratio (%)', fontsize=10)
ax4.set_title('Least Squares Fit (5 Modes)', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)
ax4.set_ylim([0, 15])

# 子图5: 所有方法对比
ax5 = fig.add_subplot(3, 3, 5)
ax5.plot(omega_continuous/(2*np.pi), xi_continuous1*100, 'b-', linewidth=1.5, alpha=0.7, label='Case 1')
ax5.plot(omega_continuous/(2*np.pi), xi_continuous2*100, 'r-', linewidth=1.5, alpha=0.7, label='Case 2')
ax5.plot(omega_continuous/(2*np.pi), xi_continuous3*100, 'g-', linewidth=1.5, alpha=0.7, label='Case 3')
ax5.plot(omega_continuous/(2*np.pi), xi_continuous_ls*100, 'm-', linewidth=2, label='Least Squares')
ax5.axhline(y=xi_target*100, color='gray', linestyle='--', alpha=0.5)
ax5.set_xlabel('Frequency (Hz)', fontsize=10)
ax5.set_ylabel('Damping Ratio (%)', fontsize=10)
ax5.set_title('Comparison of All Methods', fontsize=11, fontweight='bold')
ax5.legend(fontsize=9)
ax5.grid(True, alpha=0.3)
ax5.set_ylim([0, 15])

# 子图6: 各阶模态阻尼比对比（柱状图）
ax6 = fig.add_subplot(3, 3, 6)
x_pos = np.arange(1, 6)
width = 0.2
ax6.bar(x_pos - 1.5*width, xi_rayleigh1[:5]*100, width, label='Case 1', alpha=0.8)
ax6.bar(x_pos - 0.5*width, xi_rayleigh2[:5]*100, width, label='Case 2', alpha=0.8)
ax6.bar(x_pos + 0.5*width, xi_rayleigh3[:5]*100, width, label='Case 3', alpha=0.8)
ax6.bar(x_pos + 1.5*width, xi_rayleigh_ls[:5]*100, width, label='LS', alpha=0.8)
ax6.axhline(y=xi_target*100, color='gray', linestyle='--', alpha=0.5, label='Target')
ax6.set_xlabel('Mode Number', fontsize=10)
ax6.set_ylabel('Damping Ratio (%)', fontsize=10)
ax6.set_title('Damping Ratio by Mode', fontsize=11, fontweight='bold')
ax6.set_xticks(x_pos)
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3, axis='y')

# 子图7: 阻尼矩阵可视化（情况1）
ax7 = fig.add_subplot(3, 3, 7)
C1 = alpha1 * M + beta1 * K
im1 = ax7.imshow(C1/1e6, cmap='coolwarm', aspect='auto')
ax7.set_title('Damping Matrix C (Case 1)', fontsize=11, fontweight='bold')
ax7.set_xlabel('DOF', fontsize=10)
ax7.set_ylabel('DOF', fontsize=10)
plt.colorbar(im1, ax=ax7, label='MN.s/m')

# 子图8: 瑞利阻尼系数对比
ax8 = fig.add_subplot(3, 3, 8)
methods_names = ['Case 1', 'Case 2', 'Case 3', 'LS']
alphas = [alpha1, alpha2, alpha3, alpha_ls]
betas = [beta1*1e4, beta2*1e4, beta3*1e4, beta_ls*1e4]  # 缩放便于显示

x_pos = np.arange(len(methods_names))
width = 0.35
ax8_twin = ax8.twinx()
bars1 = ax8.bar(x_pos - width/2, alphas, width, label='Alpha', color='steelblue', alpha=0.8)
bars2 = ax8_twin.bar(x_pos + width/2, betas, width, label='Beta (x1e-4)', color='coral', alpha=0.8)
ax8.set_xlabel('Method', fontsize=10)
ax8.set_ylabel('Alpha', fontsize=10, color='steelblue')
ax8_twin.set_ylabel('Beta (x1e-4)', fontsize=10, color='coral')
ax8.set_title('Rayleigh Coefficients', fontsize=11, fontweight='bold')
ax8.set_xticks(x_pos)
ax8.set_xticklabels(methods_names)
ax8.tick_params(axis='y', labelcolor='steelblue')
ax8_twin.tick_params(axis='y', labelcolor='coral')
ax8.grid(True, alpha=0.3, axis='y')

# 子图9: 拟合误差分析
ax9 = fig.add_subplot(3, 3, 9)
errors1 = np.abs(xi_rayleigh1[:5] - xi_target) / xi_target * 100
errors2 = np.abs(xi_rayleigh2[:5] - xi_target) / xi_target * 100
errors3 = np.abs(xi_rayleigh3[:5] - xi_target) / xi_target * 100
errors_ls = np.abs(xi_rayleigh_ls[:5] - xi_target) / xi_target * 100

x_pos_modes = np.arange(1, 6)
ax9.semilogy(x_pos_modes, errors1, 'bo-', linewidth=2, markersize=8, label='Case 1')
ax9.semilogy(x_pos_modes, errors2, 'rs-', linewidth=2, markersize=8, label='Case 2')
ax9.semilogy(x_pos_modes, errors3, 'g^-', linewidth=2, markersize=8, label='Case 3')
ax9.semilogy(x_pos_modes, errors_ls, 'md-', linewidth=2, markersize=8, label='LS')
ax9.set_xlabel('Mode Number', fontsize=10)
ax9.set_ylabel('Relative Error (%)', fontsize=10)
ax9.set_title('Fitting Error', fontsize=11, fontweight='bold')
ax9.set_xticks(x_pos_modes)
ax9.legend(fontsize=9)
ax9.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rayleigh_damping_identification.png', dpi=150, bbox_inches='tight')
print("✓ 瑞利阻尼参数识别图已保存")
plt.close()

# ============================================================================
# 6. 自由振动衰减验证
# ============================================================================
print("\n\n【6. 自由振动衰减验证】")
print("=" * 70)

def free_vibration_sdof(m, k, xi, u0, v0, dt, T):
    """
    单自由度系统自由振动响应
    """
    n_steps = int(T / dt)
    t = np.linspace(0, T, n_steps)
    
    omega_n = np.sqrt(k / m)
    omega_d = omega_n * np.sqrt(1 - xi**2)
    
    # 初始条件响应
    A = np.sqrt(u0**2 + ((v0 + xi*omega_n*u0)/omega_d)**2)
    phi = np.arctan2(u0*omega_d, v0 + xi*omega_n*u0)
    
    u = A * np.exp(-xi*omega_n*t) * np.sin(omega_d*t + phi)
    
    return t, u

# 使用识别得到的阻尼比进行验证
m_test = 1000
k_test = 4 * np.pi**2 * m_test  # 周期为1s

# 不同方法的阻尼比（取第1阶模态）
xis_test = [xi_rayleigh1[0], xi_rayleigh2[0], xi_rayleigh3[0], xi_rayleigh_ls[0]]
labels_test = ['Case 1', 'Case 2', 'Case 3', 'LS']

print(f"\n自由振动衰减验证 (第1阶模态):")
print(f"{'方法':<12}{'阻尼比(%)':<12}{'理论衰减率':<15}")
print("-" * 45)

fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes = axes.flatten()

for idx, (xi, label) in enumerate(zip(xis_test, labels_test)):
    t_resp, u_resp = free_vibration_sdof(m_test, k_test, xi, 0.1, 0, 0.01, 10)
    
    # 计算对数衰减率
    peaks_idx = []
    for i in range(1, len(u_resp)-1):
        if u_resp[i] > u_resp[i-1] and u_resp[i] > u_resp[i+1] and u_resp[i] > 0:
            peaks_idx.append(i)
    
    if len(peaks_idx) >= 2:
        delta = np.log(u_resp[peaks_idx[0]] / u_resp[peaks_idx[1]])
        delta_theory = 2 * np.pi * xi / np.sqrt(1 - xi**2)
    else:
        delta = 0
        delta_theory = 0
    
    print(f"{label:<12}{xi*100:<12.2f}{delta_theory:<15.4f}")
    
    axes[idx].plot(t_resp, u_resp*1000, 'b-', linewidth=1.5)
    axes[idx].set_xlabel('Time (s)', fontsize=10)
    axes[idx].set_ylabel('Displacement (mm)', fontsize=10)
    axes[idx].set_title(f'{label}: ξ={xi*100:.2f}%', fontsize=11, fontweight='bold')
    axes[idx].grid(True, alpha=0.3)
    
    # 标记峰值
    if len(peaks_idx) > 0:
        axes[idx].scatter(t_resp[peaks_idx[:5]], u_resp[peaks_idx[:5]]*1000, 
                         c='red', s=50, zorder=5)

plt.tight_layout()
plt.savefig('free_vibration_decay.png', dpi=150, bbox_inches='tight')
print("\n✓ 自由振动衰减验证图已保存")
plt.close()

# ============================================================================
# 7. 结果总结
# ============================================================================
print("\n\n" + "=" * 70)
print("案例1总结")
print("=" * 70)

print(f"""
【瑞利阻尼参数识别总结】

1. 结构模型
   - 10层剪切型框架结构
   - 基本周期: {T[0]:.3f} s
   - 基本频率: {f[0]:.3f} Hz

2. 瑞利阻尼参数识别结果
   
   方法1(基于第1,3阶模态):
   - alpha = {alpha1:.6f}
   - beta  = {beta1:.8f}
   
   方法2(基于第1,5阶模态):
   - alpha = {alpha2:.6f}
   - beta  = {beta2:.8f}
   
   方法3(基于第2,4阶模态):
   - alpha = {alpha3:.6f}
   - beta  = {beta3:.8f}
   
   最小二乘法(拟合前5阶):
   - alpha = {alpha_ls:.6f}
   - beta  = {beta_ls:.8f}

3. 关键发现
   - 选择不同的参考模态，得到的瑞利阻尼系数差异显著
   - 基于低频模态的识别，高频阻尼比偏大
   - 基于高频模态的识别，低频阻尼比偏大
   - 最小二乘法可以平衡各阶模态的拟合误差

4. 瑞利阻尼的特性
   - 阻尼比随频率呈"V"字形变化
   - 在参考频率处阻尼比等于目标值
   - 远离参考频率时，阻尼比偏离目标值

5. 工程应用建议
   - 选择对响应贡献最大的模态作为参考
   - 对于地震分析，通常选择第1阶和第3阶
   - 对于风振分析，可能需要考虑更高阶模态
   - 当需要拟合多阶模态时，使用最小二乘法

6. 注意事项
   - 瑞利阻尼是简化模型，不能精确拟合所有模态
   - 高频模态的阻尼比可能过大，导致响应被过度抑制
   - 应根据实际结构的阻尼特性测试结果进行校准

【关键结论】
瑞利阻尼模型通过质量比例和刚度比例的线性组合来近似结构的阻尼特性。
参数识别时应根据分析目的选择合适的参考模态，并认识到模型的局限性。
""")

print("\n" + "=" * 70)
print("案例1完成！")
print("=" * 70)

代码说明

1. 参数设置区：定义系统的质量、刚度、阻尼等基本参数
2. 核心计算模块：实现振动方程的求解算法
3. 可视化模块：生成分析图表和动画
4. 结果输出：保存图片文件供文档使用

运行上述代码将生成本主题所需的所有可视化素材。

总结

阻尼模型是结构动力分析中的关键组成部分。本主题系统介绍了：

1. 瑞利阻尼模型：简单实用，适用于大多数工程分析
2. 模态阻尼模型：物理意义明确，便于模态叠加分析
3. 非比例阻尼系统：处理复杂阻尼机制的先进方法
4. 参数识别方法：从实验数据获取阻尼参数的实用技术

在实际工程应用中，应根据结构类型、分析目的和精度要求选择合适的阻尼模型，并通过实验或经验数据验证阻尼参数的合理性。