主题093：高级仿真技术与综合案例分析

摘要

本主题作为结构动力学仿真系列的高级综合篇，系统性地整合前期所学的理论与方法，深入探讨复杂工程结构的高级仿真技术。内容涵盖复杂结构的模态分析技术、非线性动力响应分析方法、以及多物理场耦合仿真原理。通过三个典型工程案例——大跨度空间结构模态分析、高层建筑非线性地震响应、以及热-结构耦合动力分析——展示如何将基础理论应用于解决实际工程中的复杂动力学问题。本主题强调仿真流程的标准化、模型验证的重要性，以及工程判断在仿真分析中的关键作用，为读者提供从理论到实践的完整技术路径。

关键词

复杂结构模态分析；非线性动力响应；多物理场耦合；大跨度结构；热-结构耦合；仿真验证

---

1. 引言

结构动力学仿真技术经过近百年的发展，已经从简单的单自由度系统分析，演变为能够处理千万自由度复杂结构的强大工具。随着计算能力的指数级增长和算法的不断创新，现代结构动力学仿真面临的挑战已经从"能否计算"转变为"如何准确、高效、可靠地计算"。

在实际工程应用中，结构动力学问题往往具有以下特征：

复杂性：真实工程结构通常具有复杂的几何形状、非均匀的材料分布、以及多种连接形式。例如，大跨度空间结构可能包含数以万计的杆件和节点，高层建筑需要考虑楼板、剪力墙、框架的协同工作，核反应堆结构则涉及压力容器、管道、支吊架等多种部件的相互作用。

非线性：结构的非线性行为是工程实践中无法回避的问题。材料非线性表现为混凝土开裂、钢材屈服、复合材料损伤；几何非线性体现在大变形、大转动、屈曲失稳；边界非线性则包括接触、碰撞、摩擦等现象。这些非线性因素往往耦合出现，使得分析难度大幅增加。

多物理场耦合：现代工程结构常常工作在复杂的环境中，需要考虑多种物理场的相互作用。热-结构耦合在航天器再入、火灾安全评估中至关重要；流-固耦合是桥梁风振、海洋平台响应分析的核心；电磁-结构耦合则在电机、变压器等设备设计中不可忽视。

不确定性：结构参数、荷载条件、边界状态都存在不确定性。如何量化这些不确定性对结构响应的影响，进行可靠性评估和鲁棒性设计，是现代结构动力学仿真的重要课题。

本主题将通过三个综合案例，系统展示如何运用高级仿真技术解决上述复杂问题。每个案例都包含完整的建模流程、分析方法、结果验证和工程应用讨论，力求为读者提供可借鉴的技术路线和工程经验。

---

2. 复杂结构模态分析技术

2.1 模态分析的理论基础

模态分析是结构动力学的基础，其核心是求解特征值问题：

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

其中，$\mathbf{K}$ 是刚度矩阵，$\mathbf{M}$ 是质量矩阵，$\omega$ 是圆频率，$\mathbf{\varphi }$ 是振型向量。

对于复杂结构，直接求解上述特征值问题面临以下挑战：

计算规模：大型结构的自由度可能达到百万甚至千万量级，存储和操作如此大规模的矩阵对计算资源提出极高要求。

频谱密集：复杂结构往往具有密集的模态分布，特别是在高频段，相邻模态的频率间隔可能非常小，这对特征值求解算法的精度和稳定性提出挑战。

局部模态：结构中可能存在局部刚度突变或质量集中，导致出现局部振动模态。这些模态对整体响应贡献很小，但会影响模态分析的收敛性。

刚体模态：对于未完全约束的结构（如航天器、浮式平台），存在刚体模态（频率为零）。这些模态需要特殊处理，否则会干扰弹性模态的求解。

2.2 子空间迭代法

子空间迭代法是求解大型特征值问题的经典算法，其基本思想是在一个低维子空间中进行迭代，逐步逼近真实的特征向量。

算法步骤：

1. 初始化：选择 $q$ 个初始向量（$q>p$，$p$ 为所需模态数），构成初始子空间基矩阵 ${\mathbf{X}}_0$。

2. 迭代循环（$k=0,1,2,...$）：

   a. 逆迭代：求解线性方程组
   $$\mathbf{K}{\bar{\mathbf{X}}}_{k+1}=\mathbf{M}{\mathbf{X}}_k$$

   b. Rayleigh-Ritz分析：在子空间中求解投影特征值问题
   $$\tilde{\mathbf{K}}\mathbf{Q}=\tilde{\mathbf{M}}\mathbf{Q}\mathbf{\Lambda }$$
   其中，$\tilde{\mathbf{K}}={\bar{\mathbf{X}}}_{k+1}^T\mathbf{K}{\bar{\mathbf{X}}}_{k+1}$，$\tilde{\mathbf{M}}={\bar{\mathbf{X}}}_{k+1}^T\mathbf{M}{\bar{\mathbf{X}}}_{k+1}$。

   c. 更新近似特征向量
   $${\mathbf{X}}_{k+1}={\bar{\mathbf{X}}}_{k+1}\mathbf{Q}$$

   d. 收敛检验：检查特征值的相对变化是否小于容差。

3. 输出结果：收敛后的特征值和特征向量即为所求模态。

算法特点：

• 收敛速度与特征值间隔有关，间隔越大收敛越快
• 适合求解低频密集模态
• 计算成本主要来自每步的线性方程组求解

