1 引言

结合面作为机械系统中一种固有的结构形式，在机械系统正常运行过程中起着传递运动、载荷和能量的作用，对机械系统整体动态特性产生显著的影响^[1]。作为一种复杂的机械系统，航空发动机静子机匣是分段机匣通过螺栓连接而成，机匣安装边结合面决定和影响着整体机匣的刚度、阻尼和稳定性。同时，在工作中的结合面摩擦、间隙接触、变载等非线性因素给机匣整体动力学建模带来了很大的困难^[2]。目前，国内外在工程设计时，通常忽略掉机匣螺栓连接处的接触非线性对其动力特性的影响^[3]。但是，实际发动机机匣结构中采用螺栓法兰连接的数量又特别多，若忽略其对动力特性的影响，则计算结果会出现较大的误差。随着人们对结构分析精度要求的日益提高，传统的建模方式难以满足现代航空发动机设计需要^[4]。因此，研究结合面的动态特性，在设计阶段对分析和预测航空发动机机匣的动态特性就显得十分重要。

近年来，国内外学者已经开始对机械结合面动力特性进行研究。李小彭等^{[5, 6]}利用结合面接触刚度的分形模型与材料应变能等效的方法验证了结合面接触刚度分形模型的正确性和结合面等效方法的合理性，并进行了考虑结合面非线性特性的非线性预应力模态分析，证明了非线性预应力模态分析的可行性。傅卫平等^[7]考虑微凸体在加卸载及动态载荷下的变形特征，建立了结合面动态接触模型，分析了结合面面压、动态位移幅值及振动频率对动态接触刚度的影响规律，研究表明：动态接触刚度增量随法向面压及动态位移幅值的增大而非线性增大。在研究结合面对螺栓法兰连接结构的固有特性影响时，一些学者提出使用薄层单元来模拟法兰连接结构中结合面的接触特性^[8~10]。但薄层单元的参数一般需要通过有限元方法优化求得，且优化过程复杂，效率低。艾延廷等^[11]基于有限元建模方法分析了装配体结合面法向接触刚度对装配体振动模态的影响，研究结果表明，考虑法向接触刚度的有限元等效模拟计算结果与解析解非常接近，不考虑接触刚度则误差较大。

Schwingshackl等^[12]建立了螺栓连接非线性模型，研究了结合面载荷分布对其结构动态特性的影响。Ma^[13]提出通过对比结构有无螺栓连接的动力学特性，从而辨识出结合面的非线性参数。Dhupia等^[14]以立柱-主轴箱为对象并建立了非线性模型，采用非线性响应耦合法进行分析，发现结合面的非线性特性使其结构固有频率和振幅峰值发生显著变化。

这些研究大多都是对结合面建模、载荷分布进行了研究，没有全面考虑结合面非线性特性对结构动态特性的影响，尤其是对薄壁机匣类型的法兰螺栓连接结构。本文依据有限元法处理非线性接触的基本理论，用拉格朗日乘子法来限制由于结合面接触引起的位移协调约束。首先以“L”型螺栓连接梁为研究对象，利用有限元分析软件ANSYS建立了考虑非线性因素的结合面接触非线性模型，进行了动力学分析，通过Newmark-β法求得位移响应。然后以薄壁机匣法兰螺栓连接模型为研究对象，分析了摩擦非线性和预紧力大小对其激励方向和非激励方向动力学响应的影响。最后利用单变量法和预应力模态法对薄壁机匣在不同摩擦系数下进行了模态分析，模态试验结果与预应力模态计算结果相符，验证了预应力模态法处理结合面非线性因素的有效性。

2 处理非线性接触问题的基本原理

本文利用有限元软件ANSYS进行有限元建模及动力学特性模拟计算。利用ANSYS求解非线性摩擦接触问题时，先将法兰螺栓连接结构划分网格后，法兰之间的相互结合面可以分别定义为接触面和目标面，从而对应生成接触单元和目标单元。当两个结合面接触之后，结合面节点的位移变量必须与目标面的位移变量同步，以便满足结合面之间的粘着接触和滑动接触条件。这种同步性体现在对应的接触节点位置上，必须要满足几何约束，可以通过迭代求解来完成。

本文用拉格朗日乘子法来限制由于结合面接触引起的位移协调约束，在时间t+Δt时的第i次迭代的所得平衡方程^[15]如下

