考虑材料温度相关性的二维轮轨弹塑性滑动接触温升分析1)

考虑材料温度相关性的二维轮轨弹塑性滑动接触温升分析 ¹⁾

伏培林^*,丁立^*,赵吉中^*,张旭^*,阚前华^,^*^,²⁾,王平^†

^*西南交通大学力学与工程学院, 成都 610031

^†西南交通大学土木工程学院, 成都 610031

FRICTIONAL TEMPERATURE ANALYSIS OF TWO-DIMENSIONAL ELASTO-PLASTIC WHEEL-RAIL SLIDING CONTACT WITH TEMPERATURE-DEPENDENT MATERIAL PROPERTIES ¹⁾

Fu Peilin^*,Ding Li^*,Zhao Jizhong^*,Zhang Xu^*,Kan Qianhua^,^*^,²⁾,Wang Ping^†

^*School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu 610031, China

^†School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China

通讯作者:2)阚前华, 教授, 主要研究方向: 材料循环本构关系及其疲劳失效. E-mail:qianhuakan@foxmail.com

收稿日期:2019-04-16接受日期:2019-04-16网络出版日期:2020-09-18

基金资助:

1)高铁联合基金 . U1734207
四川省杰出青年基金 . 2017JQ0019

Received:2019-04-16Accepted:2019-04-16Online:2020-09-18

作者简介 About authors

摘要

轮轨滑动接触温升的准确预测对于轮轨的磨耗和疲劳研究均具有重要意义. 目前的轮轨温升解析或半解析模型通常考虑Hertz弹性接触压力分布和单一材料属性的温度相关性, 与实际的轮轨传热状态尚有一定偏差, 因此在轮轨滑动温升计算模型中考虑接触压力的塑性修正和多种材料属性的温度相关性, 有望提高温升预测结果的准确性. 基于弹塑性接触理论, 同时考虑热导率、比热容和摩擦系数的温度相关性, 通过基尔霍夫变换方法以热导率温度相关性函数的积分作为待求量, 将复杂的非线性Fourier导热方程转化成含单个变系数的简单偏微分方程形式, 从而构建了一种不限制材料温度相关性函数形式的统一隐式差分求解格式, 分别讨论了对流换热系数、法向载荷、蠕滑率以及行车速度对钢轨表面滑动温升的影响. 结果表明, 当列车高速行驶时, 对流换热系数对轮轨滑动温升的影响甚微; 蠕滑率和行车速度的增大, 均会引起摩擦功率的增大, 进而导致钢轨表面温度的升高; 钢轨表面滑动温升的峰值随法向载荷的增大而近似线性上升. 此外, 在轮轨滑动温升计算模型中考虑材料属性的温度相关性可有效避免对滑动温升的过分高估, 且摩擦系数的温度相关性对温升的影响要显著强于热导率和比热容.

关键词： 弹塑性接触 ; 隐式差分法 ; 轮轨滑动接触 ; 摩擦温升 ; 温度相关性

Abstract

The accurate prediction of frictional temperature during wheel-rail sliding contact is of great significance to the studies for wear and fatigue performance of wheel-rail system. Current analytical or semi-analytical models of wheel-rail frictional temperature usually employ the elliptical distribution of contact pressure in the Hertz elastic contact theory and a single temperature-dependent material property, which differs from the actual heat-transfer state in many wheel-rail sliding contact conditions. Therefore, introducing the plastic correction of contact pressure and the temperature dependence of various properties in the calculation model of wheel-rail sliding temperature rise simultaneously could be great helpful to improve the accuracy of prediction result. Based on the elasto-plastic contact theory, considering the temperature dependence of thermal conductivity, specific heat capacity and friction coefficient simultaneously, the integration of thermal conductivity with respect to temperature is set as the quantity that needs to be solved by applying the Kirchhoff transformation method, and the nonlinear Fourier heat conduction equation is transformed into a corresponding simple partial differential equation with single variable coefficient, a unified implicit difference scheme with arbitrary form of temperature dependence is derived, and the influences of convection coefficient, vertical load, creepage and train speed on the temperature rise over rail surface are discussed, respectively. Results show that the convection coefficient has little effect on the temperature rise in the high-speed condition; the increase of creepage and train speed can enhance the friction power, and thus causes the increase of frictional temperature; and the maximum temperature rise also increases approximately linearly with the increasing vertical load. In addition, considering the temperature dependence of various thermophysical properties in the calculation model of temperature rise induced by the wheel-rail sliding contact can effectively avoid the overestimation of temperature rise, and the temperature-dependent friction coefficient has a more remarkable effect on the prediction result than thermal conductivity and specific heat capacity.

Keywords： elasto-plastic contact ; implicit difference method ; wheel-rail sliding contact ; frictional temperature rise ; temperature dependence

