一种基于装配特征信息的矩形花键接触动力学仿真方法与流程

本发明涉及仿真技术领域,尤其是一种基于装配特征信息的矩形花键接触动力学仿真方法。

背景技术：

花键是一种重要的传动联接形式，具有很高的可靠性，广泛用于轴类、齿轮类、联轴器等机械零件间的动力传递。传统花键的设计往往基于设计者的经验和国际标准，在设计过程中设置了很大的安全裕量。设计过程中也只对花键所受的静载荷进行计算，对动载荷则采用经验系数法进行估计。这样得到的结果会使花键设计过于保守，导致产品笨重，浪费材料。对于动载荷的经验估计也可能使花键在承受大冲击载荷时失效，使产品损坏。因此，在对产品设计要求日益精确的现代，提出一种能够分析花键动态分析的算法是十分必要的。

国内外对花键动态性能做了一些研究工作。国外一些研究者将花键简化成扭转弹簧，利用有限元分析和台架试验确定了花键的扭转刚度、阻尼系数和花键工作时表面接触状态；一些研究者分析了花键误差对花键动态性能的影响，推导出了含有花键误差的多体动力学模型，国内的一些研究人员利用集中质量模型，建立具有齿侧间隙的花键联接四自由度非线性动力学方程，计算了花键联接实际接触齿对数与综合时变啮合刚度，采用四阶runge-kutta法求解花键联接非线性动力学方程。这些研究对花键整体动态性能研究有一定帮助，但没有考虑装配特征信息的影响。

目前花键设计研究的几个关键技术问题：(1)花键配合公差与花键齿形公差的选择。公差选择的优良与否，不仅影响到花键的制造成本与装配成本，也影响到产品的工作质量以及花键寿命，但目前花键设计的公差选择很大程度上仍然取决于经验。花键配合公差与花键齿形公差是影响花键装配与传动最重要的两个装配特征信息，其对花键寿命和花键动态特性的影响都有待进一步研究。(2)花键接触研究。材料接触力方程一直是科研工作者关注的课题。通过合理的假设，建立材料本构方程，利用实验测定方程参数，计算材料应力，则可得到材料应变与材料所受外力的直接关系。合理的接触力计算方程，对精确计算花键接触力、分析花键受力情况以及预测花键寿命起十分关键的作用。如何建立精确的花键键齿接触与花键定心面的接触力计算方程，目前仍需深入研究。

随着产品生产要求的不断提高，国际标准化组织提出了新一代gps。gps是针对产品全生命周期中几何特征方面建立的一个几何技术标准体系，包含产品从宏观几何特征到微观几何特征的几何技术标准。传统的gps基于几何学理论构建，分别形成了尺寸公差、形位公差和表面特征系列标准，但这些标准之间不完全协调、需要整合。除此之外，gps更主要的问题在于设计、制造、检验流程中标准依据的基础理论的不统一，无法在流程之间实现几何信息的准确表达和沟通。新一代gps标准体系着重于提供一个适宜于cax集成环境的、更加清晰明确的、系统规范的几何公差定义和数字化设计、计量规范体系。根据gps所提供的理论依据，可以将尺寸公差、形位公差转化为零件尺寸、位置上的变动量，将花键装配信息加入仿真模型中。

赫兹接触力理论研究两物体因受压接触后产生的局部应力和应变分布规律。1881年h.r.赫兹最早研究了玻璃透镜在接触力作用下发生的弹性变形。他假设：①接触区发生小变形。②接触面呈椭圆形。③相接触的物体可被看作是弹性半空间，接触面上只作用有分布的垂直压力。凡满足以上假设的接触称为赫兹接触。罚函数法视物体为不可变形的刚体，物体间可以发生微小穿透。当物体发生穿透时，则在从点和靶点之间引入一个较大的法向接触力，法向接触力等于从点相对于其靶点的穿透量和表征接触碰撞界面刚度的罚因子的乘积，从而使穿透部分自动回到平衡位置。国内外学者利用赫兹接触理论和罚函数法分析了许多接触问题，如球轴承运行过程中轴承内外环、轴承滚珠和保持架之间的接触状态分析、圆柱铰链的接触状态分析等，所得的分析结果与实验结果都具有很大的贴合性。利用赫兹接触力理论和罚函数法，可以分析花键定心面处和矩形花键齿侧处的接触问题。

现有技术一的技术方案

考虑了矩形花键齿侧间隙误差与径向间隙的花键动力学模型。

