基于半光滑牛顿法的含分层双层梁受迫振动分析：一项复合材料非线性动力学研究

一、 研究团队与发表信息

本项研究由河海大学力学与工程科学学院、江苏省风电机组结构工程研究中心的刘书恒（硕士研究生）、刘庆辉（副教授，通讯作者）与山东黄河勘测设计研究院有限公司的黄晓阳共同完成。研究成果以题为《基于半光滑牛顿法的含分层双层梁受迫振动分析》的学术论文形式，于2026年5月20日在期刊《复合材料科学与工程》（Composites Science and Engineering）上进行了网络首发。

二、 研究背景与目的

本研究隶属于复合材料力学与结构非线性动力学交叉领域。分层（Delamination）是复合材料层合结构中最常见且危害性极大的缺陷之一，它会显著改变结构的动态特性，如自振频率、阻尼和振动响应，甚至可能诱发共振或疲劳破坏。因此，精确预测含分层结构在动态载荷下的行为，对于结构健康监测、损伤识别和抗疲劳设计至关重要。

传统的含分层结构动力学分析，多在线性自由振动或基于“自由模式”与“约束模式”假设的框架下进行。然而，实际工况中，分层界面在振动过程中会反复发生接触、分离的动态开合行为，这是一种强烈的接触非线性现象。以往研究在处理这种动态接触问题时，常采用拉格朗日乘子法或罚函数法等传统光滑数值方法，但在接触状态频繁切换时易出现收敛困难甚至数值发散的问题。

针对上述挑战，本研究旨在开发一种鲁棒、高效的数值方法，以精确分析含分层复合材料梁在受迫振动下的复杂非线性动力学行为。具体目标包括：1）建立一种能够有效模拟分层界面动态接触的有限元模型；2）引入半光滑牛顿法（Semi-smooth Newton Method） 来处理接触互补条件，提升求解的稳定性和收敛性；3）系统研究激励频率对含分层双层梁动力学响应（如周期运动、倍周期分岔、混沌运动）的影响规律。

三、 研究详细工作流程

本研究的工作流程是一个集成了新型有限元建模、先进非线性求解算法和动力学仿真的系统性数值分析过程。

第一步骤：建立仅具有平动自由度的梁单元模型。
为了高效模拟含分层的多层梁结构，研究摒弃了传统的每个节点具有纵向、横向位移和转角三个自由度的Timoshenko梁单元。转而采用了一种仅具有平动自由度的梁单元。该单元有四个节点，但只有六个平动自由度（每个节点组共享一个横向位移，各有独立的纵向位移）。这种单元通过特殊的形函数构造，能够描述单元内的位移场，并推导出相应的单元刚度矩阵（由弯曲刚度矩阵和剪切刚度矩阵构成）。这种单元的优势在于，当模拟具有多个分层的复杂结构时，无需在未开裂区域引入复杂的界面连续条件，建模更简洁，计算效率更高。该单元与经典梁单元的自由度可通过数学关系等价转换。

第二步骤：引入半光滑牛顿法处理动态接触问题。
这是本研究的核心创新点。研究将分层界面的接触问题抽象为结点-结点接触模型。接触的关键是满足法向的互补条件：接触间隙（g）与接触力（λ）非负，且其乘积为零（即，有接触时无间隙，有间隙时无接触）。传统牛顿法无法直接处理这种非光滑的互补条件。

1. 转化为非光滑方程：研究引入Fischer-Burmeister互补函数，将互补条件转化为一组等价的非光滑方程。对于每个接触对k，构造方程 Φ_k(λ_k, g_k) = 0，其中g_k由上下子梁对应节点的位移通过映射矩阵B计算得到。
2. 构造广义雅可比矩阵：由于互补函数在接触状态切换点（如λ_k=0, g_k=0）不可微，传统导数不存在。研究采用广义导数的概念来构造求解所需的广义雅可比矩阵（J）。该矩阵不仅包含结构的整体刚度矩阵（K），还包含了由接触条件导出的项（如d_g, d_λ），以及为了数值稳定性引入的自适应正则化参数（ε_reg） 和针对活动接触对的虚拟刚度项（K_pen）。自适应正则化策略能根据迭代情况动态调整参数，有效抑制界面穿透并加速收敛。
3. 线性化求解与全局收敛保障：在每一迭代步，求解线性方程组 J * [Δu, Δλ]^T = -r，其中r为包含平衡残差和非光滑互补残差的总体残差向量，从而更新位移u和接触力λ。为确保算法的全局收敛性（即从任意初始点出发都能收敛），研究采用了基于Merit函数的线搜索策略。该策略通过控制迭代步长α，保证每次迭代都朝着使目标函数下降的方向进行，并通过投影函数保证接触力的非负性。

