钢轨循环滚动接触过程残余应力-应变的变化规律研究

高明昕， 杨 建， 付丽华， 明乐乐， 宋 华

(1. 辽宁科技大大学 工程训练中心 ， 辽宁 鞍山 114051； 2. 辽宁科技大大学 机械工程与自动化学院，辽宁 鞍山 114051；3. 首钢股份公司迁安钢铁公司 硅钢事业部， 河北 迁安 064400)

钢轨生产过程难免产生如微小裂纹、非金属夹渣物等缺陷，铁路的高速、重载运输趋势使钢轨的服役条件更趋恶化，尤其在曲线轨道上，轮轨接触过程使钢轨产生复杂的接触应力-应变变化，是导致裂纹加深和扩散的主要原因之一[1-2]。另外，生产过程所产生的残余应力对在役钢轨的抗疲劳强度、抗断裂能力等性能有重要影响，是决定钢轨循环疲劳平均应力的关键因素之一[3]。生产残余应力和轮轨循环滚动接触所产生的残余应力共同决定了在役钢轨的残余应力-应变分布状态和损伤形式。因此，在考虑钢轨生产残余应力的基础上，研究曲线通过时钢轨循环滚动接触过程的残余应力-应变很有必要。

国内外学者对钢轨滚动接触相关问题进行了大量的研究，但鲜有考虑钢轨生产残余应力分析钢轨循环滚动接触应力-应变状态。Jiang等[4]利用接触单元法分析了钢轨二维稳态纯滚动接触法向接触压力。Kulkarni等[5]基于有限元法和半无限半空间椭圆接触算法，对钢轨三维循环纯滚动接触进行了弹塑性分析。温泽峰等[6-7]研究了多步非稳态载荷下的钢轨滚动接触应力和弹塑性变形。Wen[8]等基于循环塑性理论并利用移动的法向和切向载荷模拟了钢轨三维弹塑性滚动接触应力。Jiang等[9]利用圆形分布的Hertz接触模型和循环塑性本构模型研究了横向、纵向切向力对钢轨循环滚动接触应力-应变的影响。郭俊等[10]通过反复施加Hertz法向和切向接触力，建立了钢轨三维弹塑性循环滚动接触有限元模型，并分析了材料屈服强度对钢轨残余应力和应变的影响规律。Ringsberg等[11]利用临界平面法研究了生产残余应力对钢轨滚动接触疲劳寿命的影响，但未进行应力-应变的相关研究。常崇义等[12]基于速度变分的虚功率原理ALE有限元方法，采用相对滑移速度描述轮轨滚动接触斑内的滑动和黏着条件，分析了轮轨稳态滚动接触过程中不同横移量对法向力、切向力和横向力分布状态的影响。Jin等[13]结合Kalker非赫兹滚动接触理论、车辆/轨道耦合动力学模型和材料磨损模型，分析了曲线轨道对轮轨磨损和滚动接触应力状态的影响。目前，对于考虑生产残余应力的钢轨循环滚动接触应力-应变的相关研究仍不完善，如何有效考虑钢轨生产残余应力仍需进一步探讨。

本文提出了一种计算钢轨生产残余应力的弹塑性有限元方法，在考虑生产残余应力和循环塑性本构模型的基础下，对钢轨循环滚动接触应力-应变进行研究，并着重分析循环滚动次数、曲线半径、摩擦系数和纵向蠕滑率对钢轨残余应力-应变分布状态的影响规律。

1 钢轨残余应力

钢轨残余应力包括生产残余应力和轮轨循环滚动接触过程所产生的残余应力。钢轨生产过程主要包括热制、冷却和矫直等工序。在热制和冷却过程，钢轨因分别受到轧制力、固态相变力和收缩热应力的作用而在各部位产生不均匀的塑性变形，产生一定的残余应力。另外，在最后的矫直工序，钢轨在矫直力的作用下发生多次弹塑性弯曲变形而使其内部的残余应力重新分布，并产生矫直残余应力。因此，钢轨生产残余应力是热轧-冷却残余应力和矫直残余应力相叠加的结果。一般情况下，钢轨的生产残余应力(纵向残余应力)呈“C”形分布，即轨头和轨底部位呈拉应力状态，而轨腰部位呈压应力状态。钢轨生产残余应力的大小和分布状态将对钢轨的使用性能产生重要影响[14]。在服役过程，钢轨因受到车轮循环滚动接触载荷和摩擦力的作用而发生一定的弯曲变形和在接触表面发生塑性变形，形成一定的残余应力。钢轨循环滚动接触所产生的残余应力与其生产残余应力相互叠加，决定了钢轨最终的残余应力状态。该残余应力的分布状态对在役钢轨疲劳裂纹的萌生、断裂方式、断裂位置具有重要影响[14]。因此，研究在役钢轨循环滚动接触过程中的残余应力-应变变化规律具有重要意义。

2 数值模型