2.3 Lanczos方法

Lanczos方法是另一种高效的特征值求解算法，特别适用于求解大规模稀疏矩阵的极端特征值（最小或最大几个）。

基本思想：通过构造Krylov子空间的一组正交基，将大规模特征值问题转化为三对角矩阵的特征值问题。

算法流程：

1. 初始化：选择初始向量 ${\mathbf{q}}_1$（通常随机生成并归一化）

2. Lanczos迭代（$j=1,2,...,m$）：

   $${\mathbf{r}}_j={\mathbf{K}}^{-1}\mathbf{M}{\mathbf{q}}_j-{\alpha }_j{\mathbf{q}}_j-{\beta }_j{\mathbf{q}}_{j-1}$$

   $${\beta }_{j+1}=\sqrt{{\mathbf{r}}_j^T\mathbf{M}{\mathbf{r}}_j}$$

   $${\mathbf{q}}_{j+1}={\mathbf{r}}_j/{\beta }_{j+1}$$

   其中，${\alpha }_j={\mathbf{q}}_j^T\mathbf{M}{\mathbf{K}}^{-1}\mathbf{M}{\mathbf{q}}_j$。

3. 求解三对角矩阵特征值：
   构建三对角矩阵 ${\mathbf{T}}_m$，其对角线元素为 ${\alpha }_j$，次对角线元素为 ${\beta }_j$。求解 ${\mathbf{T}}_m$ 的特征值问题，得到Ritz值作为原问题的近似特征值。

数值稳定性：
Lanczos迭代在有限精度算术中会出现"伪收敛"现象，即多个向量可能收敛到同一个特征向量。常用的解决方案包括：

• 完全重正交化：每步都对新向量与所有先前向量进行正交化
• 选择性重正交化：仅当正交性损失超过阈值时才进行重正交化
• 分块Lanczos：同时处理多个向量，提高数值稳定性

2.4 模态验证与误差估计

获得模态结果后，必须进行严格的验证，确保结果的物理合理性和数值准确性。

物理验证：

1. 频率合理性检查：

   • 与相似结构或经验公式对比
   • 检查频率间隔是否合理（避免虚假模态）
   • 验证刚体模态（频率应接近零）

2. 振型合理性检查：

   • 观察振型是否光滑连续
   • 检查节点位置是否符合预期
   • 验证对称结构的振型对称性

3. 能量分布分析：
   计算各部件的动能和应变能分布：
   $$T_i=\frac{1}{2}{\omega }^2{\mathbf{\varphi }}_i^T{\mathbf{M}}_i{\mathbf{\varphi }}_i$$
   $$U_i=\frac{1}{2}{\mathbf{\varphi }}_i^T{\mathbf{K}}_i{\mathbf{\varphi }}_i$$
   能量分布应反映结构的刚度-质量特性。

数值验证：

1. 收敛性检验：
   通过网格加密或增加模态数，检查结果是否收敛。收敛的解应满足：
   $$\frac{|f_h-f_{h/2}|}{f_{h/2}}<ϵ$$
   其中，$f_h$ 是网格尺寸为 $h$ 时的频率，$ϵ$ 是收敛容差（通常取 1%）。

2. 正交性检验：
   验证模态关于质量矩阵和刚度矩阵的正交性：
   $${\mathbf{\varphi }}_i^T\mathbf{M}{\mathbf{\varphi }}_j={\delta }_{ij}$$
   $${\mathbf{\varphi }}_i^T\mathbf{K}{\mathbf{\varphi }}_j={\omega }_i^2{\delta }_{ij}$$
   正交性误差应小于 $10^{-6}$。

3. 残差检验：
   计算特征值问题的残差：
   $$\mathbf{r}=\mathbf{K}\mathbf{\varphi }-{\omega }^2\mathbf{M}\mathbf{\varphi }$$
   相对残差范数应满足：
   $$\frac{\Vert \mathbf{r}\Vert }{\Vert \mathbf{K}\mathbf{\varphi }\Vert }<10^{-8}$$

---

3. 非线性动力响应分析

3.1 非线性问题的分类与特点

结构非线性可分为三大类：

材料非线性：

• 弹塑性：应力超过屈服强度后产生永久变形，卸载时沿弹性路径返回
• 损伤：材料内部微裂纹扩展导致刚度退化，常见于混凝土、复合材料
• 粘弹性：应力-应变关系与时间相关，表现为蠕变和松弛

几何非线性：

• 大位移：变形后的几何构型与初始构型差异显著，需要在变形后构型上建立平衡
• 大应变：应变分量不再是小量，需要采用有限应变理论
• 跟随力：荷载方向随结构变形而改变（如风荷载、水压力）

边界非线性：

• 接触：两个或多个物体之间的相互作用，接触状态（接触/分离）随时间变化
• 摩擦：接触面上的切向力与相对滑动相关，存在粘着-滑动转变
• 间隙：连接中存在初始间隙，导致力-位移关系呈现分段特征

3.2 非线性方程的求解方法

非线性动力分析需要求解如下形式的非线性方程：

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

其中，${\mathbf{F}}^{int}$ 是内力向量，通常是位移的非线性函数。

显式积分法：

中心差分法是最常用的显式算法：

