钢轨砂带磨削温度场建模与仿真

Modeling and Simulation of Temperature Field in Rail Belt Grinding

摘要

结合钢轨与砂带的接触几何关系，建立接触区域磨削深度和轮廓计算模型，阐明磨削工艺参数对接触区参数的影响规律.随着砂带磨削半径的增加，不同钢轨廓形处接触区域面积呈对数增加，接触区域轴长半径位置也在不同截面上变化.结合钢轨砂带磨削过程特性，求解接触区域热流密度计算模型并进行理论验证.基于接触区域磨削深度和轮廓计算模型，以及区域热流密度计算模型，应用瞬时点热源温度场、连续作用点热源温度场、连续作用移动点热源温度场对接触区域连续作用移动面热源温度场进行离散化求解.研究结果表明，在设定磨削工艺参数下，钢轨磨削表面的仿真和理论温度的变化趋势相似，且几乎在同一时间达到温度最大值，最高温度的仿真和理论计算的相对误差为6.14%，验证了本文理论模型和仿真的正确性.

Abstract

Based on the geometric relationship between the steel rail and the abrasive belt， a calculation model for the grinding depth and profile of the contact area is established to clarify the influence of grinding process parameters on the contact area parameters. As the grinding radius of the abrasive belt increases， the contact area of different steel rail profiles increases logarithmically， and the position of the axis length radius of the contact area changes on different cross-sections. Based on the characteristics of the rail abrasive belt grinding process， a calculation model for heat flux density in the contact area was solved and theoretically verified. Based on the contact area grinding depth and profile calculation， and regional heat flux density calculation model solved above， the instantaneous point heat source temperature field， continuous point heat source temperature field， and continuous action moving point heat source temperature field were applied to discretize the temperature field of the continuous action moving surface heat source in the contact area. Research shows that under the set grinding process parameters， the simulation and theoretical temperature changes on the surface of rail grinding are similar and almost reach the maximum temperature at the same time. The relative error between the simulation and theoretical maximum temperature is 6.14%， which verifies the correctness of the above theoretical model and simulation.

Graphical abstract

关键词

磨削加工 / 温度分布 / 磨削过程 / 热流密度

钢轨在服役过程中产生多种损伤与缺陷 ^［1］，严重影响钢轨的疲劳寿命^［2］，利用磨削去除并修复廓形是改善轮轨匹配关系最有效的手段^［3-7］.钢轨磨削过程中，在磨削区产生的高温导致钢轨发生不均匀且不可逆的塑性变形，冷却后会在钢轨表面和内部形成残余应力，残余应力会影响钢轨的耐磨性、耐腐蚀性以及抗疲劳强度^［8］等，因此开展钢轨磨削区域温度场分布研究对提高钢轨磨削质量和服役性能具有重大意义.

国内外研究者已经对磨削温度场进行了研究，张磊^［9］通过建立磨削热源分布的综合模型得出：在不同打磨工艺参数下，磨削热源分布的综合模型可以转化为不同热源形式. Zhang等^［10］通过研究得出：三角形热源分布更加符合实际打磨工况，热源峰值在磨削区的位置代表不同磨削工艺. 聂蒙等^［11］对钢轨打磨列车打磨时的温度分布进行了分析，求解温度场的方程，并通过实验验证方程的正确性.Zhou等^［12］通过点线面热源积分方法建立了钢轨端面磨削过程中磨削区域温度场解析热模型，并通过实验验证了其正确性. 张青等^［13］根据钢轨打磨列车打磨的实际工况，讨论得出：钢轨砂轮打磨温度升高很快，降温相对较慢. 张子舆^［14］建立了钢轨砂轮打磨过程的有限元模型，对钢轨砂轮打磨后单个砂轮和多个砂轮进行温度场仿真分析，得到不同工艺参数下温度场的分布状况，并对仿真结果进行实验验证. 上述研究者通过理论、仿真和实验对矩形区域、端面磨削区域或整个椭圆形状的磨削区域的温度场进行了分析，但在实际钢轨磨削过程中，磨削区域为半个椭圆，因此需要对其进行进一步分析.

本文通过接触区域磨削深度和轮廓计算模型、接触区域热流密度计算模型对半椭圆连续作用移动面热源温度场进行了离散化求解，并通过仿真进行了理论模型验证.

1 磨削接触区域分析

钢轨磨削过程中，磨削温度是影响残余应力的主要因素^［15］，为了研究钢轨磨削过程中的温度场分布，首先对磨削区域进行分析. 不同钢轨廓形磨削示意图如图1所示.

以60 kg/m钢轨为例，从轨顶到图1中x方向钢轨的廓形半径分别为R300 mm、R80 mm和R13 mm. 建立如图1所示的坐标系，假设接触轮半径与磨削砂带厚度之和为R₁，在钢轨定点磨削过程中，砂带与不同钢轨廓形接触处的yOz截面示意图如图2所示.

由图2可知：

r_{\mathrm{1}}=\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}{h}_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{1}}-h_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{1}}^{\mathrm{2}}}（1）