钢轨生产残余应力主要来自于热轧-冷却和矫直等工艺过程，在实际生产和应用中，计算成本巨大且很难通过实验准确测量。因此，提出了一种计算钢轨生产残余应力的有限元方法，在此基础上建立了考虑生产残余应力和循环塑性本构模型的钢轨循环滚动接触有限元模型。

2.1 钢轨生产残余应力有限元计算方法

钢轨经热轧-冷却工艺后，由于各部位冷却速率不同而产生向上翘曲并形成一定的残余应力，因热轧-冷却残余应力对钢轨矫直残余应力的影响较小[11]，故仅考虑矫直残余应力。经现场测量，钢轨冷却后的矫前初始曲线半径为63 m，取钢轨长度为5 m，基于ABAQUS建立的钢轨复合矫直有限元模型见图1，该模型由9个平矫辊和8个立矫辊矫组成，各辊压下量均按照标准规程设置，见表1。因矫直辊相对钢轨刚度大，矫直过程变形量较小，故采用离散刚性体进行描述。钢轨采用双线性随动硬化材料模型，该材料模型足以描述钢轨矫直过程中的弹塑性变化[14]，相关参数见表2。为获得较为准确的残余应力数据和便于后续的循环滚动接触分析，该模型将钢轨在纵向方向上分为5段，中间100 mm范围设置为第3段，第3段两侧各200 mm分别设置为第2和第4段，两端2 250 mm范围分别设置为第1和第5段，其中第2、3、4段钢轨的生产残余应力数据将用于循环滚动接触分析。采用R3D4(4节点、双线性、四边形、3D刚体)单元对矫辊进行网格划分，采用C3D8R(8节点、线性、六面体、减缩积分、3D应力)单元对钢轨进行网格划分，为降低计算成本并保证后续循环滚动接触应力-应变的分析精度，第1和第5段钢轨网格相对粗化，第2和第4段钢轨网格较精细，仅在第3段钢轨划分精细网格。各段钢轨模型之间均采用Tie绑定约束连接。

钢轨与各矫辊间的摩擦系数设为0.3，采用全局接触建立相互作用关系。在各辊几何中心建立参考点和局部坐标，约束各辊除沿轴线旋转方向的所有自由度，并沿矫直方向分别对平矫辊、立矫辊施加3、4.5 rad/s的恒定旋转速度，同时设置钢轨矫直方向初始速度为1 800 mm/s。为加快计算速度，设置质量缩放目标时间增量为5.98×10-5s(求解时间为5.98 s)。

表1 复合矫直机组各辊压下量

表2 钢轨材料性能参数

2.2 钢轨循环塑性本构模型

钢轨在车轮次屈服状态的非对称循环载荷作用下，将产生循环蠕变和塑性变形的积累，即棘轮效应。为准确描述钢轨循环滚动接触过程中的塑性变化，本文采用Chaboche循环塑性本构模型[16]， Mises屈服函数为

(1)

(2)

Chaboche循环塑性本构模型的背应力张量为

(3)

(4)

等效塑性应变为

(5)

引入非线性各向同性硬化准则，则屈服应力为

(6)

式中：σY,0为初始屈服应力；Q、b分别为屈服应力的最大变化量、随着应变的变化速率。

对于UIC 900A型钢轨，该循环塑性本构模型的相关参数取值[15]，见表3。

表3 Chaboche循环塑性本构模型的相关参数取值

2.3 钢轨循环滚动接触有限元模型

提取钢轨复合矫直有限元模型中的第2、3、4段，在引入生产残余应力和循环塑性本构模型的基础上建立的钢轨循环滚动接触有限元模型，见图2。其中，曲线轨道钢轨的曲线半径为2 000 m，超高h0=120 mm。

轮轨之间法向接触压力采用Hertz接触理论进行计算[17]，考虑到轮轨滚动接触过程中出现的部分滑移情况，将接触区域划分为滑移区和黏着区，见图3。

因此，为引入钢轨滚动接触滑移条件，采用Kalker三维弹性体非Hertz滚动接触理论对钢轨的纵向切向力进行计算，纵向切向力分布为[17]

(7)

式中：

(8)

(9)

其中，p0为最大法向接触压力；a为椭圆接触斑的长半轴；c、d分别为黏着区的长半轴和短半轴，可通过Kalker接触理论计算获得；ds=a-c。

钢轨滚动接触的纵向蠕滑率ξ为[18]

(10)

式中：r为车轮的名义半径；Qz为总的纵向切向力；P为总的法向接触力。

在滑移区有

(11)

则式(11)可表示为

(12)

假定曲线通过时钢轨所受的横向切向力与法向压力成线性比例关系，则横向切向力分布为

(13)

式中：T为总横向切向力。

法向压力、纵向切向力、横向切向力均可通过ABAQUS中的DLOAD用户子程序进行编制定义，并建立多组分析步，每组分析步代表一次滚动接触过程，循环滚动接触通过各分析步的时间变化控制接触斑位置的循环移动，其中接触斑的移动速度为200 km/h，接触斑的横移量为-1.642 mm(“-”表示向x轴负方向平移)，初始摩擦系数和蠕滑率分别为0.3和-0.001，进而实现对曲线曲线通过时钢轨循环滚动接触过程的模拟。该钢轨循环滚动接触有限元模型的其他相关参数初始值设置，见表4。

