基于CEL法的故障齿轮热流固耦合有限元分析

曹亚超,卢 伟,崔素晓,马雷晓

(河北科技大学机械工程学院,河北石家庄 050018)

风力发电机在实际运行过程中,齿轮传动系统在环境风载荷的影响下经常出现疲劳点蚀、齿根断裂等故障。2个齿轮在啮合时,由于齿轮摩擦产生热量及润滑油的存在,状态分析具有一定的复杂性。因此,模拟齿轮传动中的热流固耦合效应对于研究齿轮传动的性能和可靠性具有重要意义[1-2]。

随着计算机计算能力的增强,计算流体力学(computational fluid dynamics, CFD)、欧拉-拉格朗日法 (coupled Eulerian Lagrangian, CEL)等方法已经应用在流体仿真的各个领域,并且可以在流固场仿真时引入温度场,同时进行热流固耦合场仿真[3-4]。MARCHESSE等[5]以斜齿轮为研究对象,利用CFD技术研究了斜齿轮的轮旋角度与风阻产生损耗功率大小的关系。CONCLI等[6]采用CFD技术对圆柱齿轮浸入润滑油深度和搅油损耗功率之间的关系进行研究。LUO等[7]研究了齿轮压力角、齿轮半径和齿廓形状对齿轮产生温度的影响。RANSEGNOLA等[8]考虑了流体-固体的动态相互作用和边界摩擦,在外啮合齿齿轮尖、侧面和滑动轴承的润滑处引入了浮膜的数值阻抗,从而实现了齿轮的强耦合仿真。国内对复杂耦合场仿真研究起步较晚,刘培坤等[9]利用模拟软件FLUENT 14.5对液体平衡环内部流场进行了数值模拟研究。朱黎明等[10]采用CFD方法对化学气相沉积制备Co-Re合金涂层反应室进行三维数值模拟分析。冯永保等[11]基于流固热耦合和流固耦合方法,利用AMESim以及FLUENT软件对高速高压下的圆弧齿轮泵摩擦副进行摩擦仿真研究。

传统齿轮热流固耦合仿真使用CFD法,而CFD法计算量较大且对边界条件的敏感性较高,需要大量的计算资源和时间,复杂的流动问题可能需要高性能计算机或集群来完成模拟,且主要研究的是正常齿轮的热流固耦合问题,而CEL法计算量相对较小且对于流动的微观细节有较好的描述能力,齿轮箱工作时复杂的热流固耦合环境需要一种更有效的方法对齿轮运行状态进行模拟[12-13],CEL法的出现使得故障齿轮热流固耦合效应机理的研究得以实现。

本研究旨在采用一种更加有效的仿真算法对齿轮热流固耦合进行模拟仿真,从而更加深入地探究齿轮传动中热流固耦合效应的机理和规律,可为齿轮模拟仿真提供指导和方法。

1 基于CEL的有限元仿真模型

QPZZ-Ⅱ故障实验平台包含一对标准渐开线直齿圆柱齿轮,通过UG软件的GC工具箱功能绘制齿轮,随后对齿轮进行修改,使其符合实际要求,齿轮形状参数见表1。

表1 齿轮形状参数

1.1 网格划分及材料赋予

图1 齿轮箱三维模型Fig.1 3D models of gear box

表2 齿轮材料参数

表3 润滑油材料参数

图2为欧拉体定义材料过程。该方法欧拉网格保持不动,材料在网格中流动。进行分析时,材料各个时变属性参数都会映射到欧拉网格上[14]。材料在欧拉网格流动时,通过单元中的欧拉体积分数(eulerian volume fraction, EVF)描述材料在网格中的状态,其中,“1.0”表示充满网格,“0.0”表示网格内没有材料,范围为0.0～1.0的数字表示材料在网格内已填充占比量[15]。

图2 欧拉体定义材料过程Fig.2 Eulerian definition of the material process