式中：r_{\mathrm{1}}为yOz截面接触区域轴长半径；h_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{1}}为yOz截面最大磨削深度.

由图2可知，轴长半径从0增加到r₁时，磨削深度h随轴长半径的变化为：

h\left(y\right)=h_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{1}}-R_{\mathrm{1}}+\sqrt[]{R_{\mathrm{1}}^{\mathrm{2}}-y^{\mathrm{2}}}，y\in \left(-r_{\mathrm{1}},r_{\mathrm{1}}\right)（2）

砂带与钢轨接触处的xOz截面示意图如图3所示.由图3可知：

r_{\mathrm{2}}=\sqrt[]{\mathrm{2}{R}_{\mathrm{2}}{h}_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{2}}-h_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{2}}^{\mathrm{2}}}（3）

式中：r_{\mathrm{2}}为xOz截面接触区域轴长半径；h_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{2}}为xOz截面最大磨削深度，即为砂带磨削深度；R_{\mathrm{2}}为钢轨与砂带接触时的钢轨廓形半径.

由图3可知，轴长半径从0增加到r_{\mathrm{2}}时，磨削深度h随轴长半径的变化为：

h\left(x\right)=h_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{2}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}，x\in \left(-r_{\mathrm{2}},r_{\mathrm{2}}\right)（4）

假设xOz截面的原点在砂带与钢轨初始接触点处，h（x）为任一xOz截面与相应yOz截面交线处的最大磨削深度，在磨削深度为h，即h_max2=h时，xOz截面任意点最大磨削深度为：

h\left(x,y\right)=h-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}-R_{\mathrm{1}}+\sqrt[]{R_{\mathrm{1}}^{\mathrm{2}}-y^{\mathrm{2}}}（5）

x\in \left(-r_{\mathrm{2}},r_{\mathrm{2}}\right)，y\in \left(-r_{\mathrm{1}},r_{\mathrm{1}}\right)

分别取R_{\mathrm{1}}=500 mm、R_{\mathrm{2}}=300 mm、h=0.1 mm，R_{\mathrm{1}}=300 mm、R_{\mathrm{2}}=300 mm、h=0.1 mm，R_{\mathrm{1}}=100 mm、R_{\mathrm{2}}= 300 mm、h=0.1 mm，得到磨削过程中体现磨削深度的MATLAB磨削区域仿真示意图如图4所示.

在钢轨磨削过程中，磨削深度一般为h=0.1 mm，当R₁=500 mm、R₂=300 mm时，通过MATLAB对上述求解接触区域的轮廓进行提取，应用Origin对数据进行非线性隐函数曲线拟合，应用Ellipse模型迭代至收敛，结果如图5所示.

由Ellipse模型所拟合出的方程为椭圆一般式， a=9.999 5 mm，b=7.744 58 mm，x_c=3.406 57×10^-6 mm， y_c=7.661 44×10^-4 mm，R²≈1.其中：a为椭圆长半轴，b为椭圆短半轴，x_c为椭圆中心点的x值，y_c为椭圆中心点的y值.

x_c与y_c的值很小，对结果的作用很小，此处忽略. 对结果进行优化，考虑到实际接触区域轮廓点在椭圆方程中的收敛性，对椭圆轮廓进行放大，优化的椭圆标准表达式为：

\frac{x^{\mathrm{2}}}{\mathrm{60}}+\frac{y^{\mathrm{2}}}{\mathrm{100}}=\mathrm{1}（6）

通过MATLAB对轮廓点的收敛性进行验证，如图6所示. 由图6可知，优化后椭圆上的每个点包含轮廓点，最大相对误差为5.9%，平均相对误差为0.5%，说明优化后的表达式是合理的.

由上述所建立的接触区域磨削深度和廓形模型可知，接触区域与磨削深度、砂带磨削半径（R₁）和钢轨廓形半径有关.在实际钢轨磨削过程中，磨削接触区域为半个椭圆，当h=0.1 mm时，接触区域面积S=\mathrm{\pi }{r}_{\mathrm{1}}{r}_{\mathrm{2}}/\mathrm{2}.不同钢轨廓形半径下，磨削接触区域面积如图7所示.

由图7可知，随着砂带磨削半径的增加，因r₁不断增加，不同钢轨廓形上的磨削接触区域面积增大；随着钢轨廓形半径的增加，因r₂不断增加，接触区域面积的累积增加量增大，接触区域面积不断增加.通过Origin非线性曲线拟合可知，随着砂带磨削半径的增加，不同钢轨廓形处的接触区域面积呈对数增加.

在给定磨削深度且砂带宽度大于2r₂条件下， x方向的椭圆半径是固定的，在固定磨削深度为 0.1 mm时，不同钢轨廓形半径下，随着砂带磨削半径的增加，接触区域在y轴方向上的轴长半径的变化如图8所示.

由图8可知，当磨削深度h=0.1mm时，R300 mm、R80 mm和R13 mm的轴长半径r₁分别为7.74 mm、 4.00 mm和1.61mm. 随着砂带磨削半径的增加，r₁不断增加，当r₁>r₂时，接触区域如图4（a）所示；当r₁=r₂时，接触区域如图4（b）所示；当r₁<r₂时，接触区域如图4（c）所示.