$\begin{array}{l} \left( {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&0\\ 0&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\rm{\Delta }}{{\mathit{\boldsymbol{\ddot U}}}^{{\rm{(}}\mathit{i}{\rm{)}}}}}\\ 0 \end{array}} \right) + \left( {\begin{array}{*{20}{c}} \mathit{\boldsymbol{C}}&0\\ 0&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\rm{\Delta }}{{\mathit{\boldsymbol{\dot U}}}^{\left( i \right)}}}\\ 0 \end{array}} \right) + \\ \left\{ {\left( {\begin{array}{*{20}{c}} {^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{K}}^{i - 1}}}&0\\ 0&0 \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0&{^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{K}}_\lambda }^{i - 1}}\\ {^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{K}}_\lambda }^{{{\left( {i - 1} \right)}^{\rm{T}}}}}&0 \end{array}} \right)} \right\}\left( {\begin{array}{*{20}{l}} {{\rm{\Delta }}{\mathit{\boldsymbol{U}}^{\left( i \right)}}}\\ {{\rm{\Delta }}{\lambda ^{\left( i \right)}}} \end{array}} \right) = \\ \left( {\begin{array}{*{20}{c}} {^{t + {\rm{\Delta }}t}\mathit{\boldsymbol{R}}}\\ 0 \end{array}} \right) - \left( {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&0\\ 0&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {^{t + {\rm{\Delta }}t}{{\mathit{\boldsymbol{\ddot U}}}^{\left( {i - 1} \right)}}}\\ 0 \end{array}} \right) - \left( {\begin{array}{*{20}{c}} \mathit{\boldsymbol{C}}&0\\ 0&0 \end{array}} \right)\cdot\\ \left( {\begin{array}{*{20}{c}} {^{t + {\rm{\Delta }}t}{{\mathit{\boldsymbol{\dot U}}}^{\left( {i - 1} \right)}}}\\ 0 \end{array}} \right) - \left( {\begin{array}{*{20}{c}} {^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{F}}^{\left( {i - 1} \right)}}}\\ 0 \end{array}} \right) + \left( {\begin{array}{*{20}{l}} {^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{R}}_{\rm{c}}}^{\left( {i - 1} \right)}}\\ {^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{c}}}^{\left( {i - 1} \right)}} \end{array}} \right) \end{array} $

(1)

式中M，C分别代表结构的质量、阻尼矩阵；^t+ΔtK^i-1代表（i-1）次迭代后的材料和几何非线性的切向刚度矩阵；^t+ΔtK_λ^i-1代表接触矩阵（i-1）次迭代后包含相容表面位移的约束方程；^t+ΔtR代表t+Δt时外载荷矢量；^t+ΔtF^(i-1)代表（i-1）次迭代后节点应力等效于单元应力矢量；^t+ΔtR_c^i-1代表（i-1）次迭代后满足库仑摩擦定律的接触力矢量；^t+ΔtΔ_c代表（i-1）次迭代后接触节点对目标面的穿透量矢量；${\rm{\Delta }}{U^{\left( i \right)}}, {\rm{\Delta }}{\mathit{\boldsymbol{\dot U}}^{\left( i \right)}}, {\rm{\Delta }}{\mathit{\boldsymbol{\ddot U}}^{\left( i \right)}}$分别代表第i次迭代所有节点的加速度、速度、位移的增量矢量Δλ⁽ⁱ⁾代表第i次迭代接触约束的Lagrange乘子；$^{t + {\rm{\Delta }}t}{\mathit{\boldsymbol{\dot U}}^{\left( {i - 1} \right)}}{, ^{t + {\rm{\Delta }}t}}{\mathit{\boldsymbol{\ddot U}}^{\left( {i - 1} \right)}}$代表第i-1次迭代后的速度和加速度矢量。使用Newmark-β法对方程（1）求解，可以求得各时间步的位移变量，从而得到结构所有节点的时间历程位移响应和加速度响应。

3 “L”型螺栓连接梁非线性建模与动力学响应

首先，建立“L”型螺栓连接梁结构的线性动力学模型和非线性摩擦接触模型以便进行对比，模型为两个长200mm，宽14mm，厚4mm，法兰高18mm的“L”型梁通过两个螺栓连接，螺栓为M6×30的标准件，如图 1所示。材料为45号钢。

计算中分别建立三种接触方式的模型：第一种为完全Bonded（绑定）模型，即螺栓、螺母、法兰互相之间的所有结合面不存在滑移完全一体化；第二种模型为螺栓位置Bonded模型，即仅在螺栓头、螺母与法兰连接位置相互绑定；最后一种模型是摩擦接触非线性模型（以下简称全摩擦模型）。考虑螺栓、螺母和法兰结构之间存在相互的接触摩擦关系，它们通过相互的摩擦接触关系将两个“L”型结构连接起来，摩擦系数取0.12，在摩擦选项中刚度矩阵设置为非对称类型。螺栓预紧力是通过预紧单元Prets179得到的，过程如下:对模型划分网格后，在螺栓中间部位节点切开并添加预紧单元将切开的两部分连起来，最后设置螺栓预紧力即可完成对螺栓施加预紧力。边界条件为一端固定，一段自由。依据公式（2）^[16]

