主题027：动力凝聚与静力凝聚技术

1. 引言

在结构动力学有限元分析中，随着模型规模的增大，计算成本急剧上升。凝聚技术（Condensation/Reduction）是一种有效的模型降阶方法，通过减少自由度数量来降低计算复杂度，同时保持足够的精度。本章将深入探讨两种主要的凝聚技术：静力凝聚（Static Condensation/Guyan Reduction）和动力凝聚（Dynamic Condensation），分析它们的理论基础、实现方法及工程应用。

1.1 凝聚技术的工程背景

现代工程结构越来越复杂，有限元模型的自由度数可达数百万甚至上千万。凝聚技术的重要性体现在：

• 计算效率：大幅降低自由度数量，减少计算时间和内存需求
• 子结构分析：支持模块化建模和并行计算
• 试验-分析相关性：便于与实验模态分析结果对比
• 设计优化：加速迭代计算过程
  
  
  
  
  
  

1.2 凝聚技术分类

技术类型	特点	适用场景
静力凝聚	基于静力假设，无频率依赖	低频模态、静力分析
动力凝聚	考虑惯性效应，频率依赖	宽频带分析、高精度要求
模态凝聚	基于模态坐标变换	模态分析、响应计算
系统等价降阶	保持系统特性	控制系统设计

2. 静力凝聚理论

2.1 基本概念

静力凝聚由Guyan于1965年提出，又称Guyan减缩。其核心思想是将自由度分为：

• 主自由度（Master DOFs）：保留的关键自由度
• 从自由度（Slave DOFs）：被凝聚消除的自由度

2.2 数学推导

考虑无阻尼自由振动方程：

$$\mathbf{Ku}={\omega }^2\mathbf{Mu}$$

将自由度分块：

$$\left[\begin{array}{cc}{\mathbf{K}}_{mm}& {\mathbf{K}}_{ms}\\ {\mathbf{K}}_{sm}& {\mathbf{K}}_{ss}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_m\\ {\mathbf{u}}_s\end{array}\right]={\omega }^2\left[\begin{array}{cc}{\mathbf{M}}_{mm}& {\mathbf{M}}_{ms}\\ {\mathbf{M}}_{sm}& {\mathbf{M}}_{ss}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_m\\ {\mathbf{u}}_s\end{array}\right]$$

静力凝聚假设从自由度与主自由度的关系由静力平衡决定：

$${\mathbf{K}}_{sm}{\mathbf{u}}_m+{\mathbf{K}}_{ss}{\mathbf{u}}_s=0$$

解得：

$${\mathbf{u}}_s=-{\mathbf{K}}_{ss}^{-1}{\mathbf{K}}_{sm}{\mathbf{u}}_m={\mathbf{T}}_{sm}{\mathbf{u}}_m$$

其中变换矩阵：

$${\mathbf{T}}_{sm}=-{\mathbf{K}}_{ss}^{-1}{\mathbf{K}}_{sm}$$

2.3 凝聚后的矩阵

定义完整的变换矩阵：

$$\mathbf{T}=\left[\begin{array}{c}\mathbf{I}\\ {\mathbf{T}}_{sm}\end{array}\right]$$

凝聚后的刚度矩阵和质量矩阵：

$${\mathbf{K}}_R={\mathbf{T}}^T\mathbf{KT}$$

$${\mathbf{M}}_R={\mathbf{T}}^T\mathbf{MT}$$

2.4 静力凝聚的精度

静力凝聚的精度取决于：

• 主自由度选择：应包含结构的主要变形模式
• 频率范围：在低频范围内精度较高
• 模态截断：高阶模态被忽略

误差主要来源于忽略了从自由度的惯性效应。

3. 动力凝聚理论

3.1 基本思想

动力凝聚考虑了从自由度的惯性效应，通过迭代或解析方法获得更精确的凝聚结果。与静力凝聚不同，动力凝聚的变换矩阵是频率依赖的。

3.2 精确动力凝聚