2 磨削过程分析

砂带的磨削过程和磨粒的扫掠示意图如图9所示.

最大磨屑厚度h_{\mathrm{c}\mathrm{u}.\mathrm{m}\mathrm{a}\mathrm{x}}和砂带与钢轨的接触弧长l_g的计算公式分别为：

h_{\mathrm{c}\mathrm{u}.\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{2}L\frac{v_{\mathrm{w}}}{v_{\mathrm{s}}}\sqrt[]{\frac{a_{\mathrm{e}}}{d_{\mathrm{e}}}}（7）

l_{\mathrm{g}}=\sqrt[]{a{d}_{\mathrm{e}}}（8）

式中：L为相邻两个磨粒之间的距离；a_{\mathrm{e}}为磨削深度，a_{\mathrm{e}}=h_{\mathrm{m}\mathrm{a}\mathrm{x}}；d_{\mathrm{e}}为砂带等效直径；v_{\mathrm{w}}为工件移动速度；v_{\mathrm{s}}为砂带速度.

弧的连线l_{\mathrm{c}}如图9所示，由图2可知，\mathrm{0}\le r_{\mathrm{1}}\left(x\right)\le a，当r_{\mathrm{1}}\left(x\right)从0增加到a的过程中，l_{\mathrm{c}}和l_g的差值逐渐变大；当x=\mathrm{0}，r_{\mathrm{1}}=a时，l_{\mathrm{c}}和l_g的差值最大. 在上述参数的磨削条件下，通过计算可得，l_{\mathrm{c}}与l_g的相对误差为0.13%，因此，可以用l_{\mathrm{c}}表示l_g.

由图9和图2可知，\mathrm{0}\le r_{\mathrm{1}}\left(x\right)\le a，当r₁（x）从0增加到a的过程中，l_{\mathrm{c}}与水平线的夹角逐渐增加，{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{0.01}. 因此，在本研究过程中，可应用不同的r₁（x）表示l_g，随着r₁（x）的增加，磨削深度近似于线性增加^［16］.

3 热量分布规律

通过能量守恒定律可知，钢轨在磨削过程中，时间t内产生的总热量为：

Q=p_{\mathrm{m}}\times t（9）

式中：p_{\mathrm{m}}为磨削砂带电机功率；t为磨削时间.

在磨削过程中，由于砂带与钢轨之间的热传导、砂带和钢轨与气体流体的热对流和热辐射会产生热量损失，根据文献［16］可知，流入钢轨的总能量q_{\mathrm{t}}为：

q_{\mathrm{t}}=q_{\mathrm{w}}+q_{\mathrm{s}}+q_{\mathrm{c}\mathrm{h}}+q_{\mathrm{f}}（10）

式中：q_{\mathrm{w}}为流入钢轨的热量；q_{\mathrm{c}\mathrm{h}}、q_{\mathrm{s}}和q_f分别为流入磨屑、砂带和热传导的能量.

其中：

q_{\mathrm{c}\mathrm{h}}=\rho c{T}_{\mathrm{m}\mathrm{p}}\frac{a_{\mathrm{e}}\cdot v_{\mathrm{w}}}{r_{\mathrm{1}}\left(x\right)}（11）

q_{\mathrm{f}}=\frac{\mathrm{2}}{\mathrm{3}}{h}_{\mathrm{f}}{T}_{\mathrm{m}\mathrm{a}\mathrm{x}}（12）

式中：\rho为工件密度；c为比热容，在室温下，c=473 J/（kg·K）；T_{\mathrm{m}\mathrm{p}}为工件熔化温度；h_f为接触区域流体冷却的对流系数；T_{\mathrm{m}\mathrm{a}\mathrm{x}}为工件区域温度的最大值.

工件和砂带的热通量分配系数为：

R_{\mathrm{w}\mathrm{s}}={\left(\mathrm{1}+\frac{\mathrm{0.974}{k}_{\mathrm{g}}}{{\beta }_{\mathrm{w}}\sqrt[]{r_{\mathrm{0}}{v}_{\mathrm{s}}}}\cdot \frac{\mathrm{1}}{F}\right)}^{-\mathrm{1}}（13）

式中：k_g为磨料的传导率；r_{\mathrm{0}}为接触半宽度；{\beta }_{\mathrm{w}}为工件材料的热特性，{\beta }_{\mathrm{w}}=\sqrt[]{k_{\mathrm{w}}{\rho }_{}{c}_{}}，k_{\mathrm{w}}为工件传导率；F为瞬态函数，稳定状态下，F=\mathrm{1}.

钢轨磨削产生的总热量为Q，表面总平均热流密度（单位时间通过单位面积的热能）为：

q_{\mathrm{t}\mathrm{0}}=\frac{Q}{tS}=\frac{p_{\mathrm{m}}t}{tS}=\frac{p_{\mathrm{m}}}{S}（14）

式中：S为接触面积.

流入钢轨的平均热流密度为：

\begin{array}{l}{q}_{\mathrm{w}\mathrm{0}}=\left(\frac{\mathrm{2}{p}_{\mathrm{m}}}{\mathrm{\pi }ab}-\frac{\mathrm{2}}{\mathrm{3}}{h}_{\mathrm{f}}{T}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\rho c{T}_{\mathrm{m}\mathrm{p}}\frac{a_{\mathrm{e}}{v}_{\mathrm{w}}}{r_{\mathrm{1}}\left(x\right)}\right)\times \\ {\left(\mathrm{1}+\frac{\mathrm{0.974}{k}_{\mathrm{g}}}{{\beta }_{\mathrm{w}}\sqrt[]{r_{\mathrm{0}}{v}_{\mathrm{s}}}}\cdot \frac{\mathrm{1}}{F}\right)}^{-\mathrm{1}}\end{array}（15）

上述给定参数下，q_{\mathrm{t}\mathrm{0}}与q_{\mathrm{c}\mathrm{h}}的数量级为{\mathrm{10}}^{\mathrm{7}}，q_f的数量级为{\mathrm{10}}^{\mathrm{2}}，此处忽略q_f的作用.

由于每一纵截面（yOz平面）上的r_{\mathrm{1}}\left(x\right)不同，则y的取值范围不同［\mathrm{0}\le y\le r_{\mathrm{1}}\left(x\right)］，对应截面上的最大磨削深度h\left(x\right)不同.

当x=x_{\mathrm{1}}时，

r_{\mathrm{1}}\left(x_{\mathrm{1}}\right)=\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}（16）

h\left(x_{\mathrm{1}}\right)=h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{2}}}（17）

