考虑接触面摩擦的柔性套齿连接结构刚度数值建模与非线性机理分析

Numerical Modeling and Nonlinear Mechanism Analysis of Spline Coupling Stiffness Considering Contact Surface Friction

摘要

以发动机柔性套齿连接结构为对象，通过在套齿连接结构轴段微元各啮合点处引入接触-摩擦模型，并结合整体结构的几何关系、受力平衡关系和数值迭代法提出了一种套齿非线性刚度的建模与计算方法，通过与ANSYS有限元模型结果进行对比，验证了所提计算方法的高求解效率和求解准确性。利用所提求解方法详细分析了套齿连接结构刚度非线性特征及关键参数的影响规律，并结合啮合齿对接触状态的变化规律揭示了其刚度非线性变化的内在机理。结果表明：套齿连接结构刚度随线位移的增加逐渐下降，但下降速率逐渐减小，此过程中套齿连接结构啮合齿接触面积减小是刚度降低的根本原因。接触面摩擦因数、啮合刚度、套齿齿宽以及扭矩增加时，套齿连接结构刚度均有所提高，其中扭矩的影响最为明显。

Abstract

Taking the flexible spline coupling structure of an engine as the research object， a contact-friction model was introduced at the meshing points of all shaft segment in the spline coupling. Combined with the geometric relationships of the structure， force equilibrium conditions， and numerical iteration methods， a modeling and calculation approach for the nonlinear stiffness of the spline couplings was proposed. The high solving efficiency and accuracy of the proposed calculation method were validated through comparisons with ANSYS finite element model results and experimental results. Using the proposed method， detailed analyses were conducted on the nonlinear stiffness characteristics of the spline couplings and the influence patterns of key parameters. The intrinsic mechanism of the nonlinear stiffness variation was revealed by combining the changing patterns of contact states in the meshing tooth pairs. Results indicate that the stiffness of the spline coupling structures gradually decreases with the increasing linear displacement， while the decreasing rate progressively diminishes. The reduction in contact area of meshing teeth during this process was identified as the fundamental cause of the stiffness degradation. The stiffness of the spline couplings shows improvement with the increase of the friction coefficient， the meshing stiffness， the spline width， and the applied torque. Among these parameters， the torque exhibits the most significant influence.

Graphical abstract

关键词

套齿连接结构 / 刚度特性 / 接触-摩擦模型 / 数值建模

Key words

spline connection / stiffness characteristic / contact-friction model / numerical modeling

0 引言

套齿连接结构由于承载能力强、结构简单紧凑、可靠性高且易于安装拆卸等优点而被广泛应用于航空发动机低压转子的风扇轴和涡轮轴连接^［1］。在发动机工作过程中，套齿连接结构不可避免地承担涡轮轴与风扇轴之间的扭矩、剪力和弯矩等多轴载荷，导致连接结构中的接触界面发生滑移、黏滞、分离等复杂的接触状态变化，诱发套齿连接结构的刚度非线性特征，从而造成转子振动特性偏离设计状态，严重时甚至引起转子及整机系统的剧烈振动。开展套齿连接结构非线性力学特性研究，并建立相应的快速计算模型，是进一步预测含套齿连接结构转子动力学特性的重要前提，对指导发动机套齿连接结构动力学设计具有重要意义。

得益于有限元理论与计算机技术的发展，国内外学者利用有限元仿真对套齿连接结构的接触载荷与结构刚度特性开展了研究。王庆国等^［2］基于有限元法计算得到了齿面应力，并发现套齿连接结构的接触区域所受载荷不均匀，输入端接触压应力、齿根弯曲应力较大。BARROT等^［3］利用有限元方法并结合解析分析，对套齿连接结构轴线方向的载荷进行了讨论，结果表明套齿连接结构接触面压力沿轴线呈U形分布。在此基础上，HONG等^［4-6］建立了直齿和斜齿载荷求解的半解析力学模型，重点关注了扭矩载荷影响下的接触齿面应力分布特征。CURÀ等^［7］采用数值仿真方法研究了套齿联轴器的反力矩、疲劳损伤和微动磨损损伤规律。针对含多种类型接触面的套齿结构，李伦绪等^［8］重点分析了其弯曲刚度非线性特征及其产生机理，基于套齿连接结构实体有限元模型详细讨论了配合间隙、拧紧力矩、扭矩的影响规律。李英杰等^［9］基于实体有限元模型重点分析了套齿连接结构在周期循环载荷下接触刚度的变化规律以及影响因素，并建立了试验器进行验证。谭援强等^［10］基于有限元法分析了套齿啮合点接触压力和滑移的分布规律以及关键参数的影响规律。