$${\ddot{\mathbf{u}}}_n={\mathbf{M}}^{-1}({\mathbf{F}}_n^{ext}-{\mathbf{F}}_n^{int}-\mathbf{C}{\dot{\mathbf{u}}}_n)$$
$${\dot{\mathbf{u}}}_{n+1/2}={\dot{\mathbf{u}}}_{n-1/2}+\Delta t\cdot {\ddot{\mathbf{u}}}_n$$
$${\mathbf{u}}_{n+1}={\mathbf{u}}_n+\Delta t\cdot {\dot{\mathbf{u}}}_{n+1/2}$$

优点：

• 每步计算量小，无需迭代
• 适合处理高度非线性问题（如冲击、爆炸）
• 容易实现并行计算

缺点：

• 条件稳定，时间步长受CFL条件限制
• 对于低频响应问题效率较低

隐式积分法：

Newmark-β法是最常用的隐式算法：

$$\mathbf{M}{\ddot{\mathbf{u}}}_{n+1}+\mathbf{C}{\dot{\mathbf{u}}}_{n+1}+{\mathbf{F}}_{n+1}^{int}={\mathbf{F}}_{n+1}^{ext}$$

其中：
$${\dot{\mathbf{u}}}_{n+1}={\dot{\mathbf{u}}}_n+(1-\gamma )\Delta t{\ddot{\mathbf{u}}}_n+\gamma \Delta t{\ddot{\mathbf{u}}}_{n+1}$$
KaTeX parse error: Expected 'EOF', got '}' at position 30: …1} = \mathbf{u}}̲_n + \Delta t \…

求解流程：

1. 预测步：基于当前状态预测下一时刻的位移和速度
2. 迭代求解：使用Newton-Raphson法求解非线性方程
3. 校正步：更新加速度、速度、位移

优点：

• 无条件稳定（对于线性问题）
• 可采用较大的时间步长
• 适合低频主导的响应

缺点：

• 每步需要迭代求解，计算量大
• 对于强非线性问题可能不收敛

3.3 Newton-Raphson迭代

隐式积分中，Newton-Raphson法是求解非线性方程的核心算法。

基本迭代格式：

$${\mathbf{K}}_T^{(i)}\Delta {\mathbf{u}}^{(i)}={\mathbf{R}}^{(i)}$$
$${\mathbf{u}}^{(i+1)}={\mathbf{u}}^{(i)}+\Delta {\mathbf{u}}^{(i)}$$

其中，${\mathbf{K}}_T$ 是切线刚度矩阵，$\mathbf{R}$ 是残差向量。

切线刚度矩阵的构成：

$${\mathbf{K}}_T=\frac{\partial {\mathbf{F}}^{int}}{\partial \mathbf{u}}={\mathbf{K}}_M+{\mathbf{K}}_G+{\mathbf{K}}_L$$

• 材料刚度矩阵 ${\mathbf{K}}_M$：反映材料本构关系的非线性
• 几何刚度矩阵 ${\mathbf{K}}_G$：反映应力对刚度的影响（应力刚化/软化）
• 载荷刚度矩阵 ${\mathbf{K}}_L$：反映跟随力效应

收敛准则：

通常采用以下收敛判据：

1. 位移收敛：$\Vert \Delta \mathbf{u}\Vert /\Vert \mathbf{u}\Vert <{ϵ}_u$
2. 力收敛：$\Vert \mathbf{R}\Vert /\Vert {\mathbf{F}}^{ext}\Vert <{ϵ}_f$
3. 能量收敛：$\Delta {\mathbf{u}}^T\mathbf{R}/({\mathbf{u}}^T{\mathbf{F}}^{ext})<{ϵ}_E$

典型取值：${ϵ}_u=10^{-3}$，${ϵ}_f=10^{-4}$，${ϵ}_E=10^{-6}$。

收敛改进技术：

1. 线搜索：在Newton方向上进行一维搜索，确定最优步长
2. 弧长法：通过控制荷载-位移路径的弧长，追踪后屈曲行为
3. BFGS更新：近似更新切线刚度矩阵，减少重新组装的计算量

3.4 接触-碰撞问题处理

接触问题是典型的边界非线性问题，其难点在于接触区域事先未知，且接触状态随时间变化。

接触条件：

对于两个物体 ${\Omega }^A$ 和 ${\Omega }^B$，接触条件包括：

1. 非穿透条件（几何条件）：
   $$g_N=({\mathbf{x}}^B-{\mathbf{x}}^A)\cdot {\mathbf{n}}^A\ge 0$$
   其中，$g_N$ 是法向间隙，${\mathbf{n}}^A$ 是接触面法向。

2. 压应力条件（力学条件）：
   $$t_N\le 0,g_N\cdot t_N=0$$
   接触面上只能存在压应力，拉应力区域自动分离。

3. 摩擦条件（Coulomb定律）：
   $$\Vert t_T\Vert \le \mu |t_N|$$
   其中，$t_T$ 是切向牵引力，$\mu$ 是摩擦系数。

数值处理方法：

1. Lagrange乘子法：
   将接触约束通过Lagrange乘子引入泛函：
   $$\Pi ={\Pi }_{elastic}+{\int }_{{\Gamma }_c}{\lambda }_N{g}_{N}d\Gamma$$
   优点：精确满足约束；缺点：增加未知量，需要特殊求解器。