当x=x_{\mathrm{2}}时，

r_{\mathrm{1}}\left(x_{\mathrm{2}}\right)=\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}（18）

h\left(x_{\mathrm{2}}\right)=h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{2}}}（19）

随着r_{\mathrm{1}}\left(x\right)的增加，磨削深度近似于线性增加. 最大磨屑厚度h_{\mathrm{c}\mathrm{u}.\mathrm{m}\mathrm{a}\mathrm{x}}和砂带与工件的接触长度l_{\mathrm{c}}的比值为磨屑厚度沿长度方向线性增加的斜率，在不同截面下的磨削深度的斜率为：

k_{\mathrm{1}}=\frac{\mathrm{2}s\sqrt[]{\frac{h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{2}}}}{d_{\mathrm{e}}}}}{\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{1}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}}（20）

k_{\mathrm{2}}=\frac{\mathrm{2}s\sqrt[]{\frac{h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{2}}}}{d_e}}}{\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}}（21）

\frac{k_{\mathrm{1}}}{k_{\mathrm{2}}}=\frac{r_{\mathrm{1}}\left(x_{\mathrm{2}}\right)}{r_{\mathrm{1}}\left(x_{\mathrm{1}}\right)}\sqrt[]{\frac{h\left(x_{\mathrm{1}}\right)}{h\left(x_{\mathrm{2}}\right)}}=\frac{\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}-h\left(x_{\mathrm{2}}\right)}}{\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}-h\left(x_{\mathrm{1}}\right)}}（22）

由于\mathrm{0}\le h\left(x\right)\le h_{\mathrm{m}\mathrm{a}\mathrm{x}}，当h=0.1 mm， R₁=500 mm，R₂=300 mm时，\mathrm{0.999}\mathrm{995}\le k_{\mathrm{1}}/k_{\mathrm{2}}\le \mathrm{1.000}\mathrm{05}，即k_{\mathrm{1}}/k_{\mathrm{2}}\approx \mathrm{1}，因此，在不同截面下，磨屑厚度和接触长度的比值相同，即具有相同斜率.

由上述磨削特点和Rowe等^［17］研究可知，在不同截面r_{\mathrm{1}}\left(x\right)上的热通量用三角形热源公式表示最合适：

q(y)=\frac{\mathrm{2}{q}_{\mathrm{w}\mathrm{0}}}{r_{\mathrm{1}}\left(x\right)}y，\mathrm{0}\le y\le r_{\mathrm{1}}(x)，-b\le x\le b（23）

假设p_{\mathrm{m}}=2 000 W，由式（11）、式（13）~式（15）可知，q_{\mathrm{t}\mathrm{0}}=\mathrm{1.64}\times {\mathrm{10}}^{\mathrm{7}} W/m²，q_{\mathrm{c}\mathrm{h}}=\mathrm{2.9}\times {\mathrm{10}}^{\mathrm{6}} W/m²，R_{\mathrm{w}\mathrm{s}}=\mathrm{0.86}. 流入钢轨的平均热通量q_{\mathrm{w}\mathrm{0}}=\mathrm{1.16}\times {\mathrm{10}}^{\mathrm{7}} W/m²，流入钢轨总体热量Q_{\mathrm{w}}=q_{\mathrm{w}\mathrm{0}}S，即Q_{\mathrm{w}}=1 409.61 W.

对于任一x上，在y的取值范围内，每一点的热通量可以通过公式求出，每一条线上热通量可以通过积分求出. 区域热通量计算式为：