PDF (569KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex 收藏本文

本文引用格式

伏培林, 丁立, 赵吉中, 张旭, 阚前华, 王平. 考虑材料温度相关性的二维轮轨弹塑性滑动接触温升分析 ¹⁾ .力学学报[J], 2020, 52(5): 1245-1254 DOI:10.6052/0459-1879-20-122

Fu Peilin, Ding Li, Zhao Jizhong, Zhang Xu, Kan Qianhua, Wang Ping. FRICTIONAL TEMPERATURE ANALYSIS OF TWO-DIMENSIONAL ELASTO-PLASTIC WHEEL-RAIL SLIDING CONTACT WITH TEMPERATURE-DEPENDENT MATERIAL PROPERTIES ¹⁾ . Chinese Journal of Theoretical and Applied Mechanics [J], 2020, 52(5): 1245-1254 DOI:10.6052/0459-1879-20-122

引言

近年来, 随着我国轨道交通行业的迅速发展, 列车的行驶速度和运量均不断提升, 这对列车的牵引和制动性能提出了越来越高的要求. 列车的牵引和制动性能通过轮轨之间的接触关系来实现, 一旦列车的牵引力或制动力超过轮轨接触关系所能提供的粘着力时, 轮轨间便会发生相对滑动^[1]. 在滑动过程中, 轮轨接触斑内的温度会因摩擦生热而迅速上升, 并随着车轮的进一步滑动脱离接触而快速衰减, 导致热应力出现交变变化. 这种交变的温升诱发热应力, 叠加在轴重引起的机械应力上, 将大幅度降低钢轨和车轮的服役寿命^[2-3]. 同时, 当滑动温升足够高时, 轮轨表层材料发生马氏体相变, 并伴随着后续的快速冷却过程而形成脆硬的马氏体组织. 在轮轨接触载荷的作用下, 这些脆硬组织容易发生脱落, 造成轮轨表面的扁疤和材料剥离, 进而影响列车行驶的安全性^[4-5]. 另外, 过高的轮轨表面温度也会导致材料硬度等力学性能的显著下降, 从而对轮轨表面的耐磨损性能带来不利影响^[6]. 由于对轮轨滑动接触温升进行实验测量往往比较困难, 因此通过数值模拟和解析分析等非测试方式开展相应的温升研究具有重要的工程意义.

随着ABAQUS、ANASYS等商业有限元软件的普及, 轮轨滑动接触温升的有限元模拟研究受到了学者们的广泛关注. 然而, 由于轮轨接触斑非常狭小, 为了准确地模拟接触斑内的温度场, 在有限元模拟过程中必须对这一区域进行大规模网格加密处理, 从而所花费的时间成本较为高昂. 因此, 针对轮轨滑动接触温升问题的解析或半解析研究也得到了大量开展. 例如, Tanvir^[7]忽略了传热参数的温度相关性, 将Hertz接触压力的椭圆型分布函数近似成四阶多项式级数形式, 通过拉普拉斯变换给出了二维轮轨滑动摩擦温升分布的解析解. 基于相同的材料属性温度无关性假设, Knothe和Liebelt^[8]通过拉普拉斯变换和格林函数法, 给出了任意接触压力下的轮轨滑动温升问题的解析解, 讨论了法向接触压力分布形式和表面粗糙度对钢轨表面温度的影响. 基于Hertz弹性接触理论, Ertz和Knothe^[9]讨论了轮轨间换热系数等因素对轮轨摩擦温升的影响. Fischer等^[10]同时考虑了摩擦功和局部塑性变形对接触界面热源的贡献, 给出了轮轨滑动接触闪温的解析表达式. 此外, 考虑到轮轨间存在磨屑等杂物, Wu等^[11]基于三体微接触模型和接触温升理论, 讨论了杂物尺寸、杂物密度、表面粗糙度、相对速度以及法向载荷等参数对接触界面温度场的影响.

由于材料属性的温度相关性, 摩擦副材料的热力学性质在摩擦温升过程中也会发生变化, 考虑材料属性的温度相关性将进一步提升温升计算结果的准确性. 然而, 材料属性温度相关性的引入会导致热传导方程的非线性化, 大大增加了温升求解的困难性. 因此, 考虑材料属性的温度相关性, 必须首先求解非线性热传导方程. Ling和Rice^[12]考虑了热导率和比热容的温度相关性, 通过傅里叶余弦变换, 给出了表面具有移动热源半无限大体的准稳态温度场的迭代数值解法. 这种求解方法允许任意形式的温度相关性函数, 但相应的计算成本高, 在实际应用中受到一定的限制. 梁钰等^[13-14]基于特征正交分解方法, 对非线性瞬态热传导问题的有限元离散格式进行降阶处理, 从而提出了一种求解考虑传热参数温度相关性的二维非线性瞬态热传导问题的高效算法. 通过假设组份材料和孔隙填充物的热导率均随温度呈三次函数变化, 戴婷等^[15]采用微分求积法求解了一类含孔隙变厚度功能梯度圆盘的稳态温度场. 当热扩散系数随温度线性变化时, Fabre和Hristov^[16]基于热平衡积分法给出了一类一维热传导方程的解析解. 林府标等^[17]通过李群分析法, 考虑扩散率关于温度的幂函数形式, 得到了该类一维广义热传导方程的显示解析解和行波解. Ghasemi等^[18]则假定热导率是温度的线性函数, 将其斜率作为扰动小量, 也给出了一类一维非线性热传导问题的摄动解法. 此外, 微分变换法^[19]、有限差分法^[20-21]和近场动力学法^[22]也被应用于非线性热传导温度场的求解.

由上述分析可见, 已有研究存在不足: (1)接触压力通常被假定为满足Hertz弹性接触理论中的半椭圆型分布形式, 而在实际轮轨接触载荷作用下, 接触斑内将发生塑性变形, 其接触压力不再满足半椭圆型分布形式, 必须进行塑性修正^[23-25]. (2)材料参数往往被简化为温度无关的常量, 或仅考虑某一个材料参数的温度相关性. 实际上, 包含摩擦系数在内的多种材料属性均会随温度的改变而发生变化^[26], 在温升计算模型中同时考虑多种材料属性的温度相关性, 有望提升温升计算结果的准确性.

因此, 本文通过拟合实验数据, 获得了热导率、比热容以及摩擦系数的温度相关性函数, 并基于弹塑性接触理论, 构建了二维轮轨滑动温升有限差分计算模型, 分别讨论了对流换热系数、法向载荷、蠕滑率以及行车速度对钢轨表面温升的影响, 可为轮轨滑动温升的准确预测提供参考.

1 计算模型

在列车运行过程中, 经过一段时间的磨损和变形, 轮轨接触斑会变得比较狭长, 可近似为矩形形状, 从而可在二维框架下讨论轮轨滑动温升问题^[1]. 忽略车轮和钢轨的内部热源, 考虑热导率和比热容的温度相关性, 基于傅里叶导热定律, 二维情况下的轮轨摩擦导热方程可以写成

(1)$ \begin{eqnarray} \label{eq1} \frac{\partial }{\partial X}\left({\lambda \left( T \right)\frac{\partial T}{\partial X}} \right)+\frac{\partial }{\partial Z}\left( {\lambda \left( T \right)\frac{\partial T}{\partial Z}} \right)=\rho c\left( T \right)\frac{\mbox{d}T}{\mbox{d}t} \end{eqnarray}$

其中, $\lambda \left( T \right)$, $c\left( T \right)$和$\rho $分别表示热导率、比热容和密度, 温度场$T$是位置坐标$(x, z)$和时间$t$的函数, 即$T=T\left({x,z;t} \right)$. 令$a$表示接触斑纵向宽度的一半, $\left( {X, Z} \right)$为固定在钢轨表面的绝对坐标系, $\left( {x, z} \right)$是原点处于接触斑前沿的随动坐标系, 如图1所示.

图1

新窗口打开|下载原图ZIP|生成PPT

图1二维轮轨接触示意图

Fig.1Schematic of a two-dimensional wheel-rail contact system

当行车速度$v$保持恒定时, 坐标$\left( {X, Z} \right)$和$\left( {x, z} \right)$满足关系式$\left( {X,Z} \right)=\left( {x-vt,z} \right)$, 即

(2)$ \begin{eqnarray} \label{eq2} \left. \begin{array}{ll} {\dfrac{\partial T}{\partial X}=\dfrac{\partial T}{\partial x}, \dfrac{\partial T}{\partial Z}=\dfrac{\partial T}{\partial z} } \\[2mm] {\dfrac{\mbox{d}T}{\mbox{d}t}=\dfrac{\partial T}{\partial x}\dfrac{\mbox{d}x}{\mbox{d}t}+\dfrac{\partial T}{\partial t}=v\dfrac{\partial T}{\partial x}+\dfrac{\partial T}{\partial t}} \\ \end{array} \right\} \end{eqnarray}$

