基于数字化齿面的弧齿锥齿轮柔性多体接触动力学分析

弧齿锥齿轮被广泛应用于汽车和航空航天领域。弧齿锥齿轮正朝着高速、高载荷质量比和高可靠性的方向发展，其动态啮合接触传动特性对整机性能有着重要的影响，也是传动误差规律、动态设计和可靠性分析的基础[1-2]。弧齿齿面的几何形状和加工方法复杂，难以获得弧齿齿面的解析解。目前，弧齿齿面的数字化建模方法应用较为广泛，一般是根据齿轮加工原理和齿面啮合理论，建立弧齿齿面方程，运用MATLAB计算弧齿齿面的数字化坐标点，再运用NURBS曲面拟合重构精确的弧齿齿面模型。Litvin等[2-4]提出弧齿锥齿轮的局部综合法，讨论了三齿有限元接触分析方法。Cheng等[5]根据齿轮啮合原理，运用精确几何模型分析准双曲面锥齿轮传动振动特性。Zhang等[6]采用三次样条拟合弧齿齿面，提出轮齿接触的数值分析方法。Sheveleva等[7]提出了螺旋锥齿轮啮合接触的接触痕、传动误差和接触压力分布的数值算法。汪中厚等[8]提出一种基于高精度数字化真实齿面的螺旋锥齿轮齿面接触分析方法。唐进元[9]，刘光磊[10]，郭伟超[11]运用MATLAB平台计算弧齿齿面的数字化数据点，分别建立了含有过渡曲面的弧齿锥齿轮三维模型和精确轮齿有限元模型。苏进[12]，刘光磊[13]，Lin[14]基于弧齿锥齿轮齿面三坐标测量网格点数据，重构离散数值弧齿齿面，分析了弧齿锥齿轮传动误差规律。张金良等[15]在ANSYS平台下建立了弧齿锥齿轮多齿模型，分析了多齿对接触下齿根接触应力的变化规律。唐进元等[16-17]运用ABAQUS和LS-DYNA平台，以动态加载接触分析有限元方法为基础，考虑惯性载荷和动态啮合关系，建立弧齿锥齿轮有限元模型，分析了啮合接触变化规律。汪中厚等[18-19]运用ABAQUS建立真实齿面的螺旋锥齿轮有限元模型，进行了传动误差分析。侯祥颖[20]运用数值化方法计算结点坐标，提出基于高精度建模的有限元分析方法。方崇德，邓效忠[21]，高建平等[22]，侯祥颖等[23]开展考虑边缘接触的弧齿锥齿轮承载接触分析。姚廷强等[24]采用Litvin的局部综合法计算得出弧齿齿面的数字化坐标点，运用三角网格方法，重构弧齿齿面模型，进行弧齿齿面动态接触多刚体动力学分析。

多数文献以静态或准静态力学分析方法为基础，建立单对或多对弧齿啮合传动的有限元模型，进行弧齿TCA轮齿接触分析和LTCA承载接触分析。这些方法未考虑弧齿锥齿轮的全齿模型，质量和惯量特性不完整，计算结果不能完全反映在弧齿齿面动态接触状态下弧齿锥齿轮的传动特性。考虑完整弧齿齿轮弹性变形和精确弧齿齿面三维动态接触特性，计算弧齿齿面啮合传动规律和边缘接触分析，成为弧齿锥齿轮系统柔性多体动力学亟待解决的研究问题。

1 弧齿锥齿轮齿面建模方法

1.1 弧齿齿面的空间离散点

Litvin提出弧齿锥齿轮的局部综合法已被文献详细论述，由于篇幅限制，此处仅给出主要参考公式。切削刀刃所形成的圆锥产形面上的任意点在刀盘坐标系下的径向位置矢量和法向矢量分别为

(1a)

(1b)

式中:rg=Ru±Pω2/2为常数刀齿半径;Ru为刀盘公称半径;Pω2为刀顶距;αg为刀具齿形角;sg、θg为刀具切削锥面的坐标参数。“±”表示对应加工弧齿齿面的凹凸面。