一般情况下花键的齿侧间隙是存在误差的，内外花键的每一对齿i在前(受载)、后(非受载)齿面都有不同的间隙值，以bi和b′i表示，如图1所示。当花键副加载后，齿侧间隙最小的键齿对最先接触，随后各键齿对以齿侧间隙的大小顺序依次开始接触。如果传递的转矩平均值所造成的键齿变形小于最大齿侧间隙时，花键齿中只有一部分键齿接触。由此可以得出一个分段线性的形变方程，方程需满足接触的键齿对数小于等于花键齿数。

式(1)为图1所示花键的动力学方程。

式中i是内花键所在零件的转动惯性矩，ct是总体转矩粘性阻尼系数，kt是单个键齿接触所产生的扭转刚度，θ(t)是内花键相对于外花键的角位移，t(t)是施加在内花键所在零件上的转矩。非线性位移方程g[θ(t)]用来代表花键间隙误差导致的刚度非线性。引入符号ωn²＝kt/i和ζ＝ct/(2iωn)，动力学方程可以写成如下的标准格式：

为了便于理解，可以按图2所示做等效替换(rp为花键的分度圆半径)，将模型中的旋转量替换为等效的平移量。替换完成后，模型转变成了两个相互啮合的齿条，齿条上带有n个短悬臂梁齿，一个齿条代表内花键所连零件，另一个齿条代表花键轴，视花键轴为相对静止。用平移位移作为新的坐标后，式(2)可写为

式(3)中的位移方程g[y(t)]是图2中给出的n个独立的分段线性位移方程的集合，图2中gi(y)代表一个处于接触状态的键齿的变形。设外花键齿宽为p，内花键齿槽宽为p，键齿i在右侧(受载)和左侧(非受载)的间隙满足

bi+b′i＝p-p,i＝1,2,...,n.(4)

由于是一维系统，花键齿可以按照间隙大小由大到小排列，例如

bn≥bn-1≥…≥b3≥b2≥b1＝a(5)

通过式(4)和(5)，左侧的齿侧间隙满足条件

b′1≥b′2≥b′3≥…≥b′n-1≥b′n＝a(6)

根据中心极限定理，花键间隙的分布可以用正态分布近似替代，分布域的宽度是由花键的制造精度决定的，可用极值法计算出的花键最大间隙与最小间隙之差。图2所示的位置为初始相对位置，两侧均有相等的最小间隙a。定义βi＝bi-a和β′i＝b′i-a,i＝1,2,...,n，可以得到

因此，当|y(t)|＜a时，没有键齿发生接触，式(3)中的g[y(t)]＝0。同样的，当(βn+a)≤y(t)或者-(β′1+a)≥y(t)时，所有键齿发生接触，g[y(t)]呈线性变化。

当y(t)不是以上两种情况时，受载的齿数r小于花键齿数n，因此是一个分段线性的方程。式(8)为所有y(t)值g[y(t)]的表达式。

再引入无量纲时间变量τ＝ωnt和无量纲位移变量x(t)＝y(t)/a，动力学方程(3)可写成一种无量纲形式

式中，α′i＝β′i/a,i＝1,...,n.

通过式(9)，即可得到花键在存在齿侧间隙的情况下的旋转运动规律和键齿的受力情况。

花键的径向间隙会使花键存在静态不对中(e0，φ0)和动态不对中。静态不对中为花键装配产生的内外花键轴心偏移(x，y)，动态不对中是在工作过程中内外花键轴心的动态振动位移。为方便模型计算，统称二者为不对中，实际不对中为：

不对中影响了花键键齿的等效啮合距离为li，而单个键齿刚度是等效啮合距离li的函数，因此花键不对中会影响花键动力学性能。为了分析方便，先假设不对中发生在x正向，即φ0＝0，y＝0；此时内花键与外花键的啮合是不均匀的，有的键啮合紧密，有的键啮合疏松；两个圆心不重合，内花键各个键的等效啮合距离不相等，如图3所示。

设实际的不对中方向上即x轴正向上键的编号为1，逆时针依次为2、3……z；各键与x轴正向的夹角为：

各键等效啮合距离为：

因此，不对中花键各键齿啮合刚度为：

ki＝f(li)(14)

当花键传递扭矩且该扭矩很大时，各键发生变形，即扭转使各键产生一个啮合力；此时若两个花键所在零件各自随系统振动，则二者之间有相对位移，导致各键又产生一个变形，即动态振动位移也使各键产生一个啮合力。

(1)扭转产生的啮合力

当花键传递扭矩为t时，

其中，r为花键齿根圆的半径；fti为：