齿轮箱采用浸油式润滑,进行齿轮传动分析时,需要同时考虑到材料变形和流固耦合问题。考虑CEL方法具有较强的耦合特性,适用于流体与固体耦合模型,对该齿轮啮合传动系统采用CEL方法进行仿真,通过ABAQUS/Explicit软件显式动力学强大的通用接触模块,利用CEL法可有效完成齿轮啮合传动的热流固仿真。相对复杂的齿轮模型用HyperMesh软件绘制网格,润滑油属于流体,在建立时选用欧拉体。各部件网格数量见表4,其中齿轮网格选用C3D8RT,润滑油网格选用EC3D8RT,有限元仿真网格模型见图3。

图3 有限元仿真网格模型Fig.3 Finite element simulation grid model

表4 各部件网格数量

1.2 边界条件加载

为方便约束和载荷的施加,在2个齿轮中心分别创建一个参考点,将参考点与对应的齿轮进行耦合,在参考点上添加约束的齿轮边界条件和载荷。

首先是齿轮之间的接触,齿轮在相互接触时将沿切向发生相对滑动,产生摩擦热,轮齿之间以 “切向行为”为主,其接触属性选择“摩擦系数0.05”。“法向行为”影响较小,可忽略接触属性,选择“硬”接触,并在热学模块选择“生热”选项。其次是齿轮与润滑油或润滑油自身接触,该类型接触可忽略模型之间的摩擦和生热,所有模块参数选择“无”。然后指派接触面和设定热传导,在分析步中选用通用接触,指派选择齿轮间接触属性为两齿轮啮合面接触属性,指派齿轮与润滑油间的接触属性为全局接触属性,完成模型接触属性的指派。

定义热传导参数时需要对模型设置表面热交换和表面辐射2个条件,查阅相关资料后,在表面热交换条件中设置齿轮与润滑油对流换热系数为1 200 W/(m2·K),在表面辐射条件中设置齿轮表面辐射率为0.8 W/m2,完成热传导的设置。设置齿轮和润滑油初始温度均为20 ℃。同时在该分析步中为润滑油模型选择“材料指派”,设置润滑油的初始位置和流动空间,完成载荷的设定。

2 温度场影响下的流固耦合仿真结果分析

图4为齿轮啮合温度场分布图,由图4可知齿轮啮合传动时温度场主要分布在齿轮间接触部位。图5为不同载荷下齿轮接触应力云图,图中左侧为主动轮,右侧为从动轮。由图可知,载荷为5 N·m时,其最大啮合应力为109.4 MPa;载荷为2.5 N·m时,其最大啮合应力为95.08 MPa;载荷为1.25 N·m时,其最大啮合应力为93.10 MPa;载荷为0 N·m时,其最大啮合应力为87.95 MPa。

图4 齿轮啮合温度场分布图Fig.4 Distribution diagram of gear meshing temperature field

由图5仿真结果可知,施加载荷不同;齿轮啮合时齿厚方向的应力也不同;当施加载荷不断增加时,相同工况下齿轮的接触应力也会增大,且最大应力均出现在啮合面上。

今年10月，由世界中餐业联合会、亚洲食学论坛组委会联合主办的“第八届亚洲食学论坛”在北京隆重举行，北京烹饪协会承办了论道京味餐饮“亚洲食学论坛”分论坛活动。论坛会上，云程会长代表北京烹饪协会做了题为《重视传承是发展京菜的核心》的主题演讲，围绕中国改革开放后京菜的发展历史，以及为什么要抓京菜、如何抓好新时代京菜的发展等问题进行阐述，受到与会者的共鸣。会上多位知名教授、专家从不同角度论述了京味餐饮文化的历史、现在和未来。

在ABAQUS显示动力学中,显式算法直接求解切线刚度不需要迭代,但显式算法要求质量矩阵为对角矩阵,因而采用减缩积分法,但该方法容易激发沙漏问题。为确保计算结果的准确性,在进行仿真时对网格单元选择完全积分。在判定计算结果是否受到沙漏效应影响时,可对比模型伪应变能与内能的比值,当伪应变能占内能的比例低于1%时仿真结果准确度较高,超过10%时仿真结果无法使用。