表4 钢轨循环滚动载荷的相关参数取值

3 结果分析

3.1 残余应力模拟结果

选取第3段钢轨中间横向对称面处的节点位置定义节点路径，其纵向残余应力曲线及相应的应力分布云见图4。由图4可见，残余应力呈“C”形分布，且整体应力状态呈拉-压-拉变化，与文献[15]中的残余应力结果基本一致，且轨底最大值小于250 MPa。

3.2 钢轨循环滚动接触残余应力-应变分析

以能较好反映钢轨循环滚动接触应力-应变状态的纵向残余应力σz、残余切应变γxy、γyz进行应力-应变分析(若无特殊说明，文中所述残余应力-应变均指纵向)。以20次循环滚动为例，轮轨循环滚动接触过程钢轨的残余应力σz、残余切应变γxy、γyz的变化规律见图5。钢轨踏面以下0.00 mm(最大残余压应力出现位置)、3.52 mm(最大残余切应变出现位置)、12.32 mm(最大残余拉应力出现位置)处节点的残余应力σz、残余切应变γxy、γyz的变化规律见图6。

由图5(a)可见，钢轨生产残余应力在第1次滚动接触后变化最大，尤其在轨头上部存在明显的拉-压变化，而其他部位仅在数值上变小而分布规律基本保持不变。随滚动次数增加，钢轨残应力主要在轨头上部发生变化，在其他部位基本保持不变。可以得到，循环滚动作用很快改变钢轨各部位的残余应力状态并使其重新分布。另外，由图5(a)、6(a)可见，随滚动次数增加，钢轨各部位残余应力均逐渐趋于稳定状态。由图5(b)、5(c)、6(b)、6(c)可见，钢轨残余切应变γxy和γyz仅在轨头上部发生变化，且随着钢轨滚动次数的增加仍近似成线性增加。

3.3 各参数对钢轨循环滚动接触残余应力-应变影响

由图5可见，钢轨在20次循环滚动后的最大残余应力和最大残余切应变均出现在轨头踏面处。因此，以钢轨踏面的一点为例进行分析。经20次循环滚动接触后，不同曲线半径对钢轨残余应力-应变的影响规律分别见图7。由图7可见，钢轨残余应力σz和残余等效塑性应变εp均随曲线半径的增加而不断减小，其中，当曲线半径由500 m增加至1 000 m时，变化速率相对较小。因此，相同工况下，较大的曲线半径将更有利于减小钢轨踏面的残余应力和残余等效塑性应变。

经20次循环滚动接触后，不同摩擦系数对钢轨残余应力-应变的影响规律见图8。由图8可见，钢轨残余应力σz和等效塑性应变εp均随摩擦系数的增加而逐渐增大。其中，纵向残余应力基本成线性变化，而残余等效塑性应变的变化速率则逐渐增加。因此，相同工况下，较小的摩擦系数将更有利于减小轨头表面的残余应力和残余等效塑性应变。

经20次循环滚动接触后，不同蠕滑率对钢轨残余应力-应变的影响规律分别见图9。由图9可见，随蠕滑率的增加，钢轨残余应力σz和残余等效塑性应变εp均先增大后减小，且最大值均出现在蠕滑率为-0.002时(蠕滑率为0.000时代表纯滚动，蠕滑率为-0.001和-0.002时代表部分滑移，蠕滑率为-0.004时代表全滑动)。同时可以看出，纯滚动较之于全滑动的情况，钢轨残余应力和残余等效塑性应变均相对较小；另外，处于纯滚动状态附近的滚动接触将更有利于减小钢轨踏面的残余应力和残余等效塑性应变。

4 结论

(1) 提出了一种引入钢轨生产残余应力的弹塑性有限元法，采用双线性随动硬化材料模型建立了钢轨生产残余应力有限元模型。所求得的生产残余应力呈“C”形分布，且整体应力状态呈拉-压-拉变化，与相关计算结果基本一致。

(2) 通过引入钢轨生产残余应力，并基于Chaboche循环塑性本构模型和 Kalker三维非赫兹接触理论，建立了用于分析钢轨循环滚动接触过程残余应力-应变的有限元模型。

(3) 随着循环滚动次数的增加，钢轨生产残余应力很快重新分布，并逐渐趋于稳定状态，且在轨头上部存在明显的拉-压变化，而残余切应变则仍近似呈线性增加。

(4) 随曲线半径的增加，钢轨残余应力和残余等效塑性应变均在数值上逐渐减小；随摩擦系数的增加，钢轨残余应力和等效塑性应变均在数值上逐渐增大，其中残余应力基本成线性变化，而残余等效塑性应变的变化速率逐渐增加；随纵向蠕滑率的增加，钢轨残余应力和残余等效塑性应变则均在数值上先增大后减小，且当蠕滑率为-0.002时数值最大。