CN117634057A - 一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 - Google Patents
一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 Download PDFInfo
- Publication number
- CN117634057A CN117634057A CN202311350092.1A CN202311350092A CN117634057A CN 117634057 A CN117634057 A CN 117634057A CN 202311350092 A CN202311350092 A CN 202311350092A CN 117634057 A CN117634057 A CN 117634057A
- Authority
- CN
- China
- Prior art keywords
- spiral bevel
- bevel gear
- gear
- contact
- point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 16
- 238000000034 method Methods 0.000 claims abstract description 35
- 238000004299 exfoliation Methods 0.000 claims abstract description 16
- 239000011159 matrix material Substances 0.000 claims description 28
- 239000013598 vector Substances 0.000 claims description 25
- 238000012546 transfer Methods 0.000 claims description 9
- 230000005540 biological transmission Effects 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 7
- 238000006073 displacement reaction Methods 0.000 claims description 6
- 238000005520 cutting process Methods 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000003754 machining Methods 0.000 description 3
- 238000005096 rolling process Methods 0.000 description 3
- 244000309464 bull Species 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000005251 gamma ray Effects 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 101100182136 Neurospora crassa (strain ATCC 24698 / 74-OR23-1A / CBS 708.71 / DSM 1257 / FGSC 987) loc-1 gene Proteins 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 230000009347 mechanical transmission Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000004901 spalling Methods 0.000 description 1
- 230000036346 tooth eruption Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Gears, Cams (AREA)
Abstract
本发明公开了一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法，包括采集弧齿锥齿轮的属性参数和剥落故障参数、轮齿接触分析、承载轮齿接触分析的步骤，并输出时变啮合刚度的计算结果，本发明所提供的方法仅需提供弧齿锥齿轮的属性参数和剥落故障参数即可运行，不受限于操作人员的技术基础，且耗时短，在保证准确性的前提下，大幅提高了计算效率。
Description
技术领域
本发明涉及齿轮传动技术领域，尤其是一种齿轮时变啮合刚度计算方法。
背景技术
弧齿锥齿轮被广泛应用于重要的机械传动系统(如直升机、航空发动机、汽车等)，具有承载能力大、效率高、功率密度大、振动和噪声小等优点。齿轮的结构参数复杂，齿面啮合力分布在三维空间中。因此，对齿轮副的啮合特性进行分析是困难的，特别是对齿轮副的故障(如剥落、裂纹等)。齿轮副时变啮合刚度(TVMS)是引起齿轮传动系统非线性振动的主要来源。TVMS也是弧齿锥齿轮传动系统动力学模型的重要组成部分，是弧齿锥齿轮传动系统动态特性分析的基础。随着直升机功率、转速等的增大，接触力增大，导致齿面失效概率增大。因此，计算弧齿锥齿轮的时变啮合刚度非常重要。目前，时变啮合刚度的计算方法主要分为四种，(a)专业或商业软件分析，虽然比较方便，但无法对故障弧齿锥齿轮的啮合特性进行分析；(b)通用商业有限元软件，如ANSYS、ABAQUS等，由于接触非线性，有限元法计算效率低，而且需要操作人员有一定的技术基础。(c)半解析法，相比于有限元法，提高了效率，但耗时仍是非常久；(d)解析法，虽然效率很高，但是大多需要通过有限元模型进行修正以提高精度。
故，需要一种新的技术方案以解决上述技术问题。
发明内容
本发明的目的针对现有技术中的不足，提出了一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法，提供一种含剥落故障的弧齿锥齿轮时变啮合刚度的计算方法，在保证准确度的同时，提高计算效率。
为达到上述目的，本发明可采用如下技术方案：
一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法，包括以下步骤：
(1)采集弧齿锥齿轮的属性参数和剥落故障参数；
(2)同时联立求解啮合方程和旋转投影面方程得到弧齿锥齿轮大齿轮齿面上任意一点的坐标以及弧齿锥齿轮小齿轮齿面上任意一点的坐标；
(3)将弧齿锥齿轮小齿轮齿面方程和法向量从轮坯坐标系转化到装配坐标系；将弧齿锥齿轮大齿轮齿面方程和法向量从轮坯坐标系转化到装配坐标系；
(4)根据弧齿锥齿轮大、小齿轮连续相切方程组计算出齿轮副在装配坐标系下第一次啮合时各自所需要的旋转的角度φg0、φp0；其中φg0为弧齿锥齿轮大齿轮第一次啮合时所需要的旋转的角度；φp0为弧齿锥齿轮小齿轮第一次啮合时所需要的旋转的角度；
(5)以所述的φp0为中间值，上下波动一定量φrange假定小齿轮转角的区间[φp0-φrange,φp0+φrange]，在该区间内等间距均分得到不同啮合时刻下弧齿锥齿轮小齿轮的转角φp；通过转角φp解得初始接触点的未知变量值，代入弧齿锥齿轮大、小齿轮齿面方程得到初始接触点的坐标，每个初始接触点代表一个啮合时刻和齿轮啮合位置；
(6)弧齿锥齿轮大、小齿轮齿面上任意一点的坐标为离散点坐标，根据该离散点坐标通过线性插值生成齿轮内部网格；采用三维双线性等参单元建立有限元模型，在齿基底部添加除绕轴旋转外的约束；在不同啮合时刻下大小齿轮接触面齿面点分别施加接触点法向的单位载荷，提取在加载点齿面以下位置的变形，组集成弧齿锥齿轮大齿轮3维整体柔度矩阵及弧齿锥齿轮小齿轮3维整体柔度矩阵/>和，n表示啮合时刻的个数，i表示齿面施力点的个数，j表示齿面提取位移点的个数；
(7)当受载时，弧齿锥齿轮的啮合状态由接触点变为接触椭圆，将接触椭圆表示为沿接触椭圆长轴的接触曲线；在接触点附近建立极坐标系(Lp，θ)，将初始接触点法向量沿着切平面移动Lp，分别与大小齿轮相交于两点，两点之间距离即为齿面间隙dp，计算不同θ下的齿面间隙，寻找最小间隙的θmin值，此方向即为接触椭圆的长轴方向，最大间隙的θmax值为接触椭圆的短轴方向；沿着长轴方向，计算不同距离Lp下的齿面间隙εn×1；
(8)根据将网格点以样条插值的方法加密，然后用距离潜在接触点最近的插值点或网格点的柔度作为潜在接触点的柔度，得到弧齿锥齿轮小齿轮潜在接触点柔度矩阵λp和弧齿锥齿轮大齿轮潜在接触点柔度矩阵λg，然后整合为潜在接触点柔度矩阵λb；
(9)根据步骤(1)中采集的剥落顶点，在投影平面内寻找剥落区域穿过的潜在接触点，其为剥落影响的范围点，在剥落影响范围内的潜在接触点间隙增加剥落深度值，得到新的潜在接触点间隙εn×1；
(10)计算时变啮合刚度Kmesh：
式中，εi表示第i个潜在接触点的间隙，Fn_all表示接触力的合力；Ste表示承载传递误差。
进一步的，步骤(1)中，属性参数包括齿轮参数、材料参数和部分加工参数；剥落故障参数包括剥落多边形的各个顶点坐标，剥落深度。
进一步的，步骤(2)中，轮坯坐标系下的弧齿锥齿轮大齿轮齿面方程
表示大齿轮轮坯的旋转角度，sg表示直线刀刃长度，即切削深度，ɑg表示刀具齿形角，θg表示刀尖旋转的角度，M2g表示Sg到S2的坐标转化矩阵；Sg为刀具坐标系；S2为轮坯坐标系；rg(sg,θg)为刀具旋转所得到的曲面矢量；
同时联立求解下述啮合方程和旋转投影面方程得到弧齿锥齿轮大齿轮齿面上任意一点的坐标；
式中，Xm为旋转投影面x轴坐标，Ym为旋转投影面y轴坐标；
同理，轮坯坐标系下的弧齿锥齿轮小齿轮齿面齿面方程
表示小齿轮轮坯的旋转角度，sp表示直线刀刃长度，即切削深度，θp表示刀尖旋转的角度，M1p表示Sp到S1的坐标转化矩阵，Sp为刀具坐标系，S1为轮坯坐标系，rp(sp,θp)为刀具旋转所得到的曲面矢量；
同时联立求解下述啮合方程和旋转投影面方程得到弧齿锥齿轮小齿轮齿面上任意一点的坐标；
进一步的，步骤(4)中，根据大小齿轮连续相切方程组计算出齿轮副在装配坐标系下第一次啮合时各自所需要的旋转的角度φg0、φp0；
和/>分别为大小齿轮装配坐标系下齿面方程，/>和/>分别为大小齿轮装配坐标系下齿面点的法向量方程。
进一步的，步骤(5)中，以所述的φp0为中间值，上下波动一定量φrange假定小齿轮转角的区间[φp0-φrange,φp0+φrange]，在该区间内等间距均分得到不同啮合时刻下小齿轮的转角φp，带入权利要求4中方程组解得初始接触点的未知变量值，代入齿面方程r1，r2得到初始接触点的坐标，每个初始接触点代表一个啮合时刻和齿轮啮合位置；并计算出空载传递误差Δφg；
式中，zg和zp分别为大小齿轮的齿数。
进一步的，步骤(8)中，
和/>分别为大小齿轮潜在接触点柔度矩阵对应元素值；i表示提取位移的接触点，j表示施加力的接触点，n为此接触时刻接触长轴的潜在接触点个数。
进一步的，步骤(10)中，首先计算接触柔度λc，
E为弹性模量；L是潜在接触点之间的距离；Fi是第i个潜在接触点上的法向接触力；
由进行变形协调迭代，
Fn×1 表示接触点法向方向上的分布力；ε n×1表示潜在接触点的间隙；Ste表示承载传递误差；Fn all表示接触力的合力。
进一步的， 步骤(10)中，δ1表示弧齿锥齿轮的分锥角；ɑn表示弧齿锥齿轮的法向压力角；β表示弧齿锥齿轮的平均螺旋角。
本发明还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述方法的步骤。
本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。
有益效果：现有LTCA技术，大多数在计算长轴时，常常参考Litvin方法，假设齿轮间隙ε为0.00635mm，此为经验数值。而本发明在计算长轴时，通过接触椭圆的空间属性， 以寻找接触点附近间隙最小的方向θmin和间隙最大的方向θmax确定长短轴，并且求取齿面内所有潜在接触点，通过变形协调迭代出实际啮合点。
同时，现有技术中基于ANSYS计算时变啮合刚度时，需要使用布尔运算来生成复杂的多边形槽，需要操作人员具有一定的技术基础，需要花费9小时左右；而本发明所提供的方法仅需提供弧齿锥齿轮的属性参数和剥落故障参数即可运行，不受限于操作人员的技术基础，且耗时短，在保证准确度的前提下，大幅提高了计算效率。
附图说明
图1为本发明含剥落故障的弧齿锥齿轮时变啮合刚度计算方法的流程图；
图2为本发明中TCA流程图；
图3为本发明中齿面加工过程和系统坐标系示意图；
图4为本发明中齿轮旋转投影面示意图；
图5为本发明中啮合状态和潜在接触点示意图：其中(a)为共轭齿面的啮合状态；(b)为极坐标系(LP，θ)下共轭齿面之间对应点的间隙；(c)为潜在接触点；
图6为本发明中剥落示意图：其中(a)为剥落区域；(b)为剥落影响范围；(c)为剥落影响下的间隙；(d)为剥落故障后，齿面的分布接触力；
图7为本发明中啮合平面及承载传递误差示意图；
图8为本发明中接触轨迹和空载传递误差示意图：其中(a)为接触轨迹；(b)为空载传递误差；
图9为本发明中潜在接触点间隙示意图；
图10为本发明中整体柔度矩阵插值示意图；
图11为有限元模型和本发明方法在不同扭矩下健康齿轮的时变啮合刚度对比示意图；
图12为本发明中不同剥落长度的剥落区域和间隙示意图：其中(a)为不同剥落宽度的剥落区域；(b)为不同剥落长度的间隙图；
图13为不同剥落长度ANSYS和本发明方法时变刚度对比图：其中(a)为ANSYS的时变刚度；(b)本发明的时变刚度；
图14为本发明中不同剥落宽度的剥落区域示意图；
图15为不同剥落宽度ANSYS和本发明方法时变刚度对比图：其中(a)为ANSYS的时变刚度；(b)本发明的时变刚度；
图16为本发明中不同剥落位置的剥落区域示意图；
图17为不同剥落位置ANSYS和本发明方法时变刚度对比图：其中(a)为ANSYS的时变刚度；(b)为本发明的时变刚度。
具体实施方式
下面结合附图和具体实施例对本发明作进一步详细说明。
结合图1所示，本发明提供的含剥落故障的弧齿锥齿轮时变啮合刚度计算方法，应用于相互啮合的弧齿锥齿轮大、小齿轮的时变啮合刚度计算。该方法的具体步骤如下。
一、采集参数
采集弧齿锥齿轮的属性参数和剥落故障参数。属性参数包括齿轮参数、材料参数和部分加工参数；剥落故障参数包括剥落多边形的各个顶点坐标，剥落深度。
二、轮齿接触分析(TCA)
轮齿接触分析(TCA)流程图如图2所示，具体包括:
(1)大齿轮齿面方程
根据采集的齿轮参数和加工参数，刀具旋转所得到的曲面Σg由矢量rg(sg,θg)表示为：
式中，Rg表示刀盘半径，sg表示直线刀刃长度，即切削深度，ɑg表示刀具齿形角，θg表示刀尖旋转的角度，±符号分别对应大齿轮刀具旋转所得到曲面的凹面和凸面。
曲面Σg的法向量ng(θg)表示为：
式中，±符号分别对应大齿轮刀具旋转所得到曲面的凹面法向量和凸面法向量。
通过坐标转换(公式(3))即可得到轮坯坐标系下齿面包络线方程(即齿面方程)
式中，表示大齿轮轮坯的旋转角度，M2g表示Sg到S2的坐标转化矩阵，由公式(4)～(9)所得，
M2g＝M2b2Mb2a2Ma2m2Mm2c2Mc2g (9)
式中，Mij角标代表从Sj(xj，yj，zj)坐标系转化到Si(xi，yi，zi)坐标系的坐标转换矩阵；Sg(xg，yg，zg)为刀具坐标系；Sm2(xm2，ym2，zm2)为机床坐标系，固定在机床上；Sa2(xa2，ya2，za2)，Sc2(xc2，yc2，zc2)，Sb2(xb2，yb2，zb2)为辅助坐标系；S2(x2，y2，z2)为轮坯坐标系，绕zb2轴旋转，旋转角度为其中Sr2为径向刀位，刀盘中心到摇台中心的距离；q2为角向刀位，表示刀盘中心与水平面夹角；/>表示机床摇台机构的转角，其与大齿轮旋转角度/>的关系如式(10)所示；XE2为垂直轮位，表示被切齿轮的中心线相对于摇台水平轴线的垂直偏置量；XB2为床位，表示工件箱沿摇台中心线方向自目标位置向前进或后退的距离，控制切齿时的深度；XD2为水平轮位，表示机床摇台中心到工件箱主轴端面的距离；γm2表示轮坯安装角。加工过程和系统坐标系请结合图3所示。
式中m2c2表示大齿轮的滚动比。
同时联立求解啮合方程(11)和旋转投影面方程(12)即可得到弧齿锥齿轮大齿轮齿面上任意一点的坐标。
式中，结合图4所示，Xm为旋转投影面x轴坐标，Ym为旋转投影面y轴坐标。
(2)小齿轮齿面方程
根据4.2.1中的齿轮参数和加工参数，刀具旋转所得到的曲面Σp由矢量rp(sp,θp)表示为：
式中，Rp表示刀盘半径，sp表示直线刀刃长度，即切削深度，ɑp表示刀具齿形角，θp表示刀尖旋转的角度，±符号分别对应小齿轮刀具旋转所得到曲面的凹面和凸面。
曲面Σp的法向量np(θp)表示为：
式中，±符号分别对应小齿轮刀具旋转所得到曲面的凹面法向量和凸面法向量。
通过坐标转换(公式)即可得到轮坯坐标系下齿面包络线方程(即齿面方程)
式中，表示小齿轮轮坯的旋转角度，M1p表示Sp到S1的坐标转化矩阵，由公式所得，
M1p＝M1b1Mb1a1Ma1m1Mm1c1Mc1p (21)
式中，Mij角标代表从Sj(xj，yj，zj)坐标系转化到Si(xi，yi，zi)坐标系的坐标转换矩阵；Sp(xp，yp，zp)为刀具坐标系；Sm1(xm1，ym1，zm1)为机床坐标系，固定在机床上；Sa1(xa1，ya1，za1)，Sc1(xc1，yc1，zc1)，Sb1(xb1，yb1，zb1)为辅助坐标系；S1(x1，y1，z1)为轮坯坐标系，绕zb1轴旋转，旋转角度为其中Sr1为径向刀位，刀盘中心到摇台中心的距离；q1为角向刀位，表示刀盘中心与水平面夹角；/>表示机床摇台机构的转角，其与小齿轮旋转角度/>的关系如式(22)所示；XE1为垂直轮位，表示被切齿轮的中心线相对于摇台水平轴线的垂直偏置量；XB1为床位，表示工件箱沿摇台中心线方向自目标位置向前进或后退的距离，控制切齿时的深度；XD1为水平轮位，表示机床摇台中心到工件箱主轴端面的距离；γm1表示轮坯安装角。
式中m1c表示小齿轮的滚动比，C、D分别为切削滚比系数。
联立求解啮合方程(23)和投影面方程(24)即可得到弧齿锥齿轮小齿轮齿面上任意一点的坐标。
(3)通过公式(25)(26)(29)，可以将小齿轮齿面方程和法向量从轮坯坐标系转化到装配坐标系。通过公式(27)(28)(30)，可以将大齿轮齿面方程和法向量从轮坯坐标系转化到装配坐标系。
式中，ΔAG和ΔAP分别为大小齿轮轴向安装误差；Σ为大小齿轮轴交角；ΔΣ为大小齿轮轴交角误差；ΔE为大小齿轮垂直偏置距，即轴交平面垂直方向上的安装误差。和分别为大小齿轮装配坐标系下齿面方程，/>和/>分别为大小齿轮装配坐标系下齿面点的法向量方程。Lij为对应Mij矩阵去掉最后一行和最后一列之后组成的3×3矩阵。
(4)大小齿轮连续相切方程组(31)中有8个未知参数，分别为sg，θg， sp，θp，φp。但因为/>和/>都是单位向量，所以方程组只能列写出7个独立方程。假定旋转投影面的中心点为接触点，计算出大小齿轮对应的sg，θg，/>sp，θp，/>便可根据大小齿轮连续相切方程组(公式(31))计算出齿轮副在装配坐标系下第一次啮合时各自所需要的旋转的角度φg0、φp0。
(5)以所述的φp0为中间值，上下波动一定量φrange假定小齿轮转角的区间[φp0-φrange,φp0+φrange]，在该区间内等间距均分得到不同啮合时刻下小齿轮的转角φp，带入方程组(31)可解得初始接触点的未知变量值，代入齿面方程r1，r2可得到初始接触点的坐标，每个初始接触点代表一个啮合时刻和齿轮啮合位置，通过式(32)计算出空载传递误差Δφg。
/>
式中，zg和zp分别为大小齿轮的齿数。
三、承载轮齿接触分析(LTCA)
(1)TCA可以得到大小齿轮齿面离散点坐标，通过线性插值生成齿轮内部网格。采用三维双线性等参单元建立有限元模型，在齿基底部添加除绕轴旋转外的约束；在不同啮合时刻下大小齿轮接触面齿面点分别施加接触点法向的单位载荷，提取在加载点齿面0.2m(m为齿轮的模数)以下位置的变形，组集成3维整体柔度矩阵和/>n表示啮合时刻的个数，i表示齿面施力点的个数，j表示齿面提取位移点的个数。
(2)当受载时，弧齿锥齿轮的啮合状态由接触点变为接触椭圆(如图5(a)所示)，本发明将接触椭圆表示为沿接触椭圆长轴的接触曲线。在接触点附近建立极坐标系(LP，θ)(如图5(b)所示)，将初始接触点法向量沿着切平面移动Lp，分别与大小齿轮相交于两点，两点之间距离即为齿面间隙dp，计算不同θ下的齿面间隙，寻找最小间隙的θmin值，此方向即为接触椭圆的长轴方向，而最大间隙的θmax值即为接触椭圆的短轴方向。沿着长轴方向，计算不同距离L(旋转投影面下要能超出齿面边界)下的齿面间隙ε(即为潜在接触点的间隙)，将旋转投影面外的点删去。考虑边缘接触的影响，还计算了超出齿轮面边界的初始接触点的潜在接触点(如图5(c)所示)。
(3)根据得到的大小齿轮的整体柔度矩阵，将网格点以样条插值的方法加密，然后用距离潜在接触点最近的插值点或网格点的柔度作为潜在接触点的柔度，得到潜在接触点的柔度矩阵λp和λg，然后根据式(33)整合为潜在接触点柔度矩阵λb。
式中，和/>分别为大小齿轮潜在接触点柔度矩阵对应元素值；i表示提取位移的接触点，j表示施加力的接触点。n为此接触时刻接触长轴的长轴点数。
(4)根据采集的剥落顶点(如图6(a))，在投影平面内寻找剥落区域穿过的潜在接触点，其为剥落影响的范围点(如图6(b))，在剥落影响范围内的潜在接触点间隙增加剥落深度值(根据采集的剥落深度)，得到新的潜在接触点间隙εn×1(如图6(c))。
(5)计算时变啮合刚度：
由公式(34)计算接触柔度λc。
式中，E为弹性模量；L是在接触点之间的距离；Fi是第i个潜在接触点上的法向接触力。
根据公式(35)进行变形协调迭代。
式中，Fn×1表示接触点法向方向上的分布力；εn×1表示潜在接触点的间隙；Ste表示承载传递误差；Fn_all表示接触力的合力，由式(36)计算而来。
式中，δ1表示弧齿锥齿轮的分锥角；ɑn表示弧齿锥齿轮的法向压力角；β表示弧齿锥齿轮的平均螺旋角(如图7所示)。
时变啮合刚度可以表示为：
式中，εi表示第i个潜在接触点的间隙。
基于上述技术方案，以下通过一个具体的弧齿锥齿轮实例作为验证案例。
表1实例弧齿锥齿轮几何参数
表2实例弧齿锥齿轮加工参数
(1)根据表1和表2，可以大齿轮得到刀具旋转所得到的曲面rg(sg,θg)和ng(θg)表示为：
大齿轮轮坯坐标系下齿面方程为：
r2(sg,θg,ψ2)＝M2g(ψ2)rg(sg,θg) (40)
大齿轮啮合方程为：
旋转投影面坐标转化公式如下：
联立式(41)和式(42)即可得到大齿轮齿面上任意一点的坐标。
根据表1和表2，可以得到小齿轮刀具旋转所得到的曲面rp(sp,θp)和np(θp)表示为：
小齿轮轮坯坐标系下齿面方程为：
小齿轮啮合方程为：
旋转投影面坐标转化公式如下：
联立式(41)和式(42)即可得到小齿轮齿面上任意一点的坐标。
(2)齿面坐标从轮坯坐标系转到装配坐标系。根据大小齿轮相切方程(如公式(31))计算出大小齿轮的初始转角φg0、φp0。在φrange＝2Π/zp的范围内改变小齿轮转角φp。得到齿面接触轨迹点。接触轨迹如图8(a)所示。根据式(32)计算计算空载传递误差空载传递误差曲线如图8(b)所示。
(3)TCA得到齿面离散点坐标，插值生成齿轮网格。采用三维双线性等参单元建立有限元模型，计算总体刚度矩阵。在齿基底部添加除绕轴旋转外的约束；在不同啮合时刻下大小齿轮接触面齿面点分别施加接触点法向的单位载荷，提取在加载点齿面0.2m(m为齿轮的模数)以下位置的变形，组集成3维整体柔度矩阵和/>
(4)在接触点附近建立极坐标系(Lp，θ)(如图5(b))通过计算初始接触点附近Lp的齿面间隙dp，寻找最小间隙的θmin值，此方向即为接触椭圆的长轴方向，而最大间隙的θmax值即为接触椭圆的短轴方向。根据长轴方向计算齿面内潜在接触点的坐标，同时计算出大小齿轮对应潜在接触点之间的间隙c，如图9所示。考虑边缘接触的影响，还计算了超出齿轮面边界的初始接触点的位置。
(5)根据大小齿轮的整体柔度矩阵，先将网格点以样条插值的方法加密，如图10所示，然后用距离潜在接触点最近的插值或网格点的柔度作为潜在接触点的柔度，进而得到大小齿轮潜在接触点的柔度矩阵，根据公式整合为整体柔度矩阵λb。
(6)确定接触柔度λc。根据(35)进行变形协调迭代，每次迭代时，寻找Fn×1中小于0的行，将(35)左边两个矩阵对应的行列都置0，直至Fn×1中所有的数都大于0。联立式(36)和式(37)计算得到不同扭矩下健康齿轮的时变啮合刚度(如图11所示)。
如图12(a)所示，Len_1，Len_2，Len_3分别为不同长度的剥落区域，A、B、C为三个啮合位置，剥落厚度为0.1mm，潜在接触点改变后的间隙如图12(b)所示。
将剥落后的间隙进行变形协调即可计算出不同剥落长度下时变啮合刚度，如图13(b)所示。图13(a)为不同剥落长度下ANSYS的时变啮合刚度，和本发明对比来看，刚度影响范围一致，都只在AC范围内下降，且下降幅度误差很小，同时可以发现，剥落长度影响刚度的下降幅度。
如图14所示，Wid_1，Wid_2，Wid_3分别为不同宽度的剥落区域，D、E、F为三个啮合位置，剥落厚度为0.1mm。
将剥落后的间隙进行变形协调即可计算出不同剥落宽度下时变啮合刚度，如图15(b)所示，(a)为不同剥落宽度下ANSYS的时变啮合刚度，和本发明对比来看，刚度变化区域一致，不同宽度下的剥落影响范围也很吻合，从中可以得到，剥落宽度主要影响刚度变化范围。
如图16所示，Loc_1，Loc_2，Loc_3分别为不同长度的剥落区域，G、H、I为三个啮合位置，剥落厚度为0.1mm。
将剥落后的间隙进行变形协调即可计算出不同剥落位置下时变啮合刚度，如图17(b)所示。图17(a)为不同剥落位置下ANSYS的时变啮合刚度，和剥落长度和剥落宽度的结论相佐证。
本发明的具体应用途径很多，以上所述仅是本发明的优选实施方式。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (10)
1.一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法，其特征在于，包括以下步骤：
(1)采集弧齿锥齿轮的属性参数和剥落故障参数；
(2)同时联立求解啮合方程和旋转投影面方程得到弧齿锥齿轮大齿轮齿面上任意一点的坐标以及弧齿锥齿轮小齿轮齿面上任意一点的坐标；
(3)将弧齿锥齿轮小齿轮齿面方程和法向量从轮坯坐标系转化到装配坐标系；将弧齿锥齿轮大齿轮齿面方程和法向量从轮坯坐标系转化到装配坐标系；
(4)根据弧齿锥齿轮大、小齿轮连续相切方程组计算出齿轮副在装配坐标系下第一次啮合时各自所需要的旋转的角度φg0、φp0；其中φg0为弧齿锥齿轮大齿轮第一次啮合时所需要的旋转的角度；φp0为弧齿锥齿轮小齿轮第一次啮合时所需要的旋转的角度；
(5)以所述的φp0为中间值，上下波动一定量φrange假定小齿轮转角的区间[φp0-φrange,φp0+φrange]，在该区间内等间距均分得到不同啮合时刻下弧齿锥齿轮小齿轮的转角φp；通过转角φp解得初始接触点的未知变量值，代入弧齿锥齿轮大、小齿轮齿面方程得到初始接触点的坐标，每个初始接触点代表一个啮合时刻和齿轮啮合位置；
(6)弧齿锥齿轮大、小齿轮齿面上任意一点的坐标为离散点坐标，根据该齿面离散点坐标通过线性插值生成齿轮内部网格；采用三维双线性等参单元建立有限元模型，在齿基底部添加除绕轴旋转外的约束；在不同啮合时刻下大小齿轮接触面齿面点分别施加接触点法向的单位载荷，提取在加载点齿面以下位置的变形，组集成弧齿锥齿轮大齿轮3维整体柔度矩阵[λg]n*i*j及弧齿锥齿轮小齿轮3维整体柔度矩阵[λp]n*i*j和，n表示啮合时刻的个数，i表示齿面施力点的个数，j表示齿面提取位移点的个数；
(7)当受载时，弧齿锥齿轮的啮合状态由接触点变为接触椭圆，将接触椭圆表示为沿接触椭圆长轴的接触曲线；在接触点附近建立极坐标系(Lp，θ)，将初始接触点法向量沿着切平面移动Lp，分别与大小齿轮相交于两点，两点之间距离即为齿面间隙dp，计算不同θ下的齿面间隙，寻找最小间隙的θmin值，此方向即为接触椭圆的长轴方向，最大间隙的θmax值为接触椭圆的短轴方向；沿着长轴方向，计算不同距离Lp下的齿面间隙ε；
(8)根据[λg]n*i*j、[λp]n*i*j将网格点以样条插值的方法加密，然后用距离潜在接触点最近的插值点或网格点的柔度作为潜在接触点的柔度，得到弧齿锥齿轮小齿轮潜在接触点柔度矩阵λp和弧齿锥齿轮大齿轮潜在接触点柔度矩阵λg，然后整合为潜在接触点柔度矩阵λb；
(9)根据步骤(1)中采集的剥落顶点，在投影平面内寻找剥落区域穿过的潜在接触点，其为剥落影响的范围点，在剥落影响范围内的潜在接触点间隙增加剥落深度值，得到新的潜在接触点间隙εn×1；
(10)计算时变啮合刚度Kmesh：
式中，εi表示第i个潜在接触点的间隙，Fn_all表示接触力的合力；Ste表示承载传递误差。
2.根据权利要求1所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(1)中，属性参数包括齿轮参数、材料参数和部分加工参数；剥落故障参数包括剥落多边形的各个顶点坐标，剥落深度。
3.根据权利要求2所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(2)中，轮坯坐标系下的弧齿锥齿轮大齿轮齿面方程
表示大齿轮轮坯的旋转角度，sg表示直线刀刃长度，即切削深度，ɑg表示刀具齿形角，θg表示刀尖旋转的角度，M2g表示Sg到S2的坐标转化矩阵；Sg为刀具坐标系；S2为轮坯坐标系；rg(sg,θg)为刀具旋转所得到的曲面矢量；
同时联立求解下述啮合方程和旋转投影面方程得到弧齿锥齿轮大齿轮齿面上任意一点的坐标；
式中，Xm为旋转投影面x轴坐标，Ym为旋转投影面y轴坐标；
同理，轮坯坐标系下的弧齿锥齿轮小齿轮齿面方程
表示小齿轮轮坯的旋转角度，sp表示直线刀刃长度，即切削深度，θp表示刀尖旋转的角度，M1p表示Sp到S1的坐标转化矩阵，Sp为刀具坐标系，S1为轮坯坐标系，rp(sp,θp)为刀具旋转所得到的曲面矢量；
同时联立求解下述啮合方程和旋转投影面方程得到弧齿锥齿轮小齿轮齿面上任意一点的坐标；
4.根据权利要求3所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(4)中，根据大小齿轮连续相切方程组计算出齿轮副在装配坐标系下第一次啮合时各自所需要的旋转的角度φg0、φp0；
和/>分别为大、小齿轮装配坐标系下齿面方程，/>和/>分别为大小齿轮装配坐标系下齿面点的法向量方程。
5.根据权利要求4所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(5)中，以所述的φp0为中间值，上下波动一定量φrange假定小齿轮转角的区间[φp0-φrange,φp0+φrange]，在该区间内等间距均分得到不同啮合时刻下小齿轮的转角φp，带入权利要求4中方程组解得初始接触点的未知变量值，代入齿面方程r1，r2得到初始接触点的坐标，每个初始接触点代表一个啮合时刻和齿轮啮合位置；并计算出空载传递误差Δφg；
式中，zg和zp分别为大小齿轮的齿数。
6.根据权利要求5所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(8)中，
和/>分别为大小齿轮潜在接触点柔度矩阵对应元素值；i表示提取位移的接触点，j表示施加力的接触点，n为此接触时刻接触长轴的潜在接触点个数。
7.根据权利要求6所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(10)中，首先计算接触柔度λc，
λc＝diag(λc1,λc2,…,λci,…,λcn),
E为弹性模量；L是潜在接触点之间的距离；Fi是第i个潜在接触点上的法向接触力；
由进行变形协调迭代，
Fn×1表示接触点法向方向上的分布力；εn×1表示潜在接触点的间隙；Ste表示承载传递误差；Fn_all表示接触力的合力。
8.根据权利要求7所述的弧齿锥齿轮时变啮合刚度计算方法，其特征在于：步骤(10)中，δ1表示弧齿锥齿轮的分锥角；ɑn表示弧齿锥齿轮的法向压力角；β表示弧齿锥齿轮的平均螺旋角。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至8中任一项所述方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至8中任一项所述的方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311350092.1A CN117634057A (zh) | 2023-10-18 | 2023-10-18 | 一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311350092.1A CN117634057A (zh) | 2023-10-18 | 2023-10-18 | 一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117634057A true CN117634057A (zh) | 2024-03-01 |
Family
ID=90036650
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311350092.1A Pending CN117634057A (zh) | 2023-10-18 | 2023-10-18 | 一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117634057A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117892459A (zh) * | 2024-03-15 | 2024-04-16 | 南京航空航天大学 | 一种弧齿锥齿轮副磨损故障预测及其啮合特性分析方法 |
| CN117892459B (zh) * | 2024-03-15 | 2024-07-02 | 南京航空航天大学 | 一种弧齿锥齿轮副啮合特性分析方法 |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107436982A (zh) * | 2017-07-27 | 2017-12-05 | 东北大学 | 考虑基体刚度修正的剥落斜齿轮副的啮合特性分析方法 |
| WO2018086160A1 (zh) * | 2016-11-09 | 2018-05-17 | 北京工业大学 | 基于粗糙表面的直齿轮三维接触刚度计算方法 |
| CN109145484A (zh) * | 2018-09-04 | 2019-01-04 | 中南大学 | 基于双曲面壳单元模型的数值载荷齿面接触分析方法 |
| CN110334460A (zh) * | 2019-07-11 | 2019-10-15 | 西北工业大学 | 圆柱齿轮啮合刚度计算方法 |
| CN110826273A (zh) * | 2019-10-28 | 2020-02-21 | 长安大学 | 一种考虑浮动特性的行星传动多体齿轮承载接触特性分析方法 |
| CN114912256A (zh) * | 2022-04-22 | 2022-08-16 | 西安交通大学 | 含裂纹故障的弧齿锥齿轮振动响应分析方法 |
| CN114925465A (zh) * | 2022-04-22 | 2022-08-19 | 西安交通大学 | 含剥落故障的弧齿锥齿轮振动响应分析方法 |
| CN114943122A (zh) * | 2022-04-22 | 2022-08-26 | 西安交通大学 | 含点蚀故障的弧齿锥齿轮振动响应分析方法 |
| CN114969616A (zh) * | 2022-04-22 | 2022-08-30 | 西安交通大学 | 基于切片法的弧齿锥齿轮啮合刚度计算方法 |
| CN115795703A (zh) * | 2022-09-26 | 2023-03-14 | 东北大学 | 一种断齿故障下平行轴齿轮系统动力学仿真模型建立方法 |
-
2023
- 2023-10-18 CN CN202311350092.1A patent/CN117634057A/zh active Pending
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2018086160A1 (zh) * | 2016-11-09 | 2018-05-17 | 北京工业大学 | 基于粗糙表面的直齿轮三维接触刚度计算方法 |
| CN107436982A (zh) * | 2017-07-27 | 2017-12-05 | 东北大学 | 考虑基体刚度修正的剥落斜齿轮副的啮合特性分析方法 |
| CN109145484A (zh) * | 2018-09-04 | 2019-01-04 | 中南大学 | 基于双曲面壳单元模型的数值载荷齿面接触分析方法 |
| CN110334460A (zh) * | 2019-07-11 | 2019-10-15 | 西北工业大学 | 圆柱齿轮啮合刚度计算方法 |
| CN110826273A (zh) * | 2019-10-28 | 2020-02-21 | 长安大学 | 一种考虑浮动特性的行星传动多体齿轮承载接触特性分析方法 |
| CN114912256A (zh) * | 2022-04-22 | 2022-08-16 | 西安交通大学 | 含裂纹故障的弧齿锥齿轮振动响应分析方法 |
| CN114925465A (zh) * | 2022-04-22 | 2022-08-19 | 西安交通大学 | 含剥落故障的弧齿锥齿轮振动响应分析方法 |
| CN114943122A (zh) * | 2022-04-22 | 2022-08-26 | 西安交通大学 | 含点蚀故障的弧齿锥齿轮振动响应分析方法 |
| CN114969616A (zh) * | 2022-04-22 | 2022-08-30 | 西安交通大学 | 基于切片法的弧齿锥齿轮啮合刚度计算方法 |
| CN115795703A (zh) * | 2022-09-26 | 2023-03-14 | 东北大学 | 一种断齿故障下平行轴齿轮系统动力学仿真模型建立方法 |
Non-Patent Citations (3)
| Title |
|---|
| ZHANWEI LI 等: "Time-varying mesh stiffness calculation of spiral bevel gear with spalling defect", MECHANISM AND MACHINE THEORY, vol. 193, 27 December 2023 (2023-12-27), XP087447794, DOI: 10.1016/j.mechmachtheory.2023.105571 * |
| 安春雷;韩振南;: "点蚀与剥落对齿轮扭转啮合刚度影响的分析", 振动、测试与诊断, no. 04, 15 December 2008 (2008-12-15) * |
| 蔡守宇;王三民;袁茹;: "基于规划法的弧齿锥齿轮扭转啮合刚度特性研究", 机械科学与技术, no. 07, 15 July 2010 (2010-07-15) * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117892459A (zh) * | 2024-03-15 | 2024-04-16 | 南京航空航天大学 | 一种弧齿锥齿轮副磨损故障预测及其啮合特性分析方法 |
| CN117892459B (zh) * | 2024-03-15 | 2024-07-02 | 南京航空航天大学 | 一种弧齿锥齿轮副啮合特性分析方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Krol et al. | Parametric modeling of gear cutting tools | |
| Zhou et al. | A new closed-form calculation of envelope surface for modeling face gears | |
| Shih et al. | Free-form flank correction in helical gear grinding using a five-axis computer numerical control gear profile grinding machine | |
| CN110587038B (zh) | 一种刮齿加工齿廓误差补偿方法 | |
| Feng et al. | Geometric design and analysis of face-gear drive with involute helical pinion | |
| Sobolewski et al. | Method of spiral bevel gear tooth contact analysis performed in CAD environment | |
| CN115841548B (zh) | 一种叶片模型的计算机辅助生成方法及系统 | |
| Shih et al. | A novel method for producing a conical skiving tool with error-free flank faces for internal gear manufacture | |
| CN111176209B (zh) | 型腔螺旋铣削加工进给率与转速离线规划方法 | |
| Fan et al. | New developments in tooth contact analysis (TCA) and loaded TCA for spiral bevel and hypoid gear drives | |
| CN113486475B (zh) | 一种圆柱齿轮滚齿加工切削力的预测方法 | |
| Katz et al. | Virtual model of gear shaping—part ii: Elastic deformations and virtual gear metrology | |
| CN117634057A (zh) | 一种含剥落故障的弧齿锥齿轮时变啮合刚度计算方法 | |
| JP4763611B2 (ja) | 研ぎ直しピニオンカッタの刃形輪郭の評価方法 | |
| CN109033669B (zh) | 基于万能运动参数驱动的螺旋锥齿轮仿真加工建模方法 | |
| CN113192180B (zh) | 一种基于插齿加工原理的椭圆齿轮参数化精确建模方法 | |
| Niu et al. | Geometrical design of variable ratio tooth profile based on Boolean subtraction operation and a novel modification method | |
| CN113010978A (zh) | 一种基于动态仿真的航空直齿轮修形方法 | |
| CN109492307B (zh) | 一种弧齿锥齿轮齿面载荷接触性能参数的数值计算方法 | |
| CN115081143B (zh) | 一种蜗杆砂轮磨齿加工磨削力的仿真计算方法及系统 | |
| Tsai | Gear module adjustment method for shrinkage compensation of injection molded plastic gears using single dedicated hob | |
| Chlost et al. | A new method of the positioning and analysis of the roughness deviation in five-axis milling of external cylindrical gear | |
| Pisula et al. | Numerical model of bevel gears cutting by duplex helical method | |
| Xue et al. | An analytical prediction methodology of the correlation between transmission error and tooth pitch error for face-hobbed hypoid gears | |
| Vivet et al. | Advanced modelling of straight bevel gear contact for multibody simulations |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |