一种间隙转动副接触应力分析的方法与流程

本发明涉及计算机辅助工程分析，更具体地说，涉及一种用于在有限元分析中模拟间隙转动副接触应力的方法，该有限元分析方法用于时间推进的工程模拟，以辅助用户在改进工程产品(例如，汽车、飞机)设计时做出决策。

背景技术：

由于设计、制造、加工、磨损等因素，机械系统(例如，汽车、飞机、机床、机器人)中的转动副不可避免存在间隙。间隙使得轴承与转轴在相对运动过程中频繁发生接触碰撞，在系统中引入振动、噪声、磨损等现象，严重影响系统运行精度，缩短系统寿命。因此，分析计算含间隙运动副机构的动力学响应具有重要的学术意义和经济价值。

有限元方法是目前最为精确的接触碰撞数值模拟方法。在采用有限元方法计算接触力虚功时，需要使用数值方法计算接触力虚功的积分。使用数值积分方法一般需要在有限元参数空间分布若干数值积分点。从施加接触约束的角度，过少的积分点容易造成接触漏检，不能保证积分精度，因此一般需要布置多个数值积分点。然而由于接触边界有限元离散后节点数目有限，过多的积分点会导致系统方程过约束，使得接触边界出现应力波动现象，影响数值分析的精度。

技术实现要素：

本发明针对间隙转动副有限元模型，提出一种消除由积分点数量增加导致系统过约束的方法，可降低间隙转动副内接触应力波动。本发明通过以下技术方案实现：

一种间隙转动副接触应力分析的方法，所述方法包括以下步骤：

s1:提取轴承接触边界有限元单元，为每个单元设置高斯积分点；

s2:搜索每一轴承接触边界高斯积分点在转轴接触边界上的投影点及其参数坐标；

s3:计算所有转轴投影点处法向间距及几何度量；

s4:为转轴接触边界每个节点计算等效投影点参数坐标、法向间距及其变分，几何度量，法向及切向牵引力；

s5:计算接触力虚功。

具体的，所述s1步骤中提取轴承接触边界有限元单元还包括：对间隙运动副的转轴和轴承进行有限元网格划分。

具体的，所述s1步骤中提取轴承接触边界有限元单元为分别提取转轴和轴承接触体接触边界上的一维单元。

具体的，所述s1步骤中设置高斯积分点包括根据高斯积分准则在轴承接触边界每个单元上分布若干高斯积分点，所述高斯积分点取插值函数阶次的2倍。

具体的，所述s2步骤中的投影点对应的参数坐标满足如下非线性方程：

其中，f(·)表示函数，τ1为转轴接触边界上任一点处的协变基矢量，x¹和x²分别为转轴和轴承接触边界上任一点的位置矢量，ξ为转轴接触边界的参数坐标。

具体的，所述s3步骤中的投影点处法向间距及几何度量根据下式计算：

高斯积分点与投影点的法向间距：

其中n为投影点法向单位矢量；

投影点处几何度量用于度量空间几何的长度、面积和角度等。

具体的，所述s4步骤中的轴承接触边界每个节点a等效投影点参数坐标及其变分满足如下方程：

其中，γ为接触边界，ra为与节点a对应的nurbs插值函数。

具体的，所述接触力虚功δw通过下式计算：

其中，tna为节点a等效投影点法向牵引力，δgna为节点a等效投影点处法向间距的变分，tna为节点a等效投影点切向牵引力。

具体的，等效投影点处法向间距gna及其变分δgna可按照下式计算：

具体的，等效投影点处法向牵引力tna可按下式计算：

tna＝∈ngna

其中，∈n为法向惩罚参数；

等效投影点切向牵引力tna可以按下式计算：

其中，μ为摩擦系数，为试验节点切向牵引力，所述试验节点切向牵引力定义为：

其中下标n指上一时间步，∈t为切向惩罚参数，等效投影点处几何度量其中m11为高斯积分点对应投影点m处几何度量，即在步骤s3中已分别算出。

相对于现有技术，本发明的有益效果为：

按本发明所公开的算法计算间隙转动副内接触应力，能够获得与接触边界节点数目相同的接触力约束数目，有效避免过多积分点导致的过约束问题，有效降低接触边界的接触应力波动，提高接触应力计算精度。

附图说明

图1为发动机活塞常用的曲柄滑块机构示意图；

图2为曲柄滑块机构网格划分示意图；

图3为间隙转动副轴承模型示意图；

图4为间隙转动副转轴模型示意图；

图5为轴承与转轴接触搜索及局部运动学几何关系示意图；

图6为未使用本发明方法的接触边界区域应力情况分析结果图；

图7为使用本发明方法的接触边界区域应力情况分析结果图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

本发明所要阐述的是一种用于在间隙转动副接触有限元分析中降低由数值积分点过约束导致的接触应力波动的计算方法。

曲柄滑块机构在机械系统中极为常见，如发动机的活塞连杆机构。图1为本实施例的曲柄滑块机构示意图。曲柄滑块机构包括曲柄10、连杆20、滑块30、轴承40和转轴50。曲柄10两转轴间距1.2m，宽度0.2m。连杆20两转轴间距1.6m，宽度0.3m。滑块30长宽均为0.3m。轴承40的内径为0.06m。转轴50的直径为0.057m。所有零部件的密度均为7800kg/m³。转轴50和滑块30的弹性模量为207gpa。曲柄10和连杆20的弹性模量为20.7mpa，泊松比为0.29。