$F = \frac{M}{{0.17d}} $

(2)

施加1N·m预紧力矩，其次，在分析模块中选择预应力计算开启状态，进行静力分析。完成静力分析以后，对前两种模型进行模态计算，表 1为前4阶模态频率，通过对比可以看出螺栓位置Bonded模型各阶模态频率都比完全Bonded模型低，差值不同。由于L型螺栓连接梁模型简单尺寸小，模态分析时不易看出区别，从图 2可以看出两种模型的振型没有明显的差异。

因为全摩擦模型中存在非线性摩擦接触关系，不能直接进行模态计算，所以先通过瞬态时间历程分析中的完全法进行分析计算，因为完全法不忽略非线性摩擦接触的设置。激励位置和数据提取点如图 1所示。分别对三种模型的Y方向施加脉冲激励，施加脉冲载荷如图 3所示，并进行结构瞬态动力学时间历程响应分析。

三种模型在该脉冲载荷激励下结构自由端中间位置沿激励方向的位移响应曲线如图 4所示。全摩擦模型与螺栓位置Bonded模型的位移响应轨迹比较接近，而完全Bonded模型的位移响应轨迹与其他两种模型有明显的差异。所以选择前两者其中一个（全摩擦模型）与完全Bonded模型进行对比分析，如图 5所示，全摩擦模型与完全Boned模型相比，前者在施加冲击后加速度呈指数形式衰减，表现出阻尼非线性特征。这是由于全摩擦模型结合面存在粘着、滑移运动，结构动能以摩擦热的形式耗散并且较大；而完全Bonded模型的能量耗散有限。

4 薄壁机匣法兰螺栓连接非线性建模与动力学响应

为进一步研究结合面非线性特性对法兰螺栓连接结构动态性能的影响，建立了复杂的薄壁机匣法兰螺栓连接的有限元模型如图 6所示，两段机匣上的法兰通过螺栓连接，机匣的两个筒体尺寸完全一样，筒体内径249mm，外径255mm即壁厚为3mm。筒体底部法兰环内径249mm，外径321mm，厚度8mm，单段机匣筒体高304mm，螺栓直径10mm。机匣模型的材料为45号钢，其材料参数：弹性模量E= 200GPa，泊松比v=0.3，密度ρ=7850kg/m³。单元选择20节点的186单元，网格为六面体。两个机匣法兰和螺栓与上下法兰的结合面节点上建立接触对，接触单元分别是Contact174和Targe170。设置常用的钢摩擦系数为0.12，在摩擦选项中刚度矩阵设置为非对称类型。在计算中分别采用三种连接方式：螺栓位置Bonded模型、全摩擦模型和完全Bonded模型。

考虑螺栓、螺母以及法兰结构之间的接触关系，根据相互之间存在的摩擦接触关系构成接触对，薄壁机匣共有12个螺栓，每个螺栓头和螺母分别与法兰形成一个接触对（12×2=24），两个法兰安装边之间形成一个接触对，所以共生成25个接触对，螺栓具体位置如图 7所示。完成非线性接触静力分析以后，选择瞬态分析模块，选择完全法求解，设置预应力开启，将前面考虑摩擦接触特性得到静力学分析结果直接导入到瞬态分析过程中。分别对螺栓位置Bonded模型、完全Bonded模型和全摩擦模型结构接近法兰部位施加沿X方向的脉冲激励，对结构进行瞬态动力学时间历程响应分析。

表 2给出了预紧力矩为20N·m的条件下三种模型自由端最大位移及相应的误差比。可以看出：螺栓位置Bonded模型和全摩擦模型之间的最大位移误差相对较小，最大误差比不超过14.2%，表明螺栓位置Bonded模型和全摩擦模型沿各个方向的位移响应比较接近；但全摩擦模型与完全Bonded模型位移响应误差较大，沿X方向的最大误差为14.1%，但是沿Y方向最大误差则达到30.1%，沿Z方向也非常明显为27.4%，表明完全Bonded模型与全摩擦模型相比，摩擦非线性因素对于顶部X方向（激励方向）响应影响较小；但对其非激励方向（Y方向和Z方向）位移响应影响较大，全摩擦模型和完全Bonded模型之间的位移响应存在较大差异。