第三步骤：集成Newmark-β时域积分进行动力学分析。
为了求解含接触非线性的动力学方程，研究将上述半光滑牛顿法与Newmark-β时域逐步积分法相结合。在每个时间步（t_i → t_i+1）内：

1. 等效静力化：利用Newmark-β法的参数（取γ=0.65, β=(γ+0.5)^2/4），根据t_i时刻的位移、速度、加速度，预测t_i+1时刻的值，并将动力学方程转化为一个等效的静力平衡方程：K_eff * u_i+1 = f_eff。其中，K_eff和f_eff是包含了质量、阻尼和上一时间步历史贡献的有效刚度矩阵和有效荷载向量。
2. 接触迭代求解：对这个等效静力问题，调用第二步骤中的半光滑牛顿法进行求解。由于接触状态未知，这是一个需要迭代的过程。算法会根据当前计算的间隙g判断接触状态，更新广义雅可比矩阵，直至残差收敛，得到t_i+1时刻满足接触条件的位移解u_i+1。
3. 状态量更新：根据Newmark-β法的公式，用求得的u_i+1回代，更新t_i+1时刻的速度v_i+1和加速度a_i+1。
4. 时间步推进：将u_i+1, v_i+1, a_i+1作为下一时间步的初始值，循环往复，直至模拟完整个动力响应历程。

整个算法的有限元实现流程如图5所示，清晰地展示了在每一个时间步内，静力接触迭代与动力时间积分是如何紧密耦合的。

第四步骤：模型验证与算例分析。
为验证所提出模型和算法的有效性与精度，研究以一个具体算例展开分析。

1. 研究对象与参数：研究对象为两端固支的含分层双层梁，上层为铝合金，下层为CFRP（碳纤维增强复合材料），具体几何尺寸和材料参数见表1。在Ansys商业软件中建立对比模型，采用Plane182单元模拟梁，并用接触单元模拟分层界面。
2. 验证内容一：自由振动分析。在“自由模式”（即假设分层界面上下梁可自由分离）下，计算了该梁的前四阶固有频率。将本文方法（使用较少单元）的结果与Ansys的精细网格模型结果对比（表2）。结果显示，两者吻合良好，最大误差为5.7%（二阶频率），验证了本文单元模型在模态分析上的准确性。
3. 验证内容二：受迫振动响应分析。在梁跨中施加正弦载荷P=F sin(2πft)。选取了两个激励频率：f=5.2 Hz（远低于一阶固有频率）和f=4805 Hz（接近一阶固有频率的0.9倍）。对比了本文方法与Ansys得到的时间-位移曲线（图6）。结果表明，在两种频率下，本文方法都能精确捕捉上下子梁的位移响应，并且能够有效反映分层界面的动态开闭状态（曲线中反映出的非线性特征）。这验证了本文所提出的半光滑牛顿法结合Newmark-β法在求解含接触非线性动力问题上的有效性和精度。

第五步骤：非线性动力学特性系统性分析。
在验证模型有效的基础上，研究系统性地探索了激励频率对系统动力学行为的影响。以一阶固有频率f1为基准，选取了f=0.1f1, 0.3f1, 0.5f1, 0.7f1, f1, 2f1, 3f1等多个频率点，并设置了较小的阻尼比（ξ=0.005）以突出接触非线性的影响。对于每个频率，研究输出了：

• 时间-位移曲线：直观显示位移随时间的变化。
• 相轨迹图：以位移为横坐标，速度为纵坐标，描绘系统在相空间中的运动轨迹。
• 庞加莱截面图：在每个激励周期起始点对相轨迹进行采样，将采样点绘制在相空间中，用于判断运动的周期性或混沌性。

四、 主要研究结果

研究通过上述流程，获得了以下关键结果：

1. 模型验证结果：如表2和图6所示，本文建立的基于半光滑牛顿法和仅具有平动自由度梁单元的模型，在计算效率和精度上取得了良好平衡。即使在单元数量较少的情况下，其计算的固有频率和动力响应时程曲线与Ansys商业软件的解高度一致，证明了该方法的可靠性。