曲柄滑块机构的转轴与轴承间隙转动副接触应力分析的方法包括以下步骤：

s1:提取轴承接触边界有限元模型，为每个单元设置高斯积分点；

s2:搜索每一轴承接触边界高斯积分点在转轴接触边界上的投影点及其参数坐标；

s3:计算所有转轴投影点处法向间距及几何度量；

s4:为转轴接触边界每个节点计算等效投影点参数坐标、法向间距及其变分，几何度量，法向及切向牵引力；

s5:计算接触力虚功，提交系统动力学方程。

具体的，首先，对间隙运动副及其他结构进行有限元网格划分。本实施例采用二阶四边形单元对机构进行网格划分，如图2-4所示。其中图2为对整个曲柄滑块机构各部件进行网格划分后的结果。图3-4分别单独给出了轴承与转轴部分的网格。本实施例仅为对曲柄滑块机构的一种网格划分示例，所提出的应力计算方法并不局限于图示的网格划分。

为了计算间隙运动副内的接触力虚功，分别指定转轴和轴承为主、从接触体。分别提取主、从接触体接触边界上的一维单元，本实施例中各8个单元，并根据高斯积分准则在轴承接触边界每个单元上分布若干高斯积分点。本实施例中取插值函数阶次的2倍，即每个接触单元内部分布4个高斯积分点。如图5所示，设s为轴承接触边界上任一高斯积分点，为了判断点s与转轴是否发生接触，首先需要搜索转轴接触边界上点s的投影点m。点m对应的参数坐标满足如下非线性方程：

其中，f(·)表示函数，ξ为转轴接触边界的参数坐标，为投影点m的参数坐标，τ1为转轴接触边界上任一点处的协变基矢量，x¹和x²分别为转轴和轴承接触边界上任一点的位置矢量，

该非线性方程使用牛顿隐式迭代求解，其切线刚度矩阵为：

其中，几何度量m11：＝τ1·τ1。

确定投影点m参数坐标后，可根据式(3)计算积分点s与投影点m的法向间距：

其中gn为积分点与投影点间的法向间距，n为点m法向单位矢量，见图3。计算投影点m处几何度量用于接触边界有限元节点等效投影点处几何度量计算，度量空间几何的长度、面积和角度等。

循环对轴承接触边界上8个单元所有积分点的投影点进行搜索并计算其法向间距和几何度量后，针对轴承接触边界每一有限元节点分别计算等效投影点参数坐标，法向间距及其变分，几何度量，法向和切向牵引力。

具体步骤如下：

对于轴承边界上任一节点a，等效投影点参数坐标及其变分可按照式(4)-(5)计算：

其中，γ为接触边界，ra为与节点a对应之nurbs插值函数。

等效投影点处法向间距gna及其变分δgna可按照式(6)-(7)计算：

等效投影点处几何度量m11a可按照式(8)计算：

其中m11为高斯积分点对应投影点m处几何度量，即前面已分别算出，避免了重复计算，计算效率高。

等效投影点处法向牵引力tna可按式(9)计算：

tna＝∈ngna(9)

其中，∈n为法向惩罚参数。

为了计算节点等效投影点切向牵引力，需要区分切向接触状态(静摩擦或滑动摩擦)。定义式(10)所示的试验节点切向牵引力：

其中下标n指上一时间步，∈t为切向惩罚参数。根据试验节点切向牵引力的值，等效投影点切向牵引力可以按式(11)区分为静摩擦或滑动摩擦两种情况：

其中，μ为摩擦系数。

循环完成对每一节点的等效投影点参数变量计算后，计算接触力虚功δw：

将接触力虚功提交系统方程组计算，供系统动力学分析备用。

图6-7对比了采用本发明的方法对某一时刻间隙运动副接触区域应力波动的影响。图6为未使用本方法的结果，从图中可以看到，当惩罚因子∈n分别取10⁴和10⁶时，接触边界区域应力产生了明显的波动。图7为使用本方法后接触边界区域应力的情况，图中黑点为未采用本方法的接触边界区域应力情况。从图中可以看到，应力波动相对于图6中的结果(即图7中黑点表示的结果)明显降低，说明了本方法对提高接触应力计算精度的效果。

本发明针对间隙转动副有限元模型，对接触边界进行网格划分，并在轴承接触边界设置数值积分点。搜索每一数值积分点在转轴接触边界上的投影点及其参数坐标，并计算该投影点处法向间距及几何度量。循环完成所有投影点的搜索及相关变量计算后，针对轴承接触边界每个节点分别计算等效投影点参数坐标、法向间距、几何度量、法向及切向牵引力。循环完成对所有节点等效变量的计算后，计算接触力虚功积分，提交系统方程。

本发明提供的间隙转动副接触应力分析方法，能够获得与接触边界节点数目相同的接触力约束数目，有效避免过多积分点导致的过约束问题，有效降低接触边界的接触应力波动。