fti＝λi·ki＝φli·ki(16)

其中，λi为各个键的变形，φ为各个键变形产生的扭转角位移，由于各个键构成的花键整体在静态常扭矩t的作用下一起扭转，引起各键的扭转角位移相同。

整理以上推导过程，可以得到：

因此，扭转角位移φ为：

则每个键传递的力：

fti＝φliki(19)

(2)动态振动位移产生的啮合力

其中，

综合以上分析，花键由于不对中而导致的啮合力为：

将各键产生的啮合力沿坐标方向分解：

其中，θi为每个键作用力方向与x轴正向的夹角，根据以上推导可知：

则花键不对中而产生的x、y向合力为：

当实际不对中发生在任意角度时，相当于上述不对中绕原点逆时针旋转φ角，φ角定义同前，如图4所示。

此时，花键不对中啮合力为：

将(25)代入花键动力学方程中，即可得到花键在不对中啮合力的作用下，花键内外轴心的相对平移运动规律及各键齿的啮合力变化。

现有技术一的缺点

1.两种算法只考虑了一个方向的间隙(齿侧方向或径向方向)，

没有将两个方向的间隙统一在一个模型下。

2.两种算法只是针对间隙的尺寸及其误差进行分析，没有上升到公差标准层面，无法直接对花键设计进行指导。

因此，对于上述问题有必要提出一种基于装配特征信息的矩形花键接触动力学仿真方法。

技术实现要素：

针对上述现有技术中存在的不足，本发明的目的在于提供一种基于装配特征信息的矩形花键接触动力学仿真方法，目的在于解决在花键设计过程中，传统仿真无法考虑公差信息等装配特征信息的问题。

一种基于装配特征信息的矩形花键接触动力学仿真方法，包括如下步骤：

第一步：分析花键装配特征信息；

第二步：建立理想花键几何外形模型；

第三步：添加装配特征信息；

第四步：建立接触动力学模型；

第五步：仿真结果求解。

优选地，所述的分析花键装配特征信息步骤，根据花键装配过程中花键定位、花键配合信息，确定建模过程中所需的花键装配特征信息。

优选地，所述的建立理想花键几何外形模型步骤，根据矩形花键的形状外形特点，建立描述花键外形特征点与特征面理想位置的二维数学模型。

优选地，所述的添加装配特征信息步骤，将第一步所分析的装配特征信息，根据gps理论中对于误差与公差的关系、几何空间中误差分布规律将公差信息添加到花键对应的几何尺寸上，根据中心极限定理和正态分布模型得到误差值。

优选地，所述的建立接触动力学模型步骤，根据赫兹接触力学理论和罚函数法，计算接触刚度和阻尼系数，以外花键特征点穿透内花键特征面的垂直距离为齿侧穿透深度，以内外花键小径定心面从点相对于靶点的穿透深度作为径向穿透深度，建立接触动力学模型。

优选地，所述的仿真结果求解步骤，利用牛顿-欧拉迭代算法，在每一迭代步中，根据当前内外花键的状态计算花键齿间与小径定心面间的接触力，将接触力代入动力学方程中，解算出下一迭代步中各物体的加速度、速度和位移，重复迭代步骤，直到系统进入平衡状态或系统失稳，记录迭代过程中产生的键齿接触力，花键角位移、角速度与角加速度。

优选地，建立基于牛顿-欧拉动力学方程的差分动力学方程，利用穿透深度、接触刚度和阻尼系数来计算接触力并用不动点迭代法解算动力学方程，其每一步解算结果即为当前时刻物体的运动状态，建立差分动力学方程如下式所示

式中，t与t+1代表当前步和下一步，δt代表仿真步的时间步长，ax，ay和αz是沿x轴、y轴的加速度以及沿z轴的角加速度。vx，vy和ωz是沿x轴、y轴的速度以及沿z轴的角速度，px，py和θz是沿x轴、y轴的位移以及沿z轴的角位移。fx，fy和mz是沿x轴、y轴的合力以及沿z轴的合力矩，由下式计算

式中，fax、fay和maz代表x、y方向的主动力和z方向的主动力矩，fcontactx、fcontacty和mcontactz代表因接触约束产生的在x、y方向的约束力和z方向的约束转矩，通过不断迭代求解动力学方程，算出每一步的花键位移、速度、加速度和接触力，即可得到花键在初始条件下每一时刻的运动状态。