下面计算分析预紧力矩大小对该机匣法兰螺栓连接结构动力学响应的影响，螺栓预紧力矩分别增加到25，30，35，40，45N·m，表 3给出了完全bonded模型和加大预紧力矩的全摩擦模型的最大位移响应及对应的误差。可以看出，在预紧力矩加大到45N·m后结构全摩擦模型和完全bonded模型沿X方向最大位移误差从13.4%降至10.8%；Y方向和Z方向位移响应误差变化更加明显，分别从27.1%，25.4%降至为13.7%和12.9%。由此说明预紧力矩大小对结构沿激励方向即X方向的动力响应影响较小，但对Y方向以及Z方向的影响较大，说明预紧力大小对结构非激励方向的动力学响应有显著影响。

5 结合面摩擦系数对薄壁机匣模态影响 5.1 不同摩擦系数下薄壁机匣的模态分析

有限元里模态计算本身是忽略非线性因素的，但是非线性因素可以通过预应力模态法以边界条件的形式进行模态计算。先对薄壁机匣进行网格划分并施加预紧力，过程同上述动力学分析的步骤。完成非线性接触静力分析以后，选择模态分析模块，选择非对称求解器，设置预应力模态分析开启，将前面考虑摩擦接触特性得到的薄壁机匣静力学分析结果作为边界条件直接导入到其预应力模态分析过程中。

本文采用单变量法和预应力模态法研究结合面摩擦系数大小对薄壁机匣振动特性的影响。依据螺栓扭矩的国家标准范围设定螺栓拧紧力矩为40N·m保持不变，在建立接触对过程中根据国家标准范围^[17]分别设定4种不同的摩擦系数（μ=0.10，0.12，0.20，0.30）进行预应力模态分析，求得薄壁机匣的各阶固有频率如图 8所示。可以看出，结合面非线性特征对薄壁机匣各阶固有频率的影响十分显著。尤其在高阶，摩擦系数为0.30的固有频率比摩擦系数为0.10的固有频率低了100多Hz。主要原因为结合面有摩擦，实际接触面积减小，接触点塑性变形的比例变大，从而结合面阻尼增加刚度减小，与参考文献[18]中接触压力随着摩擦系数变大而减小的趋势相同，接触刚度与接触压力成正比关系，所以薄壁机匣法兰螺栓连接结构的固有频率随结合面刚度的减小而减小。

5.2 薄壁机匣模态试验

采用锤击法进行薄壁机匣试验件多点激励多点响应的模态试验。试验采用LMS振动测量和分析系统及其配套的Test Lab软件、脉冲力锤、加速度传感器等。试验所用的薄壁机匣材料为45号钢，本文根据试验件实际情况取摩擦系数为0.12，与预应力模态分析时设置相同。试验时设定扭矩扳手力矩为40N·m，试验件被划分为520个测试点，每圈20个点均匀分布，一共26圈。试验模型与模态试验台如图 9和图 10所示。

由图 11和表 4看可看出，模态试验与预应力模态计算的振型吻合，固有频率与有限元模拟的固有频率有一定的误差，每一阶误差各不相同，最大为5.1%。主要原因为：（1）建模时对模型作了相应的简化，且不计算其转动惯量和剪切变形的影响。（2）有限元里材料的设置与实际试验模型有一定偏差，如弹性模量、泊松比和密度等参数。（3）实验本身有不可避免的人为误差，如力锤的真实敲击点与试验前设置的敲击点有一定的偏差。因此，有限元模拟结果与试验结果有一定的误差，但误差在可接受范围，验证了预应力模态法处理结合面非线性因素的可行性，同时也验证了该摩擦接触非线性模型的有效性。

6 结论

通过本文研究，得出如下结论：

（1）研究发现螺栓位置Bonded和全摩擦两种模型的位移响应轨迹比较一致。全摩擦模型与完全Boned模型相比，前者在施加冲击后加速度呈指数形式衰减较快，表现出阻尼非线性特征。

（2）摩擦非线性因素对本文的薄壁机匣沿激励方向响应的影响较小，最大误差比不超过14.2%；但对非激励方向响应的影响较大，沿Y方向最大误差则达到30.1%，沿Z方向也非常明显为27.4%；加大预紧力矩可使全摩擦模型和完全bonded模型之间的位移响应沿X方向最大位移误差从13.4%降至10.8%；Y方向和Z方向位移响应误差变化更加明显，分别从27.1%，25.4%降至为13.7%和12.9%。

（3）摩擦系数对薄壁机匣固有频率有重要影响，尤其是高阶固有频率。模态试验与预应力模态法计算的各阶振型吻合、模态固有频率与有限元计算的固有频率误差在5.1%以内，验证了预应力模态法处理结合面非线性因素的可行性，同时也验证了该摩擦接触非线性模型的有效性。