2. 罚函数法：
   通过罚参数将约束转化为能量项：
   $$\Pi ={\Pi }_{elastic}+\frac{{ϵ}_N}{2}{\int }_{{\Gamma }_c}⟨g_N{⟩}^{2}d\Gamma$$
   其中，$⟨\cdot ⟩$ 是Macauley括号。优点：实现简单；缺点：约束近似满足，罚参数选择影响精度和收敛。

3. 增广Lagrange法：
   结合前两者的优点，通过迭代更新Lagrange乘子，同时保持较小的罚参数。

---

4. 多物理场耦合仿真

4.1 耦合问题的分类

多物理场耦合根据物理场之间的相互作用强度，可分为：

弱耦合（单向耦合）：

• 一个物理场的输出作为另一个物理场的输入
• 两个场之间没有反馈
• 示例：稳态温度场引起的热应力分析

强耦合（双向耦合）：

• 物理场之间存在相互依赖关系
• 需要同时求解或迭代求解
• 示例：流-固耦合问题中，流体压力改变结构变形，结构变形又影响流场

耦合的数学描述：

对于两个物理场 A 和 B，耦合问题可表示为：

$${\mathcal{L}}_A({\mathbf{u}}_A,{\mathbf{u}}_B)=0$$
$${\mathcal{L}}_B({\mathbf{u}}_B,{\mathbf{u}}_A)=0$$

其中，$\mathcal{L}$ 是微分算子，$\mathbf{u}$ 是场变量。

4.2 耦合求解策略

分区求解法（Partitioned Approach）：

每个物理场使用独立的求解器，通过界面传递数据。

流程：

1. 在时间步 $t_n$，已知场 A 和场 B 的状态
2. 求解场 A，将场 B 的状态作为边界条件
3. 将场 A 的结果传递给场 B
4. 求解场 B
5. 检查收敛性，若不收敛则返回步骤2

优点：

• 可复用现有的单物理场求解器
• 各场可采用最适合的网格和时间步长
• 模块化程度高，易于维护

缺点：

• 强耦合问题收敛慢或发散
• 数据传递可能引入误差
• 稳定性依赖于耦合迭代算法

耦合迭代算法：

1. Gauss-Seidel迭代：
   依次更新各场，使用最新可用的场信息。

2. Jacobi迭代：
   各场同时更新，使用上一迭代的场信息。

3. 松弛迭代：
   引入松弛因子防止振荡：
   $${\mathbf{u}}^{new}=\alpha {\mathbf{u}}^{calculated}+(1-\alpha ){\mathbf{u}}^{old}$$

4. 准Newton法：
   利用历史迭代信息加速收敛，如IQN-ILS算法。

整体求解法（Monolithic Approach）：

将所有物理场的方程组装成一个大型方程组，同时求解。

形式：
$$\left[\begin{array}{cc}{\mathbf{K}}_{AA}& {\mathbf{K}}_{AB}\\ {\mathbf{K}}_{BA}& {\mathbf{K}}_{BB}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_A\\ {\mathbf{u}}_B\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_A\\ {\mathbf{F}}_B\end{array}\right]$$

优点：

• 强耦合问题收敛性好
• 无条件稳定（对于线性问题）
• 时间精度高

缺点：

• 需要开发专用的多物理场求解器
• 计算资源需求大
• 各场必须使用相同的网格和时间步长

4.3 热-结构耦合分析

热-结构耦合是最常见的多物理场问题之一，在航天器热防护、发动机热应力、火灾安全评估等领域有广泛应用。

控制方程：

热传导方程：
$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+Q$$

结构运动方程：
$$\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}=\nabla \cdot \mathbf{\sigma }+\mathbf{f}$$

热-力耦合本构关系：
$$\mathbf{\sigma }=\mathbf{D}:(\mathbf{\epsilon }-{\mathbf{\epsilon }}^{th})$$
$${\mathbf{\epsilon }}^{th}=\alpha (T-T_0)\mathbf{I}$$

其中，$\alpha$ 是热膨胀系数，$T_0$ 是参考温度。

耦合机制：

1. 热→结构（单向）：温度变化引起热应变，产生热应力
2. 结构→热（双向）：
   • 变形改变热边界条件（如接触热阻）
   • 塑性功转化为热量
   • 辐射换热与表面几何相关

时间尺度问题：

热传导和结构振动通常具有不同的时间尺度：

• 热扩散时间尺度：${\tau }_{th}\sim L^2/{\alpha }_{th}$
• 结构振动周期：$T_{str}\sim 1/f$

对于金属结构，通常 $T_{str}\ll {\tau }_{th}$，可采用准静态假设，在每个热时间步求解静态结构问题。

4.4 流-固耦合分析

流-固耦合（FSI）在桥梁风振、海洋平台、心血管流动等领域至关重要。

控制方程：