将刀盘坐标系下的式(1a)和式(1b)，通过一系列坐标变换至齿轮体坐标系下，可得弧齿齿轮的齿面方程rg和法向矢量ng。

(2a)

(2b)

式中:Aru为刀盘坐标系到摇台坐标系的变换矩阵;Amr为摇台坐标系到机床坐标系的变换矩阵;Aam为机床坐标系到辅助坐标系的变换矩阵;Aba为辅助坐标系到节锥坐标系的变换矩阵;Agb为节锥坐标系到齿轮体坐标系的变换矩阵。

结合齿面啮合方程可得

ng·vg=0

(3)

式中:vg为被加工齿轮与切削刃在切削点的相对速度。

由弧齿齿面方程式(2a)和啮合方程式(3)可求得以θg、ψg描述的齿面方程。

(4)

式中:xg、yg和zg为在齿轮体坐标系下弧齿齿面离散点的坐标分量。

运用MATLAB编程计算弧齿齿面的离散点的计算结果,如图1所示。

图1 弧齿齿面空间离散点

Fig.1 Coordinate points of gear tooth surface

1.2 弧齿齿面的法向矢量

以Matlab计算的数字化坐标点为弧齿齿面的表面结点，运用Hypermesh重构弧齿齿面，以六面体单元划分弧齿锥齿轮，建立真实柔性体有限元模型。将弧齿锥齿轮的单元和结点数据输出至ADAMS/MaxFlex，基于浮动参考坐标法和有限单元法的小变形柔性体建模方法[25]，运用浮动坐标系下的相对结点坐标来直接描述弧齿锥齿轮的弹性变形。选取弧齿齿面的单元表面和相邻4个结点，构成四边形面片单元，进而由四边形面片单元集合构成弧齿齿面的接触块，实现多柔性体、多弧齿面的动态接触搜索计算。图2为弧齿齿面的四边形面片网格模型。

图3为四边形面片单元的外法向矢量为在四边形面片单元上建立局部等参坐标系，采用一组等参坐标值ξ,η和双线性插值函数来描述四边形面片单元上任意一点的坐标矢量[26]，从而计算四边形面的单位外法向矢量和单位切向矢量。

图2 弧齿齿面的四边形面片单元模型

Fig.2 Quadrilateral patches of gear tooth

图3 四边形面片单元的外法向矢量

Fig.3 Normal vector of quadrilateral patch

(5)

式中:x″(ξ,η)表示局部坐标系下四边形面片单元上任意一点的参数坐标值;表示四边形面片单元上第k个结点的形函数。(ξk,ηk)为单元的局部坐标;x″k表示四边形面片单元上第k个结点坐标。

弧齿齿面的四边形面片单元的外法向矢量为

n″=x″(ξ,η)ξ×x″(ξ,η)η

(6)

式中:x″(ξ,η)ξ，x″(ξ,η)η为位置矢量x″(ξ,η)关于局部参数ξ,η的偏导数，其物理意义是四边形面片单元的切向矢量。

运用罚函数法计算弧齿齿面的法向接触力，同时考虑啮合齿面的阻尼作用和摩擦作用，建立弧齿齿面的迟滞法向接触力模型。由式(15)的接触弹性变形可得弧齿齿面间的法向接触力显示表达式为

(7a)

(7b)

(8)

式中:为单位外法向量;为正交单位切向量。

基于弧齿齿面的四边形面片接触块，编制自定义程序计算四边形面片单元中心的单位外法向矢量，提高弧齿锥齿轮柔性多体动态接触的计算效率。图4为弧齿齿面的四边形面片单元的中心外法向矢量示意图。

图4 四边形面片单元的中心外法向矢量

Fig.4 Center normal vectors of quadrilateral patches

2 弧齿齿面的动态接触力学模型

2.1 弧齿齿面的动态接触搜索方法

由参考系移动坐标、转动坐标和结点坐标来描述的柔性锥齿轮i的单元j中任一点的位置矢量的一般形式为[25]

(9)

式中:为未变形时由柔性锥齿轮的结点坐标矢量;为在浮动体坐标系下弧齿齿面上的结点坐标矢量;为惯性坐标系到体坐标系的方向变换矩阵，i=1,2。