Q_{\mathrm{w}\mathrm{1}}={\displaystyle \underset{-b}{\overset{b}{\int }}}{\displaystyle \underset{\mathrm{0}}{\overset{r_{\mathrm{1}}(x)}{\int }}}q(y)\mathrm{d}y\mathrm{d}x（24）

由式（24）可得，Q_{\mathrm{w}\mathrm{1}}=1 411.24 J，相对误差为0.12%. 因此应用三角形热源公式对椭圆热源进行计算具有较高的准确性.

4 砂带磨削温度计算

4.1 瞬时点热源温度场

假设钢轨为无限大且无内热源，在物体内某一点突然产生一个热量Q，随之热源停止继续发热，已产生的热量Q继续在物体内进行传播，初始时间后，通过解析法求解的瞬时点热源温度场可得物体内部任意一点M（U，V，W）处经过时间t后的温度T为^［18］：

T(R,t)=\frac{Q}{c\rho {(\mathrm{4}\mathrm{\pi }\alpha t)}^{\mathrm{3}/\mathrm{2}}}{\mathrm{e}}^{-R^{\mathrm{2}}/\mathrm{4}\alpha t}（25）

假设热源在（x，y，z）处，则

R^{\mathrm{2}}={(x-U)}^{\mathrm{2}}+{(y-V)}^{\mathrm{2}}+{(z-W)}^{\mathrm{2}}

式中：α为导温系数，α=k/cρ，k为热导率，ρ为工件密度.以材料U71Mn为例，取k=49 W/（m·K），ρ= 7 810 kg/m³.

4.2 连续作用点热源温度场

连续作用点热源如图10所示.

如果点热源发生后一直持续下去，单位时间内产生的热量为一定值，在任意时间t处的热流密度为q_t，从t=\mathrm{0}到t=t_{\mathrm{2}}内取任意微元\mathrm{d}t，则\mathrm{d}t=t_{\mathrm{1}}-t_{\mathrm{0}}，在时间\mathrm{d}t内产生的热量\mathrm{d}Q=q_t\mathrm{d}t，点热源持续时间为(\mathrm{0}~t_{\mathrm{2}})，对其进行积分，可得持续点热源的温度场分布：

T=\frac{q_t}{\mathrm{4}\mathrm{\pi }kR}\left[\mathrm{1}-\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{R}{\sqrt[]{\mathrm{4}\alpha t}}\right)\right]（26）

当t趋于无穷时，\mathrm{e}\mathrm{r}\mathrm{f}(\mathrm{0})=\mathrm{0}，此时温度为稳定状态.

T=\frac{q_t}{\mathrm{4}\mathrm{\pi }kR}（27）

4.3 连续作用移动点热源温度场

连续作用移动点热源在发热的同时仍移动一定的距离，R随时间发生改变，假设点热源在运动过程中x、y和z方向的点热源向目标点移动的速度分别为v_{\mathrm{1}}、v_{\mathrm{2}}和v_{\mathrm{3}}，则经过时间t后的温度为：

T={\displaystyle \underset{\mathrm{0}}{\overset{t_{\mathrm{2}}}{\int }}}\frac{q_{\mathrm{t}}}{c\rho {\left[\mathrm{4}\mathrm{\pi }\alpha (t_{\mathrm{2}}-t)\right]}^{\mathrm{3}/\mathrm{2}}}{\mathrm{e}}^{-\frac{(x-U-v_{\mathrm{1}}{t)}^{\mathrm{2}}+(y-V-v_{\mathrm{2}}{t)}^{\mathrm{2}}+(z-W-v_{\mathrm{3}}{t)}^{\mathrm{2}}}{\mathrm{4}\alpha (t_{\mathrm{2}}-t)}}\mathrm{d}t（28）

4.4 连续作用移动面热源温度场

假设面热源在xOy面上并在此平面移动，将整个面热源离散成多个微元面积为dxdy的点热源，微元面积上的平均热流密度为q_{\mathrm{w}\mathrm{0}}，则单位时间内微元面积产生的热量q_t=q_{\mathrm{w}\mathrm{0}}\mathrm{d}x\mathrm{d}y，在微元时间dt后产生的热量\mathrm{d}Q=q_{\mathrm{w}\mathrm{0}}\mathrm{d}x\mathrm{d}y\mathrm{d}t，连续作用移动面热源对物体内的某一个点引起的温度增加量为：

T={\displaystyle \underset{\mathrm{0}}{\overset{t_{\mathrm{2}}}{\int }}}{\displaystyle \underset{x_{\mathrm{1}}}{\overset{x_{\mathrm{2}}}{\int }}}{\displaystyle \underset{y_{\mathrm{1}}}{\overset{y_{\mathrm{2}}}{\int }}}\frac{q_{\mathrm{w}\mathrm{0}}}{c\rho {\left[\mathrm{4}\mathrm{\pi }\alpha (t_{\mathrm{2}}-t)\right]}^{\mathrm{3}/\mathrm{2}}}{\mathrm{e}}^{-\frac{(x-U-v_{\mathrm{1}}{t)}^{\mathrm{2}}+(y-V-v_{\mathrm{2}}{t)}^{\mathrm{2}}+W^{\mathrm{2}}}{\mathrm{4}\alpha (t_{\mathrm{2}}-t)}}\mathrm{d}x\mathrm{d}y\mathrm{d}t（29）

