0 引言我国在“十四五”期间风电装机容量预期新增250 GW[1]38，发展前景巨大，然而随着国家风电补贴政策的全面取消，对风电机组的度电成本与可靠性提出了严苛要求。风电齿轮箱是风电机组传递力与运动的关键传动装置，成本约占整机的16%，是风电机组提质增效的关键环节[1]38-39。风电齿轮箱设计寿命通常为20年（陆地）或25年（海上），全寿命周期内平均风速与湍流强度的随机变化会使风电齿轮箱面临频繁的输入载荷波动，容易造成多级齿轮早期接触疲劳失效，实际服役寿命远低于设计值[2]，如图1所示。因此，开展考虑环境参数的风电齿轮箱传动系统疲劳性能优化，对于提高其对环境工况的适应性、保障风电机组安全运行具有重要意义。10.16578/j.issn.1004.2539.2023.08.008.F001图1风电齿轮箱内齿圈齿面剥落Fig. 1Early fatigue peeling of ring gears in the wind turbine gearbox近年来，国内外学者针对风电齿轮箱传动系统的疲劳性能分析及优化开展了大量研究。建立高保真的齿轮传动系统动力学模型，是准确计算齿面动载荷、实现疲劳损伤分析及性能优化的重要前提。Helsen等[3]建立了3种保真度风电齿轮箱传动系统动力学模型，即纯扭转模型、刚性6自由度模型和柔性多体动力学模型，分析了结构柔性化程度对系统振动特性的影响。Guo等[4]对比不同保真度仿真模型的求解精度与效率，给出了风电齿轮箱传动系统动力学建模时部件刚/柔、齿轮修形以及制造安装误差等建模参考标准。许华超等[5-6]考虑内齿圈柔性，分析了齿轮传动系统振动模态与动态响应。翟洪飞等[7]采用弯曲Timoshenko梁单元，建立了箱体子结构模型，分析了其对齿轮传动系统振动特性的影响。Tan等[8-9]将主轴和齿轮箱传动轴柔性纳入传动系统动力学模型中，分析了其动态特性。孙秋云等[10]308-315[11]分析了行星轮销轴误差对风电齿轮箱传动系统动态响应和疲劳寿命的影响。风电齿轮箱传动系统常在随机气动载荷驱动下运行，属于受迫振动，其运行状态不仅与构件柔性变形和内部误差激励等相关，还取决于风速工况。在风电机组全寿命周期内，长期风速概率分布会使风电齿轮箱传动系统出现复杂的载荷波动与结构变形，容易造成齿面载荷分布不均[12]11，加剧齿轮疲劳失效风险。Xiang等[13]分析了平均风速和湍流强度对风电齿轮箱各级齿轮短期疲劳损伤的影响，指出平均风速变化会显著改变齿轮疲劳损伤值。Jorgensen等[14]结合平均风速威布尔分布，分析了风电齿轮箱各级齿轮长期疲劳损伤情况，探讨了湍流强度对风电齿轮长期疲劳损伤的影响。Wang等[15]通过建立10 MW级风电齿轮箱传动系统动力学模型，结合平均风速威布尔分布，发现行星轮系太阳轮和定轴轮系小齿轮的长期疲劳损伤值较大。为了改善风电齿轮箱传动系统动态性能，Park等[16]147-151分析了齿轮修形参数对齿面载荷分布的影响，得到了额定工况下的最佳齿轮修形方案。Wei等[17]314-335分析了齿轮修形参数对啮合刚度与系统振动响应的影响，并通过台架试验对齿轮修形结果进行了验证。Shu等[18]797-816采用代理模型技术，获取了适用于多种稳定工况的齿轮传动系统最优轴承支承刚度。刘华朝等[19]158-163以齿轮接触应力最小为优化目标，基于载荷谱对风电齿轮箱各级齿轮修形参数进行了优化。陈岩松等[12]1-14建立了考虑不同风速工况组合的风电齿轮箱行星级均载优化模型，优化后的风电齿轮箱行星级在多种风速工况下均有良好的均载性能。虽然上述研究[12]13-14[16]151[17]334[18]815[19]163对风电齿轮箱传动系统高可靠设计提供了一定的理论指导，但相关研究多基于单一或多种恒定输入转矩工况开展齿轮传动系统减振优化设计，忽略了因长期风速概率分布造成的风电齿轮箱输入转矩“短期随机波动”和“长期概率分布”的风电特征，难以保证优化后的风电齿轮箱传动系统在全寿命周期内性能达到最佳；同时，全寿命周期内工况数量多，开展全工况仿真计算疲劳损伤耗时长，计算效率低。本文以某型5 MW级海上风电齿轮箱为研究对象，考虑随机风速与传动系统拓扑结构，建立了计入全局载荷的风电齿轮箱传动系统动力学模型，基于代理模型重构了“平均风速、湍流强度-齿轮修形参数-齿轮长期接触疲劳损伤”映射关系，建立了考虑风速概率分布的多级齿轮修形参数优化函数，对比了风电齿轮箱传动系统疲劳性能优化效果。1 风电齿轮箱传动系统运行原理如图2所示，风电齿轮箱主要由低速级、中间级和高速级组成。其中，低速级和中间级均为斜齿行星轮系，主要由太阳轮（s1、s2）、行星轮（p1、p2）和行星架（c1、c2）等组成；高速级为斜齿定轴轮系，主要由大齿轮（g1）和小齿轮（g2）组成。叶轮的气动转矩、弯矩、径向载荷和推力等气动载荷通过主轴传递到低速级行星架，进而依次通过低速级行星轮-太阳轮啮合副、低速级太阳轮-中间级行星架内，花键连接将力与运动传递到中间级行星架；同理，驱动中间级太阳轮将功率传递到高速级大齿轮；最后高速级大齿轮-小齿轮啮合副使高速级输出轴驱动发电机进行发电。10.16578/j.issn.1004.2539.2023.08.008.F002图2风电齿轮箱传动系统Fig. 2Wind power gearbox transmission system2 计入全局载荷的风电齿轮箱传动系统动力学模型图3所示为计入全局载荷的风电齿轮箱传动系统动力学模型，主要由两部分构成：OpenFast风电机组整机系统动力学模型和Simpack风电齿轮箱传动系统动力学模型。采用OpenFast建立风电机组整机系统动力学模型，主要由空气动力学模块、结构动力学模块和控制模块组成。空气动力学模块根据平均风速和湍流强度等环境参数，结合叶片气动参数和桨距角等计算作用在叶轮的转矩、弯矩等气动载荷；结构动力学模块包括叶片、塔筒、齿轮箱和发电机等，其中齿轮箱被简化为传动比；控制模块包括偏航控制、桨距角控制和发电机控制，详细模型介绍见文献[20]。本文通过改变空气动力学模块的平均风速和湍流强度，计算叶轮的转矩、弯矩等气动载荷。采用Simpack建立风电齿轮箱传动系统动力学模型，各齿轮参数如表1所示。其中，传动轴、行星架和箱体采用有限元获取全局质量与刚度矩阵，然后利用模态综合法进行缩聚建模[21]；齿轮啮合采用可以考虑修形的FE225号力元模拟，轴承支撑采用FE41号力元模拟，其中，轴承支承刚度利用Romax软件计算[10]308-309，齿轮啮合刚度则根据ISO 6336系列标准计算。10.16578/j.issn.1004.2539.2023.08.008.F003图3计入全局载荷的风电齿轮箱传动系统动力学模型Fig. 3Dynamic model of the wind turbine gearbox transmission system considering global loads10.16578/j.issn.1004.2539.2023.08.008.T001表1某型5 MW级风电齿轮箱齿轮基本设计参数Tab. 1Basic design parameters of a 5 MW wind turbine gearbox参数低速级中间级高速级内齿圈行星轮太阳轮内齿圈行星轮太阳轮大齿轮小齿轮齿数932932118472312124模数/mm2424241717171212螺旋角/（°）58.59压力角/（°）202020如图4所示，为了将OpenFast风电机组整机系统动力学模型和Simpack风电齿轮箱传动系统动力学模型进行联合仿真，本文通过编制工况批处理Matlab程序，动态修改OpenFast空气动力学模块的输入工况.inp文件，计算给定平均风速和湍流强度条件下叶轮的转矩和弯矩等6自由度气动载荷；然后利用Matlab将气动载荷文件格式转为Simpack风电齿轮箱传动系统动力学模型的载荷输入.afs文件，再通过调用宏命令.sjs文件修改风电齿轮箱传动系统动力学模型中齿轮修形参数并控制模型仿真，得到风电齿轮箱传动系统齿轮副动态啮合力；最后，计算了各级齿轮接触疲劳损伤。10.16578/j.issn.1004.2539.2023.08.008.F004图4OpenFast风电机组整机-Simpack风电齿轮箱联合仿真流程Fig. 4OpenFast wind turbine model-Simpack wind turbine gearbox model cosimulation process3 考虑风速概率分布的齿轮长期疲劳损伤计算3.1　平均风速-湍流强度联合概率分布利用二元Copula函数c（·，·）构建平均风速-湍流强度联合概率密度函数fvI（u，ϑ）[22]，其表达式为fvI(u,ϑ)=c[Fvu,FIϑ]fv(u)fI(ϑ) （1）式中，fv(u)为10 min时序风速的平均风速概率密度函数，通常服从威布尔分布，fvu=θλ(uλ)θ-1e-(uλ)θ；u为平均风速；θ为形状参数；λ为尺度参数；fIϑ为10 min时序风速变化程度的概率密度函数，通常服从伽马分布，fI(ϑ)=βαΓ(α)ϑα-1e-βϑ；ϑ为湍流强度；α为尺度参数；β为形状参数；Γ（α）为伽马函数；Fv(u)、FI(ϑ)分别为fvu和fIϑ的累积分布函数；c（·，·）为二元Copula函数，本文选用Gumbel Copula函数。3.2　齿轮修形方案齿轮修形通常包括齿廓修形和螺旋线修形两种形式。齿廓修形[图5（a）]主要由修形长度、最大修形量和修形曲线确定[23]9。齿顶最大修形量ek为ek=(4+0.04Ft/bca)±4 （2）式中，Ft为与齿轮分度圆相切的圆周力；bca为有效齿宽。采用长修形方式，可得修形长度为λlong=Pb(εα-1) （3）式中，λlong为修形长度；εα、Pb分别为齿轮端面重合度和基圆节距。连接修形起始点到最大修形量ek的曲线为修形曲线，如式（4）所示。e=ek（xl）b （4）式中，l为修形起始点与单对齿啮合上界点之间的距离；x为啮合位置相对坐标值；b为指数，取b=2。螺旋线修形采用修鼓形方式，以齿宽中线为对称中心，采用二次曲线进行修形[23]11-12，齿宽中间修形量最小，齿宽两端修形量最大，呈对称分布，如图5（b）所示。螺旋线修鼓量δ表达式为δ=0.25bca×10-3+0.5fg （5）fg=A(0.1bca+10) （6）式中，fg为齿轮螺旋线误差；A为系数。图5齿轮修形Fig. 5Gear modification10.16578/j.issn.1004.2539.2023.08.008.F5a1（a）齿廓修形10.16578/j.issn.1004.2539.2023.08.008.F5a2（b）螺旋线修形3.3　基于代理模型的齿轮长期接触疲劳损伤计算采用拉丁超立方抽样法与最大差异法[24]抽取关键环境参数组合和齿形设计参数组合，作为第2节计入全局载荷的风电齿轮箱传动系统动力学模型的参数输入，可以计算出任意环境工况(u，ϑ)与齿形设计参数（ek，λ，δ）组合下的齿轮副时序动态啮合力，并利用疲劳损伤线性累加，计算齿轮接触疲劳损伤[25]，进而可得在给定u和ϑ时具有确定齿形设计参数的齿轮短期接触疲劳损伤，为Dc（u,ϑ,eks1,λs1,δs1,…,ekg2,λg2,δg2）=∑ni(u,ϑ,eks1,λs1,δs1,…,ekg2,λg2,δg2)Cσi-m=NbC∫0∞σimf(σ)dσ （7）式中，σi为第i个应力块幅值；ni为应力水平σi作用在齿轮的实际次数；C、m均为齿轮接触疲劳S-N曲线参数；Nb为应力块数量；f(σ)为应力块分布函数。在风电齿轮箱全寿命周期内，考虑u和ϑ出现的概率及持续时间，基于式（1）可得齿轮长期接触疲劳损伤为Dtol（eks1,λs1,δs1,…,ekg2,λg2,δg2）=∫u1u2∫ϑ1ϑ2TfvI（u,ϑ）Dc（u,ϑ,eks1,λs1,δs1,…,ekg2,λg2,δg2）dudϑ=∑u1u2∑ϑ1ϑ2TfvI（u,ϑ）Dc（u,ϑ,eks1,λs1,δs1,…,ekg2,λg2,δg2）ΔuΔϑ （8）式中，T为风电齿轮箱设计寿命。由于式（7）中同时存在环境参数随机变量（u˜，ϑ˜）与齿轮修形区间变量（pc¯）（⋅~表示随机变量，⋅¯表示区间变量），难以直接求解，因此，本文基于代理模型技术分别在“内层”和“外层”构建“齿轮修形变量（pc¯）-齿轮短期接触疲劳损伤Dc（u，ϑ，pc¯）”“随机环境变量（u˜，ϑ˜）-齿轮短期接触疲劳损伤Dc（u˜，ϑ˜，pc）”近似函数，具体步骤（图6）为：10.16578/j.issn.1004.2539.2023.08.008.F006图6基于代理模型的齿轮长期接触疲劳损伤计算流程Fig. 6Calculation process of gear long-term contact fatigue damage based on a surrogate modelStep 1：初始化确定随机环境变量（u˜，ϑ˜）和齿轮修形区间变量（pc¯）范围。Step 2：在“外层”采用拉丁超立方抽样法和最大差异法对u˜和ϑ˜进行抽样，生成环境参数组合工况样本矩阵X2×M。其中，M为外层样本数量。Step 3：计算齿轮修形参数（eks1，λs1，δs1，…，ekg2，λg2，δg2）分别在M组环境参数组合工况下的齿轮短期接触疲劳损伤Dc（u，ϑ，pc）1×M，具体计算过程如下：① 针对X2×M中每个样本，在“内层”采用拉丁超立方抽样法和最大差异法对[pc¯]抽样，得到齿轮修形变量样本矩阵Xdim(pc¯)×N。其中，N为内层样本数量；dim(·)为获取维度。② 调用风电齿轮箱传动系统动力学模型（图3），计算“外层”环境参数组合工况样本和“内层”齿轮修形变量样本组合下的齿轮副动态啮合力，并根据式（7）计算Dc(u，ϑ，pc)|N×M。③ 将②中第i组环境参数组合工况对应的Dc（u，ϑ，pc）|N×i分别利用克里金模型[26]、多项式响应面模型[27]和径向基函数[28]构建“齿轮修形变量-齿轮短期接触疲劳损伤”近似函数Dc(u，ϑ， pc¯)|N↔×i。然后采用均方根误差（Root Mean Square Error，RMSE）检测拟合精度[29]，选择最佳代理模型。其中，·↔表示对该变量进行拟合。最后通过近似函数Dc（u，ϑ，pc¯)|N⃡×i分别计算在M组环境参数工况下齿轮修形参数（eks1，λs1，δs1，…，ekg2，λg2，δg2）对应的Dc(u，ϑ，pc)|1×M。DcKRG（u,ϑ,pc¯）=∑i=1nκipi[eks1,λs1,…,δg2,(u,ϑ)|f]+Z[eks1,λs1,…,δg2,(u,ϑ)|f]#DcPRS（u,ϑ,pc¯）=(β0+∑i=1nβi（eks1i+…+δg2i,(u,ϑ)|f）+…+∑i=1n…∑j=1nβi…j（βi…jeks1i…δg2i,(u,ϑ)|f）#DcRBF（u,ϑ,pc¯）=∑i=1Ns[wiϕ(||（eks1,λs1,…,δg2,(u,ϑ)|f）-(eks1i,λs1i,…,δg2i,(u,ϑ)|f)||)] （9）式中，|f为第f组环境参数变量组合；∑i=1nκipi(·)为多项式全局近似模型；Z(·)为均值为0、方差为σ2、协方差非零的局部偏差模型；β0、βi、βii、βij均为多项式函数待定系数；||·||为训练样本之间的欧氏距离；w为权重系数；ϕ(·)为径向基函数；Ns为训练样本个数。则eRMSE=∑i=1Nt[Dci（u,ϑ,pc¯）-Dcis(u,ϑ,pc¯)]2Nt （10）式中，Nt为测试样本数目；Dci（u，ϑ， pc¯）为第i个测试样本中齿轮短期接触疲劳损伤真实值；Dcis（u，ϑ，pc¯）为3种代理模型技术预测第i个测试样本齿轮短期接触疲劳损伤。Step 4：获得Dc(u，ϑ，pc)|1×M后，同理，分别使用克里金模型、多项式响应面模型和径向基函数构建“关键环境参数组合工况-短期接触疲劳损伤”近似函数Dc（u˜，ϑ˜，pc），检测拟合精度，选择最佳代理模型。Step 5：结合Dc（u，ϑ， pc¯）|N⃡×i和Dc（u˜，ϑ˜，pc）计算齿轮长期接触疲劳损伤Dtol（pc¯）。4 考虑风速概率分布的齿轮接触疲劳性能优化为了得到Dtol，需要先计算Dc，其与平均风速u、湍流强度ϑ和齿轮修形变量（eks1，λs1，δs1，…，ekg2，λg2，δg2）相关。由于齿轮修形变量较多，在开展齿轮疲劳性能优化之前，采用Plackett-Burman试验设计方法筛选显著性齿轮修形参数。利用线性函数拟合各齿轮修形参数与齿轮长期接触疲劳损伤的对应关系[30]，有Y=β0+∑i=1kβixi （11）式中，β0为回归函数的常系数项；Y为目标响应；βi为回归系数；xi为第i个修形设计变量。令“+”表示每个修形设计变量的高水平，“-”表示每个修形设计变量的低水平，则第i个修形设计变量对目标响应的影响程度为E(xi)=2[∑i=1kY(xi+)-∑i=1kY(xi-)]Ni （12）式中，E(xi)为第i个修形设计变量的影响程度；Y(xi+)、Y(xi-)分别为Plackett-Burman设计中测量到的响应变量Y的最大值与最小值；Ni-、Ni+和Ni分别为对应于第i个修形设计变量的低水平试验次数、高水平试验次数和总试验次数，Ni=Ni++Ni-。通过Plackett-Burman试验设计，可得对容易失效齿轮长期接触疲劳损伤影响最为明显的齿轮修形变量组合pcf=（ekh，λh，δh，…，eks，λs，δs）。因此，可将式（8）改写为Dtolf(eks1,λs1,δs1,…,ekg2,λg2,δg2)≈∑u1u2∑ϑ1ϑ2TfvIu,ϑDc(u,ϑ,pcf)ΔuΔϑ （13）为了获取最优的齿轮修形变量组合，本文以齿轮修形变量的区间范围、未修形时容易失效齿轮长期接触疲劳损伤（Dtol1，…，Dtolj）的极大值（D1，…，Dj）为约束条件，以未修形时最容易失效齿轮长期接触疲劳损伤（Dtolf）最小为优化目标，建立齿轮修形参数优化模型，如式（14）所示，最后采用遗传算法进行求解[31]。find x=[ekh,λh,δh,…,eks,λs,δs]min Dtolf(u˜,ϑ˜,x)s.t.  Dtol1(u˜,ϑ˜,x)≤D1,…,Dtolj(u˜,ϑ˜,x)≤Djekh∈[ekh̲,ekh¯],λh∈[λh̲,λh¯],δh∈[δh̲,δh¯]⋮eks∈[eks̲,eks¯],λs∈[λs̲,λs¯],δs∈[δs̲,δs¯] （14）式中，∙̲、·＿分别表示变量的极小值与极大值。5 结果讨论与分析本文以某型5 MW级风电齿轮箱为例，各级齿轮材料及热处理如表2所示。首先，根据国内某海上风电场历年风速统计数据，建立平均风速-湍流强度联合概率密度函数，选取关键环境参数组合工况；其次，利用风电齿轮箱传动系统动力学模型，计算未修形时齿轮长期接触疲劳损伤，找出容易发生接触疲劳失效的齿轮；然后，筛选出对容易失效齿轮长期接触疲劳损伤影响较为明显的齿轮修形变量，并以容易失效齿轮长期接触疲劳损伤值最小求解最优的齿轮修形参数；最后，对比优化前、后各级齿轮接触疲劳性能改善效果。10.16578/j.issn.1004.2539.2023.08.008.T002表2风电齿轮箱各级齿轮材料及热处理Tab. 2Materials and heat treatment of multi-stage gears in the wind turbine gearbox位置齿轮名称材料热处理精度等级低速级太阳轮18CrNiMo7-6渗碳5行星轮18CrNiMo7-6渗碳5内齿圈20CrMnMo调质5中间级太阳轮18CrNiMo7-6渗碳5行星轮18CrNiMo7-6渗碳5内齿圈42CrMoA氮化5高速级大齿轮20CrMnMo调质5小齿轮18CrNiMo7-6渗碳55.1　齿轮长期疲劳损伤分析图7所示为未修形时风电齿轮箱各级齿轮长期接触疲劳损伤值。由图7可以看出，低速级太阳轮（s1）、中间级太阳轮（s2）和高速级小齿轮（g2）损伤值均大于0.55，容易发生接触疲劳失效。其中，s1损伤值大于0.7，最容易发生接触疲劳失效。因此，本文选用s1、s2和g2的长期接触疲劳损伤为目标响应（Dtols1、Dtols2和Dtolg2），采用Plackett-Burman试验设计，结合回归分析和T检验，筛选出对其影响较为明显的各级齿轮修形变量（eks1，λs1，δs1，…，ekg2，λg2，δg2，共计24个设计因子）。10.16578/j.issn.1004.2539.2023.08.008.F007图7风电齿轮长期接触疲劳损伤分布情况Fig. 7Long-term contact fatigue damage distribution of wind turbine gears表3所示为Plackett-Burman试验设计对应的仿真结果。在此基础上，表4~表6分别为24个设计因子对Dtols1、Dtols2和Dtolg2的显著性影响水平（当显著性概率0.05时，则认为该设计变量为显著变量，用**表示；当显著性概率≥0.05时，则认为该设计变量为不显著变量，用NS表示）。综合表4~表6可知，s1、s2和g2的修形参数显著影响Dtols1、Dtols2和Dtolg2，故将其作为设计变量开展后续优化。10.16578/j.issn.1004.2539.2023.08.008.T003表3Plackett-Burman试验设计仿真结果Tab. 3Simulation results of Plackett-Burman experiment design序号eks1/μmλs1/mmδs1/μm…δg2/μmDtols1Dtols2Dtolg2179.9154.9524.95…49.890.7300.6280.5492239.7354.9524.95…16.630.7040.6260.537379.9118.3224.95…16.630.7480.6140.538……………49.89………4879.9154.9574.95…49.890.7270.6240.53210.16578/j.issn.1004.2539.2023.08.008.T004表4低速级太阳轮长期接触疲劳损伤显著性分析Tab. 4Significance analysis of long-term contact fatigue damage of low-speed stage sun gears序号设计变量回归系数T检验统计量显著性概率显著性1eks1-119.8-2.4670.016 7**2λs1-0.342-1.3710.176NS3δs1297.1593.7530.000**4ekp1-233.047-4.0660.000**………………24δg2-71.426-0.2560.799NS10.16578/j.issn.1004.2539.2023.08.008.T005表5中间级太阳轮长期接触疲劳损伤显著性分析Tab. 5Significance analysis of long-term contact fatigue damage of intermediate stage sun gears序号设计变量回归系数T检验统计量显著性概率显著性1eks133.2910.5870.560NS2λs11.7292.8600.005**………………10eks2-178.255-3.0790.003**11λs23.2978.2500.000**12δs2125.3861.6410.106**13ekp2-198.246-2.9480.005**………………21δg11 788.146.3670.000**………………24δg21 601.4405.9160.000**10.16578/j.issn.1004.2539.2023.08.008.T006表6高速级小齿轮长期接触疲劳损伤显著性分析Tab. 6Significance analysis of long-term contact fatigue damage of high-speed stage pinion gears序号设计变量回归系数T检验统计量显著性概率显著性1eks118.5410.3490.728NS2λs1-0.601-2.6160.011**………………12δs2-2.481-6.6410.000**………………21δg1 -1 459.22-5.5580.000**22ekg2-27.009-0.3630.718NS23λg2-1.189-2.1020.040**24δg2-1 403.890-5.5470.000**5.2　优化效果对比表7所示为低速级太阳轮（s1）、中间级太阳轮（s2）和高速级小齿轮（g2）的修形参数取值范围和最佳修形量。其中，根据式（2）~式（6）计算得到齿轮的理论修形量，然后分别选取其150%和50%作为齿轮修形参数区间范围的上、下边界。由表7可以看出，求解的最佳齿轮修形量均离上、下边界值较远，这表明其不易受加工、装配等误差影响，适用性较好。10.16578/j.issn.1004.2539.2023.08.008.T007表7最佳齿轮修形量Tab. 7Optimum gear modification parameters齿轮符号区间范围优化值低速级太阳轮eks1［79.91，239.73］180.00λs1［18.32，54.95］38.40δs1［24.95，74.95］27.00续表10.16578/j.issn.1004.2539.2023.08.008.T008齿轮符号区间范围优化值中间级太阳轮eks2［66.66，199.99］98.00λs2［12.66，37.98］20.59δs2［17.82，53.46］52.00高速级小齿轮ekg2［57.06，171.19］127.00λg2［10.12，30.36］15.13δg2［16.63，49.89］23.00为了尽可能地反映各风速段对风电齿轮箱传动系统动载荷分布的影响差异，分别选取平均风速为3 m/s（额定风速以下）和15 m/s（额定风速以上）时观察低速级和中间级太阳轮-行星轮啮合副、低速级和中间级行星轮-内齿圈啮合副、高速级大齿轮-小齿轮啮合副的齿面载荷分布情况，如图8所示。由图8可以看出，未修形的齿轮副存在严重的偏载现象，并且在不同风速工况下其偏载程度存在一定的差异性，尤其当平均风速为15 m/s时，低速级太阳轮-行星轮啮合副的齿面最大载荷达1.82×105 N，出现了较为明显的周期性边缘接触，受载情况较为恶劣。图8优化前不同风速工况下风电齿轮齿面载荷分布Fig. 8Load distribution of wind turbine gear tooth surface under different wind speed conditions before optimization10.16578/j.issn.1004.2539.2023.08.008.F8a1（a）平均风速3 m/s，湍流强度0.1110.16578/j.issn.1004.2539.2023.08.008.F8a2（b）平均风速15 m/s，湍流强度0.06 图9所示为优化后的风电齿轮箱各级齿轮齿面载荷分布。由图9可以看出，齿轮修形可以有效地改善齿面偏载现象，使最大齿面载荷向齿面中间移动，提高其承载能力，与文献[32]结论类似。但由于追求齿轮在各工况疲劳损伤加权最小，中间级太阳轮-行星轮啮合副在部分风速工况下存在齿面偏载，但相对于优化前，优化后的齿面载荷更趋近于中间位置，同时其他齿轮啮合副在不同风速工况下齿面载荷主要集中在齿面中间位置，表明齿轮修形参数可以较好地适应风速工况变化。当平均风速为15 m/s时，低速级太阳轮-行星轮啮合副的齿面最大载荷为1.71×105 N，相对于优化前降低6.04%。图9优化后不同风速工况下风电齿轮齿面载荷分布Fig. 9Load distribution of wind turbine gear tooth surface under different wind speed conditions after optimization10.16578/j.issn.1004.2539.2023.08.008.F9a110.16578/j.issn.1004.2539.2023.08.008.F9a210.16578/j.issn.1004.2539.2023.08.008.F9a310.16578/j.issn.1004.2539.2023.08.008.F9a4（a）平均风速3 m/s，湍流强度0.1110.16578/j.issn.1004.2539.2023.08.008.F9a510.16578/j.issn.1004.2539.2023.08.008.F9a610.16578/j.issn.1004.2539.2023.08.008.F9a710.16578/j.issn.1004.2539.2023.08.008.F9a8（b）平均风速15 m/s，湍流强度0.06 表8所示为优化前后风电齿轮箱各级齿轮长期接触疲劳损伤对比。由表8可以看出，优化后的低速级太阳轮、低速级行星轮、中间级太阳轮、高速级大齿轮和高速级小齿轮的长期接触疲劳损伤值均出现了明显降低，尤其是高速级大齿轮和小齿轮，最大降幅达37.42%；最容易失效的低速级太阳轮长期接触疲劳损伤值降低了11.37%，有效地提高了风电齿轮箱传动系统疲劳性能。同时，由于本文仅对低速级太阳轮、中间级太阳轮和高速级小齿轮进行修形，优化后低速级内齿圈、中间级行星轮和内齿圈的长期接触疲劳损伤略微增加，但增加后仍较小。10.16578/j.issn.1004.2539.2023.08.008.T009表8优化前后各级齿轮长期接触疲劳损伤对比Tab. 8Comparison of long-term contact fatigue damage of each stage gears before and after optimization齿轮疲劳损伤值降幅百分比/%优化前优化后低速级太阳轮0.734 50.651 0-11.37低速级行星轮0.511 30.484 3-5.28低速级内齿圈0.190 30.197 73.88中间级太阳轮0.634 40.617 3-2.7中间级行星轮0.507 30.544 97.41中间级内齿圈0.381 40.442 215.94高速级大齿轮0.524 20.332 6-36.55高速级小齿轮0.568 80.356 0-37.426 结论考虑长期风速概率分布特征，建立计入全局载荷的风电齿轮箱传动系统动力学模型，并通过代理模型重构“齿轮修形参数-齿轮长期接触疲劳损伤”映射关系，优化了多级齿轮修形参数，并对比了优化效果。主要结论：1）在风电齿轮箱全寿命周期内，低速级太阳轮、中间级太阳轮和高速级小齿轮容易发生接触疲劳失效，尤其低速级太阳轮长期接触疲劳损伤值大于0.7，是风电齿轮箱传动系统高可靠设计的薄弱环节之一。2）优化前的风电齿轮箱各级齿轮的齿面载荷会出现明显的偏载现象，在不同风速工况下其偏载程度存在一定的差异性，容易出现周期性边缘接触；优化后的风电齿轮箱各级齿轮的齿面载荷分布在不同风速工况下仍主要集中在齿面中间位置，工况适应性强。3）优化后的风电齿轮箱各级齿轮长期接触疲劳损伤值出现了大幅降低，尤其是高速级大齿轮和小齿轮，最大降幅达37.42%。同时，最容易失效的低速级太阳轮长期接触疲劳损伤值降低了11.37%，疲劳性能提升效果显著。本文提出的考虑环境参数的风电齿轮箱传动系统疲劳性能优化方法后续将考虑应用于风电齿轮箱设计中，并开展相关试验进行验证。