图5为弧齿齿面的动态位置关系。OXYZ为绝对坐标系，和为主、从动柔性弧齿轮的浮动体坐标系，和为主动轮上齿面Sl和从动轮上齿面Ss的包围盒坐标系。

图5 弧齿齿面的动态接触关系

Fig.5 Dynamic interactions of gear teeth

在绝对坐标系下柔性锥齿轮的弧齿齿面的立方体包围坐标原点的位置矢量为

(10)

式中:为绝对坐标系下柔性主、从动锥齿轮的浮动体坐标系原点的位置矢量;和为绝对坐标系下包围盒坐标系原点的位置矢量;和为浮动体坐标系下包围盒坐标系原点的位置矢量，b=l,s。

在惯性坐标系下主从弧齿齿面上的各结点的位置矢量可由包围盒坐标系下的相对位置矢量来描述。

(11)

式中:为在绝对坐标系下未发生变形时弧齿齿面上的结点位置矢量;和为在包围盒坐标系下未发生变形时弧齿齿面上四边形面片单元的结点位置矢量。

在全域接触搜索过程中，基于弧齿齿面的立方体包围盒技术，运用式(10)对立方体刚体及其立方体单元进行搜索计算，确定发生接触的立方体。基于弧齿齿轮的小变形假设，运用式(9)和式(11)计算可能发生接触的主、从锥齿轮弧齿齿面上的少量从点和目标单元面。

在局部接触搜索过程中，弧齿啮合齿面的动态接触实质上是从结点-目标四边形面片单元的动态接触，由此构成结点-目标面的主从搜索方法。因此，主从搜索方法[26]仅对可能发生接触的从点和目标单元面进行搜索计算，以进一步确定发生接触的从点和目标单元面。由式(9)和图5可知，主动轮g1的弧齿齿面的四边形面片单元上的从结点与从动轮g2的弧齿齿面的四边形面片单元上的主结点的相对位置矢量为

(12)

式中:为在绝对坐标系下从结点与主结点的相对位置矢量。

在局部接触搜索过程中，主从搜索方法运用式(12)计算主、从结点之间距离，确定弧齿齿面上距离从结点最近的主结点，进而确定从动轮的弧齿齿面上包含该主结点的所有四边形面片单元。

图6为从结点到目标单元面的法向距离。从结点到目标单元面的法向距离等于从结点与目标单元面内的任意一点(此处选择最近的主结点)所构成的向量在外法向矢量方向上的投影长度。因此，在目标单元面的局部坐标系下从结点与主结点的最近的相对位置矢量为

(13)

式中:为目标单元面的局部坐标系的方向变换矩阵。

图6 从结点到目标单元面的法向距离

Fig.6 Normal approach from slave node to master surface

将式(13)沿着目标单元面的单位外法向矢量方向投影，计算从结点到包含最近主结点的所有目标单元面的最近法向距离，即得从结点-目标单元面的法向间隙量为

(14)

为便于问题讨论，忽略上下标符号，当从结点逐渐接近或远离目标单元面时，相对位置矢量d″与外法向矢量的方向相反，此时g<0，判断从结点与目标单元面不发生接触；当从结点穿透目标单元面时，相对位置矢量d″与外法向矢量的方向相同，此时g≥0，判断从结点与目标单元面发生接触，此时接触变形为

(15)

式中:δ为从结点与目标单元面的接触变形。

2.2 弧齿齿面的接触力计算模型

但他发现戒兄弟们之间并没有因为此事产生间隙，果智也压住了自己的愤怒。菩萨低眉，人性的卑微与弱小，纳了吧。

(16)

式中:Fgn为弧齿齿面的四边形面片单元的法向接触力，Kgc和Cgc分别为弧齿齿面弧齿中点节圆处的平均Hertz接触刚度和阻尼[24]。

图7 摩擦因数与相对速度关系

Fig.7 Relationship of friction coefficient and velocity

根据库伦摩擦定律，弧齿齿面之间的摩擦力与法向接触力间的关系为

Fgμ=μ(v)Fgn

(17)

图7为弧齿齿面的接触过程中摩擦因数与相对速度的关系，则摩擦因数为

