一种考虑齿轮和轴承非线性耦合的传动系统数值模拟方法与流程

本发明涉及一种传动系统数值模拟方法，特别是一种考虑齿轮和轴承非线性耦合的传动系统数值模拟方法，属于机械传动
技术领域：
。
背景技术：
：齿轮和轴承是驱动桥、变速箱等机械传动系统的重要组成部分，齿轮的啮合状态、轴承滚子与滚道之间的接触状态均随系统载荷工况的变化呈现非线性特性。同时，齿轮和轴承之间也存在非线性耦合：一方面，齿轮啮合状态的变化会影响传动轴系和轴承之间的载荷分配，从而导致轴承刚度发生变化；另一方面，轴承刚度的变化会影响齿轮在对应载荷工况下的错位量，进而影响齿轮的啮合状态。因此，在对传动系统进行建模和分析时，应准确考虑齿轮和轴承之间的非线性耦合影响。现有研究在对传动系统进行建模分析时，通常采用以下两种数值模拟方法：1)有限元接触计算方法：基于abaqus等商用有限元软件能够准确实现有限元接触计算，但传动系统通常包含多个滚动轴承，每个轴承的滚子和滚道之间均存在接触关系，而有限元接触计算对网格精度要求较高，受收敛性和计算规模的限制，难以在接触计算模型中同时考虑齿轮接触、每个轴承滚子与滚道的接触，所以现有研究通常仅对齿轮进行接触定义，忽略轴承非线性特性的影响，或单独对轴承进行接触计算，这样均无法真实体现轴承与齿轮之间的非线性耦合。2)基于非线性轴承单元的传动系统建模分析方法：解析形式的非线性轴承单元(harrista,kotzalasmn.essentialconceptsofbearingtechnology.5thed.crcpress,2006.)能够快速有效地模拟轴承的非线性刚度特性，但为了实现包含多个轴承的传动系统模型的非线性数值求解，齿轮的啮合关系通常采用空间弹簧形式的等效啮合单元模拟，并采用理论解析公式(donleymg,limtc,steyergc.dynamicanalysisofautomotivegearingsystems.saeinternationalcongressandexposition,detroit,michigan:1992,920762.)计算齿轮等效啮合参数，没有考虑齿面摩擦的影响，且无法准确体现齿轮非线性接触特性对传动系统的影响。目前仍缺乏一种能够准确考虑齿轮和轴承非线性耦合的传动系统数值模拟方法。技术实现要素：针对上述问题，本发明的目的是提供一种能够准确考虑齿轮和轴承非线性耦合的传动系统数值模拟方法。为实现上述目的，本发明采取以下技术方案：一种考虑齿轮和轴承非线性耦合的传动系统数值模拟方法，包括以下步骤：1)建立传动系统部件的体单元有限元模型：建立输入轴、输出轴和箱体的体单元有限元模型，其中，输入轴与输出轴之间通过齿轮啮合关系连接，输入轴和输出轴分别通过轴承与箱体连接，输入转矩施加在输入轴上，输出转矩施加在输出轴上，在输入轴、输出轴和箱体的体单元有限元模型中包含上述连接关系及加载的位置处建立边界节点，具体如下：①输入轴：轴承中心节点、输入转矩加载节点、主动齿轮等效啮合节点；②输出轴：轴承中心节点、输出转矩加载节点、从动齿轮等效啮合节点；③箱体：轴承中心节点。每个边界节点均包含6个自由度，分别用刚性连接单元将边界节点与各部件的体单元有限元模型连接，并根据实际情况在箱体的体单元有限元模型施加外部约束。2)建立包含非线性轴承单元和齿轮等效啮合单元的传动系统缩聚模型：采用解析形式的非线性轴承单元对轴承进行模拟，设轴承轴向为其局部坐标系z轴方向，则非线性轴承单元的刚度矩阵kb表示为式(1)：上式中，fx和fy为轴承单元沿局部坐标系x轴和y轴方向所受的径向力；fz为轴承单元所受的轴向力；mx和my轴承单元绕局部坐标系x轴和y轴方向所受的径向弯矩；δx和δy为轴承单元沿局部坐标系x轴和y轴方向的径向变形；δz为轴承单元的轴向变形；θx和θy为轴承单元绕局部坐标系x轴和y轴方向的转角变形。为了实现对包含非线性轴承单元的传动系统模型进行非线性迭代数值求解，需要对包含大量节点自由度的体单元模型进行缩聚变换，采用guyan缩聚法实现，仅保留边界节点自由度，计算各部件的缩聚刚度矩阵：其中，输入轴的缩聚刚度矩阵ki表示为式(2)：上式中，kiaa为输入轴边界自由度对应的刚度矩阵；kibb为输入轴内部自由度对应的刚度矩阵；kiab和kiba为输入轴边界自由度和内部自由度的刚度耦合项。输出轴的缩聚刚度矩阵ko表示为式(3)：上式中，koaa为输出轴边界自由度对应的刚度矩阵；kobb为输出轴内部自由度对应的刚度矩阵；koab和koba为输出轴边界自由度和内部自由度的刚度耦合项。箱体的缩聚刚度矩阵kh表示为式(4)：上式中，khaa为箱体边界自由度对应的刚度矩阵；khbb为箱体内部自由度对应的刚度矩阵；khab和khba为箱体边界自由度和内部自由度的刚度耦合项。轴承单元的一端连接输入轴或输出轴缩聚模型的轴承中心节点，另一端连接与之对应的箱体缩聚模型的轴承中心节点。采用等效啮合单元对齿轮副的啮合状态进行模拟，齿轮等效啮合刚度矩阵km表示为式(5)：上式中，n为等效啮合力作用方向单位矢量；km为啮合刚度系数。等效啮合单元的一端连接输入轴缩聚模型的主动齿轮等效啮合节点，另一端连接输出轴缩聚模型的从动齿轮等效啮合节点。将上述输入轴的缩聚刚度矩阵ki、输出轴的缩聚刚度矩阵ko、箱体的缩聚刚度矩阵kh、非线性轴承单元的刚度矩阵kb、齿轮等效啮合刚度矩阵km按照连接关系组集，即获得完整传动系统缩聚模型的刚度矩阵k表示为式(6)：3)传动系统缩聚模型的非线性静力学求解和线性等效轴承刚度计算：在输入转矩加载节点施加输入转矩t1，并约束输出转矩加载节点的轴向转动自由度，以消除传动系统的刚体自由度，采用newton-raphson方法对包含非线性轴承单元的传动系统缩聚模型进行非线性迭代求解，可以求得传动系统在对应载荷工况下的静平衡状态，同时获得平衡状态下式(1)所示的非线性轴承单元的刚度矩阵kb，由于该刚度矩阵是基于解析形式的非线性轴承单元在对应载荷工况下的迭代计算求得的切线刚度，所以不能直接作为线性刚度用于静力学计算。为了在后续建立的传动系统有限元接触分析模型中实现对轴承刚度的模拟，本发明提出一种线性等效轴承刚度矩阵kbe，用来等效式(1)所示的非线性轴承单元的刚度矩阵kb在传动系统静平衡时的状态，线性等效轴承刚度矩阵kbe表示为式(7)：由式(7)可知，线性等效轴承刚度矩阵kbe中仅包含对角线项，各对角线项表示轴承在该方向所受载荷与变形之比，即将轴承等效为各方向独立的线性弹簧，且满足式(8)：由式(8)可知，线性等效轴承刚度矩阵kbe对应轴承变形δb产生的轴承载荷fb与非线性轴承单元迭代平衡时产生的载荷完全相同，即在对应的载荷工况下，线性等效轴承刚度矩阵kbe在传动系统模型中所起的作用与平衡时的非线性轴承单元的刚度矩阵kb相同。4)建立包含线性等效轴承刚度的传动系统有限元接触分析模型：在abaqus软件中建立传动系统有限元接触分析模型，包括：输入轴、输出轴、箱体、主动齿轮和从动齿轮。其中，主动齿轮和从动齿轮的有限元接触分析模型根据齿轮设计参数建立，输入轴、输出轴和箱体的有限元接触分析模型与步骤1)中建立的体单元有限元模型相同，并同样在输入轴、输出轴、箱体的有限元接触分析模型中包含连接关系和加载的位置建立边界节点，具体如下：①输入轴：轴承中心节点、输入转矩加载节点；②输出轴：轴承中心节点、输出转矩加载节点；③箱体：轴承中心节点。与步骤1)相同，分别用刚性连接单元将边界节点与各部件的有限元接触分析模型连接，并根据实际情况在箱体的有限元接触分析模型施加外部约束。因为有限元接触分析模型中包含真实设计参数对应的轮齿模型，齿轮的啮合关系通过齿面接触定义体现，所以不再建立步骤1)中的齿轮等效啮合节点。将建立的主动齿轮和从动齿轮的有限元接触分析模型调整到某一啮合位置，并在有啮合关系的齿面之间定义面-面接触关系和齿面摩擦系数。主动齿轮的有限元接触分析模型与输入轴的有限元接触分析模型固连，从动齿轮的有限元接触分析模型与输出轴的有限元接触分析模型固连。用步骤3)中求得的线性等效轴承刚度矩阵kbe耦合传动轴系的轴承中心节点与桥壳的轴承中心节点，即在abaqus软件中，依次在配对的轴承中心节点之间建立bushing单元，并将刚度参数定义为kbe。加载与约束和步骤2)相同，输入转矩t1施加在输入转矩加载节点，并约束输出转矩加载节点的轴向转动自由度。利用abaqus软件能够求得相啮合的各对轮齿之间的等效接触力矢量fi和等效作用点坐标ri，分别表示为式(9)和式(10)：fi＝[fxi,fyi,fzi]t(9)ri＝[xi，yi,zi]t(10)上式中，fxi,fyi,fzi为第i对轮齿之间的等效接触力在全局坐标系中的力分量；xi,yi,zi为第i对轮齿之间的等效作用点在全局坐标系中的坐标分量。总的齿轮啮合力矢量fm表示为式(11)：fm＝[fmx,fmy,fmz]t(11)上式中，n为齿轮副之间发生啮合的轮齿对数。齿轮的等效啮合力作用方向单位矢量n表示为式(12)：n＝[nmx,nmy,nmz]t(12)上式中，nmx＝fmx/|fm|；nmy＝fmy/|fm|；nmz＝fmz/|fm|。等效啮合节点在全局坐标系中的坐标位置表示为式(13)：r＝[xm，ym,zm]t(13)上式中，xm，ym,zm为坐标分量。由式(14)～式(16)所示的力的等效和力矩平衡关系可以求得：zm＝(mmy+fmzxm)/fmx(15)ym＝(mmx+fmyzm)/fmz(16)上式中，mmy和mmx分别为齿轮啮合力对全局坐标系x轴和y轴的力矩。5)齿轮接触状态与非线性轴承刚度的平衡迭代计算：因为齿轮接触状态与非线性轴承刚度之间存在非线性耦合，需要进行平衡迭代计算，过程如下：①以采用donley等提出的理论解析公式求得的齿轮等效啮合点坐标r0和等效啮合力作用方向单位矢量n0作为初值，按照步骤2)中的方法建立包含非线性轴承单元和齿轮等效啮合单元的传动系统缩聚模型，作为迭代计算的初始模型，并基于步骤3)中的式(7)计算该初始模型在对应载荷工况下对应的线性等效轴承刚度②按照步骤4)中的方法建立传动系统有限元接触分析模型，将①中求得的线性等效轴承刚度用于接触分析模型中，计算齿轮在该轴承刚度下的接触状态，并基于步骤4)中的公式(9)～式(16)，计算新的齿轮等效啮合点坐标r1和等效啮合力作用方向单位矢量n1；③将②中求得的齿轮等效啮合点坐标r1和等效啮合力作用方向单位矢量n1用于①中的齿轮等效啮合单元，再次基于步骤3)中的式(7)计算新的线性等效轴承刚度④重复①～③，当相邻两次迭代求得的等效啮合参数满足容差时，计算收敛，收敛条件表示为式(17)：||nk-nk-1||+||rk-rk-1||<ε(17)上式中，nk-1和nk分别为第k-1次和第k次迭代求得的齿轮等效啮合力作用方向单位矢量；rk和rk-1分别为第k-1次和第k次迭代求得的齿轮等效啮合节点坐标；ε为收敛容差，为一个较小的正数，实际计算时收敛容差取为0.01。上述过程可以准确求得任意齿轮啮合方位下传动系统的平衡状态，并获得整个传动系统各部件的受力与变形结果，通过改变齿轮的方位，还可以模拟完整的齿轮啮合过程，从而对传动系统进行更为全面的计算分析。本发明由于采取以上技术方案，其具有以下优点：1、本发明提供了一种非线性轴承刚度的线性等效方法，将传动系统缩聚模型静平衡时的轴承切线刚度等效为线性刚度，并用于传动系统有限元接触分析模型，能够在齿轮有限元接触计算时准确考虑轴承刚度的影响。2、本发明综合了有限元接触分析模型和包含非线性轴承单元和齿轮等效啮合单元的缩聚模型在传动系统数值模拟中的优势，利用齿轮等效啮合参数和等效轴承刚度建立起两种模型的耦合关系，通过平衡迭代，准确实现了考虑齿轮和轴承非线性耦合的传动系统数值模拟。3、本发明能够准确计算完整传动系统和其中任意部件在不同载荷工况下的受力和变形特点，为传动系统的设计建模和分析校核提供有效手段。附图说明图1是本发明方法的流程示意图；图2是驱动桥主减速器传动系统模型平面示意图；图3是主动齿轮轴有限元模型示意图；图4是差速器壳有限元模型示意图；图5是桥壳有限元模型示意图；图6是驱动桥有限元接触分析模型；图7是轴承2的平动自由度方向的线性等效刚度随齿轮转角变化曲线。具体实施方式下面结合附图和实施例对本发明进行详细的描述。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。以如图2所示的驱动桥主减速器准双曲面齿轮传动系统为例，主动齿轮轴为输入端，差速器壳为输出端，全局坐标系原点o定义为差速器中心，主动齿轮轴线与全局坐标系x轴平行，从动齿轮轴线与全局坐标系y轴平行，主动齿轮沿全局坐标系z轴偏置，准双曲面齿轮参数如表1所示。表1准双曲面齿轮参数如图1所示，本发明提供的考虑齿轮和轴承非线性耦合的传动系统数值模拟方法，其包括以下步骤：1)建立传动系统的体单元有限元模型：建立主动齿轮轴、差速器壳和桥壳的体单元有限元模型，分别如图3～图5所示，单元类型均为四节点四面体单元，在各部件有连接关系的位置建立边界节点：①主动齿轮轴模型有5个边界节点：3个轴承中心节点(bi1,bi2,bi3)、1个输入转矩加载节点(l1)、1个主动齿轮等效啮合节点(m1)；②差速器壳模型有4个边界节点：2个轴承中心节点(bi4,bi5)、1个输出转矩加载节点(l2)、1个从动齿轮等效啮合节点(m2)；③桥壳模型有5个边界节点：轴承中心节点(bo1,bo2,bo3,bo4,bo5)。先采用理论解析公式计算准双曲面齿轮的初始等效啮合参数，求得齿轮等效啮合节点的坐标位置r0为(-181.025,-18.145,24.931)，等效啮合力作用方向单位矢量n0为(-0.7062,-0.2466,-0.6637)，初始模型的主动齿轮等效啮合节点和从动齿轮等效啮合节点的坐标位置均设置为r0。分别用刚性连接单元将边界节点与各部件体单元模型上的节点耦合，并约束桥壳的板簧座上表面节点的平动自由度，模拟实际钢板弹簧对桥壳的约束。2)建立包含非线性轴承单元和齿轮等效啮合单元的驱动桥缩聚模型：采用guyan缩聚法对各部件的体单元有限元模型进行缩聚变换，仅保留边界节点自由度，计算各部件的缩聚刚度矩阵，其中，主动齿轮轴缩聚刚度矩阵ki包含30个自由度，差速器壳刚度矩阵ko包含24个自由度，桥壳刚度矩阵kh包含30个自由度。采用解析形式的非线性轴承单元模拟滚动轴承，轴承单元非线性切线刚度矩阵分别为kb1,kb2,kb3,kb4,kb5，轴承单元一端连接传动轴系的轴承中心节点，另一端连接桥壳的轴承中心节点，即如图2所示：依次用kbj耦合边界节点bij和boj，j为轴承编号，j＝1～5。采用等效啮合单元模拟齿轮的啮合关系，等效啮合单元一端连接主动齿轮等效啮合节点，另一端连接从动齿轮等效啮合节点，即在图2中，用齿轮等效啮合刚度矩阵km耦合边界节点m1和m2。采用matlab编程的方法将上述缩聚刚度矩阵ki,ko,kh、非线性轴承刚度矩阵kb1,kb2,kb3,kb4,kb5、齿轮等效啮合刚度矩阵km按照连接关系组集，获得驱动桥系统缩聚模型的刚度矩阵k，驱动桥缩聚模型共有84个自由度。3)传动系统缩聚模型的非线性静力学求解和线性等效轴承刚度计算：在主动齿轮轴的输入转矩加载节点l1施加输入转矩t1＝1616n·m，并约束差速器壳输出转矩加载节点l2的轴向转动自由度，采用newton-raphson方法进行非线性迭代求解，迭代7次计算收敛，耗时10秒，求得驱动桥系统在对应载荷工况下的静平衡状态，并求得线性等效轴承刚度矩阵kbe1,kbe2,kbe3,kbe4,kbe5，如表2～表6所示：表2轴承1的线性等效轴承刚度矩阵kbe1δx/mmδy/mmδz/mmθx/radθy/radθz/radfx/n000000fy/n04.12e+050000fz/n004.05e+05000mx/n·mm000000my/n·mm00001.31e+070mz/n·mm000005.80e+06表3轴承2的线性等效轴承刚度矩阵kbe2δx/mmδy/mmδz/mmθx/radθy/radθz/radfx/n6.50e+0500000fy/n01.03e+060000fz/n004.77e+06000mx/n·mm000000my/n·mm00003.20e+090mz/n·mm00000-1.34e+09表4轴承3的线性等效轴承刚度矩阵kbe3δx/mmδy/mmδz/mmθx/radθy/radθz/radfx/n-2.08e+0400000fy/n05.58e+040000fz/n007.79e+04000mx/n·mm000000my/n·mm00008.35e+070mz/n·mm000001.81e+08表5轴承4的线性等效轴承刚度矩阵表6轴承5的线性等效轴承刚度矩阵kbe5δx/mmδy/mmδz/mmθx/radθy/radθz/radfx/n4.93e+0500000fy/n0-3.76e+040000fz/n004.79e+05000mx/n·mm000-1.87e+0900my/n·mm000000mz/n·mm00000-6.49e+094)建立包含等效轴承刚度的齿轮传动系统有限元接触分析模型：在abaqus软件中建立齿轮传动系统有限元接触分析模型，包括：主动齿轮轴、差速器壳、桥壳、主动齿轮、从动齿轮。主动齿轮轴、差速器壳和桥壳的有限元接触分析模型与步骤1)建立的体单元有限元模型相同，同样在各部件有连接关系的位置建立边界节点，其中：①主动齿轮轴模型有4个边界节点：3个轴承中心节点(bi1,bi2,bi3)、1个输入转矩加载节点(l1)；②差速器壳模型有3个边界节点：2个轴承中心节点(bi4,bi5)、1个输出转矩加载节点(l2)；③桥壳模型有5个边界节点：轴承中心节点(bo1,bo2,bo3,bo4,bo5)。与步骤1)相同，分别用刚性连接单元将边界节点与各部件有限元接触分析模型上的节点耦合，并约束桥壳的板簧座上表面节点的平动自由度，以模拟实际钢板弹簧对桥壳的约束作用。因为有限元接触分析模型中包含真实设计参数对应的轮齿模型，所以不再建立步骤1)中的齿轮等效啮合节点。建立主动齿轮和从动齿轮的有限元接触分析模型，单元类型为六面体单元，如图6所示，将齿轮副调整到某一啮合位置，在有啮合关系的齿面之间定义面-面接触关系，并定义齿面摩擦系数为0.15。实际建模计算时，综合考虑分析要求、计算精度和计算成本，可以仅对参与接触计算的轮齿进行建模，并将齿面网格尺寸控制在1mm以内，以保证接触计算精度和收敛性。主动齿轮的有限元接触分析模型与主动齿轮轴的有限元接触分析模型之间采用tie单元固连，从动齿轮的有限元接触分析模型与差速器壳的有限元接触分析模型之间也采用tie单元固连。用步骤3)中求得的线性等效轴承刚度kbe耦合传动轴系的轴承中心节点与桥壳的轴承中心节点，即在abaqus软件中，依次在节点bij和boj之间建立bushing单元，并将刚度参数定义为kbej，j为轴承编号，j＝1～5。加载与约束和步骤2)相同，输入转矩t1施加在输入转矩加载节点，并约束输出转矩加载节点的轴向转动自由度。5)齿轮接触状态与非线性轴承刚度的平衡迭代计算：因为齿轮接触状态与非线性轴承刚度之间存在非线性耦合，需要进行平衡迭代计算，过程如下：①以采用理论解析公式求得的齿轮等效啮合参数r0、n0作为初值，建立包含非线性轴承单元和齿轮等效啮合单元的传动系统缩聚模型，作为迭代计算的初始模型，并计算该模型在对应载荷工况下对应的线性等效轴承刚度②建立传动系统有限元接触分析模型，并将①中求得的等效轴承刚度用于有限元接触分析模型中，计算齿轮在该轴承刚度下的接触状态和新的齿轮等效啮合参数r1、n1；③将②中求得的齿轮等效啮合参数r1、n1用于①中的齿轮等效啮合单元，计算新的线性等效轴承刚度④重复①～③，当相邻两次迭代求得的等效啮合参数满足容差时，计算收敛。实际计算时收敛容差取为0.01，共迭代4次收敛，每次迭代求得的齿轮等效啮合参数如表7和表8所示，其中迭代次数0对应的是采用理论解析公式求得的初始结果，由于未考虑齿面摩擦且无法体现真实齿轮接触状态的影响，理论解析公式计算结果与有限元接触计算结果存在明显差别。表7齿轮等效啮合节点坐标迭代次数x坐标/mmy坐标/mmz坐标/mm0-181.025-18.14524.9311-179.246-14.40731.2722-179.479-14.40931.223-179.491-14.40731.2224-179.492-14.40631.225表8齿轮等效啮合力作用方向单位矢量迭代次数x分量y分量z分量0-0.7062-0.2466-0.66371-0.6451-0.1685-0.74522-0.6477-0.1624-0.74433-0.6481-0.1617-0.74424-0.6481-0.1616-0.7442进一步对比迭代收敛后两种模型的主动齿轮轴和差速器壳轴承中心边界节点的位移计算结果如表9和表10所示，计算结果吻合，说明两种模型的受力和变形状态一致，验证了驱动桥系统静平衡状态求解的正确性。表9缩聚模型边界节点位移边界节点x位移/μmy坐标/μmz坐标/μmbi1-144-44bi2-13-7-57bi3-14-12-66bi4263013bi5173116表10有限元接触分析模型边界节点位移边界节点x位移/μmy坐标/μmz坐标/μmbi1-154-43bi2-13-7-58bi3-14-11-65bi4263013bi5173116上述过程可以求得任意齿轮啮合方位下的系统的平衡状态，通过改变齿轮的方位，还可以求得完整齿轮啮合过程的静平衡状态，并获得齿轮啮合过程中的轴承刚度变化曲线，以图2中的轴承2为例，求得其平动自由度等效刚度随齿轮转角变化曲线如图7所示，轴承刚度呈周期性变化趋势，其中z向刚度波动较为明显。综上所述，本发明提出的方法能够准确实现考虑齿轮和轴承非线性耦合的传动系统数值模拟，可广泛应用于驱动桥、变速箱等各类机械传动系统的设计分析。上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。当前第1页12