将式(2)代入方程(1)中, 可得

(3)$ \begin{eqnarray} \label{eq3} &&\frac{\partial }{\partial x}\left( {\lambda \left( T \right)\frac{\partial T}{\partial x}} \right)+\frac{\partial }{\partial z}\left( {\lambda \left( T \right)\frac{\partial T}{\partial z}} \right)= \\&&\qquad \rho c\left( T \right)\left( {v\frac{\partial T}{\partial x}+\frac{\partial T}{\partial t}} \right) \end{eqnarray}$

在稳态情况下, ${\partial T} / {\partial t}=0$, 方程(3)退化为

(4)$ \begin{eqnarray} \label{eq4} \frac{\partial }{\partial x}\left( {\lambda \left( T \right)\frac{\partial T}{\partial x}} \right)+\frac{\partial }{\partial z}\left( {\lambda \left( T \right)\frac{\partial T}{\partial z}} \right)=\rho c\left( T \right)v\frac{\partial T}{\partial x} \end{eqnarray}$

对方程(4)中的各物理量进行无量纲处理

(5)$ \begin{eqnarray} \left. \begin{array}{ll} \tilde{{\lambda }}=\dfrac{\lambda \left( T \right)}{\lambda \left( {T_{\max } } \right)},\ \tilde{{c}}=\dfrac{c\left( T \right)}{c\left( {T_{\max } } \right)}\\[2mm] \tilde{{T}}=\dfrac{T}{T_{\max } },\ \tilde{{x}}=\dfrac{x}{2a},\ \tilde{{z}}=\dfrac{z}{\delta } \end{array} \right\} \end{eqnarray}$

其中, $T_{\max } $为最高温度, $\delta =\sqrt {{2a\lambda \left( {T_{\max } } \right)} / {\left[ {\rho c\left( {T_{\max } } \right)v} \right]}}$为对应的热穿透深度^[8], 则方程(4)可简化为

(6)$ \begin{eqnarray} \label{eq6} \frac{\lambda \left( {T_{\max } } \right)}{2a\rho c\left( {T_{\max } } \right)v}\frac{\partial }{\partial \tilde{{x}}}\left( {\tilde{{\lambda }}\frac{\partial \tilde{{T}}}{\partial \tilde{{x}}}} \right)+\frac{\partial }{\partial \tilde{{z}}}\left( {\tilde{{\lambda }}\frac{\partial \tilde{{T}}}{\partial \tilde{{z}}}} \right)=\tilde{{c}}\frac{\partial \tilde{{T}}}{\partial \tilde{{x}}} \end{eqnarray}$

对于本文所采用的轮轨材料(见3.1节), 系数${\lambda \left( {T_{\max } } \right)} / {\left[ {2a\rho c\left( {T_{\max } } \right)v} \right]}<3.6\times 10^{-5}$, 因此方程(4)中的第一项可忽略^[8], 对应的方程(4)退化为

(7)$ \begin{eqnarray} \label{eq7} \frac{\partial }{\partial z}\left( {\lambda \left( T \right)\frac{\partial T}{\partial z}} \right)=\rho c\left( T \right)v\frac{\partial T}{\partial x} \end{eqnarray}$

