小型联合收割机齿轮传动箱刚柔耦合动力学仿真分析

张卫华，方永锋，许 勇，张 宁，余欢乐

(贵州工程应用技术学院机械工程学院，贵州毕节 551700)

以某小型联合收割机齿轮传动箱为研究对象，该收割机底盘以3个前进速度和1个倒档速度工作[1]，分别为一档(0.52 m/s)、二档(0.82 m/s)、三档(2 m/s)，倒档速度为 0.33 m/s。在不同收割作业工况下，改变发动机传到输入轴上的转矩及转速，使整机得到不同的牵引力和速度，以达到协调稳定工作的目的。笔者已对齿轮传动箱进行了齿轮动力学仿真分析[1]，但未考虑齿轮的弹性变形，且输入轴为细长轴，有必要考虑轴的弹性变形对整个齿轮变速箱的影响，因此对齿轮系统进行刚柔耦合的动力学特性进行研究，设计制造小振动、低噪声、高可靠性、高传动性能的齿轮变速系统[2]，提高和改善收割机工作的平稳性能具有重要意义。

1 创建刚柔耦合虚拟样机模型

1.1 刚柔耦合动力学建模方法

根据参考坐标系选取的不同柔性体建模方法可以归为3类[3]：惯性坐标系方法、随转坐标系方法、浮动坐标系方法。惯性坐标系方法源于大变形非线性有限元；随转坐标系方法源于计算结构动力学；浮动坐标系方法是将多刚体动力学与结构动力学相结合的一种方法，这种方法使多刚体动力学软件[机械系统动力学自动分析(automatic dynamic analysis of mechanical systems，简称ADAMS)]扩展应用于柔性多体系统成为可能，是目前柔性多体系统建模使用最广泛的方法。

在ADAMS柔性模块中，柔性体采用离散化若干个单元的有限个结点自由度来表示物体的无限多个自由度。这些单元结点的弹性变形可近似地用少量模态的线性组合来表示。若物体坐标系的位置用其在惯性坐标系中的笛卡尔坐标X′=(x,y,z)和反映刚体方位的欧拉角ψ=(φ,θ,φ)来表示，模态坐标用q={q1,q2,…,qm}M(M为模态坐标数)表示，则柔性体的广义坐标为[4-7]：

(1)

则柔性体上任一结点I的位置向量可表示为：

rI=ri+A(sI+φIq)。

(2)

式中：ri为结点I在惯性坐标下的位置；A为物体坐标系到惯性参考系的转换矩阵；sI为结点I在物体坐标系中变形时的位置；φI对应于结点I的移动自由度的模态矩阵子块；将式(2)对时间求导，便得到该点的移动速度：

(3)

ω为物体坐标系的角速度向量；B为欧拉角的时间导数与角速度向量之间的转换矩阵；“～”为向量对应的对称矩阵。

结点I的角速度可以用物体的刚体角速度与变形角速度之和来表示：

ωI=ω+φ′qω。

(4)

式中：φ′对应于结点I的转动自由度的模态矩阵子块。

1.2 输出模态中性文件

在ANSYS中生成有限元模型，通过ANSYS与ADAMS软件之间的双向数据传输接口，输出模态中性文件(.MNF)[8]。从ADAMS中的Flex模块将导入的模态中性文件调取，引入柔性体，考虑柔性体大范围运动的高度耦合作用，从而提高仿真系统的动力学性态。将有限元模型导入ADAMS中在生成模态中性文件时关键在于建立铰接点(interface point)，在导入ADAMS后这些点会自动生成Marker点，通过这些点和其他刚体或柔性体建立联接[9]。在ANSYS中输出模态中性文件(.MNF)的步骤如下：

1.2.1 建立铰接点 在ANSYS中，输出模态中性文件时，系统要求至少选择2个以上的铰接节点。因此对于齿轮，在有限元模型的2个端面中心处创建2个关键点(key point)，其编号不能与模型中己存在的关键点编号重复；对于轴在轴承联接的中心位置处，在轴的两端建立关键点。然后用三维质量单元MASS21来划分所建立关键点，生成铰接节点。在设置MASS21的实常数时应使其质量、惯性矩等足够小(本研究均设置为1×10-6)。

1.2.2 创建刚性区域(rigid region) 选取上一步创建的关键点(key point)，再选取周围的节点创建刚性区域，建立起约束方程(constraint equations)和周围的节点建立刚性联接[也可以选用蜘蛛网刚性梁(spider web beams)]，将作用于约束副上的载荷，通过刚性梁分配到柔性体的节点上，以避免柔性体上应力集中[8]。选取齿轮与轴接触的节点(齿轮轴孔上的节点)、轴与轴承接触的节点(轴与轴承接触的表面节点)创建刚性区域。

1.2.3 输出模态中性文件 输出设置时要注意单位的转换，ANSYS中没有明确的单位定义，只能按照一定的比例进行单位换算，长度的比例因子为1 000，质量的比例因子为 0.001[8]。确定提取的模态阶数，铰接点是1个有6个自由度的节点，每个自由度都相对应1个约束模态，所以在模态中性文件中所记录的模态阶数等于ANSYS中指定的模态阶数+外联节点数×6，本研究按系统默认值设置为6阶，如果要在后处理中启动ADAMS/durably模块，查看柔性体的应力应变图，输出模态中性文件时应设置为“Include stress and strain”，否则不能查看。