根据图5中的仿真结果,选择伪应变能和内能绘制能量曲线(见图6),由图可知,显示模型的伪应变能远小于内能,所占比例在1%以下时,仿真计算精确度较高。选取载荷为5 N·m时的温度场对齿轮接触应力的影响进行研究。图7为有无温度场齿轮接触应力云图。齿轮啮合过程分为3个阶段:双齿啮入阶段、单齿接触阶段和双齿啮出阶段。各阶段有、无温度场啮合过程中最大应力见表5。

图6 能量曲线Fig.6 Energy curves

表5 有、无温度场啮合过程中最大应力

由图7可知,在齿轮啮合的3个阶段,有温度场对比无温度场时,最大接触应力分别增加了16.25,26.14,3.01 MPa。相对于无温度场3个阶段的最大接触应力分别增长了20.63%,41.28%,6.13%。无温度场时和有温度场时,齿轮在3个啮合阶段中的接触应力差异较大,对齿轮传动系统进行有限元分析时,添加温度场后最大接触应力增加且增幅较大,温度场的介入对齿轮传动仿真结果有着重要影响。

3 多故障齿轮系统应力场分析

齿轮传动系统经常处于高速、连续性的工作环境,齿轮会因为过载、达到疲劳极限或出现异物等原因,出现断齿、点蚀等故障。齿轮常见故障比例见表6,采用CEL法在上述温度场条件下对最常见的断齿和点蚀故障进行有限元仿真分析。

表6 齿轮常见故障比例

3.1 点蚀故障热流固耦合仿真与分析

齿轮传动在工作中会产生大量热量,这些热量会通过齿轮表面向外传递,因此要考虑齿轮表面的热扩散对于点蚀故障的影响。以从动轮发生点蚀故障进行研究,由齿轮传动接触特性可知,点蚀孔多数出现在齿面节线附近,并且随机产生,图8为齿轮点蚀模型。

图8 齿轮点蚀模型Fig.8 Gear pitting models

对点蚀故障的研究总共涉及到3种工况:轻度点蚀,有3个点蚀位置,直径为1 mm;中度点蚀,有5个点蚀位置,直径为1 mm;重度点蚀,有10个点蚀位置,其直径分别为1,0.5 mm。

图9为从动轮不同点蚀程度应力云图,图10为不同点蚀故障下沿齿厚方向应力图。应力沿齿厚方向传递效果不明显,点蚀故障下沿齿厚方向应力曲线见图11。由图9可知轻度点蚀最大接触应力236.6 MPa,中度点蚀最大接触应力为273.3 MPa,重度点蚀故障最大接触应力为292.7 MPa。

图9 从动轮不同点蚀故障应力云图

图10 不同点蚀故障下沿齿厚方向应力图Fig.10 Stress diagrams along the tooth thickness direction under different pitting faults

图11 点蚀故障下沿齿厚方向应力曲线Fig.11 Stress curves along the tooth thickness direction under pitting failure

由图9可知,与正常齿轮相比,点蚀故障齿轮表面接触应力分布极其不规则,在齿轮啮合碰撞时点蚀孔周围受到的挤压力较强,点蚀孔边缘的接触应力会明显的大于周边正常部位,点蚀边缘长期受到大应力挤压将会被破坏掉,从而导致点蚀扩展进一步加剧齿面点蚀程度。除点蚀孔周围应力较大之外,齿轮其余部分应力值在正常齿轮应力值附近上下浮动。由图10和图11可知,点蚀孔处齿轮接触面上的应力值大小与齿轮损坏程度有关,故障齿轮下的应力远大于正常齿轮应力,且损坏程度越大应力值就越大。在齿厚方向上,离接触面越远的部分,应力越接近于正常齿轮。由此可知,齿轮在点蚀故障下所受应力主要集中在接触面点蚀孔附近。

3.2 断齿故障热流固耦合仿真与分析

