非线性振动问题周期解的空间-时间有限元路径跟踪方法

《Computer Methods in Applied Mechanics and Engineering》：Space-time finite element path-following for periodic solutions in nonlinear vibration problems

　　

提出了一种高效路径跟踪方法，用于计算具有大量自由度的非线性动力系统的周期解。该方法在结构动力学问题（如桁架和壳体）上进行了验证。该方法采用基于位移的伽辽金有限元离散，在空间和时间上均连续。通过在装配中识别第一个和最后一个时间节点来强制位移周期性，从而在弱形式下自然确保速度周期性。空间-时间平衡方程与弧长约束迭代求解，能够相对于延续参数（如激励频率）逐步追踪周期解路径，包括多解和不稳定分支。弱形式在无量纲时间域中表述，因此时间离散不依赖于加载频率。精确雅可比矩阵（包括对延续参数的导数）由有限元贡献高效组装。利用这一特性，结合空间-时间刚度矩阵的稀疏性和修正牛顿求解器，该方法对高维模型仍保持计算效率。最后，实现了缩减时间求积作为粗时间离散的稳定化技术，使得使用有限数量的二次或三次时间单元即可稳健且准确地重建解曲线。

### 论文解读：基于空间-时间有限元路径跟踪的非线性振动周期解计算方法

#### 一、研究背景与问题

在工程实践中，承受周期性载荷或边界条件的阻尼机械系统通常经过瞬态阶段后演化为周期性稳态响应。然而，控制方程中的非线性可能导致该渐近状态非唯一，产生多个周期轨道或更复杂的吸引子，且对系统参数的小扰动高度敏感。对于由空间有限元离散产生的大量自由度（DOFs）动力系统，计算此类周期解及其稳定性和敏感性，仍是计算力学中的重大挑战。传统的时间有限差分积分方法需要模拟多个周期直至响应稳定，计算成本高且难以评估解的独特性和稳定性。因此，研究者致力于开发直接求周期解的方法，并集成到路径跟踪分析工具中，以系统研究周期解随激励频率或设计参数的变化，包括追踪多解和不稳定分支。

#### 二、研究内容与核心技术方法

研究人员提出了一种基于完全连续空间-时间有限元法（FEM）的路径跟踪策略，用于计算高维非线性结构动力学中的周期解。该方法在空间和时间上均采用连续伽辽金离散，通过装配中识别首尾时间节点强制位移周期性，弱形式下自然满足速度周期性。技术方法主要包括：（1）无量纲时间变量变换，使时间离散不依赖加载频率；（2）空间-时间平衡方程与弧长约束的迭代求解，实现周期解路径的逐步追踪；（3）精确雅可比矩阵（含对延续参数的导数）由单元贡献高效组装，结合稀疏空间-时间刚度矩阵和修正牛顿求解器；（4）采用缩减时间求积（每单元p个高斯点）作为稳定化技术，抑制粗时间离散中的虚假高频解。数值测试涵盖Duffing振子、桁架结构和浅拱壳结构，样本包括不同规模桁架（4至128跨）和壳模型（1728自由度），并与改进的射击法（shooting method）结果对比。

#### 三、主要研究结果

**1. Duffing振子稳定性与收敛性验证**
通过线性（a=1, b=0）和非线性（a=1, b=0.002）Duffing方程，研究了时间离散的离散inf-sup稳定性。采用离散inf-sup参数β（式14）评估，结果表明：全积分（每单元p+1个高斯点）在粗网格下β在特定频率趋近于零，而缩减积分（每单元p个高斯点）使β在关注频率范围内均匀有界，有效稳定离散模型。收敛性测试显示，p=1,2,3阶时间单元的峰值位移相对误差均达到p+1阶渐近收敛率，全积分精度通常更高，但缩减积分在稳定性问题出现时为首选。

**2. 桁架结构收敛性与计算效率**
以32跨平面桁架为例（几何尺寸L=50m，H=1m，材料为钢），施加两个垂直谐波力，采用线性、二次和三次时间单元进行收敛性研究。线性单元需约32个单元达到收敛；二次和三次单元在全积分下出现虚假解，采用缩减积分后仅需8个单元即可获得准确结果。与射击法（使用梯形规则，32步）对比，频率响应曲线完全一致。计算时间对比显示：对于124自由度模型，射击法耗时418.4秒，而三次缩减积分单元仅需1.60秒，效率提升超过两个数量级。

**3. 壳结构性能验证**
浅拱壳结构（L=2m，H=1m，t=0.01m，w=0.1m，256个实体壳单元，1728自由度）承受集中谐波载荷，其频率响应曲线呈现典型的软化行为。线性时间单元需64个单元收敛，二次和三次单元分别需8和4个单元。射击法（64步）耗时5小时31分43秒，而三次单元时间FEM仅需4分18秒，效率提升约77倍。对于具有1:2内共振的浅拱壳（w=0.041m），线性单元需128个单元（约3小时），三次缩减积分单元仅需8个单元（约20分钟），而射击法（128步）耗时22小时40分52秒。采用修正牛顿法（mNR）进一步将三次缩减积分单元的计算时间从20分48秒降至9分28秒。

**4. 计算规模扩展性**
通过逐步增加桁架跨数（4至128跨）和壳模型网格细化，系统评估了平均每步计算成本随空间自由度数（DOFs）的变化。射击法计算成本随DOFs增长显著，主要源于雅可比矩阵的稠密矩阵-矩阵乘积；而所提空间-时间FEM方法由于雅可比矩阵的显式局部组装和稀疏矩阵求解，总成本始终远低于射击法，在高维问题中总计算时间降低超过一个数量级。

#### 四、讨论与结论

讨论部分指出，所提方法与射击法相比，虽然系统规模更大（mN+1 vs 2N+1，m为时间节点数），但雅可比矩阵结构稀疏且易于显式组装，避免了射击法中的递归全局计算。与谐波平衡法（HBM）相比，时间近似具有局部支撑特性，便于稀疏求解器和自适应h/p细化，且无混叠问题。缩减求积规则作为稳定化技术，有效过滤了高阶时间离散中的虚假高频模式，保证了粗网格下的鲁棒性。

研究结论：开发了一种基于完全连续空间-时间有限元公式的延续策略，用于高维非线性结构动力学中的周期解计算。该方法通过空间-时间单元装配直接强制时间周期性，形成统一的空间-时间离散框架。弧长路径跟踪程序能够追踪周期解分支随激励频率或其他设计参数的演化。数值实验验证了该框架的准确性、鲁棒性和可扩展性。缩减时间求积增强了数值稳定性和收敛性，仅需每周期少量时间单元即可获得准确频率响应曲线，使该方法与基于谐波平衡的现有方法具有竞争力。未来工作将扩展到含内变量和历史依赖行为的非线性材料模型、流体相互作用，并探索并行代数多重网格（AMG）预条件广义最小残差（GMRES）迭代求解器以处理超大规模问题。该框架也有望应用于涉及电生理和流固耦合的复杂周期生物力学问题。