赵春雷，王俊龙，赵尚超，张 强

(中车齐齐哈尔车辆有限公司，齐齐哈尔161000，中国)

高速重载是我国货运机车的发展方向，轴重提高意味着需要降低摇枕自重[1-2].摇枕薄壁结构应力分布复杂且变化急剧，有限元仿真的难度相对较大.目前，在货车摇枕结构的有限元分析中，往往没有根据应力分布的变化特点进行相应的网格疏密划分[3-4]，或者忽略了实际中存在的接触问题的影响[5-6]等.同时，也缺少一定的试验结果来辅助有限元仿真时的建模和修改[7-8].这些问题都会影响到有限元仿真结果的准确性.因此，对于货车摇枕薄壁复杂结构的疲劳寿命分析，需要建立与实际情况一致的仿真计算模型从而获得更为准确的应力结果.

基于隐式边界法[9-10](IBM，Implicit Boundary Method)的全新仿真分析软件MIDAS MeshFree可解决对复杂结构进行有限元分析时网格难以处理的问题，文中运用MIDAS MeshFree软件对转K6摇枕进行结构应力有限元仿真分析.与应力测试结果相互比较，建立了基于接触非线性问题的仿真计算模型，将计算结果和测试结果进行比较分析，确定了合理的摩擦系数.提高了摇枕结构应力的仿真准确性.

1 隐式边界法原理

1.1 D函数的构造

MIDAS MeshFree用的是隐式边界法.假设图1所示为一个平面，首先在分析对象上面覆盖网格，网格单元非常规则，称为结构化网格[11].

这些单元可以分为三大类：一类是完全在分析对象内部的叫内部单元；另一类称为边界单元，边界单元又可以分为含有边界条件的和不含边界条件的；还有一类是完全在分析对象外部的叫外部单元，外部单元不参加计算.构造位移函数为

u=us+ua=DNqg+Nqa，

(1)

σ=σs+σa，

(2)

式中：ua为边界值函数；D为Dirichlet函数；N为形函数；qg为网格上分片插值/逼近的网格变量；qa为节点强制位移.这样构造函数的目的是为了满足边界条件.

D函数的表达式为

(3)

式中：D函数中的φ函数是用来描述边界的，是隐式方程.当φ≤0时，D函数的值为0，表示在边界上或分析对象外面；当0≤φ≤δ时，δ为较小的常数，一般δ≈10-5或者更小，表示在靠近边界的窄带范围内；φ≥δ时，表示在分析对象内部.

1.2 刚度矩阵的计算

虚功方程为

(4)

式中：{σ}表示应力张量；{t}表示面向量；{b}表示体力向量.

代入式(2)得

(5)

由上式可得刚度矩阵为

(6)

如果是线弹性问题，应变矩阵为

(7)

对于内部单元，因为D函数的值为1，所以应变矩阵变为

(8)

式中:Ni为形函数;m表示单元节点数.

因此，对于内部单元，采用通用的有限元方法处理.

对于不含有边界条件的边界单元，D函数的值也为1.这类单元被分为有材料和无材料两部分，通过调整高斯积分点的权重系数和高斯积分点的位置来实现单元刚度矩阵的计算.刚度矩阵的高斯积分公式为

(9)

式中:nx、ny为单元中x方向、y方向上高斯积分点的个数；(ξi,ηi)为高斯积分点的位置.

对于含有边界条件的边界单元，位移表达式中含有D函数，所以应变矩阵中含有D函数的导数.因为δ的值很小，所以D函数及其导数对刚度矩阵的影响不可忽略.将应变矩阵变成两个部分，[B1]和[B2]的表达式如下：

(10)

(11)

[B]=[B1]+[B2].

(12)

相应地，刚度矩阵变为

[Ke]=[K1]+[K2]+[K2]T+[K3],

(13)

(1)[K1]的计算.[K1]只含[B1]，只有D函数本身以及形函数的导数，而其只在单元的实体部分不为零，所以,可以采用与不含边界条件的单元矩阵相同的计算方法.

(2)[K2]的计算.[K2]同时含有D函数和形函数的导数，而D函数只在窄带区域内非零，宽度为δ,[K2]的计算可以化为线积分进行求解.

(3)[K3]含有D函数的导数，计算方法同[K2].

计算完[K1]、[K2]和[K3]后，相加就是含有边界条件的边界单元的刚度矩阵.实际上，对于MIDAS MeshFree，求解的问题都是三维的，将上述推导方法向三维空间扩展即可.

1.3 边界值函数ua的构造

边界值函数的作用是使构造的位移函数满足强制位移，用形函数来构造位移函数，能够保证解结构的完备性要求.

(14)

在结构的边界Γa上施加强制位移ua=ξ，需要在边界Γa经过的所有单元节点上施加强制位移，其他不相关的节点施加强制位移0.

2 MIDAS MeshFree软件特点

(1)该软件专门用于实体仿真分析，用MIDAS MeshFree进行分析时只需按原结构将三维实体模型建立出来倒入软件进行分析.

(2)无需模型简化和网格划分.软件采用的隐式边界法，单纯地利用结构化网格进行分析，不管模型多么复杂，都可以进行分析；而传统的有限元分析，生成网格时需要考虑模型的各种形状，为了生成高品质的网格，需要进行模型简化清理工作.网格区别如图2所示.

图2 不同网格形式