从完整的特征值问题出发：

$$(\mathbf{K}-{\omega }^2\mathbf{M})\mathbf{u}=0$$

分块形式：

$$\left[\begin{array}{cc}{\mathbf{K}}_{mm}-{\omega }^2{\mathbf{M}}_{mm}& {\mathbf{K}}_{ms}-{\omega }^2{\mathbf{M}}_{ms}\\ {\mathbf{K}}_{sm}-{\omega }^2{\mathbf{M}}_{sm}& {\mathbf{K}}_{ss}-{\omega }^2{\mathbf{M}}_{ss}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_m\\ {\mathbf{u}}_s\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

精确的动力变换矩阵：

$${\mathbf{T}}_{sm}(\omega )=-({\mathbf{K}}_{ss}-{\omega }^2{\mathbf{M}}_{ss}{)}^{-1}({\mathbf{K}}_{sm}-{\omega }^2{\mathbf{M}}_{sm})$$

3.3 迭代动力凝聚

由于精确动力凝聚需要已知频率，实际中常采用迭代方法：

步骤1：使用静力凝聚获得初始估计

$${\mathbf{T}}^{(0)}={\mathbf{T}}_{static}$$

步骤2：求解凝聚后的特征值问题

$${\mathbf{K}}_R^{(k)}{\varphi }_R={\omega }^2{\mathbf{M}}_R^{(k)}{\varphi }_R$$

步骤3：更新变换矩阵

$${\mathbf{T}}_{sm}^{(k+1)}=-({\mathbf{K}}_{ss}-({\omega }^2{)}^{(k)}{\mathbf{M}}_{ss}{)}^{-1}({\mathbf{K}}_{sm}-({\omega }^2{)}^{(k)}{\mathbf{M}}_{sm})$$

步骤4：重复步骤2-3直至收敛

3.4 IRS方法（Improved Reduced System）

IRS方法是一种改进的动力凝聚方法，通过考虑一阶惯性修正来提高精度：

$${\mathbf{T}}_{IRS}={\mathbf{T}}_{static}+{\mathbf{T}}_{inertia}$$

其中惯性修正项：

$${\mathbf{T}}_{inertia}={\mathbf{K}}_{ss}^{-1}({\mathbf{M}}_{sm}+{\mathbf{M}}_{ss}{\mathbf{T}}_{static}){\mathbf{M}}_R^{-1}{\mathbf{K}}_R$$

4. 主自由度选择策略

4.1 选择原则

主自由度的选择直接影响凝聚精度：

• 几何重要性：结构的关键节点
• 模态参与：对目标模态贡献大的自由度
• 测量可达性：便于与实验数据对比
• 边界条件：约束和加载位置

4.2 自动选择算法

能量法：
选择对角元素较大的自由度：

$$E_i=\frac{K_{ii}}{M_{ii}}$$

模态动能法：
基于模态分析结果选择：

$$M{E}_i=\sum _{j=1}^n{\varphi }_{ij}^2{M}_{ii}$$

Guyan能量法：
综合考虑刚度和质量：

$$G{E}_i=\frac{K_{ii}^2}{K_{ii}+M_{ii}{\omega }_{max}^2}$$

4.3 选择数量建议

• 保留的自由度数应为目标模态数的2-3倍
• 通常保留总自由度的5%-20%
• 在应力集中区和几何不连续处增加保留点

5. 凝聚技术的数值实现

5.1 算法流程

静力凝聚算法：

1. 选择主自由度，重新排列矩阵
2. 计算变换矩阵 T = -K_ss^(-1) * K_sm
3. 凝聚刚度矩阵 K_R = T^T * K * T
4. 凝聚质量矩阵 M_R = T^T * M * T
5. 求解凝聚后的特征值问题
6. 展开从自由度 u_s = T * u_m

动力凝聚算法：

1. 初始化变换矩阵（静力凝聚）
2. 迭代直到收敛：
   a. 凝聚矩阵
   b. 求解特征值问题
   c. 更新变换矩阵
3. 输出结果

5.2 稀疏矩阵处理

对于大规模问题，应采用稀疏矩阵技术：

• 使用稀疏存储格式（CSR、CSC）
• 采用稀疏求解器（UMFPACK、SuperLU）
• 利用矩阵分块结构

5.3 并行计算

凝聚过程可以并行化：

• 子结构并行凝聚
• 矩阵运算并行化
• 迭代过程并行加速

6. 工程应用案例

案例1：静力凝聚原理演示

通过简单结构演示静力凝聚的基本原理和计算过程。

案例2：动力凝聚方法对比

对比静力凝聚、IRS方法和迭代动力凝聚的精度和效率。

案例3：凝聚技术对模态分析的影响

分析不同凝聚方法对固有频率和振型的影响。

案例4：大型结构凝聚应用

针对实际工程结构，演示凝聚技术的应用效果。

7. 凝聚技术的优缺点

7.1 优点

• 计算效率：大幅降低自由度数量
• 内存节省：减少存储需求
• 子结构方法：支持模块化分析
• 试验相关性：便于与实验对比

7.2 缺点

• 精度损失：特别是高阶模态
• 频率限制：静力凝聚仅适用于低频
• 主自由度选择：需要经验和技巧
• 预处理开销：凝聚过程需要计算时间

7.3 改进方向

• 自适应主自由度选择
• 多级凝聚技术
• 混合凝聚方法
• 机器学习辅助选择

8. 工程实践建议

8.1 方法选择

• 低频模态分析：静力凝聚足够
• 宽频带分析：使用动力凝聚或IRS
• 高精度要求：迭代动力凝聚
• 实时应用：预计算变换矩阵

8.2 主自由度选择

• 包含所有边界条件和加载点
• 在结构的关键位置增加保留点
• 考虑模态参与因子
• 进行敏感性分析验证

8.3 精度验证

• 与完整模型结果对比
• 检查模态正交性
• 验证质量守恒
• 进行收敛性分析

8.4 软件实现

主流有限元软件的凝聚技术：

• ANSYS：Matrix Reduction功能
• ABAQUS：Substructuring和Model Reduction
• NASTRAN：Guyan Reduction和Dynamic Reduction
• MATLAB：balred函数（控制系统工具箱）

9. 总结

凝聚技术是结构动力学分析中的重要工具，通过合理选择主自由度，可以在保持足够精度的前提下大幅降低计算成本。静力凝聚简单高效，适合低频分析；动力凝聚精度更高，适合宽频带分析。在实际工程中，应根据具体问题选择合适的方法，并进行充分的精度验证。

完整Python代码实现

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

import matplotlib
matplotlib.use('Agg')
"""
案例1：静力凝聚原理演示
通过简单结构演示静力凝聚的基本原理和计算过程
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eig, inv, solve
import matplotlib.patches as patches

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

print("=" * 70)
print("案例1：静力凝聚原理演示")
print("=" * 70)

# ============================================================================
# 1. 简单框架结构模型
# ============================================================================
print("\n【结构模型】")
print("-" * 50)

# 材料参数
E = 200e9  # 弹性模量 (Pa)
A = 0.01   # 截面积 (m²)
I = 1e-4   # 惯性矩 (m⁴)
rho = 7850  # 密度 (kg/m³)

# 几何参数
L = 3.0  # 柱高 (m)
W = 4.0  # 梁跨度 (m)

print(f"材料参数: E = {E/1e9:.1f} GPa, ρ = {rho} kg/m³")
print(f"截面参数: A = {A*1e4:.2f} cm², I = {I*1e8:.2f} cm⁴")
print(f"几何尺寸: 柱高 = {L} m, 跨度 = {W} m")

# ============================================================================
# 2. 建立完整有限元模型
# ============================================================================
print("\n【完整有限元模型】")
print("-" * 50)

# 节点坐标 (x, y)
# 节点0: (0, 0) - 左下
# 节点1: (W, 0) - 右下
# 节点2: (0, L) - 左上
# 节点3: (W, L) - 右上

nodes = np.array([
    [0, 0],
    [W, 0],
    [0, L],
    [W, L]
])

n_nodes = len(nodes)
dof_per_node = 2  # x, y位移
n_dof = n_nodes * dof_per_node

print(f"节点数: {n_nodes}")
print(f"每个节点自由度数: {dof_per_node}")
print(f"总自由度数: {n_dof}")

# 单元连接
# 单元0: 柱0 (节点0-节点2)
# 单元1: 柱1 (节点1-节点3)
# 单元2: 梁 (节点2-节点3)
elements = [
    [0, 2],  # 柱0
    [1, 3],  # 柱1
    [2, 3]   # 梁
]

# ============================================================================
# 3. 组装全局刚度矩阵和质量矩阵
# ============================================================================

def beam_element_2d(node_i, node_j, E, A, I):
    """
    2D梁单元刚度矩阵（简化版，仅考虑轴向和横向位移）
    """
    xi, yi = node_i
    xj, yj = node_j
    L_elem = np.sqrt((xj - xi)**2 + (yj - yi)**2)
    
    # 轴向刚度
    k_axial = E * A / L_elem
    # 弯曲刚度（简化）
    k_bending = 12 * E * I / L_elem**3
    
    # 局部坐标系下的刚度矩阵 (4x4: ui, vi, uj, vj)
    k_local = np.array([
        [k_axial, 0, -k_axial, 0],
        [0, k_bending, 0, -k_bending],
        [-k_axial, 0, k_axial, 0],
        [0, -k_bending, 0, k_bending]
    ])
    
    # 质量矩阵（集中质量）
    m_elem = rho * A * L_elem
    m_local = np.diag([m_elem/2, m_elem/2, m_elem/2, m_elem/2])
    
    return k_local, m_local, L_elem

# 初始化全局矩阵
K_global = np.zeros((n_dof, n_dof))
M_global = np.zeros((n_dof, n_dof))

# 组装单元
for elem_idx, (i, j) in enumerate(elements):
    k_elem, m_elem, L_elem = beam_element_2d(nodes[i], nodes[j], E, A, I)
    
    # 自由度映射
    dof_i = [i * 2, i * 2 + 1]
    dof_j = [j * 2, j * 2 + 1]
    dof_map = dof_i + dof_j
    
    # 组装到全局矩阵
    for ii, gi in enumerate(dof_map):
        for jj, gj in enumerate(dof_map):
            K_global[gi, gj] += k_elem[ii, jj]
            M_global[gi, gj] += m_elem[ii, jj]

print("\n全局刚度矩阵 K (前4x4):")
print(K_global[:4, :4])
print("\n全局质量矩阵 M (对角):")
print(np.diag(M_global)[:4])

# ============================================================================
# 4. 应用边界条件
# ============================================================================
print("\n【边界条件】")
print("-" * 50)

# 固定节点0和节点1（底部）
fixed_dofs = [0, 1, 2, 3]  # 节点0和节点1的x,y位移
free_dofs = [4, 5, 6, 7]   # 节点2和节点3的x,y位移

print(f"固定自由度: {fixed_dofs}")
print(f"自由自由度: {free_dofs}")

# 提取子矩阵
K_mm = K_global[np.ix_(free_dofs, free_dofs)]  # 主自由度
K_ms = K_global[np.ix_(free_dofs, fixed_dofs)]
K_sm = K_global[np.ix_(fixed_dofs, free_dofs)]
K_ss = K_global[np.ix_(fixed_dofs, fixed_dofs)]

M_mm = M_global[np.ix_(free_dofs, free_dofs)]
M_ms = M_global[np.ix_(free_dofs, fixed_dofs)]
M_sm = M_global[np.ix_(fixed_dofs, free_dofs)]
M_ss = M_global[np.ix_(fixed_dofs, fixed_dofs)]

print(f"\nK_mm 维度: {K_mm.shape}")
print(f"K_ss 维度: {K_ss.shape}")

# ============================================================================
# 5. 完整模型模态分析
# ============================================================================
print("\n【完整模型模态分析】")
print("-" * 50)

# 求解特征值问题
eigvals_full, eigvecs_full = eig(K_mm, M_mm)
omega_full = np.sqrt(np.real(eigvals_full))
f_full = omega_full / (2 * np.pi)

# 排序
idx = np.argsort(f_full)
f_full = f_full[idx]
eigvecs_full = eigvecs_full[:, idx]

print("完整模型固有频率 (Hz):")
for i in range(len(f_full)):
    print(f"  第{i+1}阶: {f_full[i]:.4f} Hz")

# ============================================================================
# 6. 静力凝聚
# ============================================================================
print("\n【静力凝聚计算】")
print("-" * 50)

# 计算变换矩阵 T = -K_ss^(-1) * K_sm
# 注意：在我们的例子中，固定自由度是约束，所以实际上K_ss对应约束部分
# 让我们重新定义：保留自由度 = 自由自由度，凝聚自由度 = 固定自由度

# 静力凝聚变换矩阵
# 对于约束边界，从自由度（固定）= 0，所以实际上不需要凝聚
# 让我们修改问题：假设我们只想保留节点2和3的y方向位移，凝聚x方向

# 重新定义主/从自由度
master_dofs = [5, 7]  # 节点2和节点3的y方向（竖向）
slave_dofs = [4, 6]   # 节点2和节点3的x方向（水平）

print(f"主自由度 (保留): {master_dofs}")
print(f"从自由度 (凝聚): {slave_dofs}")

# 重新提取子矩阵
K_mm_new = K_global[np.ix_(master_dofs, master_dofs)]
K_ms_new = K_global[np.ix_(master_dofs, slave_dofs)]
K_sm_new = K_global[np.ix_(slave_dofs, master_dofs)]
K_ss_new = K_global[np.ix_(slave_dofs, slave_dofs)]

M_mm_new = M_global[np.ix_(master_dofs, master_dofs)]
M_ms_new = M_global[np.ix_(master_dofs, slave_dofs)]
M_sm_new = M_global[np.ix_(slave_dofs, master_dofs)]
M_ss_new = M_global[np.ix_(slave_dofs, slave_dofs)]

# 计算变换矩阵
T_sm = -np.linalg.solve(K_ss_new, K_sm_new)
print(f"\n变换矩阵 T_sm ({T_sm.shape}):")
print(T_sm)

# 完整变换矩阵
T = np.vstack([np.eye(len(master_dofs)), T_sm])
print(f"\n完整变换矩阵 T ({T.shape}):")
print(T)

# 凝聚后的矩阵
K_reduced = T.T @ K_global[np.ix_(master_dofs + slave_dofs, master_dofs + slave_dofs)] @ T
M_reduced = T.T @ M_global[np.ix_(master_dofs + slave_dofs, master_dofs + slave_dofs)] @ T

print(f"\n凝聚后刚度矩阵 K_R ({K_reduced.shape}):")
print(K_reduced)
print(f"\n凝聚后质量矩阵 M_R ({M_reduced.shape}):")
print(M_reduced)

# ============================================================================
# 7. 凝聚后模态分析
# ============================================================================
print("\n【凝聚后模态分析】")
print("-" * 50)

eigvals_reduced, eigvecs_reduced = eig(K_reduced, M_reduced)
omega_reduced = np.sqrt(np.real(eigvals_reduced))
f_reduced = omega_reduced / (2 * np.pi)

idx_r = np.argsort(f_reduced)
f_reduced = f_reduced[idx_r]
eigvecs_reduced = eigvecs_reduced[:, idx_r]

print("凝聚后固有频率 (Hz):")
for i in range(len(f_reduced)):
    print(f"  第{i+1}阶: {f_reduced[i]:.4f} Hz")

# ============================================================================
# 8. 误差分析
# ============================================================================
print("\n【误差分析】")
print("-" * 50)

print("频率对比与误差:")
print("-" * 40)
print(f"{'阶数':<6}{'完整模型(Hz)':<18}{'凝聚后(Hz)':<18}{'误差(%)':<12}")
print("-" * 40)

for i in range(min(len(f_full), len(f_reduced))):
    error = abs(f_full[i] - f_reduced[i]) / f_full[i] * 100
    print(f"{i+1:<6}{f_full[i]:<18.4f}{f_reduced[i]:<18.4f}{error:<12.2f}")

# ============================================================================
# 9. 可视化
# ============================================================================
print("\n【生成可视化图形】")
print("-" * 50)

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

# 图1: 结构模型
ax1 = fig.add_subplot(2, 3, 1)
# 绘制单元
for i, j in elements:
    ax1.plot([nodes[i, 0], nodes[j, 0]], [nodes[i, 1], nodes[j, 1]], 'b-', linewidth=3)
# 绘制节点
ax1.scatter(nodes[:, 0], nodes[:, 1], c='red', s=200, zorder=5)
# 标注节点
for i, (x, y) in enumerate(nodes):
    ax1.annotate(f'{i}', (x, y), textcoords="offset points", xytext=(0, 10), ha='center', fontsize=12)
# 标注自由度
dof_labels = ['u₀', 'v₀', 'u₁', 'v₁', 'u₂', 'v₂', 'u₃', 'v₃']
for i, (x, y) in enumerate(nodes):
    ax1.annotate(dof_labels[i*2], (x-0.3, y), fontsize=9, color='green')
    ax1.annotate(dof_labels[i*2+1], (x, y-0.3), fontsize=9, color='green')
ax1.set_xlim(-1, W+1)
ax1.set_ylim(-0.5, L+0.5)
ax1.set_aspect('equal')
ax1.set_xlabel('X (m)')
ax1.set_ylabel('Y (m)')
ax1.set_title('结构模型与自由度编号')
ax1.grid(True, alpha=0.3)

# 图2: 主/从自由度选择
ax2 = fig.add_subplot(2, 3, 2)
for i, j in elements:
    ax2.plot([nodes[i, 0], nodes[j, 0]], [nodes[i, 1], nodes[j, 1]], 'b-', linewidth=3)
# 主自由度（保留）- 绿色
ax2.scatter(nodes[2, 0], nodes[2, 1], c='green', s=300, marker='o', label='主自由度', zorder=5)
ax2.scatter(nodes[3, 0], nodes[3, 1], c='green', s=300, marker='o', zorder=5)
# 从自由度（凝聚）- 橙色
ax2.scatter(nodes[2, 0], nodes[2, 1], c='orange', s=150, marker='s', label='从自由度', zorder=6)
ax2.scatter(nodes[3, 0], nodes[3, 1], c='orange', s=150, marker='s', zorder=6)
# 固定节点
ax2.scatter(nodes[:2, 0], nodes[:2, 1], c='red', s=200, marker='^', label='固定', zorder=5)
ax2.set_xlim(-1, W+1)
ax2.set_ylim(-0.5, L+0.5)
ax2.set_aspect('equal')
ax2.set_xlabel('X (m)')
ax2.set_ylabel('Y (m)')
ax2.set_title('主/从自由度选择')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 图3: 频率对比
ax3 = fig.add_subplot(2, 3, 3)
# 只对比前2阶（因为凝聚后只有2阶）
n_compare = min(len(f_full), len(f_reduced))
x_pos = np.arange(n_compare)
width = 0.35
bars1 = ax3.bar(x_pos - width/2, f_full[:n_compare], width, label='完整模型', color='steelblue')
bars2 = ax3.bar(x_pos + width/2, f_reduced[:n_compare], width, label='静力凝聚', color='coral')
ax3.set_xlabel('模态阶数')
ax3.set_ylabel('频率 (Hz)')
ax3.set_title('固有频率对比')
ax3.set_xticks(x_pos)
ax3.set_xticklabels([f'{i+1}' for i in range(n_compare)])
ax3.legend()
ax3.grid(True, alpha=0.3, axis='y')

# 图4: 变换矩阵可视化
ax4 = fig.add_subplot(2, 3, 4)
im = ax4.imshow(T, cmap='RdBu_r', aspect='auto')
ax4.set_xlabel('主自由度')
ax4.set_ylabel('全部自由度')
ax4.set_title('静力凝聚变换矩阵')
ax4.set_xticks([0, 1])
ax4.set_xticklabels(['v₂', 'v₃'])
ax4.set_yticks([0, 1, 2, 3])
ax4.set_yticklabels(['v₂', 'v₃', 'u₂', 'u₃'])
plt.colorbar(im, ax=ax4)

# 图5: 振型对比
ax5 = fig.add_subplot(2, 3, 5)
# 绘制第一阶振型
mode_full = eigvecs_full[:, 0]
mode_reduced_expanded = np.zeros(4)
mode_reduced_expanded[0] = eigvecs_reduced[0, 0]  # v2
mode_reduced_expanded[1] = eigvecs_reduced[1, 0]  # v3
# 通过变换矩阵恢复从自由度
mode_reduced_expanded[2:] = T_sm @ eigvecs_reduced[:, 0]

x = np.arange(4)
width = 0.35
ax5.bar(x - width/2, mode_full/np.max(np.abs(mode_full)), width, label='完整模型', color='steelblue')
ax5.bar(x + width/2, mode_reduced_expanded/np.max(np.abs(mode_reduced_expanded)), width, label='凝聚后', color='coral')
ax5.set_xlabel('自由度')
ax5.set_ylabel('归一化幅值')
ax5.set_title('第一阶振型对比')
ax5.set_xticks(x)
ax5.set_xticklabels(['v₂', 'v₃', 'u₂', 'u₃'])
ax5.legend()
ax5.grid(True, alpha=0.3, axis='y')

# 图6: 误差分析
ax6 = fig.add_subplot(2, 3, 6)
errors = []
for i in range(min(len(f_full), len(f_reduced))):
    error = abs(f_full[i] - f_reduced[i]) / f_full[i] * 100
    errors.append(error)
ax6.plot(range(1, len(errors)+1), errors, 'o-', linewidth=2, markersize=8, color='darkred')
ax6.set_xlabel('模态阶数')
ax6.set_ylabel('频率误差 (%)')
ax6.set_title('静力凝聚频率误差')
ax6.grid(True, alpha=0.3)
ax6.set_xticks(range(1, len(errors)+1))

plt.tight_layout()
plt.savefig('static_condensation_demo.png', dpi=150, bbox_inches='tight')
print("✓ 静力凝聚原理演示图已保存")
plt.close()

# ============================================================================
# 10. 不同主自由度数量的影响
# ============================================================================
print("\n【不同主自由度数量的影响】")
print("-" * 50)

# 测试不同数量的主自由度
master_dof_options = [
    [5, 7],        # 2个主自由度
    [4, 5, 6, 7],  # 4个主自由度（全部）
]

results = []

for master_dofs_test in master_dof_options:
    slave_dofs_test = [d for d in free_dofs if d not in master_dofs_test]
    
    if len(slave_dofs_test) == 0:
        # 没有从自由度，就是完整模型
        K_r = K_mm
        M_r = M_mm
    else:
        # 提取子矩阵
        all_dofs = master_dofs_test + slave_dofs_test
        K_sub = K_global[np.ix_(all_dofs, all_dofs)]
        M_sub = M_global[np.ix_(all_dofs, all_dofs)]
        
        K_mm_sub = K_sub[:len(master_dofs_test), :len(master_dofs_test)]
        K_ms_sub = K_sub[:len(master_dofs_test), len(master_dofs_test):]
        K_sm_sub = K_sub[len(master_dofs_test):, :len(master_dofs_test)]
        K_ss_sub = K_sub[len(master_dofs_test):, len(master_dofs_test):]
        
        # 变换矩阵
        T_sm_sub = -np.linalg.solve(K_ss_sub, K_sm_sub)
        T_sub = np.vstack([np.eye(len(master_dofs_test)), T_sm_sub])
        
        # 凝聚矩阵
        K_r = T_sub.T @ K_sub @ T_sub
        M_r = T_sub.T @ M_sub @ T_sub
    
    # 求解特征值
    eigvals_r, _ = eig(K_r, M_r)
    f_r = np.sort(np.sqrt(np.real(eigvals_r))) / (2 * np.pi)
    
    results.append({
        'n_master': len(master_dofs_test),
        'frequencies': f_r[:2]  # 只取前两阶
    })
    
    print(f"\n主自由度数 = {len(master_dofs_test)}:")
    print(f"  第1阶频率: {f_r[0]:.4f} Hz")
    if len(f_r) > 1:
        print(f"  第2阶频率: {f_r[1]:.4f} Hz")

# 绘制对比图
fig2, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# 第一阶频率对比
ax1.bar(['完整模型', '2主DOF', '4主DOF'], 
        [f_full[0], results[0]['frequencies'][0], results[1]['frequencies'][0]],
        color=['steelblue', 'coral', 'lightgreen'])
ax1.set_ylabel('频率 (Hz)')
ax1.set_title('第一阶频率对比')
ax1.grid(True, alpha=0.3, axis='y')

# 第二阶频率对比
ax2.bar(['完整模型', '2主DOF', '4主DOF'], 
        [f_full[1], results[0]['frequencies'][1], results[1]['frequencies'][1]],
        color=['steelblue', 'coral', 'lightgreen'])
ax2.set_ylabel('频率 (Hz)')
ax2.set_title('第二阶频率对比')
ax2.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('dof_comparison.png', dpi=150, bbox_inches='tight')
print("\n✓ 主自由度数量影响图已保存")
plt.close()

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

print("""
【静力凝聚基本原理】
1. 将自由度分为主自由度（保留）和从自由度（凝聚）
2. 基于静力平衡建立主从关系: u_s = T * u_m
3. 变换矩阵: T = -K_ss^(-1) * K_sm
4. 凝聚后矩阵: K_R = T^T * K * T, M_R = T^T * M * T

【关键发现】
- 静力凝聚在低频范围内具有较高精度
- 主自由度选择直接影响凝聚精度
- 保留的自由度数应为目标模态数的2-3倍
- 静力凝聚忽略了从自由度的惯性效应

【误差来源】
- 惯性效应忽略：静力假设不考虑从自由度的质量
- 模态截断：高阶模态信息丢失
- 频率依赖性：静力变换矩阵与频率无关
""")

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


# ============================================
# 生成GIF动画
# ============================================

def create_animation():
    """创建仿真结果动画"""
    import matplotlib.pyplot as plt
    from matplotlib.animation import FuncAnimation
    import numpy as np
    
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.set_xlim(0, 10)
    ax.set_ylim(-2, 2)
    ax.set_xlabel('Time', fontsize=12)
    ax.set_ylabel('Response', fontsize=12)
    ax.set_title('Dynamic Response Animation', fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    line, = ax.plot([], [], 'b-', linewidth=2)
    
    def init():
        line.set_data([], [])
        return line,
    
    def update(frame):
        x = np.linspace(0, 10, 100)
        y = np.sin(x - frame * 0.2) * np.exp(-frame * 0.01)
        line.set_data(x, y)
        return line,
    
    anim = FuncAnimation(fig, update, init_func=init, frames=50, interval=100, blit=True)
    
    output_dir = os.path.dirname(os.path.abspath(__file__))
    anim.save(f'{output_dir}/simulation_animation.gif', writer='pillow', fps=10)
    print(f"动画已保存到: {output_dir}/simulation_animation.gif")
    plt.close()

if __name__ == '__main__':
    create_animation()

代码说明

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

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