在接触区域内, 轮轨之间因相对滑动而发生摩擦做功, 假设因轮轨滑动摩擦而损失的机械能全部转化为热能, 并以热传导的方式平均传递给车轮和钢轨, 忽略轮轨间的相互导热作用, 而接触区之外的表面则与空气发生对流换热. 因此, 边界条件可表示为

(8)$ \begin{eqnarray} \label{eq8} \left. \begin{array}{ll} T\left( {x=0,z} \right)=T\left( {x,z=\infty } \right)=T_0 \\ -\left[ {\lambda \left( T \right)\dfrac{\partial T}{\partial z}} \right]_{z=0} =q=\left\{ {{\begin{array}{ll} 0.5\mu \left( T \right)p\left( x \right)v_{\rm s}, & 0\leqslant x<2a \\ h\left( {T_{z=0} -T_0 } \right), & x\geqslant 2a \\ \end{array} }} \right. \\ \end{array} \right\}\ \ \ \end{eqnarray}$

式中, $q$为热通量; $\mu \left( T \right)$为摩擦系数; $h$为对流换热系数; $T_0 $为环境温度; $v_{\rm s}$为轮轨间的相对滑动速度, 且$v_{\rm s} =sv$, $s$为纵向蠕滑率; $p\left( x \right)$表示接触压力.

当轮轨材料满足双线性各向同性强化本构关系时, 接触压力$p\left( x \right)$可以写成^[27]

(9)$ \begin{eqnarray} \label{eq9} p\left( x \right)=\left\{\begin{array}{l} p_0^\ast \sqrt{1-\left( {\dfrac{x-a}{a}} \right)^2} ,\\\qquad a_{\rm p} \leqslant \left| {x-a} \right|

式中, $p_0^\ast$表示Hertz弹性接触下使得接触斑宽度达到$2a$时所对应的接触压力峰值; $p_{\rm y} $为临界接触压力, 且$p_{\rm y} \approx 1.6\sigma _{\rm y} $, $\sigma _{\rm y} $为屈服强度^[28]. 该模型假设当接触压力峰值超过$p_{\rm y} $时, 接触斑内开始出现局部塑性区域^[27]. $a_{\rm p} $为接触斑内塑性区域纵向宽度的一半

(10)$ \begin{eqnarray} \label{eq10} a_{\rm p} =a\sqrt {1-\left( {\frac{p_{\rm y} }{p_0^\ast }} \right)^2} \end{eqnarray}$

位于区域$0\leqslant \left| {x-a} \right|

(11)$ \begin{eqnarray} \label{eq11} \frac{1}{H^\ast }=\frac{1-\nu _1^2 }{H_1 }+\frac{1-\nu _2^2 }{H_2 }, \frac{1}{E^\ast }=\frac{1-\nu _1^2 }{E_1 }+\frac{1-\nu _2^2 }{E_2 } \end{eqnarray}$

式中, $\nu _1$和$\nu _2$分别表示车轮和钢轨的泊松比, $E_1 $和$E_2$分别为车轮和钢轨的弹性模量, $H_1 $和$H_2$分别为车轮和钢轨的塑性硬化模量.

对式(9)在接触区域进行积分, 可得

(12)$ \begin{eqnarray} \label{eq12} &&P_z -\frac{\pi E^\ast a^2}{4R}-\left( {1-\frac{H^\ast }{E^\ast }} \right)p_{\rm y} a\sqrt {1-\left( {\frac{2p_{\rm y} R}{E^\ast a}} \right)^2} + \\&&\qquad \frac{\left( {E^\ast -H^\ast } \right)a^2}{2R}\arcsin \sqrt {1-\left( {\frac{2p_{\rm y} R}{E^\ast a}} \right)^2} =0 \end{eqnarray}$

式中, $R$为车轮滚动半径; $P_z $为横向单位长度上的法向载荷, $P_z =P / {b_0 }$, $P$表示总法向载荷, $b_0$为矩形接触斑的横向等效宽度. 当$H^\ast =E^\ast $时, 方程(12)退化为

(13)$ \begin{eqnarray} \label{eq13} a=\sqrt {\frac{4RP_z}{\pi E^\ast}} \end{eqnarray}$

这与Hertz弹性接触理论在二维接触状态下的计算结果完全一致^[1,29].

由于速度$v$的影响, 车轮高速行进时的总法向载荷(动载)会高于车轮静止时的总法向载荷(静载), 二者之间满足^[30]

(14)$ \begin{eqnarray} \label{eq14} P=KP_{\rm s} ,\ \ K=1+\frac{0.4v}{100} \end{eqnarray}$

式中, $P_{\rm s}$为静载, 等于轴重之半.

2 计算模型的求解

观察式(7)和式(8)可以发现, 同时考虑热导率、比热容和摩擦系数的温度相关性, 使得方程(7)成为具有复杂边界条件(8)的含多个变系数的非线性抛物型偏微分方程, 直接进行求解十分困难. 因此, 基于基尔霍夫变换法, 对热导率的温度相关性函数进行积分

(15)$ \begin{eqnarray} \label{eq15} \vartheta =\vartheta \left( T \right)=\int_{T_0 }^T {\lambda \left( T \right)\mbox{d}T} \end{eqnarray}$

从而方程(7)可以转化成一种更简洁的形式

(16)$ \begin{eqnarray} \label{eq16} \frac{\kappa \left( \vartheta \right)}{v} {\frac{\partial ^2\vartheta }{\partial z^2}}=\frac{\partial \vartheta }{\partial x},\ \ \kappa \left( \vartheta \right)=\frac{\lambda \left( {\vartheta ^{-1}\left( T \right)} \right)}{\rho c\left( {\vartheta ^{-1}\left( T \right)} \right)} \end{eqnarray}$

其中$\vartheta ^{-1}\left( T \right)$表示$\vartheta \left( T \right)$的反函数.

类似地, 边界条件(8)也可以改写为

(17)$ \begin{eqnarray} \label{eq17} \left. \begin{array}{ll} \vartheta \left( {x=0,z} \right)=\vartheta \left( {x,z=\infty } \right)=0\\ -\left[ {\dfrac{\partial \vartheta }{\partial z}} \right]_{z=0} =q=\\ \qquad \left\{ \begin{array}{ll} 0.5sv\mu \left( {\vartheta ^{-1}\left( T \right)} \right)p\left( x \right), & 0\leqslant x<2a \\ h\left( {\left[ {\vartheta ^{-1}\left( T \right)} \right]_{z=0} -T_0 } \right), &x\geqslant 2a \\ \end{array}\right. \\ \end{array} \right\} \end{eqnarray}$

真实材料的$\kappa \left( \vartheta \right)$和$\mu \left( \vartheta \right)$需要通过拟合实验结果获得, 可能并不总是满足线性函数等简单形式, 这意味着摄动法以及拉普拉斯变换等常见方法可能不再适用于方程(16)的求解. 因此, 具有更强适用性的隐式差分法被用于该问题的求解.

根据边界条件(17), 可令计算区域为$\bar{\varOmega}=\lt[0,x_{\max} ]\times\lt[0, \bar{z}_{\max} ]$, 其中$x_{\max } $和$\bar{z}_{\max}$分别表示$x$和$z$的取值上限, $\bar{z}_{\max}=\infty$. 然而, 数学意义上的$\infty $并不能直接用于实际计算, 可用有限大值$z_{\max } $来近似替代. 因此, 实际计算区域为$\varOmega =\left[ {0,x_{\max } } \right]\times \left[ {0,z_{\max } } \right]$, 如图2所示. 用矩形网格对区域$\varOmega $进行离散, 沿$x$方向设置$m$+1个节点, 步长$\delta _1 ={x_{\max } } / m$; 沿$z$方向设置$n$+1个节点, 步长$\delta _2 ={z_{\max } } / n$. $\vartheta_j^i $表示在节点$\left( {i,j} \right)$, 即$x=\left( {i-1} \right)\delta _1 $且$z=\left( {j-1} \right)\delta _2 $处$\vartheta $的取值, $i=1$, 2, $\cdots$, $m+1$, $j=1$, 2, $\cdots$, $n+1$.

图2

新窗口打开|下载原图ZIP|生成PPT

图2计算区域的离散

Fig.2Schematic of calculation domain discretization

基于以差商近似代替微商的有限差分思想, 可将方程(16)和(17)转化成对应的差分方程形式

(18)$ \begin{eqnarray} \label{eq18} &&\left. \begin{array}{ll} \dfrac{\kappa _{j+0.5}^{i-1} }{v}\left[ {\dfrac{\vartheta _{j+1}^i -2\vartheta _j^i +\vartheta _{j-1}^i }{\left( {\delta _2 } \right)^2}} \right]=\dfrac{\vartheta _j^i -\vartheta _j^{i-1} }{\delta _1 } \\ \kappa _{j+0.5}^{i-1} =\kappa \left( {\dfrac{\vartheta _j^{i-1} +\vartheta _{j+1}^{i-1} }{2}} \right) \\ \end{array} \right\} \end{eqnarray}$