(3)软件中有丰富的CAD数据接口.常见的CAD软件,如UG、ProE、Solidworks、CATIA等生产的自身格式都可以导入到软件中，也支持内核格式Parasolid、ACIS和商业格式STEP和IGES的导入.

(4)简单直观的三步骤分析流程.只需导入 3D CAD模型，添加荷载和边界条件就可以进行分析和查看结果，与以往的有限元分析相比，缩短了大量的时间且结果很相近.

(5)丰富的分析功能.目前的功能有线性和非线性静力分析、模态分析、稳态/瞬态热传递分析、热应力分析、疲劳分析、线性动力分析和拓扑优化分析等.

3 基于试验建立分析模型

3.1 试验测试

转K6型摇枕是我国重载货车用转向架摇枕，为对称结构，长2 430 mm，心盘座部位宽度470 mm.侧壁厚20 mm，双筋板厚18 mm，上平面厚23 mm，底部厚 30 mm.车体载荷施加于心盘上，二系弹簧组支撑于摇枕两端的底面[12-13].该摇枕依据TB/T1959-2006[14]进行静载试验，两端底面支撑于两个台架上，两个台架下部各由一根滚柱支撑，在心盘上施加1 040.96 kN的垂向载荷.加载设备为5 000 kN四柱压力试验机；数据采集设备为UCAM-60静态数据采集仪1台；试验测试件经过精确地划线确定了应力测点，贴片时避开铸造缺陷和凹坑.粘贴应变片区域分为：A区域应力测点、B区域应力测点，试验测试粘贴应变片位置和测点区域在摇枕上的坐标(如图3所示).各测点传感器与测量仪器之间的导线连接以及仪器调试后，应首先进行预加载，所加载可按两级或三级缓慢加载，在各级载荷下预加载时，被测结构和测试系统均处于正常状态，方可进行正式试验.正式试验时，每种工况加(卸)载次数不少于三次，取其中两次数据比较稳定的进行平均值处理，测试摇枕数量为两件.

图3 A区域和B区域测点位置

3.2 建立分析计算模型

转K6摇枕材料为B+级钢，材料的弹性模量为172 000 MPa,泊松比0.3.在进行结构分析时，首先分配和定义材料；其次添加约束条件，为了保证试验数据和仿真计算的可比对性，合理简化试验工装进行仿真计算模型建立；考虑到摇枕两端底面与台架顶面之间的接触非线性问题，试验工装与转K6摇枕接触处建立非线性接触连接，上心盘与摇枕接触处建立绑定连接；在上心盘与试验加载工装接触面添加试验载荷；进行网格设置，摇枕和上心盘采用用户指定(数量)进行网格划分，试验工装采用自动网格(密度)进行网格划分，设置分析控制参数后进行计算分析.

在转K6型转向架摇枕试验时，为了使摇枕和支撑工装接触充分，将摇枕与工装接触面放置铅板，在仿真计算时,摇枕和支撑工装接触面处建立接触连接设置参数，其中,摩擦系数的不同直接影响摇枕仿真计算精度，因此,设置不同的摩擦系数进行计算分析研究，依据试验台架建立的仿真计算模型如图4所示.

图4 仿真计算模型

4 结果对比分析

分析摩擦系数对应力计算结果的影响，分别设置摩擦系数μ=0.1、0.15、0.2、0.3进行计算分析，提取计算结果测点位置处的应力值列于表中.去除试验过程中带来的数据波动影响，将试验测试对称测点数据进行求平均值处理，将处理后的试验数据列于表格中.计算不同摩擦系数下试验测试结果与仿真计算结果差值，绘制不同摩擦系数下对应力差值的影响曲线.应力计算结果、试验测试结果及差值见表1和表2，不同摩擦系数下应力差值影响曲线如图5所示.

表1 B区域实测应力、计算应力及差值 MPa

表2 A区域实测应力、计算应力及差值 MPa

从图5中曲线变化情况可以看出，计算中摩擦系数对测点A区域的应力值影响比测点B区域的应力值影响大.原因一是A区域处于试验件与试验工装接触区域较近位置，摩擦系数的边界条件对其影响较敏感.原因二是A区域处于几何外形突变区域，这种区域是应力相对集中区域，应力梯度大.摩擦系数对其计算应力值的影响也较大.从图中曲线可以看出，随着摩擦系数的增大，A区域和B区域应力值都将减小，A区域相对B区域减小快.当μ=0.15时，A区域和B区域所有测点应力差值最小，最大应力差值为4.9 MPa，且应力差值分布均匀.当μ<0.15或是μ>0.15时，应力差值相对较大且随着摩擦系数的变化，应力差值分布离散.因此,计算摩擦系数μ=0.15时，应力差值水平已经满足,工程上的仿真计算误差要求，可以将μ=0.15作为相对最为合理的仿真计算摩擦系数输入值.

图5 摩擦系数对差值的影响

5 结 论

通过文中对比分析结果得出，对于转K6型摇枕的计算分析，工装与试验件接触处摩擦系数对应力计算结果的影响较为显著，其计算中最为合理的接触摩擦系数值应为0.15.文中获得的摩擦系数值可以作为基于摇枕静载荷试验台架计算仿真分析时的摩擦系数输入值，其建立的基于台架的仿真计算模型，可获得合理的转向架摇枕的疲劳寿命分析结果.