基于无网格法的风力发电机主轴强度仿真研究*

聂新宇,汪小芳,刘 勇,白聪儿,孙哲杰,夏国俊

(1.浙江运达风电股份有限公司,浙江 杭州 310012;2.浙江省风力发电技术重点实验室,浙江 杭州 310012;3.浙江大学 能源工程学院,浙江 杭州 310027)

0 引 言

随着碳达峰、碳中和战略目标的提出,风电行业又将迎来快速发展的良好机遇[1]。为了降低度电成本,风力发电机组的单机容量被不断提高,风轮直径与载荷也相应增大[2-3]。

作为机组传动链的关键一环,主轴负责将源自风轮的扭矩传递至齿轮箱、发电机,将倾覆力矩、推力等其他载荷传递至主轴承等支撑结构[4-5]。这些载荷的增大会显著提高主轴不同区域的应力,因此,为了保障机组的安全运行,基于强度仿真的主轴结构设计显得非常重要。

目前,关于风力发电机主轴的强度分析主要以有限元法为主。

王斌等人[4]利用ANSYS建立了主轴疲劳强度的分析模型,探究了卸荷槽宽度及深度等尺寸对锁紧螺纹段疲劳损伤的影响。汪亚洲等人[6]建立了主轴的极限强度有限元分析模型,并借此对轴承前沿大圆角和椭圆长径尺寸进行了优化,降低了主轴的最大应力。赵佰余等人[7]采用Workbench计算了主轴的极限与疲劳强度,并对可设计区进行了结构优化,在主轴减重的同时还提高了其安全余量。白儒等人[8]通过有限元分析,从结构刚度协调性和表面粗糙度两方面对主轴进行了优化,在满足强度要求的前提下,实现了主轴的轻量化设计目标。

为准确模拟载荷的传递路径,这些主轴模型中均包含主轴承结构,通常为双列调心滚子轴承(SRB)或单列圆锥滚子轴承(TRB),包含数十至上百个滚子,一般采用杆单元、弹簧单元或间隙单元对其作等效处理[9-10],设定仅受压属性。其建模与计算过程较为复杂耗时。

为了提高仿真效率,近年来,基于无网格法的仿真方式在电子设备、车辆、医学等领域得到越来越多的应用[11]。例如,MAGALHES F C等人[12-13]对3D打印、器官等复杂结构进行了静强度评估,简化了网格划分过程,有效识别了结构集中应力。顾张祺等人[14]对电子产品组件开展了模态分析与加速度过载分析,简化了几何清理、接触设置等环节,有效缩短了产品研制初期的设计时间。

由上述研究可知,无网格仿真可以简化传统结构仿真中的前处理环节,适合对含接触的大型装配体强度进行快速评估。然而,目前该方法在含轴承的主轴结构上还少有应用。

基于上述原因,笔者将以单SRB和双TRB两类机型的主轴结构为例,分别建立其有限元和无网格仿真模型,探究适用于三维原型的无网格模型前处理方法,分析无网格法在主轴装配体上的计算可靠性,及其在建模过程与计算时间上的优势,进而提出一种利用无网格法提高迭代效率的主轴结构设计方法。

1 主轴有限元模型

1.1 模型概述

笔者利用ANSYS软件建立主轴的有限元模型,采用GL2010规定的轮毂转动坐标系[15]。

轮毂转动坐标系如图1所示。

图1 轮毂转动坐标系

图1中,坐标系随轮毂转动,原点在轮毂中心,XR为转动轴方向,ZR朝向叶片,并与XR垂直,YR垂直于XR和ZR,根据右手定则确定其方向。

两类机型的主轴装配结构如图2所示。

图2 两类主轴装配结构示意图

装配体一般包含主轴、轮毂、轴承、挡圈、胀紧套及齿轮箱接口法兰等部件:单SRB机型只含一个主轴承,配合齿轮箱扭力臂支撑主轴;双TRB机型主要通过2个轴承对主轴进行支撑[10]。

ANSYS网格划分图如图3所示。

图3 ANSYS网格划分图

主轴结构:为准确模拟其集中应力,笔者采用高阶六面体单元SOLID186划分网格;轮毂结构:做了适当简化并采用较疏的网格。

SRB结构模型单元总数为368 343,节点总数为1 481 364;TRB结构单元总数为431 014,节点总数为1 638 017。

1.2 材料属性

主轴是整体锻件,材料为合金结构钢42CrMo4,弹性模量为2.1×105MPa,泊松比为0.3;

轮毂材料为QT400,弹性模量为1.69×105MPa,泊松比为0.275;

其他部件采用钢结构,材料属性参照主轴设定。

1.3 轴承等效

轴承滚子数量较多,且与内外滚道间为非线性接触。为了减少计算量,笔者采用多组受压不受拉的杆单元LINK10等效模拟[7],通过调整杆单元的截面积,使其刚度和滚子、滚道间的接触刚度一致。