(19)$ \begin{eqnarray} \label{eq19} &&\left. \begin{array}{ll} \vartheta _j^1 =\vartheta _{n+1}^i =0 \\[2mm] -\dfrac{\vartheta _2^i -\vartheta _1^i }{\delta _2 }=q^{i-1}=\left\{ \begin{array}{l} {0.5sv\mu _1^{i-1} p^i, 0\leqslant x<2a} \\[1mm] {h\left( {T_1^{i-1} -T_0 } \right), x\geqslant 2a} \\ \end{array} \right. \\ \end{array} \right\} \end{eqnarray}$

其中, $\mu _1^{i-1} =\mu \left( {\vartheta _1^{i-1} } \right)$, $p^i=p\left({\left( {i-1} \right)\delta _1 } \right)$, $T_1^{i-1} =T\left( {\vartheta _1^{i-1} } \right)$, $i=2,3$, $\cdots$, $m+1$, $j=2,3$, $\cdots$, $n+1$.

将方程(19)代入方程(18)中, 便可以推导出非线性热传导方程(16)的二层隐式差分格式

(20)$ \begin{eqnarray} \label{eq20} A^{i-1}\bar{ \vartheta}^i+ F^{i-1}=\bar{ \vartheta}^{i-1} \end{eqnarray}$

其中, $ F^{i-1}=\left[ -\varepsilon \kappa _{2+0.5}^{i-1}\ \ q^{i-1}\delta _2 \ \ 0\ \ \cdots \ \ 0\ \ 0 \right]^{\rm T}$, $\bar{\vartheta}^{i-1}=\left[\vartheta^{i-1}_2\ \ \vartheta^{i-1}_3\ \ \cdots\ \ \vartheta^{i-1}_{n-1}\ \ \vartheta^{i-1}_n\right]^{\rm T}$, 上标T为转置符号, $ A^{i-1}$满足

