轮轨滚动接触疲劳和磨损是铁路运营中不可避免的问题,尤其在重载铁路运营中,钢轨磨耗更为严重.而轮轨接触应力状态与钢轨磨耗程度密切相关.
关于轮轨接触理论,国内外学者进行了大量研究.20世纪90年代初,Kalker等[1]发展了三维弹性体非Hertz滚动接触理论及其数值方法CONTACT,该理论被称做轮轨滚动接触分析的完全理论;随后又开发了简化理论及其数值方法Fastsim.Hedrick[2]对VJ非线性蠕滑定律做了改造,形成了小自旋情形下三维非线性蠕滑率/力计算模型.Li[3]在非Hertz滚动接触理论基础上,在四分之一弹性空间内用有限元法对轮轨接触影响系数进行了修正,针对钢轨轨距角处发生的共形接触问题进行了研究.杨新文等[4-6]利用切片投影法寻找接触区域,并以分区计算模型修正Hertz理论求出了非椭圆斑上的法向应力.Teliskivi等[7]利用有限元软件ANSYS分析了轮轨接触应力和接触区域,并与Hertz接触理论和完全接触理论进行了比较.
由于当车轮轮缘与钢轨在轨距角处发生接触时,非Hertz理论中基于弹性半空间条件下的影响系数已不适用,因此,本文利用切片投影法找出轮轨接触点,然后利用有限元法在全空间内对非Hertz理论中的影响系数进行插值计算,在保证计算效率的前提下,更加精确地求解轮轨法向应力,为轮轨滚动接触疲劳和磨损问题的研究提供理论基础,并与P_M模型求出的轮轨法向应力进行了对比.
1 轮轨法向接触应力计算方法 1.1 非Hertz接触理论1982-1992年,Kalker基于弹性力学能量原理发展了三维非Hertz滚动接触理论CONTACT,在弹性半空间的假设下,利用BossinisqeCerruti公式及牛顿迭代法,可以求出三维轮轨滚动接触的精确解.图 1为矩形区域中离散的接触斑,Δx1、Δx2为单元网格的纵向边长和横向边长.区域Ω为接触斑范围,也即离散的区域中发生接触的部分.pI为单个网格内的应力.
在该模型中,基于弹性力学,余能原理的离散形式为[8]
$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\left( {{p_{Jj}}} \right)} C = \frac{1}{2}{p_{Ii}}{A_{IiJj}}{p_{Jj}} + \left[ {\left( {{h_J} - q} \right){p_{J3}} + } \right.}\\ {\left. {\left( {{W_{J\tau }} - {{u'}_{J\tau }}} \right){p_{J\tau }}} \right]{A_0}} \end{array}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;{p_{J3}} \ge 0\;\;\;\;\left| {{p_{J\tau }}} \right| \le {b_J}\;\;\;\;\forall x \in {A_{\rm{c}}}\\ {A_0}\sum\limits_J^M {{p_{J3}} = {F_{\rm{N}}}} \end{array} \right. $ | (1) |
式中:i,j=1,2,3, 为接触斑切平面局部坐标系的坐标轴x1,x2,x3的3个方向;τ=1,2, 为接触斑局部坐标系的坐标轴x1,x2的2个方向;I, J为矩形单元的编号;pJi为作用在单元J上的应力沿i轴的作用力分量,在同一个单元上,把它看成常数;q为轮轨在x3方向的弹性压缩量;AIiJj为影响系数,表示单元I上作用沿i轴方向的单位力引起单元J中心处沿j轴方向发生的位移;FN为法向轮轨力;hJ为轮轨接触面之间的法向间隙在单元J中心处的分量,通过轮轨接触几何计算确定;WJτ为轮轨间的刚性滑动量;uJτ为轮轨间的弹性位移量;Ac为全部矩形单元;A0为矩形单元面积,A0=Δx×Δx=0.5 mm×0.5 mm=0.25 mm2;bJ为J单元中心处库仑极限摩擦力.
1.2 法向问题的求解如果不考虑轮轨材料本构模型的差异,求解法向问题时,可把式(1) 化简成下式进行计算:
$ \left\{ \begin{array}{l} \mathop {\min }\limits_{{p_{J3}}} C = \frac{1}{2}{p_{I3}}{A_{I3J3}}{p_{J3}} + \left( {{h_J} - q} \right){p_{J3}}{A_0}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;{p_{J3}} \ge 0\;\;\;\;\forall X \in {A_{\rm{c}}}\\ {A_0}\sum\limits_J^M {{p_{J3}} = {F_{\rm{N}}}} \end{array} \right. $ | (2) |
引入Lagrange乘子λJ3,法向问题的KuhnTucker条件为
$ \begin{array}{*{20}{c}} {{h_1} - q + {\lambda _{13}} + \sum\limits_I {{A_{I313}}{p_{I3}}} = 0}\\ {{h_2} - q + {\lambda _{23}} + \sum\limits_I {{A_{I323}}{p_{I3}}} = 0}\\ \vdots \\ {{h_M} - q + {\lambda _{M3}} + \sum\limits_I {{A_{I3M3}}{p_{I3}}} = 0}\\ {{\lambda _{M3}}{p_{13}} = 0}\\ {{\lambda _{M3}}{p_{23}} = 0}\\ \vdots \\ {{\lambda _{M3}}{p_{M3}} = 0}\\ {{F_n} - {A_0}\sum\limits_J {{p_{J3}}} = 0}\\ {{p_{J3}} \ge 0,J = 1,2, \cdots ,M}\\ {{\lambda _{J3}} \ge 0,J = 1,2, \cdots ,M} \end{array} $ | (3) |
方程(3) 总共有2M+1个独立变量,其中pJ3为M个,λJ3为M个,q为1个.pJ3≥0,λJ3≥0为约束条件.针对非线方程组(3),利用NewtonRaphson法即可求得法向应力pJ3.
1.3 影响系数的修正影响系数AI3J3是指在接触区某一点施加一个单位荷载,从而引起的另一点的位移.它对式(3) 求解的精度和迭代过程中的数值稳定性都有重要影响.在CONTACT中,影响系数可通过BossinisqeCerruti公式求出.但是当车辆通过曲线段时,钢轨轨距角和车轮轮缘经常发生接触,轮轨接触半径和接触区尺寸接近,半空间假设不再成立,而有限单元法则不受半空间假设的限制,因此可以用有限单元法来计算影响系数.
根据我国重载铁路标准LM型踏面和CHN75钢轨型面,建立轮轨结构三维有限元模型,如图 2和图 3所示.车轮和钢轨均用实体单元SOLID45进行离散.钢轨与车轮泊松比分别为0.29和0.30,弹性压缩模量分别为214 000 MPa和210 000 MPa.根据影响系数的对称性,为了节省计算时间,车轮有限元模型建立一半,并在可能发生的接触范围进行网格加密,施加对称约束.
由于网格尺寸微小,有限元模型单元数太大,因此计算速度会很慢.如图 4、图 5所示,为了保证计算效率,在有限元模型中只在主轮廓线上的节点施加单位力,并提取出整个可能出现的接触范围内的法向位移.在求出法向间隙后,通过三次多项式插值的方法求出法向间隙内的影响系数AI3J3.
具体步骤如下:
(1) 有限元模型中,车轮和钢轨分别在各自坐标系下建模,所以首先要统一坐标系.如图 6所示,在轮对坐标系o′z′y′中,将车轮坐标x′i、y′i、z′i按公式(4) 转过一个侧滚角φw和一个摇头角θw,并偏移一个横移量yw后,转到钢轨坐标系OZY中,则分别得到车轮和钢轨的所有可能发生的接触范围Ω.
(2) 如图 7所示,在确定轮轨的法向间隙范围Ω′后,如果要求在点A′作用法向单位力时,Ω′范围内其他点的法向位移,则要找到点A′距离范围Ω内的最近点A,并将点A′移动到位于主轮廓线上的,与A点同在一个ZY平面上的B点,同时范围Ω′内的其他点移动相同的距离.
$ \left\{ \begin{array}{l} {x_i} = {{x'}_i}\cos {\varphi _{\rm{w}}} - {{y'}_i}\cos {\theta _{\rm{w}}}\sin {\varphi _{\rm{w}}} + {{z'}_i}\sin {\theta _{\rm{w}}}\sin {\varphi _{\rm{w}}}\\ {y_i} = {{x'}_i}\sin {\varphi _{\rm{w}}} + {{y'}_i}\cos {\theta _{\rm{w}}}\cos {\varphi _{\rm{w}}} - \\ \;\;\;\;\;\;{{z'}_i}\sin {\theta _{\rm{w}}}\cos {\varphi _{\rm{w}}} + {y_{\rm{w}}}\\ {z_i} = {{y'}_i}\sin {\theta _{\rm{w}}} + {{z'}_i}\cos {\theta _{\rm{w}}} \end{array} \right. $ | (4) |
(3) 得到移动后的法向间隙范围Ω′内的坐标x′、y′后,则可以通过MATLAB中的griddata函数对范围Ω内的坐标x、y及已经求得的B点的影响系数进行插值计算, 进而可以求得A′点的范围Ω′内的影响系数.有限元法得到钢轨的接触影响系数为AI3J3r,车轮的接触影响系数为AI3J3w.则最终需要代入方程(3) 获得综合的影响系数为
$ {A_{I3J3}} = \left( {A_{I3J3}^{\rm{r}} + A_{I3J3}^{\rm{w}}} \right)/2 $ | (5) |
有限元法计算得到的影响系数及BossinisqeCerruti公式求出轮轨接触的影响系数见图 8.
图 8为在轨顶处45 mm×43 mm,边长0.5 mm的网格范围中心施加单位力得到的影响系数,可以看出2种方法得到的影响系数形式基本一致.BossinisqeCerruti公式求出的影响系数最大值为9.731×10-9 m,有限元法得到的钢轨影响系数最大值为1.61×10-8 m,车轮影响系数最大值为1.38×10-8 m.
图 9为2种方法在钢轨不同位置处施加单位力得到的主轮廓线上的影响系数对比.在轨顶处,Bossinisqe-Cerruti公式求出的影响系数最大值为9.731×10-9 m,有限元法得到的钢轨影响系数最大值为1.54×10-8 m.在轨距角处,Bossinisqe-Cerruti公式求出的影响系数最大值为9.731×10-9 m,有限元法得到的钢轨影响系数最大值为1.61×10-8 m.由图 8可知,BossinisqeCerruti公式求出的影响系数轨顶与轨距角都一样,而利用有限元法得到的钢轨影响系数轨顶与轨距角处都不一样,在轨距角处的影响系数大于轨顶处,这是因为轨距角处型面变化明显,已不再符合半空间的条件.
P_M(partition model)接触模型[6]是一种基于Hertz接触理论,在车轮和钢轨接触轮廓型面不同弧段进行分区后,先分开分析,然后通过连续性条件将各块接触斑连接起来,进而求出法向应力的接触模型.
对于两个椭圆弹性物体接触,Hertz接触理论可得到应力呈椭圆形分布的接触斑.在已知法向接触力N的条件下,接触斑的弹性压缩量δ、长半轴b、短半轴a及最大应力pO可由下式求出:
$ a = m{\left[ {\frac{{3{G^ * }N}}{{4\left( {A + B} \right)}}} \right]^{\frac{1}{3}}} $ | (6) |
$ b = n{\left[ {\frac{{3{G^ * }N}}{{4\left( {A + B} \right)}}} \right]^{\frac{1}{3}}} $ | (7) |
$ \delta = \frac{{K\left( e \right)}}{{{\rm{\pi }}m}}{\left[ {\frac{9}{2}\left( {A + B} \right){G^{ * 2}}{N^2}} \right]^{\frac{1}{3}}} $ | (8) |
$ {p_O} = \frac{1}{{{\rm{\pi }}mn}}{\left[ {\frac{{6{{\left( {A + B} \right)}^2}N}}{{{G^{ * 2}}}}} \right]^{1/3}} $ | (9) |
式(6)~(9) 中:A和B为接触参数;K(e)表示第一类椭圆积分,其他接触参数如下:
$ {G^ * } = \frac{{1 - {\nu _1}}}{{2{G_1}}} + \frac{{1 - {\nu _2}}}{{2{G_2}}} = \frac{{1 - \nu _1^2}}{{{E_1}}} + \frac{{1 - \nu _2^2}}{{{E_2}}} $ | (10) |
$ m = {\left[ {\frac{{2E\left( e \right)}}{{{\rm{\pi }}\left( {1 - {e^2}} \right)}}} \right]^{1/3}} $ | (11) |
$ n = m\sqrt {1 - {e^2}} = {\left[ {\frac{{2E\left( e \right)}}{{\rm{\pi }}}} \right]^{1/3}}{\left( {1 - {e^2}} \right)^{1/6}} $ | (12) |
式(10)~(12) 中:ν1和ν2分别为车轮和钢轨的泊松比;G1和G2分别为车轮和钢轨的剪切模量;E1和E2分别为车轮和钢轨的弹性模量;E(e)表示第二类椭圆积分;e为
在轮轨发生接触的区域,钢轨及车轮的曲率半径会有变化,因此若只用轮轨刚性接触点处的曲率半径去计算接触参数A、B,进而得到的接触斑和接触应力必然是不准确的.P_M模型首先根据车轮和钢轨的曲率变化对接触区域进行分区,再根据变化后的接触形状对接触应力进行修正.
如图 10所示,根据曲率半径的不同将接触区域分为3部分,Zone1、Zone2、Zone3,保持中心区域Zone1不变,修正其他区域,修正的原则是确保接触边界上的法向间隙相等.接触区域修正后,法向接触应力的合力也跟着发生改变,不再等于轮轨法向接触应力N.刚性接触点O处的接触应力pO根据Hertz理论可以得到,中心接触区域的边界点J和K处的接触应力也可以得到,分别为pJ和pK,进而可以得到3个接触区域各自的应力分布.通过求和即可得到新的法向力N′, 若N′与N不满足误差要求,则根据式(12) 对最大应力pO进行修正,从而得到新的接触斑.重复以上步骤,直到N′与N满足误差要求为止.
$ {{p'}_O} = {p_O}{\left( {\frac{N}{{N'}}} \right)^{1/3}} $ | (13) |
采用重载铁路标准CHN75轨型面和LM型车轮踏面,如图 11和图 12所示.
30 t轴重重载铁路轨道在摇头角等于0°的情况下,用修正的非Hertz法与P_M分区法得到在轮对横移量为0、4和8 mm时的法向应力p如图 13~15所示.
图 16为轮轨法向接触应力最大值随轮对横移量变化曲线.由图 16可知,修正的理论和P_M模型计算得到的最大接触应力的变化趋势基本一致.横移量为负时,两种计算方法差别较小,最大相差8%.横移量为正时,两种计算方法差别较大,最大相差22%.随横移量变化,修正的理论始终比P_M模型计算得到的最大接触应力要大,最大可达3 659.6 MPa.当轮对横移量为-8~0 mm时,由于左轮左轨踏面不断贴合,最大接触应力有先增大后减小的趋势,但变化幅度不大,接触应力在1 500~2 500 MPa之间变化;然而当横移量从0 mm变化为8 mm时,最大接触应力有先减小而后急剧增大的变化趋势,这是由于当轮对横移量为0~4 mm时,轮轨型面较为匹配,接触斑面积较大;当横移量持续增大时,由于车轮与钢轨轨距角接触,接触面积急剧降低,造成接触应力急剧增大.
图 17为轮轨接触斑面积随轮对横移量变化曲线.由图 17可知,修正的理论和P_M模型计算的接触斑面积的变化趋势基本一致.随横移量变化,修正的理论始终比P_M模型计算得到的接触斑面积要大.当轮对横移量从0 mm变化为-8 mm时,由于左轮左轨踏面不断贴合,接触斑面积有增大的趋势,最大接触斑面积可达182.75 mm2;然而右轮右轨接触斑面积有先增大后急剧减小的变化趋势,这是由于当轮对横移量为0~4 mm时,轮轨型面较为匹配,最大接触斑面积可达173.75 mm2;当横移量持续增大时,由于车轮与钢轨轨距角接触,接触面积急剧降低.
(1) 用有限元法计算出的影响系数大于BossinisqeCerruti公式求出的影响系数.利用有限元法得到的钢轨影响系数在轨顶与轨距角处不同,在轨距角处的影响系数大于轨顶处.这是由于钢轨轨距角型面变化已不符合半空间的假设,因此用有限元法计算的影响系数更适用于车轮与钢轨发生两点接触时的应力计算.
(2) 修正的理论和P_M模型计算出的法向应力,随横移量变化而变化的趋势一致,且修正的理论计算的法向接触应力始终比P_M模型计算的法向应力略大.
(3) 修正的理论和P_M模型计算出的接触斑面积,随横移量变化而变化的趋势一致,且修正的理论计算的接触斑面积始终要比P_M模型计算的更大.当横移量为0~4 mm时,轮轨型面较为匹配,最大接触斑面积可达173.75 mm2,当横移量持续增大时,由于车轮与钢轨轨距角接触,接触面积急剧降低,同时法向接触应力急剧增大.
[1] |
KALKER J J, JOHNSON K L. Three-dimensional elastic bodies in rolling contact[M]. Delft: Kluwer Academic Publishers, 1990
|
[2] |
HEDRICK J K. A Comparison of alternative creep force models for rail vehicle dynamic analysis[J]. Vehicle System Dynamics, 1983, 12(1/2/3): 79 |
[3] |
LI Z. Wheel-rail rolling contact and its application to wear simulation[J]. Fatigue & Fracture of Engineering Materials & Structures, 2002, 26(10): 999 |
[4] |
杨新文, 顾少杰, 练松良. 30 t轴重重载列车轮轨法向接触应力分析[J]. 铁道学报, 2015, 37(6): 19 YANG Xinwen, GU Shaojie, LIAN Songliang. Normal stress of wheel/rail contact under 30 t axle's heavy haul wagon[J]. Journal of Chinese Rail Society, 2015, 37(6): 19 |
[5] |
YANG X W, GU S J, ZHOU S H, et al. A method for improved accuracy in three dimensions for determining wheel/rail contact points[J]. Vehicle System Dynamics, 2015, 53(11): 1620 DOI:10.1080/00423114.2015.1066508 |
[6] |
GU S J, YANG X W, ZHOU S H, et al. An innovative contact partition model for wheel/rail normal contact[J]. Wear, 2016, 366(1): 38 |
[7] |
TELLISKIVI T, OLOFSSON U. Contact mechanics analysis of measured wheel-rail profiles using the finite element method[J]. Proceedings of the Institution of Mechanical Engineers Part F: Journal of Rail & Rapid Transit, 2001, 215(2): 65 |
[8] |
金学松, 刘启跃. 轮轨摩擦学[M]. 北京: 中国铁道出版社, 2004 JIN Xuesong, LIU Qiyue. Wheel/rail tribology[M]. Beijing: China Railway Publishing House, 2004 |