行星齿轮传动具有体积小、质量轻、传动比大、效率高等特点，在机械行业广泛应用．关于行星齿轮传动研究已取得许多成果．杨望灿等[1]提出了一种基于自适应学习率的深度学习的行星齿轮箱故障诊断模型．Chaari等[2]建立了行星齿轮系统平面模型，分析了制造误差对动力学性能的影响．Ericson等[3]针对两种直齿行星齿轮装置进行了模态振动试验测试研究，测量了动响应、固有频率和振型．Xue等[4]提出了单支太阳轮-行星轮-内齿圈以及整个行星轮系扭转刚度模型．Ren等[5]考虑制造偏心误差、齿廓误差、支撑变形和时变啮合刚度等，提出了弯-扭-轴人字齿轮动力学模型，研究了浮动特性．Tatar等[6]提出了六自由度行星齿轮转子系统动力学模型，该模型采用集中参数齿轮箱模型和转轴Timoshenko梁有限元模型．Bahk等[7]提出了一种可以获得太阳轮-行星轮和内齿圈-行星轮啮合处齿廓修形激励的模型．张俊等[8]为抑制系统的振动和噪声，以动态传动误差波动量为指标，研究了斜齿行星轮系的齿轮修形策略．马辉等[9]研究了齿顶修形对啮合刚度的影响．Wang等[10]以电动汽车行星齿轮动力总成为研究对象，提出了一种基于齿面和齿背三维修形的抑制瞬态扭转振动方法．Öztürk等[11]采用承载静态传输误差(承载传动误差)计算时变啮合刚度函数，分析了线性和抛物线齿廓修形对行星齿轮系统动态响应特性的影响．齿轮修形在减小误差敏感性、减振降噪、改善齿面载荷分布等方面具有重要意义，一直为研究者广泛关注[12]，但是目前对行星轮系统修形优化技术的研究文献很少．考虑到直齿行星传动是一种常见行星传动形式，随着技术发展，对其传动性能要求越来越高．为了提高直齿行星轮系的性能，须要对其修形优化技术进行深入研究．本研究提出了直齿轮承载接触分析方法及修形优化方法，并对修形效果进行动力学分析验证，可为行星轮工程修形设计提供理论依据．1 直齿行星轮系齿轮参数本研究的直齿行星轮系输入轴为太阳轮轴，输出轴为行星架，内齿圈固定．行星轮系齿轮参数如表1所示(分度圆模数为2.0 mm，分度圆压力角为20°)，行星轮数目为3．10.13245/j.hust.240604.T001表1某行星轮系齿轮参数参数太阳轮行星轮内齿圈齿数172567齿顶高/mm2.000 02.000 01.774 6齿根高/mm2.52.52.5内齿圈齿顶高计算公式[13]为har=(1-7.55/zr)m，(1)式中har，zr和m分别为内齿圈齿顶高、齿数和分度圆模数．2 直齿行星轮系承载接触分析方法2.1　直齿行星轮系承载接触分析模型直齿轮承载接触分析模型如图1所示．10.13245/j.hust.240604.F001图1直齿轮承载接触分析模型采用数学规划法求解加载轮齿接触问题，min∑j=12n+1Xj;-Fp+Z+d+X=w;P=eTp+X2n+1;s.t.    pj≥0,dj≥0,Xj≥0,Z≥0;pj=0或dj=0, (2)式中：Xj( j=1，2，…，2n+1)为人工变量，X=[X1，X2，…，X2n]T；e为单位2n维列向量；F为齿轮副的齿面法向柔度矩阵；p为离散点载荷向量；pj为齿对k的瞬时接触椭圆长轴离散点j处的法向载荷；Z=Z [1，1，…，1]T，Z为轮齿变形后的法向位移(实际上就是啮合时的变形，因此式(2)可以求解得到啮合变形Z，求解时将一个啮合周期离散为若干啮合位置，可以得到按照啮合周期连续变化的啮合变形量)；w为啮合点处初始法向齿面间隙向量；d为啮合点处变形后的法向齿面间隙向量；dj为齿对k的瞬时接触椭圆长轴离散点j处变形后的齿面间隙；P为总法向载荷．式(2)是由已知参数F，P，w和未知参数p，d，Z组成的一个非线性规划，目标函数是变形能最小，用一种改进的正规形法进行求解，从而可以求得一个啮合周期内当前啮合位置的Z．详细的承载接触分析求解过程见文献[14-15]．2.2　直齿行星轮系承载接触分析方法验证为了验证以上方法的正确性，与Romax进行对比分析，结果如图2所示，图中：δex为外啮合分析时沿啮合线方向的承载传动误差；θs为太阳轮转角；δin为内啮合分析时沿啮合线方向的承载传动误差；θp为行星轮转角．10.13245/j.hust.240604.F002图2承载传动误差对比从图2可以看出：Romax方法和本文承载接触分析方法承载传动误差曲线及波动幅值接近，外啮合波动幅值相对误差为5.2%，内啮合波动幅值相对误差为-1.5%，说明本文方法计算承载传动误差可行．3 直齿行星轮系齿轮修形优化方法3.1　行星轮系齿轮修形方式修形模式可以采用大轮和小轮同时齿廓及齿向双向修形的方法，也可以只进行大轮或小轮双向修形(对于外啮合齿轮副，小轮为太阳轮，大轮为行星轮；对于内啮合齿轮副，小轮为行星轮，大轮为内齿圈)．这里采用小轮齿廓和齿向修形的方法，对于行星轮系中的外啮合齿轮只进行太阳轮修形，对于内啮合齿轮副只进行行星轮齿廓和齿向修形．齿廓和齿向修形曲线有直线和抛物线等，这里采用抛物线修形曲线，且齿向采用鼓形修形曲线．本研究采用沿齿高方向表示的齿廓修形图如图3所示，图中：y1为齿顶最大修形量；y2为齿高方向齿顶修形长度；y3为齿根最大修形量；y4为齿高方向齿根修形长度．10.13245/j.hust.240604.F003图3修形方式齿向采用鼓形修形，y5为最大齿向修形量，考虑到直齿轮齿向修形对承载传动误差影响较小，齿向修形参数不进行优化，而采用鼓形修形经验确定．对于高精度、高可靠度的高速齿轮，鼓形修形量为10~25 μm，加上制造误差5 μm．3.2　行星轮系齿轮修形优化变量及目标根据修形方式可知修形参数为小轮齿廓修形参数和齿向鼓形修形参数，考虑到直齿轮理想安装条件下齿向修形量对承载传动误差影响比较小，为了便于优化，齿向修形量不进行优化，直接取经验值y5=15 μm，因此这里只对齿廓修形进行优化，修形优化参数为y1，y2，y3和y4．承载传动误差反映了齿轮实际啮合与理想啮合的偏差情况，大量研究表明承载传动误差幅值变化是工作过程中振动的直接激励，是产生振动、噪声的重要因素．直齿行星轮系齿轮修形优化目标就是降低承载传动误差波动幅值，本研究采用承载接触分析方法计算承载传动误差．采用承载接触分析方法计算得到轮齿承载后的法向位移Z，将Z转化为用转角形式表示，即一个啮合周期的承载传动误差．3.3　直齿行星轮系齿轮修形优化方法及流程齿轮修形优化是一个多变量的复杂非线性优化问题，目前用的比较多的是遗传算法，当遗传算法进化代数及种群规模比较大时可以取得较好的优化结果．但是随之而来的缺点是优化时间比较长，不利于工程应用．本研究提出基于罚函数法的优化算法，该算法耗时少，为了增强优化效果，初值可以随机取多组数据．基于罚函数法的修形优化方法见下式，Y=y1,y2,y3,y4;gu(Y)=-Y(u)≤yu_min      (u=1,2,3,4);gu(Y)=Y(u-4)≤y(u-4)_max      (u=5,6,7,8);f(Y)=max(Te)-min(Te);φ(Y,r(k))=f(Y)-r(k)∑u=181/gu(Y),(3)式中：Y为修形参数组成的向量；gu(Y)为约束函数；yu_min和y(u-4)_max为修形变量的下限和上限；Te为承载传动误差；f (Y)为承载传动误差波动幅值；φ(Y，r(k))为罚函数；r(k)为惩罚因子．优化流程如图4所示．优化时分别按内啮合齿轮副和外啮合齿轮副单独优化．10.13245/j.hust.240604.F004图4修形优化流程4 直齿轮传动修形优化设计4.1　直齿轮副工况参数修形设计与工况密切相关，这里行星轮系修形工况为功率12 kW，根据行星轮系工况参数等可以计算得到内外啮合直齿轮副工况参数．4.2　内啮合齿轮副修形优化设计根据修形优化设计方法，进行内啮合齿轮副优化，修形优化后得到行星轮的齿廓修形参数．内啮合齿面修形优化和不修形仿真结果见图5，图中ε为承载传动误差，修形前承载传动误差波动幅值为4.81"，修形后波动幅值为1.15"，可以看出修形优化后承载传动误差波动幅值有明显的减小．10.13245/j.hust.240604.F005图5内啮合齿轮承载传动误差对比4.3　外啮合齿轮副修形优化设计同样根据修形优化设计方法，外啮合齿轮副修形优化后得到太阳轮齿廓修形参数，外啮合齿面修形优化和不修形仿真结果见图6，图中修形前承载传动误差波动幅值为17.36"，修形后波动幅值为4.80"，可以看出修形优化后承载传动误差波动幅值显著减小．10.13245/j.hust.240604.F006图6外啮合齿轮承载传动误差对比5 动力学建模及分析5.1　动力学模型为了分析修形对行星齿轮传动系统振动特性的影响，须要建立行星齿轮系统动力学模型．在动力学建模之前，首先须要根据承载传动误差得到时变啮合刚度．原始承载传动误差表示了承载传动误差与小轮转角的关系，根据转速条件可以得到承载传动误差与时间的关系，载荷除以承载传动误差就可以得到啮合刚度，这样就得到啮合刚度表达式ksn=Ts/(3rbs)εexn(t);krn=Tp/rbpεinn(t), (4)式中：n=1,2,3；ksn和krn为第n对外啮合和内啮合齿轮副啮合刚度；齿轮Ts和Tp为太阳轮和行星轮力矩；rbs和rbp为太阳轮和行星轮基圆半径；εexn(t)和εinn(t)为第n对外啮合和内啮合齿轮副沿啮合线方向的承载传动误差．本文算例中的行星齿轮系统物理模型见图7，参考文献[16]建立包含支撑阻尼的动力学模型．10.13245/j.hust.240604.F007图7直齿行星齿轮系统物理模型行星齿轮系统动力学模型的广义坐标为      X=[xs,ys,us;xc,yc,uc;xr,yr,ur;xn,yn,un]      (n=1,2,3).构件之间的相对位移如下    δsn=(us-xssin ϕsn+yscos ϕsn)-(-un+xnsin α+yncos α);    δrn=(ur-xrsin ϕrn+yrcos ϕrn)-(un-xnsin α+yncos α);    δpnx=-(xn-xccos ϕn-ycsin ϕn);    δpny=-(yn-uc+xcsin ϕn-yccos ϕn),式中：下标s，c，r，p分别为太阳轮、行星架、内齿圈和行星轮；x，y，u分别为三个方向的位移；us和ur为啮合线方向位移，uc为行星架沿旋转切线方向在行星轮中心处的位移；δsn为太阳轮与行星轮n之间的相对位移在啮合线方向的投影；δrn为内齿圈与行星轮n之间的相对位移在啮合线方向的投影；δpnx和δpny分别为行星架与行星轮n之间沿径向和切向的相对位移；ϕsn=ϕn-α，ϕrn=ϕn+α，ϕn为行星轮与太阳轮之间的相位角，α为齿轮副啮合角．根据直齿行星齿轮系统物理模型，分别建立太阳轮、行星架、内齿圈和行星轮动力学方程    ms(x¨s-2ωcy˙s-ωc2xs)+csxx˙s+ksxxs=∑Fsnsin ϕsn,    ms(y¨s+2ωcx˙s-ωc2ys)+csyy˙s+ksyys=-∑Fsncosϕsn,    (Is/rs2)u¨s+csuu˙s+ksuus=Ts/rs-∑Fsn;    mc(x¨c-2ωcy˙c-ωc2xc)+∑kpx(δpnxcos ϕrn-δpnysin ϕrn)+∑cpx(δ˙pnxcos ϕrn-δ˙pnysin ϕrn)+ccxx˙c+kcxxc=0,    mc(y¨c+2ωcx˙c-ωc2yc)+∑kpy(δpnxsin ϕrn+δpnycos ϕrn)+∑cpy(δ˙pnxsin ϕrn+δ˙pnycos ϕrn)+ccxy˙c+kcyyc=0,    (Ic/rc2)u¨c+∑cpuδ˙pny+∑kpuδpny+ccuu˙c+kcuuc=Tc/rc;    mr(x¨r-2ωcy˙r-ωc2xr)+crxx˙r+krxxr=∑Frnsin ϕrn,    mr(y¨r+2ωcx˙r-ωc2yr)+cryy˙r+kryyr=-∑Frncos ϕrn,    (Ir/rr2)u¨r+cruu˙r+kruur=-∑Frn;    mp(x¨n-2ωcy˙n-ωc2xn)-cpxδ˙pnx-kpxδpnx=Fsnsin αsp-Frnsin αrp,    mp(y¨n+2ωcx˙n-ωc2yn)-cpyδ˙pny-kpyδpny=Fsncos αsp+Frncos αrp,    (Ipn/rpn2)u¨n=Frn-Fsn,式中：ms，mc，mr，mp分别为太阳轮、行星架、内齿圈和行星轮的质量；rs，rc，rr，rpn分别为太阳轮、行星架、内齿圈以及第n个行星轮的基圆半径；ksx，ksy，ksu，kcx，kcy，kcu，krx，kry，kru，kpx，kpy分别为太阳轮、行星架、内齿圈和行星轮在x，y自由度上的支撑刚度和u自由度上的扭转支撑刚度(平移方向)；csx，csy，csu，ccx，ccy，ccu，crx，cry，cru，cpx，cpy分别为太阳轮、行星架、内齿圈以及行星轮在x，y平移方向上的支撑阻尼和u方向上的扭转支撑阻尼(平移方向)；Ts，Tc分别为太阳轮输入扭矩和行星架的输出扭矩；αsp，αrp分别为太阳轮、内齿圈与行星轮的啮合角；Ipn，Is，Ic，Ir分别为行星轮、太阳轮、行星架和内齿圈的转动惯量；ωc为行星架角速度．假设齿轮啮合间隙为bb(本文中bb取0.2 mm)，行星轮-太阳轮、行星轮-内齿圈啮合副之间的啮合力分别为Fsn=ksnf(δsn,bb)+csnδ˙sn;Frn=krnf(δrn,bb)+crnδ˙rn,式中：Fsn和Frn为外啮合及内啮合的啮合力；f (δsn，bb)和f (δrn，bb)为啮合时变形量函数；ksn，krn，csn，crn分别为齿轮副之间的啮合刚度和啮合阻尼；ζ为阻尼比(取值范围一般为0.03~0.17，本研究取值为0.1)，行星轮与太阳轮副、行星轮与内齿圈副的啮合阻尼csn=2ζkspnIsIp/(Isrpn2+Iprs2);crn=2ζkrpnIrIp/(Irrpn2+Iprr2),式中kspn和krpn为第n对外啮合平均啮合刚度和内啮合平均啮合刚度．5.2　动力学分析振动加速度是评价齿轮系统振动的重要指标，考虑到振动通过轴承传递到箱体，然后通过箱体传递到其他装置．这里分析轴承处的振动加速度，并只考虑太阳轮、内齿圈和行星架轴承支撑处的振动加速度．首先仿真得到不修形的结果如图8所示；然后分别仿真分析外啮合齿轮修形、内啮合齿轮修形及内外啮合齿轮同时修形的情况，仿真分析结果见表2，as-max，ar-max和ac-max分别为太阳轮支撑轴承、内齿圈支撑轴承和行星架支撑轴承最大加速度．10.13245/j.hust.240604.F008图8修形前动力学特性10.13245/j.hust.240604.T002表2振动加速度分析修形|as-max||ar-max||ac-max|xyxyxy修形前9.70313.8171.6342.3650.36147.094外修形优化4.5386.4731.6382.3710.36147.081内修形优化9.69713.8090.8031.1620.36147.083内外修形优化4.5356.4690.8031.1620.36147.069m/s2根据表2可以看出：外啮合修形后能减小太阳轮支撑处的振动加速度，内啮合修形后能减小内齿圈支撑处的振动加速度，内、外修形对行星架支撑处振动加速度影响很小．另外，太阳轮支撑处与内齿圈支撑处振动加速度相比，太阳轮支撑处振动加速度较大，外啮合修形减振效果很好，因此建议只进行太阳轮修形，不仅成本低，而且减振效果明显．6 结语a. 研究了内外啮合直齿轮承载接触分析方法．针对算例进行承载接触分析方法与Romax软件计算结果对比分析，两种方法得到的承载传动误差曲线及波动幅值接近，外啮合波动幅值相对误差为5.2%，内啮合波动幅值相对误差为-1.5%，对比分析证明了本文提出的承载接触分析方法的可行性．b. 提出了直齿行星传动内外啮合齿轮副的罚函数修形优化方法，并进行了算例直齿行星传动系统修形优化设计，使内啮合及外啮合齿轮副的承载传动误差波动幅值显著下降(分别下降76.1%和72.4%)．c. 将承载传动误差转换为啮合刚度，建立了直齿行星传动系统动力学模型，分析了内外啮合齿轮副修形对轴承支撑处振动加速度的影响，结果表明外啮合修形优化可以显著减小太阳轮支撑处振动加速度(算例中asx和asy都减少了53.2%)，内啮合修形优化可以减小内齿圈支撑处振动加速度(算例中arx和ary都减少了50.9%)，虽然内外同时修形能取得较好效果，但外啮合修形减振效果明显，因此建议只进行太阳轮修形．