(21)$ A^{i-1}=\begin{equation*} \begin{bmatrix} 1+\varepsilon \kappa _{2+0.5}^{i-1} & -\kappa _{2+0.5}^{i-1} & 0 & 0 & {\cdots } &0 & 0 & 0 \\ -\varepsilon \kappa _{3+0.5}^{i-1} & 1+2\varepsilon \kappa _{3+0.5}^{i-1} &-\varepsilon \kappa _{3+0.5}^{i-1} & 0 & {\cdots }&0 & 0 & 0 \\ 0 & -\varepsilon \kappa _{4+0.5}^{i-1} & 1+2\varepsilon \kappa _{4+0.5}^{i-1} & -\varepsilon \kappa _{4+0.5}^{i-1}& {\cdots }&0 & 0 & 0 \\ {\cdots } & {\cdots } & {\cdots } & {\cdots } & {\cdots } &{\cdots } & {\cdots } & {\cdots } \\ 0 & 0 & 0 & 0 & {\cdots } & -\varepsilon \kappa _{n-1+0.5}^{i-1} & 1+2\varepsilon \kappa _{n-1+0.5}^{i-1} & -\varepsilon \kappa _{n-1+0.5}^{i-1} \\ 0 & 0 & 0 & 0 & {\cdots } &0 & -\varepsilon \kappa _{n+0.5}^{i-1} & 1+2\varepsilon \kappa _{n+0.5}^{i-1} \end{bmatrix} \end{equation*} \\ \varepsilon ={\delta _1 } \Big/\left[ {\left( {\delta _2 } \right)^2v} \right], \ \ i=2,3, \cdots, m+1. $

在上述推导过程中, 并没有提前给定温度相关性函数$\lambda \left( T \right)$, $c\left( T \right)$和$\mu \left( T \right)$的形式, 这表明所给出的隐式差分格式(20)可以适用于具有任意温度相关性函数形式的热导率、比热容和摩擦系数. 联立求解方程(15)和(20), 便可获得任意节点位置处的温度$T_j^i $, $i=1$, 2, $\cdots$, $m+1$, $j=1$, 2, $\cdots$, $n+1$.

需要说明的是, 隐式差分法没有计算区域维度的限制, 当三维弹塑性接触压力分布函数已知时, 上述计算方法可较容易地推广至三维轮轨滑动接触温升的求解.

3 结果与讨论

3.1 材料属性的确定

考虑CHN60N型钢轨和轮径(2$R)$为890 mm的LM$_{\rm A}$型车轮. 环境温度$T_{0}=25{^\circ}$C, 材料密度$\rho =7790$ kg/m$^{3}$, 弹性模量$E_1 =E_2 = 209$ GPa, 硬化模量$H_1 =H_2 =20.9$ GPa, 泊松比$\nu _1 =\nu _2 =0.3$, 屈服强度$\sigma _{\rm y} =608$ MPa. 根据文献^[31]中的实验数据, 热导率$\lambda \left( T \right)$、比热容$c\left( T \right)$和摩擦系数$\mu \left( T\right)$的拟合曲线如图3(a)$\sim$图3(c)所示, 对应的拟合函数为

(22)$ \begin{eqnarray} \label{eq22} \left.\begin{array}{l} {\lambda \left( T \right)=a_{11} T^2+a_{12} T+a_{13} } \\ {c\left( T \right)=a_{21} T+a_{22} } \\ {\mu \left( T \right)=a_{31} \exp \left( {a_{32} T} \right)} \\ \end{array}\right\} \end{eqnarray}$

图3

新窗口打开|下载原图ZIP|生成PPT

图3热导率、比热容和摩擦系数的拟合曲线

Fig.3The fitting curves of thermal conductivity, specific heat capacity and friction coefficient

其中, 拟合系数$a_{11} =6.61\times 10^{-6}$, $a_{12} =0.01$, $a_{13} =47.71$, $a_{21} =0.13$, $a_{22} =487.50$, $a_{31}=0.35$, $a_{32}= -1.40\times 10^{-3}$.

尽管实际材料的屈服强度和弹性模量也会表现出温度相关性, 但其温度相关性函数的引入会大大增加弹塑性接触问题的复杂性, 导致接触压力分布函数的显式形式不再存在, 进而使轮轨滑动温升的计算变得更加困难. 因此, 出于简化问题的考虑, 这里忽略了屈服强度和弹性模量的温度相关性.

3.2 不同模型间的比较

不同轮轨滑动温升计算模型下的钢轨表面温升曲线如图4所示. 图中上标(ep), (e), (de)和(in)分别表示弹塑性接触、弹性接触、温度相关性和温度无关性, 后文亦同. 由于钢轨的摩擦热影响区深度很浅, 仅为2 mm左右^[32-33], 为了兼顾计算精度与工作效率, 此处可先试取计算区域深度$z_{\max }=4$ mm. 可以发现, 当忽略材料属性的温度相关性时, 所提出方法计算的温升曲线与文献^[8]得到的曲线符合较好, 表明本文所采用的隐式差分法和选取的计算区域深度是可行的. 同时, 考虑材料属性温度相关性时的温升曲线明显偏离忽略材料属性温度相关性时的曲线, 二者峰值温升之间的差值甚至可达640${^\circ}$C. 由此可见, 在轮轨温升计算模型中考虑材料属性的温度相关性是非常有必要的. 另外, 在较高温升下, 摩擦系数会出现剧烈的衰减现象(见图3(c)), 对应的热通量随之降低, 进而对钢轨表面温度的继续升高产生抑制作用, 最终造成考虑摩擦系数温度相关性时的温升计算结果远远低于忽略该特性时的计算结果. 因此, 在轮轨滑动温升模拟中十分有必要考虑摩擦系数的温度相关性.

图4

新窗口打开|下载原图ZIP|生成PPT

图4不同轮轨滑动温升计算模型下的钢轨表面温升曲线($P_{\rm s} =90$ kN, $h=100$ W$\cdot$(m$^{2}\cdot {^\circ}{\rm C})^{-1}$, $v=350$ km/h, $s=8$%)