(18)

式中:havsin(…)半正弦函数描述摩擦因数从静摩擦因数到动摩擦因数的平滑非线性变化，vgτ为弧齿齿面中发生Hertz接触作用时所对应的接触靶点处的相对切向速度，由从结点与目标单元面在接触对处切线方向的速度分量来计算。

图8 目标单元的法向接触力分配关系

Fig.8 Normal force distribution on target element

图8为目标单元的法向接触力分配关系。采用从结点在目标单元面的投影靶点p与四边形面片的4个结点构成的4个不同面积的三角形的线性插值方法，将目标单元面上的法向接触力分配给四个结点，获得结点的法向接触力分量[27]。

各结点处的法向接触力的分量力满足：

(19)

式中:fi为各个结点按照二级权重因子方法来线性分解得到的法向接触力为二级权重因子，由4个不同面积的三角形的线性插值方法确定的一级权重因子ωi描述，请参见文献[27]的详细推导过程。

2.3 弧齿锥齿轮柔性多体接触动力学模型

柔性锥齿轮系统的动力学方程为

(20)

式中:Mg为弧齿锥齿轮的质量矩为广义坐标。F为广义力矢量，分别为约束雅可比矩阵和拉格朗日乘子列向量为约束方程列向量

基于精确弧齿齿面模型和弧齿轮真实结构的柔性体有限元模型，基于ADAMS/MaxFlex多体动力学仿真平台，基于弧齿齿面的四边形面片接触块，实现弧齿齿面的多界面动态接触搜索算法和迟滞法向接触力计算，分析弧齿锥齿轮的多柔性体、多界面动态接触力学特性。

3 柔性弧齿锥齿轮动力学计算实例

3.1 弧齿锥齿轮传动特性分析

弧齿啮合齿面的椭圆接触面和接触刚度满足Hertz点接触条件，以弧齿齿面中点处的平均Hertz接触刚度参数，来近似描述单对弧齿啮合齿面的接触刚度参数。锥齿轮的材料弹性模量为2.1×1011 N/mm2，泊松比为0.3，弧齿锥齿轮在啮合齿面中点处的接触刚度参数为Kgc=7.2×105 N/mm1.5。从动轮受扭矩载荷，弧齿轮啮合传动周期T=0.02 s，啮合频率f=1 000 Hz，单齿啮合时间Tz=0.001，仿真计算时间t=0.05 s，主动轮的启动加速驱动函数为(time,0,0,0.01,3 000)r/min。

表1 弧齿锥齿轮设计参数

Tab.1 Parameters of bevel gear

图9和图10分别为刚性弧齿轮和柔性弧齿轮的啮合传动过程中，从动齿轮的角速度和转角误差的变化规律。分析可知，刚性弧齿锥齿轮传动的转速波动幅度比柔性弧齿锥齿轮的要大。考虑弧齿锥齿轮的结构弹性变形时，从动齿轮转速以3 000 r/min为中心，以12 r/min幅值波动，幅度相对较小，转角动态误差较小，且成周期性变化，传动更加平稳。图中放大区域为0.02～0.03 s半个啮合周期内(T/2=0.01 s)，由单齿啮合时间Tz=0.001知，有10对相邻弧齿啮入和啮出，受啮合频率的影响，角速度出现10个波峰波谷的周期变化规律。计算结果说明弹性变形对弧齿锥齿轮的传动特性有重要的影响，较为真实地表现了弧齿锥齿轮的传动特性。

图9 从动齿轮的角速度

Fig.9 Angular speed of driven gear

图10 从动齿轮的转角动态误差

Fig.10 Rotation angular errors of driven gear

图11分别为刚、柔性弧齿轮啮合传动过程中，相邻弧齿啮合齿面的动态接触力的变化规律。弧齿锥齿轮传动过程中，同时有相邻3对轮齿参与承载啮合传动，当有弧齿逐渐啮入或啮出时，相邻弧齿的动态接触力会表现出波峰或波谷的变化规律。由于柔性弧齿轮传动的转速波动较小，单齿啮合承载时间短于刚性轮齿，柔性弧齿啮合齿面存在明显的弧齿齿顶边缘接触(图13)，使得柔性弧齿齿面的啮合接触力比刚性轮齿的要大。图12为从动弧齿锥齿轮的约束反力。柔性弧齿锥齿轮啮合传动过程中，由于柔性弧齿存在明显的边缘接触，约束位置处的合力和各方向分量力均呈周期性的变化规律，合力(FM)、径向分量力(FX和FY)和轴向分量力(FZ)的波动幅值大于刚性齿轮的。柔性齿轮的径向分量力FX的数值是以刚性齿轮的径向分量力FX为最小值，以2 180 N成倍的幅度周期波动，而柔性齿轮的径向分量力FY的数值是以刚性齿轮的径向分量力FY为最大值，以1 850 N成倍的幅度周期波动。结果表明考虑弹性变形和弧齿齿面动态接触关系的完整弧齿锥齿轮更真实地模拟齿轮传动特性。弧齿锥齿轮的刚体假设方法忽略弧齿锥齿轮的弹性变形，会低估弧齿齿面的动态接触力和齿轮系统的约束反力，从而影响弧齿锥齿轮和滚动轴承的可靠性分析和疲劳寿命的计算结果。

图11 相邻弧齿啮合齿面的动态接触力

Fig.11 Dynamic contact forces of adjacent teeth

图12 从动齿轮的约束反力

Fig.12 Constraint forces of driven gear

图13 弧齿齿面的啮合接触力和综合应力

Fig.13 Contact force and Von Mises stress of gear tooth surface

图13为弧齿齿面啮合接触力、齿顶边缘各结点的综合应力和齿面中部各结点的综合应力。分析可知，主动轮的凹面和从动轮的凸面为弧齿啮合齿面，齿面从大端啮入，小端啮出。对比分析图11和13可知，考虑圆弧齿的结构弹性变形影响，弧齿啮合传动过程中，柔性弧齿锥齿轮存在明显的齿顶边缘接触和弧齿小端的边缘接触，弧齿齿顶边缘接触对齿面综合应力和齿面啮合接触力有着重要的影响。当齿顶边缘或弧齿小端发生接触时，由于承载区域较小，弧齿齿面存在较大的综合应力，尤其是齿顶中部发生边缘接触时的综合应力数值较大。当齿顶中部边缘和齿面中部同时发生接触时，弧齿齿面啮合接触力达到最大值，随着齿顶边缘接触逐渐退出啮合，接触力逐渐减小，此时弧齿齿面啮合接触力与刚性齿轮的齿面啮合接触力大小接近。

3.2 弧齿齿面的动态接触特性

图14为弧齿齿面的综合应力云图。图16为弧齿锥齿轮啮合齿面的最大综合应力的变化规律。在启动加速阶段，弧齿啮合齿面出现较大的应力值，达到稳定后，呈现每个单齿啮合时产生的最大应力的周期变化规律。对比图14的应力云图结果分析，图15中A类极大值为弧齿齿顶啮入时齿顶中部边缘呈现的周期性波峰应力值，D类极小值为弧齿齿顶啮出时齿顶中部边缘呈现的周期性波谷应力值，B类极大值为弧齿小端边缘的波峰应力值，C类极大值为弧齿小端齿顶角边缘的波峰应力值。图16为弧齿齿面齿高方向系列结点的动态应力规律(节点序号从弧齿大端向小端，沿齿宽方向编号)。第1条曲线为弧齿齿根位置处的动态应力。当弧齿面啮入时，齿根位置处的大端位置具有较大的动态应力，沿着齿宽方向动态应力逐渐减小，到齿宽中部后齿根应力变化较为稳定平缓。第4条曲线为弧齿齿高方向上中部位置的动态应力。弧齿大端位置具有较大的动态应力，沿着齿宽方向动态应力逐渐减小至一定数值而保持相对稳定。结果表明弧齿锥齿轮的负载扭矩主要由弧齿啮合齿面的中部承载，第4条曲线的动态应力相对稳定，有利于弧齿锥齿轮的传动平稳和疲劳寿命。第6条曲线为弧齿齿顶边缘位置的动态应力。当弧齿啮入或啮出时，弧齿齿顶发生边缘接触的大端和小端都存在较大的动态应力。随着弧齿齿面啮入，负载扭矩主要由弧齿啮合齿面的中部承载，此时弧齿齿顶边缘接触的中部也有较大的动态应力。