由于采用上述技术方案，本发明有益效果：本发明将矩形花键的装配特征信息纳入仿真模型的考虑范围，利用接触动力学对花键工作状态进行仿真，仿真结果可以得到花键配合公差和花键齿形公差对花键动态特性的影响，得到不同公差选择对花键键齿接触力的影响，仿真结果更加符合真实结果，可以对冲击载荷作用下花键寿命进行分析，通过本算法进行仿真计算，可以使设计结果更加贴近现实，减少安全裕量，避免冲击失效，为矩形花键规范化设计提供指导。

附图说明

图1是现有技术的花键前后齿面的齿侧间隙示意图；

图2是现有技术的花键等效位移模型；

图3是现有技术的不对中花键等效啮合距离示意图；

图4是现有技术的任意方向不对中花键等效啮合距离示意图；

图5-6是本发明的矩形花键键齿(槽)位置度公差(a)内花键(b)外花键；

图7是本发明的花键二维平面模型；

图8是本发明的理想花键外形的几何建模；

图9是本发明的添加了装配特征信息的花键几何外形；

图10是本发明的判断碰撞的辅助向量；

图11是本发明的内外花键小径定心面的理想位置；

图12是本发明的发生碰撞的内外花键小径定心面；

图13是本发明的仿真迭代求解的流程图；

图14是本发明的方法总流程图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。

如图5并结合图6至图13所示，一种基于装配特征信息的矩形花键接触动力学仿真方法，包括如下步骤：

第一步：分析花键装配特征信息；第二步：建立理想花键几何外形模型；

第三步：添加装配特征信息；第四步：建立接触动力学模型；第五步：仿真结果求解。

进一步的，所述的分析花键装配特征信息步骤，根据花键装配过程中花键定位、花键配合信息，确定建模过程中所需的花键装配特征信息，所述的建立理想花键几何外形模型步骤，根据矩形花键的形状外形特点，建立描述花键外形特征点与特征面理想位置的二维数学模型。

进一步的，所述的添加装配特征信息步骤，将第一步所分析的装配特征信息，根据gps理论中对于误差与公差的关系、几何空间中误差分布规律将公差信息添加到花键对应的几何尺寸上，根据中心极限定理和正态分布模型得到误差值。

进一步的，所述的建立接触动力学模型步骤，根据赫兹接触力学理论和罚函数法，计算接触刚度和阻尼系数，以外花键特征点穿透内花键特征面的垂直距离为齿侧穿透深度，以内外花键小径定心面从点相对于靶点的穿透深度作为径向穿透深度，建立接触动力学模型。

进一步的，所述的仿真结果求解步骤，利用牛顿-欧拉迭代算法，在每一迭代步中，根据当前内外花键的状态计算花键齿间与小径定心面间的接触力，将接触力代入动力学方程中，解算出下一迭代步中各物体的加速度、速度和位移，重复迭代步骤，直到系统达到平衡状态，记录迭代过程中产生的键齿接触力，花键角位移、角速度与角加速度。

第一步：分析花键装配特征信息，影响花键动态特性的装配特征信息有花键的配合公差与花键键齿的尺寸公差和形位公差。配合公差是约束花键在装配过程中内外花键松紧程度一种规范，分为间隙配合、过渡配合和过盈配合，不同的配合方式对花键动态特性都有不同的影响。对于矩形花键而言，其配合分为小径配合、大径配合和齿侧配合，一般情况下都使用小径配合。花键键齿的尺寸公差和形位公差影响内外花键之间的齿侧间隙，齿侧间隙的不均匀性也影响花键传动的动态特性，花键键齿的位置度公差如图5和图6所示。经过分析，选取矩形花键的小径配合公差和花键键齿的尺寸公差与位置度公差作为考虑的装配特征信息。

第二步：建立理想花键几何外形模型，图7所示为花键二维平面模型。a1、a2、b1、b2、b3、b4为内外花键上的特征点，其中a1、a2为外花键齿顶面两端端点，b1、b2为内花键齿槽顶部两端端点，b3、b4为内花键齿槽底部两端端点。通过获取特征点的精确位置计算花键穿透时的齿侧穿透深度，建立如图8所示的理想花键的几何模型。图中，ri'为花键第i个齿对应的位置向量，向量方向由花键质心指向第i个齿的齿顶中点，rib'为键齿(槽)侧面相对于齿顶中点的偏移矢量，其值为实际键(槽)宽b的一半，rih'为齿(槽)的偏移矢量。如此，建立了理想花键的几何外形模型。

第三步：添加装配特征信息，根据步骤1的分析，以及gps和中心极限定理，将花键的配合公差、键齿尺寸公差和位置度公差添加到理想花键几何外形模型中，如图9所示。