连续点热源温度场，定积分的原函数为误差函数，不可积，因此将三重积分进行离散化求解.

4.5 连续作用移动面热源温度场的离散化求解

4.5.1 面热源离散

将椭圆热源分成无数个具有相同大小的矩形，每一个矩形微元面热源看成是连续作用移动点热源.

已知x为(-b,b)，将x的总长度分成n个微元长度h_{\mathrm{1}}，h_{\mathrm{1}}=\mathrm{2}b/n. 椭圆热源的最大长度为a，a为(\mathrm{0},a)，将a分成的m个微元长度h_{\mathrm{2}}，h_{\mathrm{2}}=a/m. 由前述磨削区域分析可知，r_{\mathrm{1}}(x)在(-b,b)范围内取值，每一个x处的r_{\mathrm{1}}(x)大小不同，因此每一个r_{\mathrm{1}}(x)/h_{\mathrm{2}}不一定全部为整数，为了保证每一个微元面积相同，将r_{\mathrm{1}}(x)/h_{\mathrm{2}}进行四舍五入取整，每一个微元面积s=h_{\mathrm{1}}{h}_{\mathrm{2}}=\mathrm{2}ab/(nm). 面积离散示意图如图11所示.

4.5.2 时间离散

假设初始位置到运动结束位置的距离为L，面热源向目标点移动的整体速度为v，则面热源移动的总时间为t=L/v，t的取值范围为(\mathrm{0},L/v)，将时间分成g个微元长度t_{\mathrm{0}}，t_{\mathrm{0}}=L/(vg).

4.5.3 微元面积上的热量计算

每一个微元面积s在第k(\mathrm{0}\le k\le g)个位置停留时间t_{\mathrm{0}}所发出的热量为\mathrm{d}Q=q_{\mathrm{w}\mathrm{0}}\mathrm{d}x\mathrm{d}y\mathrm{d}t.

4.5.4 微元面积热通量计算

在微元面积s上，每一点的热通量不同，因此q_{\mathrm{w}\mathrm{0}}s需要进一步进行计算. 已知：

q(y)=\frac{\mathrm{2}{q}_{\mathrm{w}\mathrm{0}}}{r_{\mathrm{1}}\left(x\right)}y，\mathrm{0}\le y\le r_{\mathrm{1}}(x)（30）

r_{\mathrm{1}}(x)=\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}\right)}^{\mathrm{2}}}（31）

如前述热量分布规律可知，r_{\mathrm{1}}(x)在(-b,b)范围内取值，每一个x处的r_{\mathrm{1}}(x)斜率相同，不同r_{\mathrm{1}}(x)上的热流密度都可以用上述热流密度公式进行计算，除了y轴穿过的微元，在同一个微元s内，若x_i在上半个椭圆，则r_{\mathrm{1}}(x_i)\ge r_{\mathrm{1}}(x_{i+\mathrm{1}/\mathrm{2}})\ge r_{\mathrm{1}}(x_{i+\mathrm{1}})，r_{\mathrm{1}}(x)随着x的增加而减小；若x_i在下半个椭圆，则r_{\mathrm{1}}(x_i)\le r_{\mathrm{1}}(x_{i+\mathrm{1}/\mathrm{2}})\le r_{\mathrm{1}}(x_{i+\mathrm{1}})，r_{\mathrm{1}}(x)随着x的增加而减小. 在同一y值上，r_{\mathrm{1}}(x_i)到r_{\mathrm{1}}(x_{i+\mathrm{1}})的线热流密度的计算为：

q_{\mathrm{l}\mathrm{x}}={\displaystyle \underset{x_i}{\overset{x_{i+\mathrm{1}}}{\int }}}\frac{\mathrm{2}{q}_{\mathrm{w}\mathrm{0}}y}{\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}\right)}^{\mathrm{2}}}}\mathrm{d}x（32）

通过MATLAB对式（32）积分可知，被积分函数不可积.已知\mathrm{2}{q}_{\mathrm{w}\mathrm{0}}y为常数，假设R_{\mathrm{1}}=500 mm，R_{\mathrm{2}}= 300 mm，h_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.1 mm，则被积分函数为\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}{R}_{\mathrm{1}}\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}\right)-{\left(h_{\mathrm{m}\mathrm{a}\mathrm{x}}-R_{\mathrm{2}}+\sqrt[]{R_{\mathrm{2}}^{\mathrm{2}}-x^{\mathrm{2}}}\right)}^{\mathrm{2}}}}.被积分函数图如图12所示.

由图12可知，为了简化计算，应用梯形公式进行计算. 由计算公式可知，在同一y值上，x_i~x_{i+\mathrm{1}}的线平均热流密度可等效为该线的中点热流密度.