1.3 创建刚柔耦合仿真模型

在ADAMS/View中通过操作build→ADAMS/Flex→Rigid to flex便可读入上一步创建的模态中性文件，然后直接将柔性体替换掉刚性体，替换后相应的约束和载荷会自动转移到柔性体上，但要求柔性体与刚性体完全重合[10-13]。对于柔性体，不能修改其质量属性、几何、体积、惯性矩等，对柔性体的编辑主要是决定哪几阶模态参与运算，以节省运算时间，获得准确结果。如前所述，模态中性文件的模态为ANSYS系统默认的模态数(6阶)与2个关联节点的12阶自由模态组成。由于前6阶模态为构件的刚性平移和转动，其频率几乎为0，所以将第1～6阶模态过滤掉，仅保留余下的第7～18阶模态参与运算。轴和齿轮的柔性体替换刚性体后，刚性体上的载荷和约束以及Marker点都已自动转移到相应的柔性体上，所以不用重新添加载荷及约束关系，对于接触力(contact)，ADAMS中提供了柔性体对刚性体的接触力，即柔性体和实体的接触(flex to solid)。为测量输入轴(Ⅰ)以及中间轴(Ⅱ)两端轴承的支反力，选择柔性体和箱体，在轴承质心的位置分别创建旋转副。

为了和完全刚性体的动力学仿真结果作比较，接触力设置参数(stiffness、damping、force exponent、penetration depth)、负载、转动驱动、仿真步长和仿真时间与完全刚性时的仿真参数一致。图1为Ⅰ、Ⅱ、Ⅲ轴的柔性体以及第1、2级主齿轮的刚柔耦合模型。

2 结果与分析

2.1 各轴转速分析

图2至图4分别为Ⅰ、Ⅱ、Ⅴ轴转速随时间的变化曲线；图5至图7为各轴角加速度随时间的变化曲线。由图2至图4可知，各轴转速在加速阶段一直到稳定运行的阶段都有不同层次的波动。具体表现为输入轴在0.010～0.025 s之间有个微小的波动，之后平稳运行；中间轴在0.01～0.03 s有微小的波动(相对于稳定阶段)；而输出轴在0.00～0.04 s有较大的波动，波动幅值反方向最大值为202 °/s，正方向达287.8 °/s，与完全刚性时的速度有很大差别，但其转速波动逐渐减小，呈收敛趋势，0.04 s之后在平均值附近作微小的波动，平稳运行，并与理论值非常接近。

从加速度也可以看出，在发动机加速启动到稳定运行的一段时间内有一定波动，在加速度为0(速度最大)后，随时间很快衰减。这主要是因为柔性轴的振动以及齿轮的弹性变形，导致输出轴转速和角加速度在稳定运行之前按某一固有频率波动，但从平稳运行阶段来看，与完全刚性相比，各轴转速及加速度曲线，波动幅值减小。说明刚柔耦合中柔性体弹性变形对各轴转速以及加速度有一定的影响，更加接近齿轮转动实际的动态特性。图8为输出轴(Ⅴ)对角加速度曲线进行快速傅立叶变换得到的频域曲线，由图8可知，其主要以235 Hz的频率波动。

2.2 齿轮啮合接触力分析

图9、图10分别为第一级齿轮啮合的圆周力、径向力的时域及频域；图11、图12分别为第二级齿轮啮合的圆周力和径向力时域及频域。各图中上部曲线表示啮合力时域图，圆周力沿ADAMS坐标系统y方向，相同为正，相反为负；径向力沿坐标系统z方向，相同为正，相反为负。下部曲线表示啮合力频域。

2.2.1 从时间历程上分析 由仿真结果可知，第一级齿轮刚柔耦合接触力曲线变化与完全刚性体时的曲线变化基本一致。即在0.00～0.01 s为发动机启动加速阶段，启动瞬间有个较大的冲击力，与刚性时的冲击力比较明显减小，在0.03 s后平稳运行，在均值附近波动，圆周力和径向力波动幅值也比完全刚性时的小。产生这一结果是因为刚柔耦合分析中考虑到应力应变的问题，在啮合过程中，作为柔性的齿轮比作为刚性体的轮齿具有更大的弹性变形，这些变形在啮合过程中由于阻尼的存在会产生能量损失。

第二级齿轮啮合圆周力和径向力在加速阶段均有明显的波动且波动幅值较大，波动的频率下降，周期有所增大，这是因为此时的柔性轴可以等效为具有扭转刚度和弯曲刚度的梁，使齿轮的扭转刚度降低，啮合力振动的固有频率降低，周期增加，但随时间逐渐收敛。在0.03 s后平稳运行，圆周力波动幅值比完全刚性时的仿真值相对要小，而径向力则周期性地有个较大的波动值，且在0.048 7 s有个振动冲击，这是因为柔性体受力后产生弹性变形，使啮入冲击和啮出冲击导致啮合力存在一定的突变，但没有特别突出的极值出现，变化相对平稳，其波动范围都在均值附近，且变化趋势也体现轮齿有周期性的啮入和啮出的动态特性。由表1可以看出，刚柔耦合分析的动态接触力均值与完全刚性体分析得到的接触力值非常接近，但从图9至图12可知，刚柔耦合的结果更接近齿轮系统实际的动态特性。