Fig.4The temperature rise curves of rail surface for different calculation models of wheel-rail sliding temperature rise ($P_{\rm s} =90$ kN, $h=100$ W$\cdot$(m$^{2}\cdot {^\circ}{\rm C})^{-1}$, $v=350$ km/h, $s=8$%)

另外, 与Hertz弹性接触理论相比, 当法向载荷过大、接触斑内表面材料出现局部塑性变形时, 弹塑性接触理论所给出的峰值接触压力更小、接触斑宽度更大(见图5), 进而造成钢轨表面峰值温度所在位置朝接触斑后沿方向移动. 由于所采用的高强度钢轨CHN60N具有很高的临界接触压力, 相应的塑性变形区域很小, 因而两种接触压力形式下的最大滑动温升比较接近. 但对于屈服强度或临界接触压力较小的接触体材料, 这两种温升计算结果间的差异将更加显著.

图5

新窗口打开|下载原图ZIP|生成PPT

图5弹性接触压力与弹塑性接触压力的比较

Fig.5The comparison of elastic contact pressure and elasto-plastic contact pressure

3.3 对流换热系数对滑动温升的影响

在轮轨滑动温升计算模型中, 钢轨非接触区域表面的边界条件直接取决于对流换热系数, 然而该物理参数往往受到车速、风速以及环境温度等诸多因素的影响, 其具体取值尚无定论^[34]. Wang等^[35]通过比例模型实验测得了钢轨表面对流换热系数在列车高速行驶下的取值范围为80$\sim $160 W/(m$^{2}\cdot {^\circ}{\rm C})$. 因此,图6和图7给出了不同对流换热系数下的钢轨表面温升和热通量分布曲线.

图6

新窗口打开|下载原图ZIP|生成PPT

图6不同对流换热系数下的钢轨表面温升曲线($P_{\rm s}=90$ kN, $v=350$ km/h, $s=8$%)

Fig.6The temperature rise curves of rail surface with different convection coefficients ($P_{\rm s}=90$ kN, $v=350$ km/h, $s=8$%)

图7

新窗口打开|下载原图ZIP|生成PPT

图7不同对流换热系数下的钢轨表面热通量分布曲线($P_{\rm s}=90$ kN, $v=350$ km/h, $s=8$%)

Fig.7The heat flux distribution of rail surface with different convection coefficients ($P_{\rm s}=90$ kN, $v=350$ km/h, $s=8$%)

不难看出, 钢轨表面峰值温度一直出现在接触斑内, 并靠近接触斑后沿, 这与文献^[36,37]的结论相一致. 同时, 不同对流换热系数下的温升曲线几乎完全重合, 表明当列车高速行驶时, 所讨论的对流换热系数变化对轮轨滑动温升的影响甚微, 这与文献^[38]的结论一致. 进一步观察图7可知: 在车轮高速行进情况下, 摩擦功率对热通量的影响会显著高于对流换热系数变化的影响, 进而体现为轮轨滑动温升对对流换热系数变化的不敏感性, 因此在后续的讨论中对流换热系数将被设为定值100 W/(m$^{2}\cdot{^\circ}{\rm C})$.

3.4 法向载荷对滑动温升峰值的影响

不同法向载荷下的钢轨表面温升峰值如图8所示. 从图8中可以发现, 在不同计算模型下, 对应的温升峰值均会随法向载荷的增大近似呈线性上升, 但温升取值则存在明显的差异. 与其他计算模型相比, 忽略材料属性温度相关性时的温升峰值计算结果最大, 表明忽略材料属性的温度相关性会大大高估轮轨间的滑动温升, 考虑材料属性的温度相关性可以减小这种高估程度, 进而提高温升计算结果的准确性. 同时, 仅考虑摩擦系数温度相关性时的温升峰值计算结果明显低于仅考虑热导率和比热容的温度相关性时的结果, 与同时考虑热导率、比热容和摩擦系数的温度相关性时的温升计算结果更为接近, 表明就式(22)所表征的材料而言, 摩擦系数的温度相关性对最终温升计算结果的影响显著强于热导率和比热容. 然而, 即便已经在温升计算模型中引入了温度相关的摩擦系数, 继续考虑热导率和比热容的温度相关性仍会影响最终的温升峰值计算结果, 二者之间的差值可接近70${^\circ}$C, 这再一次凸显了同时考虑多种材料属性的温度相关性对于提高轮轨滑动温升模拟结果准确性的意义.

图8

新窗口打开|下载原图ZIP|生成PPT

图8钢轨表面最大温升随法向载荷的变化($v=350$ km/h, $s=8$%)

Fig.8The maximum temperature rises of rail surface versus vertical loads ($v=350$ km/h, $s=8$%)

3.5 蠕滑率对滑动温升峰值的影响

钢轨表面温升峰值随蠕滑率的变化曲线如图9所示. 可以看出, 温升峰值随蠕滑率的增大而上升, 这是因为在相同接触压力和行车速度的情况下, 蠕滑率的增大会引起轮轨接触点对间相对速度的增大, 提高摩擦热通量, 进而导致钢轨表面温度的上升. 当蠕滑率较低($s\leqslant 1$%)时, 轮轨滑动温升并不显著, 热导率、比热容和摩擦系数的变化量很小, 不同计算模型的温升峰值曲线接近重合; 但随着蠕滑率的进一步增大, 钢轨表面温度达到较高水平, 材料属性的温度相关性开始凸显, 因而该4组曲线互相偏离, 且偏离差异随蠕滑率的增加而增大. 忽略材料属性温度相关性时的钢轨滑动温升峰值与蠕滑率呈线性关系, 这与文献^[8]中的结论相一致; 考虑材料属性温度相关性时的温升峰值-蠕滑率曲线则体现出非线性特征, 且位于考虑材料属性温度无关性时的曲线的下方.