q_{\mathrm{w}\mathrm{0}}s=q(x_{i+\mathrm{1}/\mathrm{2}},y_{j+\mathrm{1}/\mathrm{2}})h_{\mathrm{1}}{h}_{\mathrm{2}}（33）

4.5.5 温度计算

由于在钢轨磨削过程中面热源在钢轨表面上，钢轨可以等效为半无限大的物体，通过绝热边界的镜像热源原理可得物体表面的温度是无限大物体表面温度的2倍：

T_{\mathrm{S}}=\mathrm{2}{\displaystyle \sum _{i=\mathrm{0}}^{i=n-\mathrm{1}}}{\displaystyle \sum _{j=\mathrm{0}}^{j=m-\mathrm{1}}}{\displaystyle \sum _{k=\mathrm{1}}^{k=g}}\frac{q(x_{i+\mathrm{1}/\mathrm{2}},y_{j+\mathrm{1}/\mathrm{2}})h_{\mathrm{1}}{h}_{\mathrm{2}}{t}_{\mathrm{0}}}{c\rho {\left\{\mathrm{4}\mathrm{\pi }\alpha \left[L/v-(g-k)t_{\mathrm{0}}\right]\right\}}^{\mathrm{3}/\mathrm{2}}}{\mathrm{e}}^{-\frac{(x-U-v_{\mathrm{1}}k{t}_{\mathrm{0}}{)}^{\mathrm{2}}+(y-V-v_{\mathrm{2}}k{t}_{\mathrm{0}}{)}^{\mathrm{2}}+W^{\mathrm{2}}}{\mathrm{4}\alpha \left[L/v-(g-k)t_{\mathrm{0}}\right]}}（34）

4.5.6 磨削过程温度计算

通过温度计算公式，可得目标点P在接触区域面热源移动过程中的温度随时间的变化，如图13所示.由图13可知，面热源在逐渐靠近目标点的过程中，温度先上升后下降，当t=1.382 4 s时，温度开始上升，当t=1.513 6 s时，温度达到最高，为412.46 ℃.

5 仿真分析

基于ABAQUS的用户自定义子程序DFLUX与Visual Studio进行联合仿真，根据上述热量分布规律，应用三角形热源模型在钢轨R300 mm的廓形上建立钢轨磨削的有限元温度仿真模型，应用瞬态热传递求解器对温度进行求解.以60 kg/m钢轨为例建立钢轨模型，取钢轨长度为0.05 m，环境温度为25 ℃. 有限元仿真示意图如图14所示.

通过JmatRro对材料的温变性能进行分析，可得仿真过程中钢轨材料物理性能参数，如表1所示.

通过创建输出数据文件（Output Database File， ODB）对仿真中的中心结点温度进行输出，将仿真分析与理论计算值进行对比. 温度仿真分析与理论计算对比曲线如图15所示.

由图15可知，面热源在走近、重合和走过目标点的过程中，目标点在仿真分析和理论计算过程中的温度都是先上升后下降，且几乎在同一时间点获得最高温度，由仿真可得其最高温度为437.806 °C，相对误差为6.14%. 结合砂轮和砂带的磨削特点，Zhou等^［12］、Lin等^［19］的研究结果证明了上述理论的正确性.

6 结 论

结合钢轨与砂带的接触几何关系建立的接触区域磨削深度和轮廓计算模型，阐明了磨削工艺参数对接触区参数的影响规律：随着砂带磨削半径的增加，不同钢轨廓形处接触区域面积呈对数增加. 当砂带磨削半径大于钢轨廓形半径时，椭圆接触区域轴长在yOz截面上；当砂带磨削半径等于钢轨廓形半径时，钢轨磨削接触区域为半圆；当砂带磨削半径小于钢轨廓形半径时，椭圆接触区域轴长在xOz截面上.

结合钢轨砂带磨削过程特性及磨削过程中的热量分布规律，得到了接触区域热流密度计算模型. 基于瞬时点热源温度场、连续作用点热源温度场、连续作用移动点热源温度场对连续作用移动面热源温度场进行了离散化求解，通过仿真对理论结果进行了验证. 结果得出：在给定磨削工艺参数下，仿真和理论的温度变化趋势相似且几乎在同一时间达到温度最大值，仿真和理论最高温度的相对误差为6.14%，通过已有研究验证了上述理论的正确性.

参考文献

[1] 刘启跃， 张波， 周仲荣 ．铁路钢轨损伤机理研究［J］．中国机械工程， 2002， 13（18）： 1596-1599．

[2] LIU Q Y， ZHANG B， ZHOU Z R ．The research of wear characteristics on steel rail［J］．China Mechanical Engineering，2002， 13（18）： 1596-1599．（in Chinese）

[3] 唐文龙 ．断轨原因分析及其预防措施［J］．铁道建筑，2001， 41（12）：31-32．

[4] TANG W L ．Cause analysis of broken track and its preventive measures［J］．Railway Engineering，2001， 41（12）： 31-32．（in Chinese）