图9中，点划线代表理想键齿(槽)轮廓，虚线代表只包含位置度误差ep的键齿(槽)轮廓，实线代表键齿(槽)的规范表面模型。rip'为键齿(槽)相对于理想位置的位置偏差矢量，rit'为齿(槽)侧面相对于理想轮廓的偏差矢量，b’和b”为理想键(槽)宽，b为带有尺寸误差的实际键(槽)宽，尺寸误差用eb表示。d和d’分别表示非理想和理想的小径直径，两者之差为小径的尺寸误差ed。gps中认为误差在公差带上呈正态分布，因此ep、eb、ed的取值在各自的公差带上满足正态分布。

第四步：建立接触动力学模型。为了判断特征点a2是否与内花键左侧齿侧发生碰撞，引入两个辅助向量v1、v2，如图10所示。

图10中，向量v1从特征点b4指向b2，向量v2从特征点b4指向a2。当a2在b4b2的右侧时，a2不与齿侧发生碰撞，此时向量v1和v2的叉乘积符号为负。当发生碰撞时，a2与b4b2重合或者处于b4b2左侧，v1和v2的叉乘积符号为正。相反方向的碰撞检测可以按相同的方法推导。齿侧穿透深度的计算表达式为

引入内外花键相对角速度ω'＝ωint-ωext，相对穿透速度的计算表达式为

式中，为特征点a2相对于固连在花键相对坐标系的矢径，为与向量v1垂直的向量，其方向指向内花键里侧。

因此，结合接触力计算公式，键齿i在圆周方向上产生的接触力为

式中，kc为单键啮合刚度，dc为阻尼系数。

在径向方向，内外花键的小径定心面碰撞问题可以简化为含有间隙的圆柱铰的碰撞问题。

如图11所示，内外花键在不受力的状态下处于完全同心的理想状态，此时内外花键的小径定心面之间存在间隙c，其值为内外花键小径直径之差的二分之一，即

c＝ri-re＝(di-de)/2

式中，ri与re分别为内外花键小径定心面的半径。

当花键开始受载时，内外花键产生偏心距eij，如图12所示。eij的方向从内花键小径圆心指向外花键小径圆心，其值为内外花键小径圆心的距离。当eij的值小于c时，内外花键的小径定心面没有发生碰撞。当eij大于等于c时，内外花键小径定心面发生碰撞。定义qi为外花键小径定心面上的碰撞点，qj为内花键上对应的碰撞点，和分别为碰撞点在惯性坐标系中的矢量坐标。定义径向穿透深度

通过几何关系分析，可以得到径向穿透深度的方向与eij相同，其大小等于eij与c的差值，即

径向方向的相对穿透速度的方向与eij相同，其表达式为

因此，根据hertz接触理论，内外花键的小径定心面的接触力为

式中，δr是内外花键小径定心面的穿透深度，是相对碰撞速度，是初始碰撞速度，ec是恢复系数，kr是碰撞刚度。

第5步：仿真结果求解

建立基于牛顿-欧拉动力学方程的差分动力学方程，利用穿透深度、接触刚度和阻尼系数来计算接触力并用不动点迭代法解算动力学方程，其每一步解算结果即为当前时刻物体的运动状态。建立差分动力学方程如下式所示

式中，t与t+1代表当前步和下一步。δt代表仿真步的时间步长。ax，ay和αz是沿x轴、y轴的加速度以及沿z轴的角加速度。vx，vy和ωz是沿x轴、y轴的速度以及沿z轴的角速度，px，py和θz是沿x轴、y轴的位移以及沿z轴的角位移。fx，fy和mz是沿x轴、y轴的合力以及沿z轴的合力矩，由下式计算

式中，fa和ma代表主动力和主动力矩，fcontact和mcontact代表接触约束产生的约束力和约束转矩。

通过不断迭代求解动力学方程，算出每一步的花键位移、速度、加速度和接触力，即可得到花键在初始条件下每一时刻的运动状态，仿真迭代求解的运算流程图如图13所示

本发明将矩形花键的装配特征信息纳入仿真模型的考虑范围，利用接触动力学对花键工作状态进行仿真，仿真结果可以得到花键配合公差和花键齿形公差对花键动态特性的影响，得到不同公差选择对花键键齿接触力的影响，仿真结果更加符合真实结果，可以对冲击载荷作用下花键寿命进行分析，通过本算法进行仿真计算，可以使设计结果更加贴近现实，减少安全裕量，避免冲击失效，为矩形花键规范化设计提供指导。

以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。