图14 弧齿齿面的综合应力云图

Fig.14 Von Mises stress of gear teeth surfaces

图15 弧齿齿面的最大综合应力

Fig.15 Max Von Mises stress of gear teeth surfaces

(a) 齿高方向的第1曲线

(b) 齿高方向的第4曲线

(c) 齿高方向的第6曲线

图16 弧齿齿面齿宽方向结点的动态应力规律

Fig.16 Dynamic stress of gear teeth along tooth width

由于工况载荷和加工装配误差等外部因素，弹性变形和动态接触等内部因素，弧齿锥齿轮在实际工况下通常会出现边缘接触现象。图17为相邻3对弧齿单齿啮合时间内(Tz=0.001 s)齿面的最大应力规律。在弧齿锥齿轮传动过程中，由于弧齿锥齿轮的弹性变形和动态接触特性，相邻弧齿在单齿啮合时间内齿面不同位置发生动态接触，存在不同的最大应力规律。给定参数对应的弧齿锥齿轮实例中，弧齿锥齿轮主要发生弧齿齿面中部和齿顶中部的承载接触现象，而且在弧齿单齿啮合时间内齿顶中部的边缘接触承载时间相对较长。齿顶中部和齿面小端中部有明显的边缘接触现象，发生边缘接触时存在相对平稳变化的较大应力值。相邻弧齿齿面中部承载接触时最大应力的变化规律有所不同，当最大应力未从齿顶边缘过度至齿面中部时存在一定的波动(图17(a))。图为18相邻10对弧齿齿顶边缘接触时的最大应力。不同弧齿齿顶的边缘接触规律有所不同，但总体呈现出相邻5对弧齿齿顶的周期性边缘接触规律。结果表明边缘接触改变了弧齿啮合齿面的载荷分布规律，反映出弧齿齿轮的承载啮合性能，更为真实地表现了弧齿齿轮动态接触传动的全过程。

(a) 相邻弧齿1的单齿啮合时间内齿面最大应力

(b) 相邻弧齿2的单齿啮合时间内齿面最大应力

(c) 相邻弧齿3的单齿啮合时间内齿面最大应力

图17 相邻3对弧齿齿面的最大应力规律

Fig.17 Max Von Mises stress of adjacent three teeth

图18 相邻10对弧齿齿顶边缘接触时的最大应力

Fig.18 Max Von Mises stress of adjacent ten teeth addendum edge

4 结 论

本文提出弧齿锥齿轮的多柔性体接触动力学分析方法，分析了完整弧齿锥齿轮的弹性变形和精确弧齿齿面的三维动态接触特性，对弧齿锥齿轮的动态啮合传动性能的影响规律，真实地模拟了弧齿锥齿轮的多柔性体、多界面动态接触动力学特性，对高性能弧齿锥齿轮的动力学设计具有参考意义。

(1) 从动齿轮角速度是以啮合频率为基础的小幅波动的周期变化规律，考虑弹性变形的柔性弧齿锥齿轮的角速度和转角动态误差较小，较为稳定。弹性变形对弧齿锥齿轮承载传动性能有很大的影响，刚性弧齿锥齿轮方法会低估弧齿齿面的动态接触力和齿轮系统的约束反力。

(2) 在启动加速阶段，弧齿锥齿轮的啮合弧齿出现较大的应力值，稳定后呈现每个单齿啮合时的最大应力的周期变化规律。弧齿大端的齿根应力值大于小端的，齿面中部承载应力分布较为稳定，齿顶大端、小端在啮入啮出时有较大的应力变化，齿顶中部承载应力分布较为平缓变化。

(3) 弧齿锥齿轮主要发生弧齿齿面中部和齿顶中部的承载接触区域，齿顶中部发生边缘接触现象，承载时间相对较长，较大应力的平稳变化规律，总体呈现出相邻5对弧齿齿顶的周期性边缘接触规律。