根据ISO 16281[16],滚子与滚道间为线接触,其弹性常数CL可通过下式计算得到:

(1)

式中:Lwe为滚子有效接触长度。

进而可得接触载荷Q与位移δ间的关系式为:

Q=CLδ10/9

(2)

等效后的轴承模型如图4所示。

图4 轴承有限元模型剖面图

1.4 边界条件

模型边界条件设置如下:挡圈环向面与主轴、轴承间采用标准接触,摩擦系数为0.2,其他相邻部件间均为绑定接触;轴承安装于轴承座内,外侧环向面可视为固定约束。在实际运行中,主轴后端与齿轮箱相连,其轴向转动受到限制。

因此,笔者在齿轮箱重心与胀紧套或齿轮箱接口法兰间建立刚性连接,约束重心点ROTX=UY=UZ=0;风轮受到重力和风动载荷等,利用Bladed软件计算其各向力矩达到极限时的载荷,包括Mx_max、My_max、Mz_max和Myz_max4种工况;在轮毂中心建立加载点,通过刚性梁单元与变桨轴承的安装面相连,将载荷施加于该加载点。

两类主轴的4种极限工况载荷如表1所示。

表1 极限载荷表

2 无网格模型

2.1 背景介绍

在有限元仿真中,存在网格划分复杂、大型装配体计算耗时的问题。近年来,一些简化甚至略去网格划分的数值分析方法被提了出来,工程上统称为无网格法。

无网格法的理论基础多样:如外部有限元逼近法(external finite element approximation,EFEA),采用式(3)所示近似函数,在单元边界自由度基础上附加内部自由度,可用不规则几何体离散求解域,而非四面体、六面体等标准单元,因此较容易采用软件实现自动网格划分[17];此外还有隐式边界法(implicit boundary method,IBM),采用不受分析对象结构限制的规则网格,对完全处于结构内部、部分处于结构内部的网格,采取不同算法进行分析计算[18]。

笔者将以基于外部有限元逼近法的SimSolid软件为例,进行主轴的建模分析。

外部有限元逼近法的近似函数如下:

(3)

式中:uh为近似解;ai为单元内部自由度;pi为内部自由度基函数;bk为单元边界自由度;pk为边界自由度基函数。

2.2 建模方法

SimSolid建模步骤如下:1)导入几何体;2)设定材料参数;3)创建接触关系;4)建立分析模型。

材料、载荷等参数与上文有限元模型一致;所不同的是:1)不含网格划分过程;2)导入模型为几何体原型,主轴承不做简化。

在该无网格模型中,接触是通过几何体间的接触对实现。接触对的建立由如下3个参数控制:间隙容差、穿透容差和连接分辨率。合理的参数设置要保证滚子与滚道间充分接触,同时刚度不至太高。

模型边界条件设置方法如下:

1)远程载荷施加。以轮毂中心为远程点,轮毂上的3个变桨轴承的安装面为加载面,施加表1所示载荷。

2)主轴尾端约束。在胀紧套或齿轮箱接口法兰端面上,建立“接地衬套”形式的虚拟连接,设定轴向转动刚度与垂直轴向的平动刚度为较大数值(此处设定为1020N·m/rad和1020N/m)。

3)轴承移动约束。在轴承外侧环向面上施加铰链约束,保证轴承可以绕轴向转动,同时轴向、径向平动自由度为零。

以SRB模型为例,滚子与滚道间的接触对呈长条状。轴承内部接触示意图如图5所示。

图5 轴承内部接触示意图

3 结果与讨论

3.1 刚度分析

有限元分析是业内公认的主轴强度分析方法[4,7,10]。因此,笔者以有限元结果为基准,判断无网格模型的准确度。

首先,比较结构变形。此处以Myz_max工况为例,两类结构在不同仿真方式下的计算合位移,如图6所示。

由图6可以看出:对于单SRB结构,SimSolid计算位移分布规律与ANSYS一致,其最大值均出现在轮毂端部,分别是31.06 mm和31.61 mm,相对误差为1.77%,该处远离轴承和齿轮箱的约束,因此位移较大;

对于双TRB结构,最大位移同样出现在轮毂端部,数值分别是12.06 mm和11.88 mm,相对误差为-1.49%。

在不同工况下,2种模型的最大合位移对比如表2所示。

表2 最大合位移对比表

由表2可知:最大合位移的相对误差在-6.06%～2.28%之间。在仿真模型中,轮毂与主轴绑定,主轴的变形受到轴承的影响,轴承刚度越大,对主轴的变形限制越明显,轮毂位移量越小;反之增大。

位移相近,说明无网格法对轴承刚度的模拟较为可靠。

3.2 单SRB结构集中应力分析

在刚度计算准确的基础上,笔者分析了无网格法计算主轴应力的准确度。

此处,笔者首先以单SRB结构Myz_max工况为例。

