0 引言柔性均载方法常用于解决风电齿轮箱行星传动机构因制造与安装误差引起的载荷分布不均问题，可以有效地解决风电主齿轮箱因不均载现象而产生的振动、轴承失效、齿轮断裂等一系列故障[1-4]。采用柔性行星销轴是常用的柔性均载方式之一。柔性销主要是利用其变形特性提高风电齿轮箱行星轮系的均载性能：一是其变形的挠度比普通销大，可以降低行星轮系的齿间均载系数；二是其变形前后接触位置的轴线接近平行，可以改善行星齿轮的齿面载荷分布。目前，国内外的风电企业关于柔性销主要开展的是等截面式直销的研究与应用，包括柔性销轴的均载性能影响因素、柔性销结构的有限元仿真分析与计算、柔性销结构的行星轮系均载系数计算、柔性销的有限元结构优化设计方法以及柔性销零件参数的计算等[5-9]。变截面构件使用较多的是建筑行业，因此，可以结合其他变截面构件的研究方法来研究变截面柔性销。李银山等[10]研究了阶梯形变截面梁、等高变宽变截面梁和等宽变高变截面梁变形的快速解析法；刘庆玲[11]研究了任意载荷作用下的变截面柔性构件的变形特性；周期源等[12]计算了变截面梁在考虑剪切变形影响时的挠度；刘胜来[13]利用格林法求解了变刚度梁挠度曲线的微分方程。变截面柔性销的变形计算过程相较于等截面柔性销要复杂，一些研究在计算时都忽略了剪切变形的影响，从而导致计算结果和有限元仿真与实际工程情况产生较大的误差。同时，目前关于变截面形式的柔性销研究较少。因此，本文建立了考虑剪切效应时，变截面柔性销的挠曲线微分方程，并利用格林函数法进行求解，分析其在额定工况下的变形特性；对采用变截面柔性销结构的风电齿轮箱行星传动进行仿真分析，研究销轴支承刚度和啮合错位对均载性能的影响规律，并与采用普通行星销轴的风电齿轮箱行星结构进行了比较分析。1 变截面柔性销的结构与受力分析柔性行星销轴组件由柔性销轴和套筒组成。柔性销的一端与行星架连接，另一端与套筒连接。套筒通过行星轮轴承和行星齿轮相连接。如图1所示，柔性销根据截面形状可分为等截面销以及变截面销。等截面柔性销是根据最大弯矩所在截面的工作应力来确定直径，在柔性销受载时，两端弯矩大而中间弯矩小，等截面结构使柔性销在中间段的材料未能充分发挥作用。而变截面柔性销在设计时，按弯矩沿销的轴线变化情况相应地改变截面尺寸。在两端连接直径相同的情况下，为了得到相同的最大挠度，变截面柔性销的轴向尺寸比等截面柔性销的更短，因此能减少箱体尺寸与质量，节省材料，降低吊装成本。如图2所示，变截面柔性销通常由大端、小端与变截面部分组成，变截面部分为两段圆弧面和一段锥面，锥面的角度可为0°。图 1等截面柔性销组件与变截面柔性销组件Fig. 1Constant cross-sectional flexible pin and the variable cross-sectional flexible pin10.16578/j.issn.1004.2539.2023.10.017.F1a1（a）等截面柔性销组件10.16578/j.issn.1004.2539.2023.10.017.F1a2（b）变截面柔性销组件 10.16578/j.issn.1004.2539.2023.10.017.F002图2变截面柔性销结构与尺寸Fig. 2Structure and dimension of the variable cross-sectional flexible pin如图3所示，对变截面柔性销轴及其套筒进行受力分析。将行星轮轴承作用于套筒上的力简化为一集中力F，记F的作用点与行星架的距离为K。此时，套筒右端受到柔性销的反作用力为力F1和转矩M1，且右端弯矩大小为图3变截面柔性销及其套筒的受力分析Fig. 3Force analysis of the variable cross-sectional flexible pin and its sleeve10.16578/j.issn.1004.2539.2023.10.017.F3a1（a）柔性销与套筒受力分析10.16578/j.issn.1004.2539.2023.10.017.F3a2（b）柔性销受力分析简图 M1=F1(L-K)=F(L-K) （1）视行星架为刚性结构，则柔性销左端为固定连接，柔性销右端受力为1个集中力和1个弯矩，作用点在右端接触面中点。因此，左端弯矩与右端弯矩的关系为M2=M1-F1L=-FK （2）因此，柔性销各个截面上的剪力Q=F，而各个截面的弯矩为M(x)=M1-Q(L-x)=F(x-K) （3）2 变截面柔性销的挠曲线微分方程求解2.1　考虑剪切效应的挠曲线微分方程根据材料力学中求解弯曲变形的方法，不考虑剪切变形时，变截面柔性销的挠曲线微分方程为d2wMdx2=M(x)EI(x) （4）其转角方程和挠度方程分别为θ=∫M(x)EI(x)dx+C1 （5）wM=∬M(x)EI(x)dxdx+C1x+C2 （6）式（4）适用于跨度远大于截面尺寸的梁，忽略了剪力对弯曲变形的影响。铁木辛柯提出，横力弯曲时的梁，除了在弯矩作用下产生挠度wM，还会在剪切作用下产生附加挠度wQ[14]。附加挠度wQ的斜率为截面中性轴上的剪应变，即dwQdx=γ=QksGA(x) （7）式中，ks为考虑泊松比效应的剪切修正系数，对于圆形截面，ks=6(1+ν)(7+6ν)，ν为泊松比。根据变截面柔性销的受力分析可知，各个截面的剪力Q为常数。因此，考虑剪切效应时，变截面柔性销的挠曲线微分方程为d2wdx2=M(x)EI(x)-QdA(x)ksGA2(x)dx （8）记式（8）的右端为f(x)，且记f1(x)=M(x)/EI(x)。为了计算方便，在求解挠曲线微分方程时，可以先根据式（6）求解wM，而wQ=∫γdx，则w=wM+wQ。受剪切变形的影响，截面转角的计算不再为dw/dx，而是θ=dw/dx-γ，即θ=dwMdx （9）2.2　积分法与叠加法由式（4）~式（6）可知，在使用积分法求解变截面柔性销的挠度和转角时，被积函数f1(x)为复杂解析函数，因此，不定积分的解一般不能用初等函数表示。由图2可知，变截面柔性销的母线方程需用4段函数表示，则被积函数f1(x)也需用4段函数表示。因此，在计算4个区间的不定积分时，需要求解8个待定积分常数。为了避免冗繁的定积分求解和待定积分常数的求解，常用的做法是将变截面部分的f1(x)进行多项式拟合，将其转换成有理函数，则不定积分的解可以获得初等函数表达式。拟合后被积函数减少为2段，此时需要求解的待定积分常数减少为4个。然而，使用多项式拟合被积函数的方法在计算时引入了插值误差，在挠度曲线的求解精度上会有损失，而且不方便进行误差估计。同时，为了提高拟合精度，相应地需要提高拟合阶数，从而造成高阶项系数过小，会产生振荡现象。本文案例的变截面柔性销f1(x)进行3阶多项式拟合后，最高项系数为-4.574 9×10-12。对于不同尺寸的变截面柔性销，其积分法计算过程不便于进行编程。由于上述原因，在初步设计时，变截面柔性销挠度与转角计算很少用到积分法。目前，企业中常用的变截面柔性销挠度与转角的计算方法为叠加法。叠加法是求解弹性小变形梁弯曲变形的常用方法。在已有的文献中，叠加法一般用来计算阶梯形变截面梁的变形。对于连续型变截面梁，可以将其切分成若干微段，每一段截面取平均值，从而等价为阶梯形变截面梁，然后分段计算挠度和转角并进行叠加。使用叠加法计算变截面柔性销的过程如图4所示，取每一微段的长度为l。取第1微段、第2微段为研究对象，首先刚化第2段，计算节点2的转角和挠度；然后刚化第1段，计算此时节点3的转角和挠度，再加上第1段变形产生的刚性位移，得到节点3的实际转角和挠度。图4叠加法分段计算挠度和转角Fig. 4Calculation of deflection and rotation angle in sections by the superposition method10.16578/j.issn.1004.2539.2023.10.017.F4a1（a）取柔性销前2个微段为研究对象10.16578/j.issn.1004.2539.2023.10.017.F4a2（b）刚化柔性销第2微段10.16578/j.issn.1004.2539.2023.10.017.F4a3（c）刚化柔性销第1微段 以此类推，第i个节点的挠度和转角为wi=∑k=2i(wk,1)F+(wk,1)Mk+θk-1l （10）θi=∑k=2i(θk,1)F+(θk,1)Mk （11）企业在应用叠加法时忽略了剪切变形的影响。根据式（8）可知，使用叠加法计算式（10）中的(wi)F时需要加上(wi)Q。叠加法的计算精度受微段长度l的影响。计算任意一点的挠度时需要计算前面各点的挠度和转角进行叠加。2.3　结合格林函数法的求解方程格林函数法常用于求解有初始条件的非齐次微分方程。式（4）所示为不考虑剪切效应的变截面柔性销的挠曲线微分方程，可视为以下2阶非齐次常微分方程在边界条件下的求解：wM″(x)-f1(x)=0wM(0)=w1wM(L)=w2 （12）根据2阶常微分方程格林函数的构造方法[15]，式（4）的格林公式为G(x,xa)=x(L-xa)L,0≤x≤xaxa(L-x)L,xax≤L （13）根据叠加原理，式（8）有形式解wM=u+v，其中，u″=-f1(x)；u(0)=0；u(L)=0；v″=0;v(0)=w1;v(L)=w2。利用特征方程可以求出v对应的齐次方程的通解为v(xa)=C1+C2xa，代入边值条件后求得C1=w1，C2=(w2-w1)/L。根据格林函数的性质[16]，u的解可写为u(xa)=-∫0LG(x,xa)f1(x)dx （14）因此，不考虑剪切变形时，变截面柔性销的挠度方程可写为wM(xa)=u(x)+v(x)=-∫0LG(x,xa)f1(x)dx+L-xaLw1+xaLw2=(xa-L)L∫0xaxM(x)EI(x)dx+xaL∫xaL(x-L)M(x)EI(x)dx+L-xaLw1+xaLw2 （15）变截面柔性销的截面转角方程为θ(xa)=w2-w1L+1L∫0xaxM(x)EI(x)dx+1L∫xaL(x-L)M(x)EI(x)dx （16）根据变截面柔性销的受力分析，可视其为悬臂结构。因此，其左端的挠度和转角分别为wM(0)=0，θ(0)=0，则w2=∫0L(x-L)M(x)EI(x)dx （17）结合式（7）与边界条件可以得到，考虑剪切效应时，变截面柔性销任意一点的挠度方程为w(xa)=xaL∫0Lx-LM(x)EI(x)dx+(xa-L)L∫0xaxM(x)EI(x)dx+xaL∫xaL(x-L)M(x)EI(x)dx+∫0LQksGA(x)dx （18）转角计算方程为θ(xa)=1L∫0Lx-LM(x)EI(x)dx+1L∫0xaxM(x)EI(x)dx+1L∫xaLx-LM(x)EI(x)dx （19）对比积分法，利用本文方法求解的挠度方程与转角方程将不定积分转化为定积分形式，避免了求解多个积分常数的冗繁过程，在数值计算软件中可以不计算原函数而得出计算值，便于编程实现。对比叠加法，本文方法求解各点的计算过程相互独立，在计算任意一点的挠度时，无须计算前面各点的挠度和转角，加快了计算效率，同时避免了累积误差。2.4　变截面柔性销的计算实例与变形特性以某大兆瓦风电主齿轮箱一级行星传动上的变截面柔性销轴为例，其弹性模量E=2.06 GPa，泊松比为0.3，销轴受力简化载荷F=1 455 238 N，弯矩M=253 939 031 N∙mm。柔性销的基本尺寸数据如表1所示。10.16578/j.issn.1004.2539.2023.10.017.T001表1变截面柔性销尺寸Tab. 1Dimensions of the variable cross-sectional flexible pinD1/mmD2/mmD3/mmR1/mmR2/mmα/(°)H1/mmH2/mmH3/mm2451603503502000210270175根据本文推导的挠度方程与转角方程，计算变截面柔性销各个位置的挠度如表2所示。由表2可以看出，随着长度的累加，不考虑剪切效应时，各点挠度的计算误差越来越大，末端的相对挠度误差达到了46.4%。因此，考虑剪切效应对变截面柔性销变形的影响十分有必要。10.16578/j.issn.1004.2539.2023.10.017.T002表 2变截面柔性销各个位置的挠度Tab. 2Deflection of the variable cross-sectional flexible计算位置考虑剪切的挠度不考虑剪切的挠度差值35.75-0.024-0.005-0.01871.5-0.067-0.024-0.042107.25-0.132-0.060-0.072143-0.218-0.111-0.107178.75-0.317-0.173-0.144214.5-0.417-0.237-0.180250.25-0.507-0.291-0.216286-0.577-0.331-0.246321.75-0.628-0.353-0.275357.5-0.655-0.351-0.304pin at each position 单位：mm比较本文方法、工程方法与有限元法在柔性销与套筒接触中点（x=L）处的挠度如表3所示，有限元计算的挠度值为-0.660 9 mm，本文方法计算的挠度值为-0.655 1 mm（与有限元结果相差-0.88%），结合考虑剪切效应的挠曲线方程，改进后的工程叠加法计算的挠度值为-0.694 8 mm（与有限元结果相差5.13%），3种计算方法的结果十分接近。同时，3种方法计算得到的挠度曲线如图5所示，说明本文推导求解的考虑剪切变形得到的挠度方程计算有效。10.16578/j.issn.1004.2539.2023.10.017.T003表3变截面柔性销末端挠度Tab. 3Deflection at the end of the variable cross-sectional flexible pin方法挠度/mm与有限元结果的相对误差/%本文方法-0.655 1-0.88工程方法-0.694 85.13有限元法-0.660 9/10.16578/j.issn.1004.2539.2023.10.017.F005图5变截面柔性销的挠度曲线Fig. 5Deflection curves of the variable cross-sectional flexible pin由图5~图7可以看出，变截面柔性销受载后的变形曲线呈“Z”字形，柔性销的两端轴线变形后依然接近平行。由图8可以看出，柔性销末端的转角接近0 rad。由于行星齿轮是通过套筒与柔性销末端相连接的，因此，使用柔性销时其在受载后可以保持齿面平行、减少偏载，从而改善齿面载荷分布及齿间均载性能。10.16578/j.issn.1004.2539.2023.10.017.F006图6变截面柔性销的总变形Fig. 6Total deformation of the variable cross-sectional flexible pin10.16578/j.issn.1004.2539.2023.10.017.F007图7变截面柔性销中性轴的变形Fig. 7Deformation of the variable cross-sectional flexible pin neutral axis10.16578/j.issn.1004.2539.2023.10.017.F008图8变截面柔性销的转角曲线Fig. 8Rotation angle curve of the variable cross-sectional flexible pin3 变截面柔性销的均载作用及影响因素为了探究变截面柔性销的变形特性对风电齿轮箱行星轮系的均载作用，在Romax软件中建立采用变截面柔性销的风电齿轮箱仿真模型，分析其支承刚度对均载性能的影响，同时建立采用普通行星销轴的模型，对比不同啮合错位下的行星轮齿面载荷分布情况。大兆瓦风电齿轮箱及其变截面柔性销结构的仿真模型分别如图9、图10所示，其输入转速为11 r/min，输入转矩为3.82×106 N∙m，一级行星传动部分的主要参数如表4所示。10.16578/j.issn.1004.2539.2023.10.017.F009图9风电齿轮箱仿真模型Fig. 9Simulation model of the wind turbine gearbox10.16578/j.issn.1004.2539.2023.10.017.F010图10变截面柔性销仿真模型Fig. 10Simulation model of the variable cross-sectional flexible pin10.16578/j.issn.1004.2539.2023.10.017.T004表4一级行星传动的主要参数Tab. 4Main parameters of the first planetary gear train构件太阳轮s行星轮p内齿圈r齿数263189模数/mm1818-18压力角/（°）22.522.522.5螺旋角/（°）000齿宽/mm460449-0.144变位系数0.50.209 3-0.378 73.1　均载系数的定义柔性销在行星传动中最主要的作用是改善行星轮间的载荷分布以及齿面载荷分布。本文使用均载系数Kγ衡量行星轮间载荷分布的均匀程度，而齿面载荷分布通过啮合齿轮的接触应力云图直观地表示。均载系数Kγ定义为行星传动中最大啮合力与平均啮合力的比值，即Kγ=(Fspi)maxFspi¯=(Fspi)max∑(Fspi/n) （20）Kγ的取值范围为1≤Kγ≤n，且数值越大，载荷分布越不均匀。根据行星传动系统的均载机制分析，减小等效啮合刚度可以有效地降低Kγ的值，即在一定条件下，系统的等效啮合刚度越低，行星轮间的载荷分布越均匀。而对于齿面载荷分布，决定其分布情况的重要依据是初始啮合螺旋线误差，这个误差是啮合齿轮沿齿宽的综合变形、位移和制造误差的合成量，通过改变制造啮合错位来研究齿面载荷分布情况。3.2　支承刚度对改善齿间载荷分布的影响在风电齿轮箱行星传动结构中，等效啮合刚度的大小对行星轮间的均载性能影响很大。使用支承刚度较低的行星销轴可以改变系统的等效啮合刚度。本文通过修改变截面销轴最小截面的尺寸d来体现变截面柔性销支承刚度的改变，销轴的最小截面尺寸越小，其支承刚度越低。在Romax软件中将变截面柔性销的最小截面尺寸分别修改为180、170、160、150、140、130 mm，根据太阳轮与行星轮外啮合的啮合力求得系统均载系数Kγ随最小截面尺寸d的变化规律如图11所示。由图11可知，随着变截面柔性销最小截面尺寸的减小，柔性销支承刚度降低，行星传动的均载系数Kγ先减小后增大。这是由于减小支承刚度可以增加销轴受载后产生的弹性变形，从而使得各构件在一定程度内自由浮动，提高风电齿轮箱行星传动机构的均载性能。但是，随着销轴支承刚度的进一步降低，变截面柔性销的变形程度过大，行星轮齿面载荷分布变得不均匀，系统的均载性能会变差。10.16578/j.issn.1004.2539.2023.10.017.F011图11柔性销轴支承刚度对均载系数的影响Fig. 11Influence of the support stiffness on the load-sharing coefficient因此，在设计变截面柔性销时，可以通过改变最小截面尺寸来改变行星轮间的均载性能。与等截面柔性销通过修改销轴长度的方式相比，变截面柔性销的使用可以缩短齿轮箱的轴向尺寸，减轻齿轮箱质量。但是，变截面柔性销最小截面的尺寸并不是越小越好，还应该考虑齿轮的承载能力、销轴的强度以及系统振动等因素。3.3　啮合错位对改善齿面载荷分布的影响制造加工误差以及构件变形等因素会引起啮合错位，啮合错位会使得齿面载荷发生偏载，即影响齿轮的齿面载荷分布[17]。本文通过改变制造加工引起的啮合错位来探究齿面载荷分布情况，并对比两种不同的行星支承，即变截面柔性销结构与普通行星销轴结构。在Romax软件中将制造加工引起的啮合误差f分别设置为0、20、40、60 μm，并给定相同的齿轮修形条件，分析行星轮1与太阳轮的接触应力，如图12所示。图12风电齿轮箱一级行星齿轮接触应力云图Fig. 12Nephogram of contact stress for the first stage planetary gear10.16578/j.issn.1004.2539.2023.10.017.F12a1（a）行星轮1变截面柔性销系统——啮合误差0 μm10.16578/j.issn.1004.2539.2023.10.017.F12a2（b）行星轮1变截面柔性销系统——啮合误差20 μm10.16578/j.issn.1004.2539.2023.10.017.F12a3（c）行星轮1变截面柔性销系统——啮合误差40 μm10.16578/j.issn.1004.2539.2023.10.017.F12a4（d）行星轮1变截面柔性销系统——啮合误差60 μm10.16578/j.issn.1004.2539.2023.10.017.F12a5（e）太阳轮变截面柔性销系统——啮合误差0 μm10.16578/j.issn.1004.2539.2023.10.017.F12a6（f）太阳轮变截面柔性销系统——啮合误差20 μm10.16578/j.issn.1004.2539.2023.10.017.F12a7（g）太阳轮变截面柔性销系统——啮合误差40 μm10.16578/j.issn.1004.2539.2023.10.017.F12a8（h）太阳轮变截面柔性销系统——啮合误差60 μm 由图12可以发现，随着啮合误差的增大，接触斑点沿着齿宽往一侧偏移，接触区域面积不断减小。对比相同啮合误差作用下两种支承结构的行星齿轮接触应力云图，可以直观地发现，使用变截面柔性销结构可以有效改善齿面载荷分布情况。根据变截面柔性销的变形特性可知，其受载后可以利用自身的柔性浮动来调整，在变形后依然与行星轴初始轴线保持相对平行，使得行星齿轮啮合后的倾斜度减小，从而改善啮合齿轮的齿面载荷分布。4 结论1）推导了变截面柔性销的挠曲线微分方程，既考虑了各个截面惯性矩的不同，也考虑了剪切效应对变截面柔性销弯曲变形的影响。同时，基于格林函数法求解了其挠度计算方程和转角计算方程。本文方法计算结果与有限元分析的误差为0.88%，变截面柔性销变形计算方法为变截面柔性销的设计提供了理论计算参考。2）变截面柔性销的变形特性使其在受载后的变形曲线呈“Z”字形，左右两端轴线趋向平行。利用变截面柔性销的变形特性可以改善行星轮间的载荷分布以及齿面载荷分布情况。3）随着变截面柔性销支承刚度的减小，风电齿轮箱一级行星传动系统的均载系数先减小后增大。因此，在设计变截面柔性销时，不能一味地减小支承刚度，需综合考虑齿轮承载能力、销轴强度条件及系统振动等因素，才能更好地改善系统的均载性能。4）随着啮合误差的增加，行星轮的接触斑点逐渐沿着齿宽往一侧偏移，接触面积逐渐减小，啮合齿轮的齿面载荷分布越不均匀。对比两种行星销轴结构的接触应力云图，采用变截面柔性销结构系统可以改善因制造加工引起的啮合错位对齿面载荷分布的影响，从而可以降低齿轮修形难度。