有限元方法能够精细地考虑套齿连接结构复杂结构，但通常需要划分较多网格，因而计算耗时，且在涉及接触非线性计算时还通常存在计算发散问题，为此部分学者从解析/数值分析角度对套齿连接结构力学特性进行分析。MARMOL等^［11］从套齿的变形出发推导出了套齿连接结构的刚度，但该研究中假设接触面载荷沿轴线方向均匀分布，存在较大的局限性。朱彬等^［12］使用悬臂梁法计算得到了套齿连接结构的刚度，并与有限元仿真对比，研究发现啮合齿对套齿连接结构整体刚度存在不可忽视的影响。ZHAO等^［13］计算了套齿连接结构中一个齿的刚度，研究表明，在小载荷下套齿连接结构的刚度是关于齿数、单齿刚度及压力角的常数，而当载荷较大时，其刚度表现出明显的非线性特征。VAUJANY等^［14］通过啮合齿之间的间隙判断潜在的接触点并划分局部网格，结合有限元计算和插值法，提出了一种分析套齿连接结构瞬时载荷接触的模型，并指出轴段的几何参数会影响套齿连接结构的刚度。HONG等^［15］基于齿的几何构型以及荷载分布建立了半解析的套齿连接结构刚度矩阵，通过矩阵的分块反演方法减少了计算时间，但由于计算时需要多次迭代柔度矩阵，使得计算时间较长。上述工作主要是利用短梁理论计算齿的刚度，并进一步综合所有齿的刚度得到套齿连接结构整体刚度，未关注齿对啮合状态变化带来的影响。YU等^［16］针对现有模型不足，采用切片方法建立了柔性套齿连接结构刚度计算的解析模型，该模型能够考虑随外载荷变化时各齿啮合状态的变化。在YU等模型基础上，ZHANG 等^［17］进一步发展了可考虑初始不对中的套齿连接结构刚度计算模型，但其模型中均假设接触界面光滑无摩擦，与真实情况具有一定差异。

总地来说，现有套齿连接结构模型中，有限元法虽然能够获得准确的套齿连接结构非线性刚度，但受限于模型规模和计算效率，不能直接用于套齿连接结构转子动力学计算。与之相对，数值方法在套齿连接结构转子动力学计算时显现出较大应用潜力，然而目前的套齿连接结构数值模型对齿啮合状态的考虑较少，更未纳入接触面摩擦的影响。有鉴于此，本文以发动机柔性套齿连接结构为对象，结合轴段微元的啮合接触-摩擦模型和整体结构受力平衡关系，提出一套针对套齿连接结构刚度的建模方法，不仅能考虑套齿连接结构各啮合齿对的啮合状态及其沿轴线方向的变化，还能反映接触面摩擦特征。通过与有限元模型结果和试验结果对比，验证所提建模方法的精度和效率，在此基础上对套齿连接结构的刚度非线性机理及关键参数影响规律进行了讨论，所提模型易于修改关键参数以预测不同结构尺寸参数的套齿连接结构刚度特性。

1 含接触面摩擦的套齿刚度建模

目前航空发动机中采用的套齿连接结构分为柔性套齿连接结构和刚性套齿连接结构两类。其中，柔性套齿连接结构具有结构简单、装配方便及不同心适应能力强等优势，在俄系发动机及我国部分主力战机发动机中均有应用。某型航空发动机低压转子的柔性套齿连接结构如图1所示，主要由含内齿的风扇轴、含外齿的涡轮轴和锁紧螺母三部分组成。涡轮在工作过程中产生的巨大扭矩使外齿和内齿的齿面接触并压紧，内外齿啮合传递扭矩，锁紧螺母施加轴向力实现轴向压紧，提高结构整体连接的稳定程度。

1.1 套齿连接结构啮合接触-摩擦模型

考虑套齿连接结构的结构与载荷特征，构建模型时作以下假设：

1）弧形垫片允许涡轮轴和风扇轴之间发生一定的相对移动，它对套齿连接结构整体的刚度贡献不大，而套齿连接结构受载变形过程中，轴段变形相对于啮合齿变形较小，因此，建模时仅考虑啮合齿刚度，忽略轴段刚度^［18］。

2）考虑到连接结构小变形特点，且啮合齿通常存在齿侧间隙^［6］，建模时仅考虑单侧齿面接触。

3）每对啮合齿在接触面上保持直线接触，当套齿连接结构仅受扭矩时，其接触位置为节径^［19］（下称“啮合圆”）。

4）弧形垫片仅传递轴向载荷，对套齿连接结构角向变形的约束作用很小，刚度建模时其影响不予考虑。

基于上述假设，建立了套齿连接结构刚度计算的简化力学模型，如图2所示。模型中，风扇轴左端固定约束，涡轮轴右侧受到横向剪力/弯矩和扭矩，在不失一般性的前提下，在y方向上施加横向载荷。在外部荷载作用下的变形如图2b所示，可看成涡轮轴绕某一中心点（x_a，0）的转动，绕z轴的转角为θ_z。

为了得到外部载荷与变形的关系，将套齿连接结构划分为厚度为dx的轴段微元（图2b）。首先确定单个轴段微元的受力，然后沿轴向积分得到整个套齿连接结构的载荷。各齿的编号规律如图3所示，则各齿中线所在的角度如下：