2种模型计算得到的主轴Mises应力分布如图7所示。

图7 单SRB结构Myz_max工况下主轴Mises应力云图

由图7可以看出:主轴应力均呈对称分布,对称轴为合弯矩Myz方向,主轴首尾两端应力较高。前端圆角位于轴承与轮毂之间,受弯矩最大,因此应力水平普遍较高。叠加该区域的截面积变化梯度大,应力集中现象较为明显。主轴尾端半径最小,受扭矩引起的切应力最高。

笔者统计这两个区域不同工况下的最大Mises应力结果,如图8所示。

图8 单SRB结构主轴不同工况下2种模型应力对比图

由图8可以看出:2种模型的最大Mises应力相对误差在-2.36%～2.94%之间,最大值均出现在Myz_max工况下的尾端,分别是295 MPa和289 MPa;

主轴所用合金钢的屈服强度σs=560 MPa,根据GL2010规定,材料安全系数γm取1.1[15],由式(4)可得主轴许用应力[σ]=509 MPa,2种模型计算主轴应力均满足极限强度要求。

许用应力计算表达式如下:

[σ]=σs/γm

(4)

3.3 双TRB结构集中应力分析

笔者接下来分析双TRB结构主轴的应力。

同样以Myz_max工况为例,ANSYS与SimSolid模型计算得到的主轴Mises应力分布,如图9所示。

由图9可以清楚地看到:主轴最大应力均位于前端圆角处,分别是147 MPa和148 MPa,相对误差为0.68%,尾端最大应力分别是90 MPa与92 MPa,相对误差为2.22%(如表1所示);该工况下,扭矩在各向载荷中的占比较低,导致尾端应力要明显小于前端圆角。

笔者统计了该结构的2种模型在不同工况下的应力计算结果,如图10所示。

图10 双TRB结构主轴不同工况下2种模型应力对比图

由图10可以看出:2种模型最大Mises应力相对误差在-6.71%～4.76%之间,最大值均出现在My_max工况下的前端圆角处,分别是168 MPa和176 MPa,同样满足极限强度设计要求。

3.4 计算效率对比

在建模步骤、计算时间与结果方面,笔者对2种模型的差异进行了统计,结果如表3所示。

表3 2种模型对比表

由表3可知:相较于ANSYS,基于无网格法的SimSolid建模省去了网格划分和轴承等效等步骤,大大缩短了建模时间,降低了建模难度;

在计算时间上,针对单SRB和双TRB结构,ANSYS计算时间平均为13 min 44 s和29 min 24 s,SimSolid为3 min 12 s和3 min 36 s,相对用时缩短了77%与88%,有限元法计算时间受模型规模、结构复杂度影响较大,无网格法则相对较小,展现了其在风力发电机组结构越来越复杂背景下的应用潜力。

3.5 无网格辅助主轴结构设计方法

鉴于无网格法在计算效率上的优势,笔者提出一种无网格辅助主轴结构设计方法。

该主轴结构设计方法的流程图如图11所示。

图11 无网格辅助主轴结构设计流程图

在新型尤其是大功率主轴的研发过程中,基于经验设计的初始结构往往和最终结构差异较大,需经历多次迭代过程,以得到满足极限、疲劳强度要求的结构。

由于有限元分析的迭代成本较高,因此,笔者在传统的设计流程中添加图11虚线框所示无网格法,利用无网格法计算高效的特点,预先优化主轴结构,减少后续迭代次数,从而缩短整体设计时间[19]。

4 结束语

无网格仿真可以简化传统结构仿真中的前处理环节,然而,目前该方法在含轴承的主轴结构上还少有应用。

为了提高风力发电机中传递载荷的关键部件—主轴的设计效率,笔者以单SRB和双TRB两类机型的主轴结构为例,分别建立了有限元仿真模型和无网格仿真模型,分析无网格法在主轴装配体上的计算可靠性,及其在建模过程与计算时间上的优势,进而提出了一种利用无网格法提高迭代效率的主轴结构设计方法。

研究结论如下:

1)对于两类机型,无网格法计算主轴变形结果与有限元分析结果相吻合,最大位移相对误差在-6.06%～2.28%之间,说明无网格法在主轴及相关轴承的刚度模拟上的准确度较高;

2)比较2种模型可知,在主轴前端圆角和尾端等应力集中区域,两类机型的最大应力相对误差分别在-2.36%～2.94%及-6.71%～4.76%之间,说明无网格法应力计算结果是可靠的;

3)在保证计算精度的同时,无网格法的建模过程更加简便,同时计算时间显著减低,由此,提出了一种无网格辅助主轴结构设计的方法。

在后续的研究中,笔者拟将无网格法应用于轮毂、机架等含复杂曲面,以及装配关系更为复杂的结构件中,以扩展该方法的应用范围,进一步提升风力发电机组的结构设计效率。