表1 齿轮啮合力刚柔耦合与完全刚性时仿真均值比较

2.2.2 从频域上分析 第一级齿轮圆周力和径向力啮合频率最大值均为397.1 Hz，相应力的幅值分别为62.9、26.0 N，且更好地符合齿轮的理论激振频率(396.7 Hz)。第二级齿轮圆周力和径向力的啮合频率最大值为235.1 Hz，其相应力的幅值分别为37.0、13.6 N，径向力振动较明显，其啮合频率逐级降低。从整体上来看，刚柔耦合仿真得到的齿轮啮合力在频域上基本以啮合频率及低次倍频为主的频谱结构，与完全刚性体时的振动情况相比，振幅有所减小，运行状态良好。综上所述，当考虑啮合轮齿的弹性变形时，啮合传动振动减小，并且具有较好的平稳性，更加符合齿轮传动设计的要求。

2.3 强度分析

由于轮齿在受载后产生变形会使其在理论啮合线外发生啮合，造成啮入和啮出冲击。齿轮在键联接处易造成应力集中对齿轮啮合传动、使用寿命等有很大的影响，所以有必要对齿轮进行强度校核。在ADAMS/view中加载Durability插件，在后处理中即可查看柔性体的变形及应力应变云图。由图13可知，轴的变形主要在轴承联接附近，而由图13-a可见，齿轮的变形分布在1、2、3、4等4个区域上；齿轮的最大应力(von mises stress)分布在轴孔周围以及键槽的位置上，并没有分布在轮齿上(图13-b)，图14为齿轮受力综合最大的10个节点列表。

由图14可以看出，齿轮的柔性体最大应力值为 102.339 MPa，发生在时间为0.006 7 s时对应的5 122节点上，图15为仿真过程中最大应力节点5 122的应力变化曲线。由图15可以得出，柔性体在整个运动周期内的动态等效应力的均值为22.400 MPa左右，而最大应力为102.339 MPa。本研究选用材料为45号钢，该材料热处理后的屈服极限为360.000 MPa，所以最大节点应力102.339 MPa小于材料屈服极限，齿轮结构强度符合要求。同理可对第二级主齿轮进行强度分析(图13-c)。

3 结论

基于柔性多体理论，通过柔性体直接替换刚性体的方法，建立齿轮传动箱的刚柔耦合动力学模型，设置与完全刚性的动力学模型相同的仿真参数、求解、仿真步长进行动力学仿真，仿真分析结果表明，刚柔耦合所得的结果基本吻合，但刚柔耦合更加符合齿轮传动的实际动态特性，且其运行平稳性相对提高。对齿轮的强度分析结果表明，其结构设计符合要求。刚柔耦合仿真所得到的载荷谱可以作为进一步有限元分析的边界条件，也可为齿轮系统零部件的寿命计算、故障诊断和预测提供依据。

[1]张卫华，尹 健，蔡 群.小型联合收割机变速器虚拟样机建模及动力学仿真[J].数字技术与应用，2011(12)：245-247.

[2]李润方，王建军.齿轮系统动力学——振动、冲击、噪声[M].北京：科学出版社，1997.

[3]Wasfy T M，Noor A K.Computational strategies for flexible multibody systems[J].Applied Mechanics Reviews，2003，56(6)：553-613.

[4]Shabana A A.Dynamic of multibody system[M].New York：Cambridge University Press，1989.

[5]Garcia-Vallejo D，Sugiyama H，Shabana A A.Finite element analysis of the geometric stiffening effect.Part 1：a correction in the floating frame of reference formulation[J].Proceedings of the Institution of Mechanical Engineers，Part K：Journal of Multi-Body Dynamics，2005，219(2)：187-202.

[6]ADAMS/FLEX user’s manual[Z].Mechanical Dynamics Inc.，2000.

[7]刘 俊，林砺宗，刘小平，等.ADAMS柔性体运动仿真分析研究及运用[J].现代制造工程，2004(5)：53-55.

[8]崔玉玲.履带式挖掘机行走减速机构动力学仿真研究[D].长春：吉林大学，2011.

[9]王 炎，马吉胜，蒙 刚.齿轮系统刚柔耦合动力学建模与仿真研究[J].机械传动，2009，33(4)：32-35.

[10]王建军，李润方.齿轮系统动力学理论体系[J].中国机械工程，1998(12)：55-58.

[11]陈立平.机械系统动力学分析及ADAMS应用教程[M].北京：清华大学出版社，2005.

[12]曾 攀，雷丽萍，方 刚.基于ANSYS平台有限元分析手册：结构的建模与分析[M].北京：机械工业出版社，2011.

[13]梁 君，赵登峰.模态分析方法综述[J].现代制造工程，2006(8)：139-141.