SOLARENERGY技术产品与工程太阳能一种风电机组轴承简化建模方法的研究特变电工新疆新能源股份有限公司■黄文摘要：针对风电机组变桨轴承及偏航轴承结构与受力的复杂性，提出了一种新的滚动轴承简化建模方法。通过利用有限元分析软件MSC.Marc/Mentat中的GAP单元模拟滚动体与内、外圈之间的接触来实现对轴承建模的简化，同时对简化模型进行加载计算，得到内、外圈的应力及应变结果；并用Hertz理论及Newton迭代法计算同种工况下内、外圈的应力值及位移值，来验证此种建模方法的可行性和正确性。结果表明：滚动轴承的内、外圈接触可采用GAP单元进行模拟，从而为风电机组轴承有限元建模提供了一种简便且有效的方法。关键词：滚动轴承；GAP单元；MSC.Marc/Mentat；有限元0引言滚动轴承的应用非常广泛，不仅涉及到航天、船舶等领域，也涉及到风力发电设备[1]。对滚动轴承的分析起初采用的是动力学方法，对轴承的整个运动过程建立动力学方程，求得任意瞬间轴承各部件的运动参数和接触应力。目前，轴承的分析主要采用计算机建模，通过仿真的形式分析轴承的各相关性能指标[2]，但滚动轴承的仿真建模主要还是以实体建模为主，所需计算量非常大。因此，本文针对风电机组上的变桨轴承与偏航轴承，采用GAP单元简化轴承建模的方法，并用理论方法验证其正确性。结果表明，这种简化方法可有效模拟滚动轴承受力特性，对结构尺寸较大的轴承建模可起到一定的借鉴作用。递形式是直线形式，如图1所示。图1荷载传递路线根据荷载的传递路径特征，在对轴承进行简化建模时，需充分考虑MSC.Marc/Mentat中的GAP单元特性，通过对单元特性的设置，完成最终模型的建立。GAP单元支持MSC.Marc/Mentat软件中的任意其他单元连接，分为左右2个连接点，分别连接在模型的不同节点上，用以描述模型中任意两1采用MSC.Marc/Mentat软件中的GAP单元对轴承滚动体进行分析四点接触轴承有3个关键点，分别是滚动体中心点、滚动体与外圈接触点、滚动体与内圈接触点。四点接触轴承在受力的过程中，其荷载传收稿日期：2017-01-18通信作者：黄文(1983—)，男，硕士，主要从事风力发电机组研发方面的工作。88177347@qq.com55SOLAR ENERGY 10/2018Copyright©博看网 www.bookan.com.cn. All Rights Reserved.SOLARENERGY太阳能技术产品与工程节点的相对关系(包含摩擦、滑动或热影响等)[3]。入，GAP单元的法向刚度将变得足够大；而在未接触区，它对分析对象的运动状态不产生影响[5]。当接触点对的相对位移大于初始间隙时，法向力Fn的值就会是零，而接触点对就会处于分离状态，这时GAP单元不会发生力的传递，从而不会影响分析对象的运动状态；而当相对位移小于初始间隙时，表明接触点对处于接触状态，此时GAP单元就像一个线性弹簧，其法向接触刚度为KA，在GAP单元的法线方向就会存在一个法向力Fn(负值)，内外圈就会通过GAP单元在B、C两点之间传递荷载。1.1实体滚动轴承参数及三维模型采用三维软件对滚动轴承(以某一单排四点接触球轴承为例)进行建模，包含轴承内圈、轴承外圈、滚动体(模型中省略保持架，在后期的有限元分析软件中采用一个微小弹簧代替保持架)来实现滚动体的固定。此模型建立比较简单，具体如图2所示。轴承各参数如表1所示。1.2.2GAP单元对轴承的简化建模轴承三维模型建好之后，要对模型进行网格划分。网格划分是否合理将关系到后期模型计算的准确性和效率，所以网格要过渡合理、均匀；同时网格类型的选择也至关重要，因网格类型有数值2.07×105840.78539726.5150.801205.131010.32.200.77781187.601224.40图2四点接触球轴承三维模型表1四点接触球轴承参数表参数弹性模量E/MPa滚动体数量Z初始接触角α0/(°)滚道圆弧半径ri/mm滚动体直径Dw/mm滚动体回转直径dm/mm轴承列数n柏松比μ加载前两对角曲率中心距离A/mm曲率中心参数a/mm轴承内圈沟底直径F/mm轴承外圈沟底直径F′/mm四面体、六面体、八面体等，每种类型适用于不同的模型，此处是对单个轴承进行分析，为了达到使用较少的网格重划分次数而得到较高计算精度的目的，最终选择六面体网格。因轴承整圈具有左右对称性，只取半个轴承模型进行分析，滚动体也只取半个进行实体建模(只考虑此滚动体的应力应变情况，与理论计算进行对比)，其他的滚动体只考虑力的传递，不再具体分析每个滚动体的应力应变，在这里将其他滚动体用“GAP单元”来替代定义。刚度是滚动轴承的一项重要指标，而GAP单元就像一个线性弹簧，将其长度设置为滚动体直径后，在软件中再将其初始轴向刚度和径向刚度设置为零(初始状态不受力)，一旦轴承承受外界作用力，GAP单元可出现荷载－位移关系，轴向刚度和径向刚度就会出现一定的数值，接触刚度这时就可近似模拟滚动体刚度[6]。接触角是滚动轴承的一个主要设计参数，在分析轴承受力、变形、运动关系和确定轴承承载能力时，都需要首先确定接触角的大小。GAP单元由2个节点连接组成，一定数量的GAP单元561.2滚动轴承的简化建模分析1.2.1GAP单元的特性假设B点代表外圈接触点，C点代表内圈接触点，那么B、C两点可组成GAP单元，B、C两点之间的距离可设定为滚动体的直径Dw。当轴承受外界作用力时，GAP单元可出现荷载－位移关系，其关系可通过接触面的切向分量和法向分量来表述[4]。在接触区域，为了防止接触体相互侵SOLAR ENERGY 10/2018Copyright©博看网 www.bookan.com.cn. All Rights Reserved.SOLARENERGY技术产品与工程太阳能可以近似代替一个滚动体，也就是在建立之初就已经指定了力的传递点，所以GAP单元无法精确描述接触角的具体数值。这里针对接触角不做具体分析，但GAP单元越多，越能更准确地反映滚动体和滚道之间的接触情况。GAP单元可以模拟任意直径滚动体建模，只需在软件中将GAP单元长度设置为滚动体直径即可。所谓一定数量的GAP单元代替一个滚动体，并不是说GAP单元越少越好，也不是越多越好。GAP单元太少时，不足以模拟滚动体与滚道的接触状态；而GAP单元太多时，就会增加软件计算时间，降低效率。笔者前期已做过大量实验验证，为了使简化模型更好地与理论吻合，GAP单元需要均匀布置到接触角范围内，相邻GAP单元之间的角度差应控制在1°～3°之间。本模型采用14对GAP单元代替1个滚动体，相邻GAP单元之间的角度差约为2.14°(本模型轴承接触角为30°)。网格划分单元的大小根据位置有不同的取值，其在内外圈接触的地方大小约为2mm，其他地方单元格可适当放大，如图3所示。1.2.4接触设置接触设置主要是滚动体与内、外圈之间的状态设定。首先，针对GAP单元，根据软件中的要求，不需要对其进行接触设置，只需将其与内、外圈节点连接即可；其次，对半个滚动体与内、外圈的接触设置，这里设置为接触，即“touching”状态；最后，对其他参数进行设置，因为内圈、外圈和滚动体都假定采用轴承钢，材料为GCr15，所以将它们设置为可变形体，相应的摩擦系数设置为0.15。1.2.5边界条件、荷载的设置1)边界条件的设置：滚动轴承受力时，其运动状态是不固定的，可以在轴承外圈上设置约束条件，同时对轴承内圈施加荷载；也可以在轴承内圈上设置约束条件，同时对轴承外圈施加荷载。在这里采取对外圈进行固定的方式，即对外圈采用全约束方式，同时约束轴承剖面沿法向的平动自由度，以及约束内圈沿轴向的转动自由度。2)荷载的设置：以上边界条件设置好之后，开始在轴承内圈上加载来计算轴承接触部位各应力变化。根据MSC.Marc/Mentat软件对MPC点的介绍，这种点可以作为一个荷载施加中心，通过与模型相连，将荷载传递到模型上。此例中，将MPC点设置在轴承中心点处，并且将此点与内圈上表面所有点相连，荷载的传递路径就是通过这些连线进行传递。MPC点设置好之后，在其上面施加混合荷载(包含翻转力矩、径向力和轴向力)。边界条件的设置及荷载的施加如图5及表2所示。)MPC*RBE3BCGAP图3基于GAP单元的滚动轴承简化建模1.2.3保持架的模拟在MSC.Marc/Mentat软件中，选用刚度很小的弹簧单元对滚动体进行固定，小弹簧的作用类似于保持架，将各滚动体固定，并且均匀隔开。为了准确模拟保持架的作用，弹簧的一端与滚动体上节点相连，另一端设置为接地不动，如图4所示。YZX图4弹簧单元固定滚动体57SOLAR ENERGY 10/2018图5边界条件、荷载设置示意图Copyright©博看网 www.bookan.com.cn. All Rights Reserved.SOLARENERGY太阳能技术产品与工程表2荷载值大小参数翻转力矩M/N·mm翻转力矩与径向力的夹角Φ/(°)轴向力Fa/N径向力Fr/N数值1.80×10990800001963742滚动轴承理论计算理论计算采用的轴承模型与有限元计算模型一致，参数表和荷载表也采取表1、表2中的数据。2.1径向位移δr、轴向位移δa和倾斜角θ数值的确定四点接触轴承在受力时，会发生一定的位移。要通过理论方法计算轴承在荷载作用下的相关量，首先需要确定3个未知量：径向位移δr、轴向位移δa和倾斜角θ，具体如图8所示。δrθδa1.2.6滚动轴承简化建模有限元计算结果根据表2中的相关参数，可以得到滚动轴承的应力和应变结果，其中，内圈与滚动体之间的最大接触应力值约为4603.2MPa，最大位移值约为0.1511mm；外圈与滚动体之间的最大接触应力值约为4521.6MPa，最大位移值约为0.1465mm。内、外圈与滚动体接触的应力云图如图6、图7所示。/MPa4.603×1034.373×1034.143×1033.913×1033.683×1033.452×1033.222×1032.992×1032.762×1032.532×1032.302×1032.071×1031.841×1031.611×1031.381×1031.151×1039.206×1026.905×1024.603×1022.302×1020.000dm图83个未知量示意图针对以上3个未知量，需要根据文献[7-9]对其进行如下推导求解。图9是在混合荷载作用下，对轴承外圈固定，轴承内圈曲率中心轨迹的变化曲线。其中，S为受载时内、外圈滚道曲率中心之间的最短距离。图10为滚动体与滚道未ZYX加载与受载后的接触示意图，初始接触角α0的取值为π/4。ZZ′θZ′′Q′图6内圈与滚动体接触的应力云图/MPa4.522×1034.296×1034.069×1033.843×1033.617×1033.391×1033.165×1032.939×1032.713×1032.487×1032.261×1032.035×1031.809×1031.583×1031.356×1031.130×1039.043×1026.782×1024.522×1022.261×1020.000δRoSRiΨO′δrδaXX′YZXP′YY′Y′′图7外圈与滚动体接触的应力云图58图9内、外圈滚道曲率中心改变后轨迹示意图SOLAR ENERGY 10/2018Copyright©博看网 www.bookan.com.cn. All Rights Reserved.SOLARENERGY技术产品与工程太阳能图9、图10中的其他主要参数可通过下文中S=A+δi+δo(1)式中，δi为受载后内圈曲率中心位移量；δo为受载后外圈曲率中心位移量。所以，受载后内、外圈曲率中心变化总量δn为：δn=δi+δo=S-Aα0R的公式进行确定：Dw2Dw2iF2(2)F(ρi)=图11内圈与滚动体接触示意图|ρ1-ρ2|+|ρ3-ρ4|式中，ρ1为滚动体与内圈接触时的主曲率，ρ2为滚动体与外圈接触时的主曲率，且Arori∑ρi(4)为轴承内圈沟底曲率，ρ4=2。Fρ1=ρ2=2；ρ3为内圈滚道圆弧曲率，ρ3=1；ρDwri4根据式(3)、式(4)的推导，同理可得出轴承外圈辅助变量F(ρo)及外圈与滚动体之间的主曲率之和∑ρo。a.滚动体与滚道未加载时的接触示意图αQδi内圈滚道曲率中心点半径Ri为：Ri=1dm+(ri-0.5Dw)cosα02外圈滚道曲率中心点半径Ro：Ro=Ri－Acosα0(5)(6)滚动体总共有84个，假设其中某个滚动体A的中心在X轴上，设其与X轴夹角为ψ，则初始么绕其逆时针方向分布的第i个滚动体，其其他滚动体在滚道圆周上均匀分布，那值ψ0=0，δoψi的范围在[0,2π]之间，如图12所示。ψi=2πi，84ZQb.滚动体与滚道受载后的接触示意图图10滚动体与滚道未加载与受载后的接触示意图Yψi为：内圈与滚动体之间的主曲率之和∑ρi可表示i12∑ρ=ρ+ρX+ρ3+ρ4(3)图12滚动体在滚道上的分布轴承内圈辅助变量F(ρi)可表示为：59根据参考文献[9]，图9中S值可表示为：SOLAR ENERGY 10/2018Copyright©博看网 www.bookan.com.cn. All Rights Reserved.SOLARENERGY太阳能技术产品与工程S=(Asinα0+δa+Riθcosψi)+(Acosα0+δrcosψi)[2212]Ko=(7)212将式(7)带入式(2)可得：δn=(Asinα0+δa+Riθcosψi)+(Acosα0+δrcosψi)2式中，eδ和eδ可通过公式(4)计算出F(ρi)的oi(eδ)o32(∑ρ)o112(12)[]-A[9]值后，查《球轴承的设计计算》中的赫兹接触系数表得到。对每个滚动体上的力进行分解，使每个滚动体上的受力均分解为与初始径向力Fr、轴向力Fa00(8)某个滚动体上所受荷载Qψ为：iQψ=K(δnn)=Ki{[(Asinα+δ+Rθcosψ)+0aii212以及力矩M0相同方向的分力，然后对各分力在各自方向上分别求和，具体如下：(9)Fa=∑Qψsinαψi=2πψi=0ψi=0ii2(Acosα0+δrcosψi)]-A}nψi=2π(13)(14)(15)式中，K为滚动体与滚道接触的总刚度值；对于球轴承而言，n取1.5。K可由式(10)求得：éùêúêúêú1úK=ê11êêæ1önæ1önúúêú+ç÷ç÷êKúëèiøèKoøûnFr=∑Qψcosψicosαψi=2πiidM=m∑Qψcosψisinα2ψ=0(10)式中，α为轴承受载后的接触角，则有：sinα=式中，Ki为轴承内圈与滚动体之间的相对刚度；Ko为轴承外圈与滚动体之间的相对刚度。其中：Ki=(eδ)i[(Asinα0+δa+Riθcosψi)2+(Acosα0+δrcosψi)2]Acosα0+δrcosψiii0Asinα0+δa+Riθcosψi12(16)212cosα=32(∑ρ)i1[(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)]a20ri12(11)由此可推导出式(18)~式(20)：0aii20ri2121.5(17)[(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)]-A}{F-K∑ψi=2πψi=0a0aii20ψi=2πψi=00aii20[(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)][(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)]-A}×(Acosα+δcosψ)cosψ{F-K∑=0(19)[(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)][(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)]-A}×(Asinα+δ+Rθcosψ)cosψ{dM-K∑=02[(Asinα+δ+Rθcosψ)+(Acosα+δcosψ)]ri212×(Asinα0+δa+Riθcosψi)=0(18)ri2121.50riir0aii20ri212mψi=2πψi=00aii20ri2121.50aiii0aii20ri212(20)轴承所受荷载如表2所示，根据式(16)~式径向力时)、初始轴向位移δa′(仅受轴向力时)、(18)，可分别计算轴承的初始径向位移δr′(仅受计算出各初始参数后，可根据牛顿迭代法进行收敛计算，具体如下[10]：X(k+1)初始倾斜角θ′(仅受弯矩时)。其中：60SOLAR ENERGY 10/2018=X-J(Xk[k)]f(X)-1k(21)Copyright©博看网 www.bookan.com.cn. All Rights Reserved.技术产品与工程SOLARENERGY太阳能é∂f1∂f1ù...êúê∂x1∂xnúúêúJ(X)=ê⋮⋮êúê∂f∂fúnnúêê...úë∂x1∂xnû通过计算所得的轴承3个未知量(轴向位移δa、(22)径向位移δr和倾斜角θ)，结合公式(9)得到。针对以上有限元的计算结果，理论计算为了能与有限元模型更好的匹配，选择同一位置处的滚动体进行荷载理论计算，即取ψ=π，则计算所得滚动体的荷载Q=58708.58383。2.2.1轴承内圈与滚动体接触的最大应力及位移值0.061145401，F(ρi)=0.6309204。pi=1.5×πeaebmaxii根据以上计算所得的3个初始值，将式(16)~式(18)联立，采用牛顿迭代法进行求解。求解的过程是迭代收敛的过程，当最终收敛后，就会得轴到轴承在外荷载(见表2)作用下的径向位移δr、向位移δa和倾斜角θ。根据表1、表2中的数值及以上相关公式，可求得轴承其他参数值，具体如表3所示。表3计算所得参数值大小根据式(3)、式(4)及表3可知，∑ρi=内圈与滚动体接触的最大接触应力pi为：3参数轴承外圈与滚动体之间相对刚度Ko轴承内圈与滚动体之间相对刚度Ki总相对刚度K滚动体与内圈接触的主曲率ρ1/mm轴承内圈滚道圆弧曲率ρ3/mm轴承内圈沟底曲率ρ5/mm滚动体与外圈接触的主曲率ρ2/mm轴承外圈滚道圆弧曲率ρ4/mm轴承外圈沟底曲率ρ6/mm初始径向位移δr′/mm径向位移δr/mm初始轴向位移δa′/mm轴向位移δa/mm倾斜角θ/(°)初始倾斜角θ′/(°)数值1049636.0301049537.791371085.010.0490196080.049019608-0.037735849-0.0377358490.000842034-0.0008167271205.9088170.0500.1562310320.033-0.0936757380.00020.000291258[9]根据F(ρi)的值，查《球轴承的设计计算》(∑ρiQ2)max(23)中的赫兹接触系数表，并利用插值法可得：πeaeb=1.96082836×10-3、eδ=2.48944154×10-4，则pi=4614.5MPa。maxiii滚动体与内圈的最大位移δi为：将相关数值代入式(24)可得δi=0.1481mm。δi=eδ3i(∑ρ)Qi2(24)2.2.2轴承外圈与滚动体接触的最大应力及位移值参照轴承内圈与滚动体的接触计算，同理可计算出oo内圈滚道曲率中心点半径Ri/mmπeaeb=1.95156592×10-3、eδ=2.50127688×10-4，o∑ρo=0.05948664、F(ρo)=0.6206288、则轴承外圈与滚动体接触的最大接触应力值最大位移δo=0.1475mm。po=4552.2MPa，max3对比分析基于理论计算的结果与GAP单元对轴承的简化建模计算结果进行对比，具体如表4所示。最大应力/MPa4603.24521.64614.54552.2612.2接触应力、位移求解首先计算84个滚动体所受的荷载大小，可计算方式GAP建模计算理论计算接触位置内圈与滚动体外圈与滚动体内圈与滚动体外圈与滚动体表4结果对比最大位移/mm0.15110.14650.14810.1475SOLAR ENERGY 10/2018Copyright©博看网 www.bookan.com.cn. All Rights Reserved.SOLARENERGY太阳能技术产品与工程根据表4中的数值可知，GAP单元建模计算的结果与理论计算的结果误差较小，基本吻合。克服了实体轴承建模难度大、复杂性高的缺点，同时也大幅提高了计算速度。6)这种简化方法可以应用于风力发电机组中变桨轴承及偏航轴承的建模，从而为风电轴承建模提供一种新的思路。参考文献[1]吴云鹏,孙立红.滚动轴承力学模型的研究及其发展[J].煤矿机械,2004，(2):5－7.[2]林腾蛟.深沟球轴承运转过程动态特性有限元分析[J].振动与冲击学报,2009,28(1):118－122.[3]陈火红,杨剑,薛小香，等.新编Marc有限元实例教程[M].北京:机械工业出版社,2007.[4]陈达亮,郝志勇.GAP单元在气门弹簧动态特性有限元分析中的应用[J].拖拉机与农用运输车,2004，(2):24－26.[5]王国林,张建,王启唐，等.基于GAP单元的车架有限元分析[J].浙江大学学报,2008,29(3):206－209.[6]皮亚南,董懿琼,戴莉莉，等.滚动轴承保持架不稳定性分析[J].江西科学,2008,10(5):770－773.[7]HarrisTA，KotzalasMN.Advancedconceptsofbearingtech-nology[M].England&Wales：Taylor&FrancisGroup，2006.[8]AmasorrainJI，SagartzazuX，DamiánJ.Loaddistributioninafourcontact-pointslewingbearing[J].MechanismandMachineTheorz，2003，(38)：479－496.[9]冈本纯三.球轴承的设计计算[M].北京：机械工业出版社，2003.[10]翟瑞彩,谢伟松.数值分析[M].天津：天津大学出版社，2000.4结论本文采用有限元分析方法，用GAP单元模拟轴承内、外圈的接触，然后在添加外荷载的情况下对模型进行静强度分析计算；同时，采用理论方法，在同种工况下对该轴承分析计算，针对两种不同方法的计算结果分析对比，可得到以下结论：1)GAP单元建模计算结果与理论计算结果误差较小，吻合度高。2)GAP单元模拟轴承滚动体进行建模，可实现轴承内、外圈力的传递。3)当不针对滚动体进行分析时，可用GAP单元代替所有滚动体进行简化建模，不会影响其他连接部件的计算结果。4)GAP单元无法精确模拟轴承接触角的具体数值，只能近似模拟滚动体与滚道间的接触情况。要研究接触角时，需采取实体建模进行分析。5)GAP单元模拟滚动体对轴承的简化建模，(接第44页)在外加电场的作用下，可使较少的进行封装，然后分析PID现象产生的原因，并得出以下结论：1)使用EVA封装的 p型PERC双面双玻光伏组件易出现PID现象；2)即使使用POE材料封装的p型PERC双面双玻光伏组件，其背面出现PID现象的风险也较大，这与电池片本身的结构有关。参考文献[1] 葛华云. 基于光伏组件的电位诱发功率衰减的研究[D].长春: 吉林大学, 2013.[2] 王征, 郭海玲. 光伏组件PID现象抑制方法研究[J].太阳能, 2016, (5): 41－44.[3] 吴翠姑, 于波, 韩帅, 等. 多晶硅光伏组件功率衰减的原因分析及优化措施[J]. 电气技术, 2009, (8): 113－114.Na+通过封装材料到达电池片表面。使用POE封装的光伏组件背面更易出现PID现象是因为双面PERC电池片正面为化学钝化，其氮化硅中含有高密度的固定正电荷，对Na+有一定的排斥作用，会减弱一部分Na+的富集；但是其背面为场钝化，Al2O3/Si接触面具有较高的固定负电荷密度，背面玻璃中析出的Na+使氧化铝内的电荷发生再分布，导致钝化效果恶化。同时，双面PERC电池片正面含有一层氧化硅减反射层，可以起到抗PID效应，而背面没有。4 结论本文分别使用EVA和POE材料对光伏组件62SOLAR ENERGY 10/2018Copyright©博看网 www.bookan.com.cn. All Rights Reserved.