{\phi }_i=\mathrm{2}\mathrm{\pi }(i-\mathrm{1})/z（1）

式中：{\phi }_i为第i个齿的中线所在的角度；i为齿的编号；z为齿数。

此时由于涡轮轴的偏转，该轴段微元中的涡轮轴将在y轴方向上发生偏移\delta (x)，表示为

\delta (x)=(x-x_{\mathrm{a}}){\theta }_z（2）

式中：x为微元的轴向坐标；x_a为转动中心的轴向坐标；θ_z 为绕z轴的转角。

此时，该微元对应的各齿的啮合距离将发生变化（图3），涡轮轴未偏移前各齿啮合半径假设为节圆半径r_{\mathrm{m}}，当涡轮轴套齿对应截面产生偏移时，该截面内各个齿的啮合距离r_{\mathrm{m}i}近似为

r_{\mathrm{m}i}=r_{\mathrm{m}}+\delta \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_i（3）

对于任一微元，当截面发生偏移\delta时，该截面内各个齿将沿啮合力方向产生附加的侵入量以及切向（与啮合力方向垂直）的相对滑移量，如图4a所示，对于第i个齿，侵入量和滑移量分别为

\left.\begin{array}{l}{\Delta }_{\mathrm{b}i}^n=\delta \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_i\\ {\Delta }_{\mathrm{b}i}^t=\delta \mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_i\end{array}\right\}（4）

式中：{\Delta }_{\mathrm{b}i}^n为涡轮轴套齿偏移时第i个齿在接触点法线方向（n）上的侵入量；{\Delta }_{\mathrm{b}i}^t为涡轮轴套齿偏移时第i个齿在接触点切线方向（t）上的相对滑移量，其正方向按图4a所示的局部坐标系定义。

图5所示为齿的几何参数，{\alpha }_i为第i个齿啮合点处啮合点外法线与竖直y向的夹角，且{\alpha }_i={\phi }_i+{\beta }_i，其中{\beta }_i为第i个齿啮合点处的啮合外法线与齿的中线的夹角（锐角），可由齿的几何参数获得。

在涡轮扭矩作用下，涡轮轴还将产生绕中心线的扭转变形，所产生的扭转角记为\varphi，取绕+x轴旋转为正，如图4b所示。由扭转产生的沿啮合方向的侵入量和切向（与啮合力方向垂直）的相对滑移量为

\left.\begin{array}{l}{\Delta }_{ti}^n=\varphi r_{\mathrm{m}i}\mathrm{s}\mathrm{i}\mathrm{n}{p}_i\\ {\Delta }_{ti}^t=-\varphi r_{\mathrm{m}i}\mathrm{c}\mathrm{o}\mathrm{s}{p}_i\end{array}\right\}（5）

式中：{\Delta }_{ti}^n为涡轮轴套齿扭转时第i个齿在接触点法线方向上的侵入量；{\Delta }_{ti}^t为涡轮轴套齿扭转时第i个齿在接触点切线方向上的相对滑移量；p_i 为第i个齿的啮合点位置处的压力角。

综上，涡轮轴套齿各个齿在横向载荷和扭矩共同作用下，各齿在啮合点处沿啮合力方向的总侵入量{\Delta }_i^n和总滑移量{\Delta }_i^t为

\left.\begin{array}{l}{\Delta }_i^n={\Delta }_{\mathrm{b}i}^n+{\Delta }_{ti}^n\\ {\Delta }_i^t={\Delta }_{\mathrm{b}i}^t+{\Delta }_{ti}^t\end{array}\right\}（6）

对于任一齿的微元，由于上述法向侵入和切向滑动，在接触点处产生法向正压力和切向摩擦力，结合Jenkins单元引入一种可考虑法向力变化的接触-摩擦单元，建立套齿微元啮合的接触-摩擦模型，如图6所示。n向的弹簧表示法向啮合刚度k_{\mathrm{g}i}^n，它是内外啮合齿刚度及接触刚度的串联刚度，具体计算公式见文献［16］；t向弹簧称为“迟滞弹簧”，其刚度表示为k_{\mathrm{g}i}^t，用于模拟接触面滑移行为并计算摩擦力，其值可认为与法向啮合刚度成正比，即k_{\mathrm{g}i}^t=\beta k_{\mathrm{g}i}^n，β为比例系数。

根据上述力学模型，当接触点处存在相对变形时，沿啮合点法向方向的力为

\mathrm{d}{F}_i^n=-H({\Delta }_i^n)k_{\mathrm{g}i}^n{\Delta }_i^n\mathrm{d}x（7）

其中，H({\Delta }_i^n)为Heaviside函数，形式如下：

H({\Delta }_i^n)=\left\{\begin{array}{c}\mathrm{0}{\Delta }_i^n\le \mathrm{0}\\ \mathrm{1}{\Delta }_i^n>\mathrm{0}\end{array}\right.（8）

当法向力为0时，切向载荷亦为0。对于第i个齿，当法向力不为0时，沿啮合点切向方向存在两种状态：若啮合点切向位移{\Delta }_i^t与滑块位移{\epsilon }_i之差较小时，“迟滞弹簧”拉力小于临界摩擦力\mu \mathrm{d}{F}_i^n（\mu为接触面摩擦因数），此时滑块不发生滑动，即“阻滞状态”，摩擦力大小与“迟滞弹簧”的拉力大小相等，方向相反；若{\Delta }_i^t与{\epsilon }_i之差增大到一定程度，“迟滞弹簧”拉力达到\mu \mathrm{d}{F}_i^n时，滑块发生滑动，即“滑移状态”，此时摩擦力等于临界摩擦力。据此可得到摩擦力的表达式：

\mathrm{d}{F}_i^t=\left\{\begin{array}{l}-H({\Delta }_i^n)k_{\mathrm{g}i}^t({\Delta }_i^t-{\epsilon }_i)\mathrm{d}x{k}_{\mathrm{g}i}^t({\Delta }_i^t-{\epsilon }_i)\mathrm{d}x\le \mu \mathrm{d}{F}_i^n\\ -H({\Delta }_i^n)\mu \mathrm{d}{F}_i^n\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}({\dot{\epsilon }}_i)k_{\mathrm{g}i}^t({\Delta }_i^t-{\epsilon }_i)\mathrm{d}x>\mu \mathrm{d}{F}_i^n\end{array}\right.（9）

基于时间推进法的基本思想求解啮合点摩擦力，基本过程如下：

1）获取第m步的啮合点切向滑移量{\Delta }_i^t(m)以及第m-\mathrm{1}步的滑块位移{\epsilon }_i(m-\mathrm{1})；若是程序的第一次迭代计算，则滑块位移{\epsilon }_i的初值设为0。

2）假设计算中第m步为“阻滞状态”，根据下式计算得到预估的摩擦力\mathrm{d}{F}_{i,\mathrm{p}\mathrm{r}\mathrm{e}}^t：

\mathrm{d}{F}_{i,\mathrm{p}\mathrm{r}\mathrm{e}}^t=-k_{\mathrm{g}i}^t({\Delta }_i^t(m)-{\epsilon }_i(m-\mathrm{1}))\mathrm{d}x（10）

3）将预估的摩擦力与此时的接触面临界摩擦力相比较，若小于临界摩擦力，则第m步中啮合点为“阻滞状态”，此时实际摩擦力与预估的摩擦力相等；反之则表示啮合点为“滑移状态”，此时实际摩擦力等于临界摩擦力，即

（11）

4）时间步m时的滑块位移的计算式如下：

{\epsilon }_i(m)={\Delta }_i^t(m)+\mathrm{d}{F}_i^t/k_{\mathrm{g}i}^t（12）

5）重复步骤2）~4），求解下一步的摩擦力，直至完成预先规划的载荷步范围。

根据啮合点处法向力和切向摩擦力，进一步可得到该微元所有齿的合力：

\left.\begin{array}{l}\mathrm{d}{F}_y={\displaystyle \sum _{i=\mathrm{1}}^z}\mathrm{d}{F}_i^n\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_i+{\displaystyle \sum _{i=\mathrm{1}}^z}\mathrm{d}{F}_i^t\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_i\\ \mathrm{d}{F}_z={\displaystyle \sum _{i=\mathrm{1}}^z}\mathrm{d}{F}_i^n\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_i+{\displaystyle \sum _{i=\mathrm{1}}^z}-\mathrm{d}{F}_i^t\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_i\\ \mathrm{d}{M}_x={\displaystyle \sum _{i=\mathrm{1}}^z}\mathrm{d}{F}_i^n\mathrm{c}\mathrm{o}\mathrm{s}(p_i)r_{\mathrm{m}i}+{\displaystyle \sum _{i=\mathrm{1}}^z}-\mathrm{d}{F}_i^t\mathrm{s}\mathrm{i}\mathrm{n}(p_i)r_{\mathrm{m}i}\end{array}\right\}（13）

根据微元受力结果，沿轴向进行积分可得到所有啮合齿形成的合力与合力矩：

\left.\begin{array}{l}{F}_y={\int }_{\mathrm{0}}^L\mathrm{d}{F}_y\hspace{1em}\hspace{1em}{F}_z={\int }_{\mathrm{0}}^L\mathrm{d}{F}_z\\ M_x={\int }_{\mathrm{0}}^L\mathrm{d}{M}_x\hspace{1em}\hspace{1em}{M}_y={\int }_{\mathrm{0}}^L-x\mathrm{d}{F}_z\hspace{1em}\hspace{1em}{M}_z={\int }_{\mathrm{0}}^{L}x\mathrm{d}{F}_y\end{array}\right\}（14）