流体（Navier-Stokes方程）：
$${\rho }_f\left(\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla p+\mu {\nabla }^2\mathbf{v}+\mathbf{f}$$
$$\nabla \cdot \mathbf{v}=0$$

固体（结构动力学方程）：
$${\rho }_s\frac{{\partial }^2\mathbf{u}}{\partial t^2}=\nabla \cdot \mathbf{\sigma }+\mathbf{f}$$

界面条件：

1. 运动学条件（位移/速度连续）：
   $${\mathbf{v}}_f={\dot{\mathbf{u}}}_s\text{on}{\Gamma }_{fsi}$$

2. 动力学条件（力平衡）：
   $${\mathbf{\sigma }}_f\cdot \mathbf{n}={\mathbf{\sigma }}_s\cdot \mathbf{n}\text{on}{\Gamma }_{fsi}$$

数值挑战：

1. 附加质量效应：
   对于不可压缩流体，界面处的压力波动会产生"附加质量"效应，降低系统的有效频率：
   $${\omega }_{eff}^2=\frac{k}{m+m_{added}}$$
   这可能导致隐式耦合迭代的收敛困难。

2. 网格变形：
   结构变形会改变流体域的几何，需要：

   • 动网格技术（ALE描述）
   • 网格重生成
   • 浸没边界法（IBM）

3. 时间步长限制：
   需要同时满足流体CFL条件和结构稳定性条件。

---

5. 案例1：大跨度空间结构模态分析

5.1 工程背景

大跨度空间结构（如体育馆、机场航站楼、展览馆）通常采用网架、网壳、弦支穹顶等形式。这类结构具有跨度大、自重轻、刚度相对较小的特点，对风荷载和地震作用敏感。

本案例以某体育馆弦支穹顶结构为对象，展示复杂空间结构的模态分析流程。

5.2 结构描述

几何参数：

• 平面尺寸：120m × 100m（椭圆形平面）
• 矢高：15m
• 网格尺寸：约5m × 5m
• 总节点数：约2500个
• 总杆件数：约8000根

构件规格：

• 上弦杆：圆钢管，直径300mm，壁厚12mm
• 下弦杆：圆钢管，直径250mm，壁厚10mm
• 腹杆：圆钢管，直径150mm，壁厚8mm
• 拉索：钢绞线，直径50mm

材料参数：

• 钢材：Q345，弹性模量 $E=206$ GPa，密度 $\rho =7850$ kg/m³
• 索材：高强度钢丝，弹性模量 $E=195$ GPa，密度 $\rho =7850$ kg/m³

边界条件：

• 周边支座：固定铰支座（约束竖向和水平位移，允许转动）
• 支座数量：沿周边均匀布置，共48个

5.3 有限元建模

建模策略：

1. 杆件建模：
   采用空间杆单元（ truss element），每个节点3个平动自由度。

   单元刚度矩阵：
   $${\mathbf{k}}^e=\frac{EA}{L}\left[\begin{array}{c}\mathbf{n}\\ -\mathbf{n}\end{array}\right]\left[\begin{array}{cc}{\mathbf{n}}^T& -{\mathbf{n}}^T\end{array}\right]$$
   其中，$\mathbf{n}$ 是单元方向余弦向量。

2. 拉索建模：
   拉索需要考虑几何非线性（悬链线效应），采用两节点索单元或等效杆单元。

   索的等效弹性模量（Ernst公式）：
   $$E_{eq}=\frac{E}{1+\frac{(\gamma L{)}^{2}EA}{12T^3}}$$
   其中，$\gamma$ 是索的容重，$T$ 是索张力。

3. 预应力处理：
   弦支穹顶的拉索需要施加预张力，通过几何刚度矩阵考虑预应力的刚化效应：
   $${\mathbf{K}}_G=\sum _e\frac{T^e}{L^e}{\mathbf{G}}^{eT}{\mathbf{G}}^e$$

4. 质量矩阵：
   采用一致质量矩阵或集中质量矩阵。对于动力分析，集中质量矩阵计算效率更高：
   $$\mathbf{M}=\text{diag}(m_1,m_1,m_1,m_2,m_2,m_2,...)$$

Python实现：

"""
主题093 - 案例1: 大跨度空间结构模态分析
工程背景: 弦支穹顶体育馆结构
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation, PillowWriter
from scipy import linalg
import matplotlib
matplotlib.use('Agg')

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

print("=" * 70)
print("案例1: 大跨度空间结构模态分析")
print("工程背景: 弦支穹顶体育馆结构")
print("=" * 70)

# =============================================================================
# 1. 结构几何参数
# =============================================================================
print("\n【1. 结构几何参数】")

# 穹顶几何
L_x = 120.0      # 长轴 (m)
L_y = 100.0      # 短轴 (m)
f_rise = 15.0    # 矢高 (m)
n_x = 25         # 长轴方向网格数
n_y = 21         # 短轴方向网格数

# 生成节点坐标
def generate_dome_nodes():
    """生成弦支穹顶节点坐标"""
    nodes = []
    node_id = 0
    
    # 上弦节点 (椭圆抛物面)
    for i in range(n_x):
        for j in range(n_y):
            x = (i / (n_x - 1) - 0.5) * L_x
            y = (j / (n_y - 1) - 0.5) * L_y
            # 椭圆抛物面方程: z = f * (1 - (x/a)^2 - (y/b)^2)
            z = f_rise * (1 - (2*x/L_x)**2 - (2*y/L_y)**2)
            if z >= 0:  # 只保留穹顶部分
                nodes.append([node_id, x, y, z, 'top'])
                node_id += 1
    
    # 下弦节点 (位于上弦下方一定距离)
    bottom_offset = 3.0  # m
    for i in range(n_x):
        for j in range(n_y):
            x = (i / (n_x - 1) - 0.5) * L_x
            y = (j / (n_y - 1) - 0.5) * L_y
            z = f_rise * (1 - (2*x/L_x)**2 - (2*y/L_y)**2) - bottom_offset
            if z >= -bottom_offset + 1:  # 只保留有效区域
                nodes.append([node_id, x, y, z, 'bottom'])
                node_id += 1
    
    return np.array(nodes, dtype=object)

nodes = generate_dome_nodes()
n_nodes = len(nodes)
print(f"总节点数: {n_nodes}")

# 提取坐标
node_coords = np.array([[n[1], n[2], n[3]] for n in nodes])
node_types = [n[4] for n in nodes]

# =============================================================================
# 2. 单元连接
# =============================================================================
print("\n【2. 单元连接】")

def generate_elements(nodes):
    """生成杆件和索单元"""
    elements = []
    elem_id = 0
    
    # 创建节点位置字典用于快速查找
    node_dict = {}
    for i, n in enumerate(nodes):
        key = (round(n[1], 3), round(n[2], 3))
        node_dict[key] = i
    
    # 上弦杆 (连接相邻上弦节点)
    top_nodes = [i for i, n in enumerate(nodes) if n[4] == 'top']
    for i in top_nodes:
        xi, yi, zi = nodes[i][1], nodes[i][2], nodes[i][3]
        for j in top_nodes:
            if i >= j:
                continue
            xj, yj, zj = nodes[j][1], nodes[j][2], nodes[j][3]
            dist = np.sqrt((xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2)
            if 4.5 < dist < 6.0:  # 网格尺寸约5m
                elements.append([elem_id, i, j, 'top_chord', 0.3, 0.012])
                elem_id += 1
    
    # 下弦杆
    bottom_nodes = [i for i, n in enumerate(nodes) if n[4] == 'bottom']
    for i in bottom_nodes:
        xi, yi, zi = nodes[i][1], nodes[i][2], nodes[i][3]
        for j in bottom_nodes:
            if i >= j:
                continue
            xj, yj, zj = nodes[j][1], nodes[j][2], nodes[j][3]
            dist = np.sqrt((xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2)
            if 4.5 < dist < 6.0:
                elements.append([elem_id, i, j, 'bottom_chord', 0.25, 0.010])
                elem_id += 1
    
    # 腹杆 (连接上下弦)
    for i_top in top_nodes:
        x_t, y_t = nodes[i_top][1], nodes[i_top][2]
        # 找最近的下弦节点
        min_dist = float('inf')
        nearest_bottom = None
        for i_bot in bottom_nodes:
            x_b, y_b = nodes[i_bot][1], nodes[i_bot][2]
            dist = np.sqrt((x_t-x_b)**2 + (y_t-y_b)**2)
            if dist < min_dist:
                min_dist = dist
                nearest_bottom = i_bot
        if nearest_bottom is not None and min_dist < 5.0:
            elements.append([elem_id, i_top, nearest_bottom, 'web', 0.15, 0.008])
            elem_id += 1
    
    # 拉索 (连接下弦节点到支座)
    # 简化: 下弦周边节点连接到地面环梁
    support_z = -5.0
    for i_bot in bottom_nodes:
        x, y, z = nodes[i_bot][1], nodes[i_bot][2], nodes[i_bot][3]
        r = np.sqrt(x**2 + y**2)
        if r > 45:  # 周边节点
            elements.append([elem_id, i_bot, None, 'cable', 0.05, 0.0, support_z])
            elem_id += 1
    
    return elements

elements = generate_elements(nodes)
n_elements = len(elements)
print(f"总单元数: {n_elements}")

# 统计各类单元
n_top = sum(1 for e in elements if e[3] == 'top_chord')
n_bottom = sum(1 for e in elements if e[3] == 'bottom_chord')
n_web = sum(1 for e in elements if e[3] == 'web')
n_cable = sum(1 for e in elements if e[3] == 'cable')
print(f"  上弦杆: {n_top}")
print(f"  下弦杆: {n_bottom}")
print(f"  腹杆: {n_web}")
print(f"  拉索: {n_cable}")

# =============================================================================
# 3. 材料参数与截面特性
# =============================================================================
print("\n【3. 材料参数】")

E_steel = 206e9      # 钢材弹性模量 (Pa)
rho_steel = 7850     # 钢材密度 (kg/m³)
E_cable = 195e9      # 索材弹性模量 (Pa)
cable_tension = 500e3  # 索预张力 (N)

print(f"钢材弹性模量: {E_steel/1e9:.0f} GPa")
print(f"索材弹性模量: {E_cable/1e9:.0f} GPa")
print(f"索预张力: {cable_tension/1e3:.0f} kN")

# =============================================================================
# 4. 组装刚度矩阵和质量矩阵
# =============================================================================
print("\n【4. 组装系统矩阵】")

n_dof = n_nodes * 3  # 每个节点3个平动自由度

# 初始化矩阵
K = np.zeros((n_dof, n_dof))
M = np.zeros((n_dof, n_dof))

# 组装单元矩阵
for elem in elements:
    elem_id, node_i, node_j, elem_type, d, t = elem[0], elem[1], elem[2], elem[3], elem[4], elem[5]
    
    if elem_type == 'cable':
        # 拉索单元简化处理
        support_z = elem[6]
        xi, yi, zi = nodes[node_i][1], nodes[node_i][2], nodes[node_i][3]
        xj, yj, zj = xi, yi, support_z  # 支座点
        A = np.pi * (d/2)**2
        E = E_cable
        # 几何刚度贡献
        L = np.sqrt((xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2)
        if L > 0:
            k_geom = cable_tension / L
            # 简化为对节点的刚度贡献
            dof_i = [node_i*3, node_i*3+1, node_i*3+2]
            for ii in range(3):
                M[dof_i[ii], dof_i[ii]] += rho_steel * A * L / 3
    else:
        # 普通杆件
        xi, yi, zi = nodes[node_i][1], nodes[node_i][2], nodes[node_i][3]
        xj, yj, zj = nodes[node_j][1], nodes[node_j][2], nodes[node_j][3]
        
        L = np.sqrt((xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2)
        A = np.pi * ((d/2)**2 - (d/2-t)**2)  # 环形截面
        
        # 方向余弦
        cx, cy, cz = (xj-xi)/L, (yj-yi)/L, (zj-zi)/L
        
        # 单元刚度矩阵 (6x6)
        k_local = (E_steel * A / L) * np.array([
            [ cx*cx,  cx*cy,  cx*cz, -cx*cx, -cx*cy, -cx*cz],
            [ cy*cx,  cy*cy,  cy*cz, -cy*cx, -cy*cy, -cy*cz],
            [ cz*cx,  cz*cy,  cz*cz, -cz*cx, -cz*cy, -cz*cz],
            [-cx*cx, -cx*cy, -cx*cz,  cx*cx,  cx*cy,  cx*cz],
            [-cy*cx, -cy*cy, -cy*cz,  cy*cx,  cy*cy,  cy*cz],
            [-cz*cx, -cz*cy, -cz*cz,  cz*cx,  cz*cy,  cz*cz]
        ])
        
        # 组装到全局矩阵
        dof_map = [node_i*3, node_i*3+1, node_i*3+2, node_j*3, node_j*3+1, node_j*3+2]
        for ii in range(6):
            for jj in range(6):
                K[dof_map[ii], dof_map[jj]] += k_local[ii, jj]
        
        # 一致质量矩阵 (简化为集中质量)
        m_elem = rho_steel * A * L
        for ii in range(6):
            M[dof_map[ii], dof_map[ii]] += m_elem / 2

print(f"系统自由度: {n_dof}")

# =============================================================================
# 5. 边界条件处理
# =============================================================================
print("\n【5. 边界条件处理】")

# 识别边界节点 (周边节点)
boundary_nodes = []
for i, n in enumerate(nodes):
    x, y = n[1], n[2]
    r = np.sqrt(x**2 + y**2)
    if r > 50:  # 周边区域
        boundary_nodes.append(i)

print(f"边界节点数: {len(boundary_nodes)}")

# 约束边界节点的所有平动自由度
constrained_dofs = []
for node_id in boundary_nodes:
    constrained_dofs.extend([node_id*3, node_id*3+1, node_id*3+2])

# 创建降阶矩阵
free_dofs = [i for i in range(n_dof) if i not in constrained_dofs]
n_free = len(free_dofs)

K_reduced = K[np.ix_(free_dofs, free_dofs)]
M_reduced = M[np.ix_(free_dofs, free_dofs)]

print(f"约束后自由度: {n_free}")

# =============================================================================
# 6. 求解特征值问题
# =============================================================================
print("\n【6. 模态分析】")

# 求解广义特征值问题
n_modes = 20
eigvals, eigvecs = linalg.eigh(K_reduced, M_reduced, eigvals=(0, n_modes-1))

# 计算频率和周期
omega = np.sqrt(np.abs(eigvals))
freqs = omega / (2 * np.pi)
periods = 1 / freqs

print(f"前{n_modes}阶固有频率:")
for i in range(min(10, n_modes)):
    print(f"  模态{i+1}: {freqs[i]:.3f} Hz, 周期 {periods[i]:.3f} s")

# 将振型扩展回完整自由度
phi_full = np.zeros((n_dof, n_modes))
for i in range(n_modes):
    phi_full[free_dofs, i] = eigvecs[:, i]

# =============================================================================
# 7. 模态验证
# =============================================================================
print("\n【7. 模态验证】")

# 正交性检验
orth_error = np.zeros((n_modes, n_modes))
for i in range(n_modes):
    for j in range(n_modes):
        orth_error[i, j] = np.abs(phi_full[:, i].T @ M @ phi_full[:, j] - (1 if i==j else 0))

max_orth_error = np.max(orth_error)
print(f"正交性最大误差: {max_orth_error:.2e}")

# 参与质量计算
participation_mass = np.zeros(n_modes)
total_mass = np.sum(M.diagonal())
for i in range(n_modes):
    # Z方向参与系数
    gamma_z = phi_full[2::3, i].T @ M[2::3, 2::3].diagonal()
    participation_mass[i] = gamma_z**2 / (phi_full[:, i].T @ M @ phi_full[:, i])

cumulative_mass = np.cumsum(participation_mass) / total_mass * 100
print(f"\n前5阶模态Z向参与质量比:")
for i in range(5):
    print(f"  模态{i+1}: {participation_mass[i]/total_mass*100:.2f}%, 累计: {cumulative_mass[i]:.2f}%")

# =============================================================================
# 8. 生成可视化
# =============================================================================
print("\n【8. 生成可视化】")

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

# 1. 结构三维模型
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
# 绘制杆件
for elem in elements:
    if elem[3] != 'cable':
        node_i, node_j = elem[1], elem[2]
        xi, yi, zi = nodes[node_i][1], nodes[node_i][2], nodes[node_i][3]
        xj, yj, zj = nodes[node_j][1], nodes[node_j][2], nodes[node_j][3]
        color = 'steelblue' if elem[3] == 'top_chord' else 'coral' if elem[3] == 'bottom_chord' else 'gray'
        ax1.plot([xi, xj], [yi, yj], [zi, zj], color=color, linewidth=0.8, alpha=0.7)

ax1.scatter(node_coords[:, 0], node_coords[:, 1], node_coords[:, 2], 
           c='red', s=10, alpha=0.5)
ax1.set_xlabel('X (m)')
ax1.set_ylabel('Y (m)')
ax1.set_zlabel('Z (m)')
ax1.set_title('弦支穹顶结构模型')

# 2. 频率分布
ax2 = fig.add_subplot(2, 3, 2)
mode_nums = np.arange(1, n_modes + 1)
ax2.bar(mode_nums, freqs, color='steelblue', alpha=0.7)
ax2.set_xlabel('模态阶数')
ax2.set_ylabel('频率 (Hz)')
ax2.set_title('固有频率分布')
ax2.grid(True, alpha=0.3)

# 3. 振型展示 (前4阶)
for i in range(4):
    ax = fig.add_subplot(2, 3, 3+i, projection='3d')
    
    # 变形后的坐标
    scale = 20  # 放大系数
    disp = phi_full[:, i].reshape(-1, 3)
    coords_deformed = node_coords + disp * scale
    
    # 绘制变形结构
    for elem in elements[:min(200, len(elements))]:  # 限制绘制数量
        if elem[3] != 'cable':
            node_i, node_j = elem[1], elem[2]
            xi, yi, zi = coords_deformed[node_i]
            xj, yj, zj = coords_deformed[node_j]
            ax.plot([xi, xj], [yi, yj], [zi, zj], 'b-', linewidth=0.5, alpha=0.6)
    
    ax.scatter(coords_deformed[:, 0], coords_deformed[:, 1], coords_deformed[:, 2],
              c='red', s=5, alpha=0.5)
    ax.set_title(f'模态{i+1}: {freqs[i]:.3f} Hz')
    ax.set_xlabel('X (m)')
    ax.set_ylabel('Y (m)')
    ax.set_zlabel('Z (m)')

plt.tight_layout()
plt.savefig('案例1_大跨度空间结构模态分析.png', dpi=150, bbox_inches='tight')
print("静态图已保存: 案例1_大跨度空间结构模态分析.png")

# =============================================================================
# 9. 生成动画 - 振型动画
# =============================================================================
print("\n【9. 生成动画】")

fig_anim, axes = plt.subplots(2, 2, figsize=(12, 10), subplot_kw={'projection': '3d'})
axes = axes.flatten()

# 准备动画数据
n_frames = 50
scale_anim = 15

def update(frame):
    t = frame / n_frames * 2 * np.pi
    
    for i, ax in enumerate(axes):
        ax.clear()
        
        # 变形坐标 (随时间变化)
        disp = phi_full[:, i].reshape(-1, 3)
        coords_anim = node_coords + disp * scale_anim * np.sin(t)
        
        # 绘制结构
        for elem in elements[:min(150, len(elements))]:
            if elem[3] != 'cable':
                node_i, node_j = elem[1], elem[2]
                xi, yi, zi = coords_anim[node_i]
                xj, yj, zj = coords_anim[node_j]
                ax.plot([xi, xj], [yi, yj], [zi, zj], 'b-', linewidth=0.5, alpha=0.6)
        
        ax.scatter(coords_anim[:, 0], coords_anim[:, 1], coords_anim[:, 2],
                  c='red', s=5, alpha=0.5)
        ax.set_title(f'模态{i+1}: {freqs[i]:.3f} Hz')
        ax.set_xlabel('X (m)')
        ax.set_ylabel('Y (m)')
        ax.set_zlabel('Z (m)')
        
        # 固定视角范围
        ax.set_xlim([-70, 70])
        ax.set_ylim([-60, 60])
        ax.set_zlim([-10, 20])
    
    return axes

anim = FuncAnimation(fig_anim, update, frames=n_frames, interval=100, blit=False)
writer = PillowWriter(fps=10)
anim.save('案例1_穹顶振型动画.gif', writer=writer)
print("动画已保存: 案例1_穹顶振型动画.gif")

plt.close('all')

# =============================================================================
# 10. 结果汇总
# =============================================================================
print("\n" + "=" * 70)
print("【评估结果汇总】")
print("=" * 70)

print(f"\n1. 结构参数:")
print(f"   - 平面尺寸: {L_x}m × {L_y}m")
print(f"   - 矢高: {f_rise}m")
print(f"   - 节点数: {n_nodes}")
print(f"   - 单元数: {n_elements}")

print(f"\n2. 模态特性:")
print(f"   - 基频: {freqs[0]:.3f} Hz (周期 {periods[0]:.3f}s)")
print(f"   - 前5阶频率范围: {freqs[0]:.3f} - {freqs[4]:.3f} Hz")
print(f"   - 正交性误差: {max_orth_error:.2e}")

print(f"\n3. 质量参与:")
print(f"   - 前5阶累计Z向参与质量: {cumulative_mass[4]:.1f}%")

print(f"\n4. 工程评估:")
if freqs[0] > 1.0:
    print(f"   - 基频较高，结构刚度良好")
else:
    print(f"   - 基频较低，对风荷载敏感")

print("\n" + "=" * 70)
print("案例1分析完成!")
print("=" * 70)