图9

新窗口打开|下载原图ZIP|生成PPT

图9钢轨表面最大温升随蠕滑率的变化($P_{\rm s}=90$ kN, $v=350$ km/h)

Fig.9The maximum temperature rises of rail surface versus creepages ($P_{\rm s}=90$ kN, $v=350$ km/h)

3.6 行车速度对滑动温升峰值的影响

图10给出了钢轨表面温升峰值随行车速度的变化曲线. 与蠕滑率对钢轨温升峰值的影响相类似, 当法向载荷和蠕滑率恒定时, 行车速度的增大会引起轮轨接触点对间相对速度的增大, 提高摩擦热通量, 进而导致钢轨表面温升峰值的增加. 因此, 热导率、比热容和摩擦系数在高速下的变化量比低速情况更为显著, 从而导致不同温升模型计算结果之间的差异随行车速度的增大而加大. 在轮轨滑动温升计算模型中引入材料属性的温度相关性函数, 可以避免做出过分高估的温升评价, 这一特征在图10中表现为考虑材料属性温度相关性时的温升峰值计算结果明显低于忽略温度相关性时的计算结果. 同时, 仅考虑摩擦系数温度相关性时的温升峰值曲线远远低于仅考虑热导率和比热容温度相关性时的结果, 并接近于同时考虑这三者的温度相关性时的曲线, 表明图3(c)所呈现出的摩擦系数热衰减特性对温升计算结果的影响显著高于热导率和比热容. 因此, 在计算模型中增加摩擦系数的温度相关性函数可大幅提高温升预测的准确性.

图10

新窗口打开|下载原图ZIP|生成PPT

图10钢轨表面最大温升随行车速度的变化($P_{\rm s}=90$ kN, $s=8$%)

Fig.10The maximum temperature rises of rail surface versus train speeds ($P_{\rm s}=90$ kN, $s=8$%)

尽管所构建的二维轮轨滑动接触温升计算模型同时考虑了接触压力的塑性修正和热导率、比热容以及摩擦系数的温度相关性, 但忽略了屈服强度和弹性模量的温度相关性等因素对温升结果的影响, 与真实的轮轨传热状态尚有一定的偏差; 同时, 为了便于获得接触压力的解析解, 材料模型采用了简化的双线性各向同性强化模型来考虑车轮和钢轨的弹塑性响应, 而实际的车轮和钢轨的弹塑性响应是非线性的, 这将对计算结果产生一定的影响. 上述不足之处在目前提出的计算模型中尚无法被考虑, 需要在进一步的工作中予以完善.

4 结论

在二维弹塑性接触模型框架下, 同时考虑热导率、比热容和摩擦系数的温度相关性, 基于基尔霍夫变换方法, 将复杂的非线性Fourier导热方程转化成一种以热导率关于温度的积分函数为待求量的简单方程形式, 进而构建了相应的轮轨滑动接触温升计算模型. 基于隐式差分法, 推导出了一种不限材料温度相关性函数形式的统一差分求解格式, 并分别讨论了对流换热系数、法向载荷、蠕滑率以及行车速度对钢轨表面滑动温升的影响, 主要结论如下:

(1) 由于钢轨温升峰值总是位于接触斑内的后沿附近, 并且当接触斑内出现塑性变形区域时, 弹塑性接触理论所给出的接触斑宽度较弹性接触理论更大, 相应的温升峰值所在位置距接触斑前沿更远.

(2) 当行车速度处于较高水平时, 对流换热系数对轮轨滑动温升的影响甚微, 在温升计算模型中将对流换热系数设为定值是可行的. 蠕滑率和行车速度的增加, 均会引起轮轨接触点对间相对速度的增大, 摩擦热通量随之上升, 进而导致钢轨表面温升峰值的增加. 另外, 钢轨的温升峰值也会随法向载荷的增加而近似线性上升.

(3) 当实际的轮轨滑动温升处于较低水平时, 材料属性的温度相关性对温升计算结果的影响并不明显; 当实际温升较显著时, 在计算模型中同时考虑多种材料属性的温度相关性可以避免过分高估的温升评价, 且摩擦系数的温度相关性对温升计算结果的影响会明显高于热导率和比热容.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

金学松 , 刘启跃 . 轮轨摩擦学. 北京 : 中国铁道出版社 , 2004

[本文引用: 3]

( Jin

Xuesong

, Liu

Qiyue

. Tribology of Wheel and Rail . Beijing : China Railway Publishing House , 2004 (in Chinese))

[本文引用: 3]

[2]

, Wen

, Li

, et al .

Thermo-elastic-plastic finite element analysis of wheel/rail sliding contact

Wear, 2011 , 271 : 437 - 443

DOIURL[本文引用: 1]

[3]

Srivastava

, Sarkar

, Ranjan

Effects of thermal load on wheel-rail contacts: A review

Journal of Thermal Stresses, 2016 , 39 : 1389 - 1418

DOIURL[本文引用: 1]

[4]

Ahlström

, Karlsson

Modelling of heat conduction and phase transformations during sliding of railway wheels

Wear, 2002 , 253 : 291 - 300

DOIURL[本文引用: 1]

[5]

Chen

, He

, Zhao

, et al .

The influence of wheel flats formed from different braking conditions on rolling contact fatigue of railway wheel

Engineering Failure Analysis, 2018 , 93 : 183 - 199