结合断齿故障实物零件与相关研究,将从动轮断齿长度分别为5,10,15 mm和中间断齿的齿轮断齿模型(见图12)作为研究对象进行有限元仿真,齿轮中间断裂的有效长度在5 mm左右。图13为不同断齿故障接触应力云图,由图可知正常情况下齿面上最大接触应力为73.91 MPa,中间断齿轮齿最大接触应力为118.6 MPa,断齿长度为5 mm轮齿最大接触应力为127.8 MPa,断齿长度为10 mm最大接触应力为133.3 MPa,断齿长度为15 mm最大接触应力为168.6 MPa。图14为不同断齿故障沿齿宽方向接触应力分布,由图13和图14可知在断齿长度较为接近的条件下,中间断齿齿轮的接触应力要小于一侧断齿齿轮的接触应力。正常齿轮沿齿宽方向的接触应力分布较为均匀,基本呈直线分布,断齿故障齿轮在断裂位置最大接触应力会出现明显的增大,并且越靠近断裂端处接触应力越大且断裂深度越大,最大接触应力也越大。

图12 齿轮断齿模型

图13 不同断齿故障接触应力云图

图14 不同断齿故障沿齿宽方向接触应力分布Fig.14 Distribution of contact stress along the tooth width direction of broken tooth fault

图15为沿齿轮齿厚方向应力路径云图。沿齿厚方向应力曲线见图16,由齿厚方向应力仿真结果可知,沿齿厚方向应力路径是由齿轮接触面上接触应力最大位置走向齿轮另一侧齿根处,在不同断齿工况下,应力沿齿厚方向的应力趋势基本保持一致,断齿故障越严重齿厚方向上的整体应力越大,齿轮受到应力的影响就越大。

图15 沿齿厚方向应力路径云图Fig.15 Nephograms of stress path along the tooth thickness direction

4 结 语

本文主要针对风力发电机齿轮传动副在风载荷影响下经常出现的疲劳点蚀、断齿机械故障,采用CEL法模拟缺陷齿轮热流固耦合传动,研究缺陷齿轮接触应力演化规律。通过对比分析不同损坏程度下的应力仿真结果,得出以下结论。

1)施加载荷不同,齿轮啮合时齿厚方向的最大接触应力也不同。当施加载荷不断增加,同一工况下齿轮的接触应力增大,且最大接触应力均出现在啮合面上。

2)啮合前、啮合中和啮合后3个阶段在未添加热力场时的最大接触应力分别为78.78,63.33和49.07 MPa,添加热力场后3个阶段最大接触应力分别增长了20.63%,41.28%,6.13%,说明热力场对齿轮流固耦合分析结果的影响不容忽视,应在齿轮传动的有限元仿真中予以考虑。

3)从齿宽和齿厚方面来看,点蚀孔处接触应力明显大于齿面正常部分且点蚀损坏程度与其呈正相关;点蚀孔处被长期挤压会导致点蚀扩展,从而进一步加剧齿面破坏程度;从齿厚方面看,离接触面越远应力值越接近正常齿轮,说明点蚀故障对齿面的影响较大。

4)断齿故障齿轮在断裂位置接触应力会激增,并且越靠近断裂端处接触应力越大,且随断裂程度增加而进一步增大;在不同断齿工况下,应力沿齿厚方向的变化趋势保持一致,断齿程度越高整体应力值越大。

采用CEL热流固耦合法模拟齿轮啮合传动过程,不仅可以考虑流热场在齿轮啮合时的大变形行为,还能有效避免网格畸变导致的计算不收敛问题,可为探究齿轮传动中热流固耦合效应的机理和规律提供指导和方法。

本文主要采用CEL法对风机传动齿轮啮合过程进行了数值仿真研究,仿真过程中忽略了齿轮结构变形对流场空间的改变,模拟环境较为理想,未来还需要进一步针对实际工况对研究进行完善,如考虑摩擦生热、双向耦合等因素,并深入开展相关实验研究。