[5] 周清跃， 田常海， 张银花，等 ．高速铁路钢轨打磨关键技术研究［J］．中国铁道科学， 2012， 33（2）： 66-70．

[6] ZHOU Q Y， TIAN C H， ZHANG Y H， et al ．Research on key rail grinding technology of high-speed railway［J］．China Railway Science， 2012， 33（2）： 66-70．（in Chinese）

[7] 王文健， 陈明韬， 郭俊， 等. 高速铁路钢轨打磨技术及其应 用［J］. 西南交通大学学报， 2007， 42（5）： 574-577．

[8] WANG W J， CHEN M T， GUO J， et al ．Rail grinding technique and its application in high-speed railway［J］. Journal of Southwest Jiaotong University， 2007， 42（5）： 574-577．（in Chinese）

[9] 贺振中 ．国外钢轨打磨技术的应用与思考［J］．中国铁路，2000（10）： 38-40．

[10] HE Z Z ．Application and thinking of rail grinding technology abroad［J］．Chinese Railways，2000（10）：38-40．（in Chinese）

[11] 欧阳明 ．高速铁路钢轨打磨及轮轨接触应力分析［D］．成都：西南交通大学，2015：10-15．

[12] OUYANG M. Rail grinding and wheel-rail contact stress analysis of high-speed railway［D］. Chengdu：Southwest Jiaotong University， 2015： 10-15．（in Chinese）

[13] 郭战伟 ．普速铁路钢轨打磨对轮轨接触关系的影响［J］．中国铁道科学， 2020， 41（6）： 109-116．

[14] GUO Z W ．Influence of rail grinding on wheel-rail contact relationship in general speed railway［J］．China Railway Science，2020， 41（6）： 109-116．（in Chinese）

[15] SEO J，GOO B， CHOI J，et al ．Effects of metal removal and residual stress on the contact fatigue life of railway wheels［J］．International Journal of Fatigue，2008，30（10/11）：2021-2029．

[16] 张磊 ．单程平面磨削淬硬技术的理论分析和试验研究［D］．济南： 山东大学， 2006： 28-109．

[17] ZHANG L ．Theoretical analysis and experimental study on hardening technology of single-pass plane grinding［D］．Jinan：Shandong University， 2006： 28-109．（in Chinese）

[18] ZHANG L C， MAHDI M ．Applied mechanics in grinding-IV．the mechanism of grinding induced phase transformation［J］．International Journal of Machine Tools and Manufacture， 1995， 35（10）： 1397-1409．

[19] 聂蒙， 李建勇， 沈海阔， 等 ．钢轨打磨作业过程中磨削温度建模与仿真［J］．铁道学报， 2013， 35（10）： 89-93．

[20] NIE M， LI J Y， SHEN H K，et al ．Modeling and simulation of temperature in process of rail grinding operation［J］．Journal of the China Railway Society， 2013， 35（10）： 89-93．（in Chinese）

[21] ZHOU K， DING H H， STEENBERGEN M， et al ．Temperature field and material response as a function of rail grinding parameters［J］．International Journal of Heat and Mass Transfer，2021， 175： 121366．

[22] 张青， 崔晓璐， 陈光雄， 等 ．钢轨打磨热力耦合分析［J］．铁道学报， 2015， 37（5）： 78-82．

[23] ZHANG Q， CUI X L， CHEN G X， et al ．Thermal-mechanical coupling analysis in rail grinding［J］. Journal of the China Railway Society， 2015， 37（5）： 78-82．（in Chinese）

[24] 张子舆 ．钢轨打磨过程中钢轨温度场及热机应力研究［D］．成都： 西南交通大学， 2017： 16-65．

[25] ZHANG Z Y. Study on temperature field and thermal stress of rail during rail grinding［D］. Chengdu： Southwest Jiaotong University， 2017： 16-65．（in Chinese）

[26] 杨天勇 ．钢轨砂带磨削的残余应力及其影响因素研究［D］．北京：北京交通大学，2019：35-65．

[27] YANG T Y. Study on residual stress and its influencing factors in rail abrasive belt grinding［D］. Beijing： Beijing Jiaotong University，2019：35-65．（in Chinese）

[28] MARINESCU I D， ROWE W B， DIMITROV B，et al ．Tribosystems of abrasive machining processes［M］//Tribology of Abrasive Machining Processes．Amsterdam：Elsevier，2004：13-40．

[29] ROWE W B， BLACK S C E， MILLS B，et al ．Grinding temperatures and energy partitioning［J］．Proceedings of the Royal Society of London Series A： Mathematical， Physical and Engineering Sciences， 1997， 453（1960）： 1083-1104．

[30] 蔡在亶 ．金属切削原理［M］．上海： 同济大学出版社， 1994：67-80．

[31] CAI Z D ．Metal cutting principle［M］．Shanghai：Tongji University Press， 1994： 67-80．（in Chinese）

[32] LIN B， ZHOU K， GUO J， et al ．Influence of grinding parameters on surface temperature and burn behaviors of grinding rail［J］．Tribology International， 2018， 122： 151-162．

基金资助

国家自然科学基金资助项目（52275399）