2. 非线性动力学行为图谱：通过对不同激励频率下的响应分析，研究清晰地揭示了系统随激励频率变化而经历的丰富动力学行为（图7至图13）：

   • 低频区（f=0.1f1）：系统响应接近准静态，时间-位移曲线为规则的正弦波，相轨迹为单一闭合曲线，庞加莱截面为一个孤立点。这表明系统处于单周期运动状态。此时激励频率低，接触非线性效应较弱。
   • 中低频区（f=0.3f1）：系统行为发生质变。时间-位移曲线显示每两个激励周期完成一次完整循环，相轨迹出现两条闭合曲线，庞加莱截面出现两个孤立点。这表明系统发生了倍周期分岔，进入2倍周期运动。这是系统从周期运动向更复杂运动演化的临界标志。
   • 中频区（f=0.5f1, 0.7f1）：系统响应变得极其复杂。时间-位移曲线幅值无规律波动，相轨迹杂乱且充满折叠，庞加莱截面上的点呈散乱分布，形成具有某种结构的点集。这些是混沌运动的典型特征。此时，激励频率使得接触状态的切换变得频繁且无固定规律，非线性效应占据主导，系统对初始条件极度敏感。
   • 共振区（f=f1）：当激励频率等于系统一阶固有频率时，线性共振效应起主导作用。时间-位移曲线恢复为大幅值的规则周期振动，相轨迹为单一闭合曲线，庞加莱截面回归为一个点。系统从混沌状态返回单周期运动。这是因为强烈的共振响应抑制了接触状态的频繁切换。
   • 高频区（f=2f1, 3f1）：当激励频率为固有频率的整数倍时，系统再次进入混沌运动状态。此时激发了亚谐波共振（即响应中出现频率为激励频率分数倍的成分），接触非线性效应与亚谐波共振耦合，导致运动再度失稳。

3. 相轨迹的非光滑特征：在所有频率的相轨迹图中，均观察到一个共同现象：轨迹由光滑曲线段和非光滑的折线段交替构成。光滑段对应分层界面分离（无接触）的振动阶段，此时上下梁独立运动；非光滑的折线段则精确对应了接触发生的瞬间，此时速度发生突变。这从几何上直观证实了接触非线性的本质，即系统动力学在接触点处是非光滑的。

五、 研究结论与价值

本研究成功发展了一套基于半光滑牛顿法的鲁棒数值框架，用于分析含分层双层梁的受迫振动接触非线性问题。主要结论如下：

1. 方法有效性：所提出的方法能够稳定、精确地求解分层界面动态开闭引起的强非线性动力问题，克服了传统方法在接触状态切换处易发散的缺点。
2. 动力学行为规律：即使不考虑几何非线性，仅由分层界面接触条件改变引起的非线性，就足以导致含分层双层梁产生极其丰富的动力学行为，包括单周期运动、倍周期分岔和混沌运动。激励频率是控制这些行为转换的关键参数。
3. 物理机制：系统在特定频率区间进入混沌，是由于接触状态的切换与系统固有振动模态发生了复杂的相互作用。而在共振频率附近，线性共振效应压制了接触非线性，使系统回归周期运动。

本研究的科学价值在于，为复合材料层合结构的接触非线性动力学分析提供了一种新颖、高效的数值工具。其应用价值体现在为工程中复合材料结构的损伤检测（通过识别非线性振动特征）、振动控制以及抗疲劳设计提供了重要的理论依据和预测手段。

六、 研究亮点

1. 算法创新：首次将半光滑牛顿法系统性地应用于含分层复合材料梁的动态接触振动分析。通过Fischer-Burmeister函数和广义雅可比矩阵的构造，实现了对非光滑接触问题的超线性收敛求解，显著提升了计算的鲁棒性。
2. 模型特色：采用仅具有平动自由度的梁单元进行建模，简化了多层含分层结构的有限元离散过程，提高了计算效率，为复杂分层缺陷的模拟提供了便利。
3. 发现系统性：全面、系统地揭示了含分层梁在宽频带激励下的非线性动力学行为演化图谱（周期→倍周期→混沌→周期→混沌），明确了激励频率的核心调控作用，深化了对这类结构非线性机理的认识。
4. 验证充分：通过与前四阶固有频率和动力时程响应的Ansys解进行定量对比，有力地证明了所开发方法的精度，增强了研究结论的可信度。

七、 其他有价值内容

研究在算法实现细节上提供了宝贵信息，例如自适应正则化策略和基于Merit函数的线搜索全局化策略的具体实施方法，这对于其他学者复现或发展类似算法具有参考价值。此外，研究明确指出了该模型目前仅考虑了法向接触，未引入切向摩擦约束，这为未来的研究指明了可能的扩展方向（如考虑摩擦效应）。