式中：F_y 、F_z 分别为套齿连接结构在y和z方向的反力；M_x 、M_y 、M_z 分别为套齿连接结构上绕x轴、y轴和z轴的力矩。

1.2 平衡方程及其迭代求解

涡轮轴套齿所有齿面上的合力和合力矩与涡轮轴右侧所受外载荷相平衡，于是得到套齿连接结构的受力平衡方程，外载荷仅为扭矩和横向剪力时的受力平衡方程为

\mathit{F}\left(\mathit{X}\right)=\left\{\begin{array}{l}{F}_y+{\widehat{F}}_y\\ F_{y}L+M_z\\ T+M_x\end{array}\right.（15）

\mathit{X}=[x_{\mathrm{a}}\hspace{1em}{\theta }_z{\hspace{1em}\varphi ]}^{\mathrm{T}}

式中：{\widehat{F}}_y为给定的横向剪力；{\widehat{M}}_z为给定的横向弯矩。

求解式（15）可获得套齿连接结构线刚度及齿面接触载荷。外载荷仅为扭矩和横向弯矩时的受力平衡方程为

\mathit{F}\left(\mathit{X}\right)=\left\{\begin{array}{l}{F}_y\\ {\widehat{M}}_z+M_z\\ T+M_x\end{array}\right.（16）

对其进行求解可得到套齿连接结构角刚度及齿面接触载荷。上述平衡方程为非线性方程，采用牛顿迭代法进行求解，迭代格式如下：

{\mathit{X}}^{(k+\mathrm{1})}={\mathit{X}}^{(k)}-{(\frac{\partial \mathit{F}({\mathit{X}}^{(k)})}{\partial \mathit{X}})}^{-\mathrm{1}}\mathit{F}({\mathit{X}}^{(k)})（17）

其中，{\mathit{X}}^{(k)}为迭代过程中第k步结束时的自变量。求解时，初始变量可取{\mathit{X}}_{}^{(\mathrm{0})}=[L/\mathrm{2}\hspace{1em}{\theta }_z^{\xi }\hspace{1em}{\varphi }_{\mathrm{0}}{]}^{\mathrm{T}}，其中，{\theta }_z^{\xi }可以取一个极小的非零实数，{\varphi }_{\mathrm{0}}为仅存在扭矩时的扭转角。代入初值至式（17）进行迭代求解，直至获得满足收敛条件的{\mathit{X}}^{(k)}，然后根据x_a和{\theta }_z^{}求解得到端面线位移{\delta }_{\mathrm{l}}，进而求解得到线刚度k_{\mathrm{l}}和角刚度k_{\theta }（迭代求解流程如图7所示）：

k_{\mathrm{l}}={\widehat{F}}_y/{\delta }_{\mathrm{l}}（18）

k_{\theta }={\widehat{M}}_z/{\theta }_z（19）

2 数值模型的对比验证

2.1 基于有限元模型的对比验证

由于线刚度与角刚度的变化趋势相似，故后续仅对线刚度进行讨论。某型套齿连接结构的结构参数见表1，各齿的编号如图8所示，根据其结构参数建立得到套齿连接结构的实体有限元模型，如图9所示，其中齿宽为图9中套齿连接结构前后端面的距离。在风扇轴左端面施加全约束，在涡轮轴的右端面上施加扭矩和横向剪力，其中扭矩以等效力偶的方式施加在节点上，提取涡轮轴的右端面在横向剪力方向上的位移作为套齿连接结构的线位移。模型采用SOLID185单元建立套齿连接结构实体结构，采用CONTA 174与TARGET 170单元模拟啮合齿的接触面，有限元模型共包含71 282个实体单元、82 282个节点。

根据套齿啮合刚度公式计算得到其单位长度下啮合点法向接触刚度64.1 kN/mm，给定比例系数β=0.4。图10所示为本文模型计算得到的刚度曲线与ANSYS模型获得的结果对比，可以看出：套齿连接结构刚度随线位移的增加呈非线性下降，小位移阶段刚度随线位移的增大迅速下降，当线位移增大到一定程度时刚度逐渐趋于稳定。总体而言，两者变化规律和数值比较吻合，证明了本文模型的正确性。

进一步对比本文模型和有限元方法的求解效率，求解本文模型和实体有限元模型中所使用的CPU型号为i7-12700，内存大小为16 GB。两类模型共计算了32个数据点，求解有限元模型需要约0.5 h，而本文模型仅需要约2 s，可以看出本文模型具有优异的求解效率。此外，使用有限元模型对不同结构尺寸参数的套齿连接进行仿真时需要重新建模、划分网格等操作，而在本文模型中能够快速修改套齿连接的齿数、基圆半径等各项参数并对其力学性能进行预测，因此，本文模型适用于工程方案阶段的快速迭代计算。

2.2 实际套齿模型试验验证

为验证本文数值方法的有效性，设计了一个与发动机低压转子套齿联轴器刚度特征相似的试验装置，试验装置现场图见图11a，套齿结构如图11b所示，该套齿的尺寸参数见表2，所使用的材料为45钢。将砝码10放置在砝码挂钩上，经过空心法兰轴8向外花键轴施加扭矩载荷，而横向剪力通过立式千斤顶9施加，如图11c所示。

试验中，扭矩为150 N·m，剪力载荷最大值为10 kN，将数值计算得到的结果与试验得到的结果进行对比，如图12所示。可以看出本文模型与实验得到的刚度变化规律基本一致，验证了本文模型的有效性，数值模型中的“接触-摩擦模型”能在一定程度上模拟接触面接触状态的变化。但由于存在一定的简化，且装配过程中存在误差，导致数值计算结果与试验结果存在一定误差。

3 接触状态与非线性刚度的机理分析

套齿非线性刚度特征与啮合齿接触状态密切相关，基于前文所建模型，从齿的啮合接触状态出发对其刚度非线性变化特性进行解释。给定扭矩T=1 kN·m，不同载荷下接触面积占比变化情况如图13a所示，“接触面积占比”指处于黏滞/滑移状态的接触区域的面积与齿的总啮合面积之比。可以看出：当外载荷小于4.7 kN时，所有啮合齿接触良好，接触面积占比100%，此时套齿连接结构刚度处于较高阶段，且相对较稳定；当载荷进一步增大，接触面积占比迅速下降，载荷增加到15 kN时约降至60%，结合图10 a，此时线位移约5 μm，对应刚度为3 MN/mm；当载荷增加至25 kN时，接触面积占比降至10%，此时线位移达到0.04 mm，对应刚度小于1 MN/mm，当载荷或线位移进一步增加时，接触面积占比降幅较小，故套齿连接结构刚度趋于稳定。本质上，啮合齿接触区域降低使参与承载的面积减小，进而造成了刚度的非线性下降。

图13b为套齿连接结构涡轮轴扭转角随外载荷变化曲线，可以看出，未受外载荷时扭转角为正，与扭矩方向相同；之后随剪力载荷的增大，扭转角逐渐减小，且减小的速率逐渐增大，当剪力载荷超过16 kN时，扭转角由正值变为负值，这表明该区域内的啮合齿接触平衡状态发生了较为显著的改变。对于传递一定扭矩的套齿连接结构，当套齿连接结构受外载荷时，啮合齿接触状态变化，导致啮合齿形成的载荷合力矩发生改变，为了维持与扭矩平衡，套齿扭矩角必然发生变化，并进一步影响啮合齿接触状态。

图14和图15所示分别为四种剪力载荷下啮合齿接触状态及侵入位移云图，各齿的编号如图8所示。图14中的x_a表示套齿涡轮轴的偏转中心位置，啮合齿的各区域接触时为粉色，分离时为蓝色；图15中，侵入为正，分离为负。可以看出：载荷10 kN时，各啮合齿大部分区域处于接触侵入状态，仅编号4~11齿的部分区域以及齿16~21的小部分区域处于分离状态，编号18的啮合齿最右侧区域处于侵入位移最大，而编号为7的啮齿最右侧区域分离位移最大。随载荷增大，出现分离区域有所扩散，单个啮合齿分离区域也明显增加。载荷进一步增大，仅齿5~10和15~22的少部分区域处于接触状态，且由于涡轮轴扭转角变负，各齿在偏转中心线及其两侧的区域均处于分离状态，啮合齿接触区域在减小，但接触区域对应的侵入位移在不断增大。

4 关键参数的影响规律

4.1 摩擦因数

不同装配标准和工况条件下，套齿连接结构摩擦因数在装配/工作过程中变化较大^［20］，基于本文考虑接触面摩擦力的数值模型，摩擦因数对刚度的影响如图16所示。可以看出，增大摩擦因数会增大套齿连接结构的刚度，在较小位移下增幅明显，而当位移较大时，对其刚度的影响很小。相同位移下，相较于无摩擦（μ=0），摩擦因数增大刚度的百分比依次为11.16%（μ=0.15），19.24%（μ=0.3）和38.19 %（μ=0.45）。这是由于各齿空间位置不同，位移较小时，大部分齿处于接触状态并传递载荷，此时部分齿上的摩擦力方向与横向载荷方向夹角很小，能在一定程度上减小横向载荷产生的变形，进而提高刚度；而位移较大时，大部分齿处于分离状态，主要承受外载荷的齿的法线方向与横向载荷方向夹角很小，相应地，这些接触面的切线方向与载荷方向几乎垂直，对平衡横向载荷几乎没有影响。

4.2 啮合刚度

在套齿连接结构中，由于齿的结构较为细长，相较于加工齿的基座（轴段），套齿更容易发生变形，故套齿的柔性（或刚度）会对套齿连接结构的刚度产生影响。而根据文献［16］中计算啮合刚度的方法，不同的材料参数和套齿构型参数会对套齿的啮合刚度产生较大影响。啮合刚度k_n对套齿连接结构刚度的影响如图17所示。结果表明，啮合刚度提高会显著增大套齿连接结构的刚度，但增幅随啮合刚度的提高而降低。较小位移下，啮合刚度对套齿连接结构刚度的影响较大，而随位移增大，影响变小。这是由于啮合刚度对套齿连接结构刚度的影响取决于处于接触状态的接触面积，当位移增大，接触面积减小且处于接触状态齿的数量减少，套齿的啮合刚度的影响减小。

4.3 扭矩

在转子工作过程中，随涡轮输出功率变化，套齿连接结构所传递的扭矩也会随之变化，扭矩对套齿连接结构刚度特性的影响如图18所示。结果表明，相同位移下，扭矩能显著增大套齿连接结构的刚度。这是由于扭矩增大时，各对啮合齿之间接触更加紧密，改善了套齿连接结构的接触状态，处于接触状态的接触面积会随之增大，导致套齿连接结构刚度增大。

4.4 齿宽

在设计/加工套齿连接结构时，通常内外齿的齿宽长度不一样，以包容一定的安装误差和工作导致的变形。套齿齿宽对套齿连接结构刚度的影响如图19所示。结果表明：套齿齿宽增大时，套齿连接结构刚度有所增大，但随位移增大，刚度增大的程度逐渐减小。这是由于齿宽增大时，处于接触状态的接触面积增大，在较小位移时能显著提高套齿连接结构刚度。但各对啮合齿处于接触状态的接触面积和齿的数量随位移的增大而减小，此时不同齿宽下接触面积基本一致，使得刚度变化不大。

5 结论

1）本文通过在套齿连接结构轴段微元各啮合点处引入接触-摩擦模型，并结合整体结构的几何关系、受力平衡关系和数值迭代法，提出了一种柔性套齿连接结构刚度的模型，该模型能够十分简便地通过修改尺寸参数获得不同尺寸的套齿连接结构刚度特性，还能够考虑啮合齿对的接触摩擦特性、接触状态变化及其沿轴向变化特征。通过与ANSYS有限元模型刚度计算结果对比，证明了本文建模方法兼具高计算效率和高计算精度。

2）基于本文刚度求解方法分析了套齿连接结构刚度非线性特征，结合啮合接触变化规律揭示了刚度非线性产生机理，结果表明：随线位移的增大，套齿连接结构刚度逐渐下降，且下降速率逐渐减小。在外载荷或线位移增大过程中，套齿涡轮轴扭转角逐渐减小，甚至扭转角方向逐渐变为与扭矩方向相反。此过程中啮合接触区域也不断减小，这是套齿刚度降低的主要原因。当载荷或线位移增加到一定程度时，啮合接触区域近似呈现对称分布特征，且基本保持稳定，此时套齿连接结构刚度亦基本稳定。

3）研究了关键参数对套齿连接结构刚度的影响规律，结果表明：扭矩、啮合刚度、接触面摩擦因数及套齿齿宽增大，套齿连接结构刚度均有所提高。其中，摩擦因数和套齿齿宽主要对小位移阶段的刚度影响明显，啮合刚度增大到一定程度后与摩擦因数和齿宽的影响相似；而扭矩增大时，整个位移区间的套齿连接结构刚度均提高较为显著。

参考文献

[1] 王永亮， 赵广， 孙绪聪， 等. 航空花键研究综述［J］. 航空制造技术， 2017， 60（3）： 91-100.

[2] WANG Yongliang， ZHAO Guang， SUN Xucong， et al. Summary of Aviation Spline Research［J］. Aeronautical Manufacturing Technology， 2017， 60（3）： 91-100.

[3] 王庆国， 陈大兵， 魏静， 等. 基于有限元法的渐开线花键联接接触分析［J］. 机械传动， 2014， 38（1）： 134-137.

[4] WANG Qingguo， CHEN Dabing， WEI Jing， et al. Contact Analysis of Involute Spline Joint Based on FEM［J］. Journal of Mechanical Transmission， 2014， 38（1）： 134-137.

[5] BARROT A， PAREDES M， SARTOR M. Extended Equations of Load Distribution in the Axial Direction in a Spline Coupling［J］. Engineering Failure Analysis， 2009， 16（1）： 200-211.

[6] HONG J， TALBOT D， KAHRAMAN A. A Semi-analytical Load Distribution Model for Side-fit Involute Splines［J］. Mechanism and Machine Theory， 2014， 76： 39-55.

[7] HONG J， TALBOT D， KAHRAMAN A. Load Distribution Analysis of Clearance-fit Spline Joints Using Finite Elements［J］. Mechanism and Machine Theory， 2014， 74： 42-57.

[8] HONG J， TALBOT D， KAHRAMAN A. A Generalized Semi-analytical Load Distribution Model for Clearance-fit， Major-fit， Minor-fit， and Mismatched Splines［J］. Proceedings of the Institution of Mechanical Engineers， Part C： Journal of Mechanical Engineering Science， 2016， 230（7/8）： 1126-1138.

[9] CURÀ F， MURA A， ADAMO F. Fatigue Damage in Spline Couplings： Numerical Simulations and Experimental Validation［J］. Procedia Structural Integrity， 2017， 5： 1326-1333.

[10] 李伦绪， 陈果， 杨默晗. 航空发动机套齿连接结构刚度特性仿真分析及试验研究［J］. 中国机械工程， 2022， 33（18）： 2249-2257.

[11] LI Lunxu， CHEN Guo， YANG Mohan. Simulation Analysis and Experimental Study of Stiffness Characteristics of Aero-engine Spline Couplings［J］. China Mechanical Engineering， 2022， 33（18）： 2249-2257.

[12] 李英杰， 赵广， 袁运博， 等. 航空花键接触刚度仿真与试验［J］. 航空动力学报， 2024， 39（12）： 20230070.

[13] LI Yingjie， ZHAO Guang， YUAN Yunbo， et al. Simulation and Experiment on Contact Stiffness of Aviation Splines［J］. Journal of Aerospace Power， 2024， 39（12）： 20230070.

[14] 谭援强， 蒋理宽， 姜胜强， 等. 渐开线花键副微动摩擦接触分析［J］. 机械工程学报，2018，54（7）： 123-130.

[15] TAN Yuanqiang， JIANG Likuan， JIANG Shengqiang， et al. Analysis of Fretting Friction Contact of Involute Spline Pair［J］. Journal of Mechanical Engineering， 2018，54（7）： 123-130.

[16] MARMOL R A， SMALLEY A J， TECZA J A. Spline Coupling Induced Nonsynchronous Rotor Vibrations［J］. Journal of Mechanical Design， 1980， 102（1）： 168-176.

[17] 朱彬， 杨诚， 刘烨辉， 等. 套齿连接结构力学模型及其刚度影响因素［J］. 机械设计与制造， 2019（）： 86-90.

[18] ZHU Bin， YANG Cheng， LIU Yehui， et al. Mechanical Model of Sleeve Tooth Connection Structure and Its Stiffness Influencing Factors［J］. Machinery Design & Manufacture， 2019（S1）： 86-90.

[19] ZHAO Guang， SU Juncong， ZHAI Jingyu， et al. Study on Nonlinear Meshing Stiffness of Spline［C］∥Proceedings of the 9th IFToMM International Conference on Rotor Dynamics. Cham： Springer， 2015： 1315-1321.

[20] de VAUJANY J P， GUINGAND M， TEIXEIRA ALVES J. Quasi-static Analysis of Spline Coupling under Load［J］. Mechanics & Industry， 2017， 18（2）： 202.

[21] HONG J， TALBOT D， KAHRAMAN A. A Stiffness Formulation for Spline Joints［J］. Journal of Mechanical Design， 2016， 138（4）： 043301.

[22] YU Pingchao， WANG Cun， LIU Yunlong， et al. Analytical Modeling of the Lateral Stiffness of a Spline Coupling Considering Teeth Engagement and Influence on Rotor Dynamics［J］. European Journal of Mechanics—A/Solids， 2022， 92： 104468.

[23] ZHANG Chao， CAO Peng， ZHU Rupeng， et al. Dynamic Modeling and Analysis of the Spline Joint-flexible Coupling-rotor System with Misalignment［J］. Journal of Sound and Vibration， 2023， 554： 117696.

[24] GUO Yi， LAMBERT S， WALLEN R， et al. Theoretical and Experimental Study on Gear-coupling Contact and Loads Considering Misalignment， Torque， and Friction Influences［J］. Mechanism and Machine Theory， 2016， 98： 242-262.

[25] CURÀ F， MURA A， GRAVINA M. Load Distribution in Spline Coupling Teeth with Parallel Offset Misalignment［J］. Proceedings of the Institution of Mechanical Engineers， Part C： Journal of Mechanical Engineering Science， 2013， 227（10）： 2195-2205.

[26] 李英杰， 赵广， 吴学深， 等. 航空花键-转子系统自激振动研究综述［J］. 航空学报， 2022， 43（8）： 625532.

[27] LI Yingjie， ZHAO Guang， WU Xueshen， et al. Review of Research on Self-excited Vibration of Aviation Spline-rotor System［J］. Acta Aeronautica et Astronautica Sinica， 2022， 43（8）： 625532.

基金资助

国家自然科学基金(52372387)

民航飞机健康监测与智能维护重点实验室自主基金专项资金(NJ2024022)

中国博士后科学基金(2022M711615)

南京航空航天大学研究生科研与实践创新计划(xcxjh20230743)

南京航空航天大学研究生科研与实践创新计划(xcxjh20240721)