CN114925498B - 考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法 - Google Patents

Abstract

本发明涉及轮轨法向载荷计算领域，具体为考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法。其包括如下步骤：分别建立车轮和钢轨接触区附近车轮和钢轨三维曲面廓形函数；计算初始接触点位置及初始轮轨间隙函数；计算每个节点处的挠度组成挠度矩阵；根据初始轮轨间隙函数得到初始轮轨间隙矩阵；计算影响系数矩阵和梯度折减矩阵卷积；计算下一步即将产生接触压力的集合；通过接触压力合力与法向外载荷之比对总体接触压力沿折减方向进行控制；计算车轮、钢轨材料的硬化模量；计算轮轨材料的等效弹性模量和等效硬化模量。本发明相比现场直接检测或实验研究而言，成本小；相比现有接触理论而言，更加接近轮轨接触的真实情形。

Description

考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法

技术领域

本发明涉及轮轨法向载荷计算领域，特别是涉及考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法。

背景技术

随着中国铁路技术的迅速提高，其铁路运营规模、车辆轴重、速度、运行频次不断增加，形成了交错密集的铁路网络。轮轨系统则是保证列车安全行进的关键，其接触法向接触载荷的计算是进一步分析轮轨服役性能的重要前提。得益于经典Hertz接触理论的提出，轮轨接触理论在此基础上有了很大发展，建立了一系列的弹性接触算法。但是，这些方法均在弹性假设下提出，并没有考虑材料的塑性变形带来的影响。同时，轮轨接触过程中极易发生多段曲率同时接触的情况，这将为得到更为精确的法向接触载荷分布带来更大的挑战。

目前，在相关理论的支持下，结合有限元计算进而考虑材料塑性属性进行计算是较为有效的方法。但是该方法涉及高度几何非线性、材料非线性、接触非线性等问题，其计算耗时长、收敛性差，不利于进一步的分析和计算。

发明内容

本发明目的是针对背景技术中存在的问题，提出考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法，相比现场直接检测或实验研究而言，成本小；相比建立车轮和钢轨接触模型进行计算而言，耗时少；相比现有接触理论而言，更加接近轮轨接触的真实情形。

本发明的技术方案，考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法，包括如下步骤：

S1、根据车轮和钢轨真实几何形貌，分别建立接触区附近车轮和钢轨三维曲面廓形函数f_车轮(x,y,z)和f_钢轨(x,y,z)，其中X轴为钢轨横向，Y轴垂直向上，Z轴沿列车行进方向，坐标原点为车轮标准轮径检测点；

S2、根据车轮和钢轨真实材料属性，获取其单调拉伸应力-应变曲线及其弹性模量E_车轮和E_钢轨、弹性泊松比v_e,车轮和v_e,钢轨、塑性泊松比v_p,车轮和v_p,钢轨、屈服强度σ_s,车轮和σ_s,钢轨、抗拉强度σ_b,车轮和σ_b,钢轨等基本材料属性；

S3、根据步骤S1建立的车轮和钢轨的曲面函数f_车轮(x,y,z)和f_钢轨(x,y,z)对接触区I_g内进行离散，车轮和钢轨接触区均在OXZ平面内离散为2ⁿ×2ⁿ的均匀矩形网格，n为正整数，即离散后车轮和钢轨上下对应节点具有相同的X、Z坐标，仅Y方向坐标不同；

S4、考虑在横移量为X₀情况下的轮轨真实接触情形，根据步骤S1建立的车轮和钢轨的曲面函数f_车轮(x,y,z)和f_钢轨(x,y,z)计算初始接触点位置M₀(x₀,y₀,z₀)及初始轮轨间隙函数h(x,y,z)；

S5、根据公式1计算车轮和钢轨每个坐标处的影响系数K(x,z)，进而根据公式2计算进一步得到影响系数矩阵K_ij：

K_ij＝∫∫K(x_i-x',z_i-z')dx'dz',(i,j)∈I_g  (公式2)

S6、在接触区I_g内给出满足公式3的任意初始接触压力分布p_ij：

其中，P₀为所考虑的真实法向外载荷，Δx为当前离散后在X方向的节点距离，Δz为当前离散后在Z方向的节点距离；

S7、根据步骤S5得到的影响系数矩阵K_ij和公式4的卷积计算过程计算每个节点处的挠度组成挠度矩阵u_ij：

S8、根据步骤S4得到的初始轮轨间隙函数h(x,y,z)得到初始轮轨间隙矩阵h_ij，进而结合步骤S7计算的挠度矩阵u_ij，得到当前情况下每个节点的间隙组成间隙矩阵g_ij：

g_ij＝u_ij+h_ij,(i,j)∈I_g  (公式5)

S9、根据步骤S8，进而可以得到接触斑内所有节点的平均间隙：

其中，I_c为接触压力p_ij非负的接触斑内区域，I_c∈I_g；N_c为接触斑内接触压力p_ij非负的节点个数；

S10、根据步骤S9得到的平均间隙

通过公式7对间隙矩阵进行更新，进而根据公式8求得此时法向间隙的平方和G：

S11、计算梯度折减矩阵t_ij：

其中，在初始条件下，t_ij ^old＝0，G_old＝1，δ＝0，其值在迭代过程中会发生变化；

S12、根据公式10将G_old的值进行更新，并根据公式11计算影响系数矩阵K_ij和梯度折减矩阵t_ij卷积：

G_old＝G  (公式10)

S13、根据步骤S12，通过公式12对修正矩阵r_ij进行更新：

其中，

为在接触斑内修正矩阵r_ij的平均值；

S14、更新后的修正矩阵r_ij ^*用于在梯度折减矩阵t_ij方向上的折减系数τ：

S15、对当前步骤下的接触压力进行更新：

p_ij ^*＝p_ij-τt_ij,(i,j)∈I_c  (公式14)

S16、计算下一步即将产生接触压力的集合I_ol：

I_ol＝{(i,j)∈I_g:p_ij＝0,g_ij ^*＜0}  (公式15)

S17、当

时，则下一步计算过程中δ＝1，否则δ＝0，并且，通过公式16对I_ol内的节点压力进行更新：

p_ij ^#＝p_ij ^*-τg_ij ^*,(i,j)∈I_ol  (公式16)

S18、通过接触压力合力与法向外载荷之比对总体接触压力沿折减方向进行控制：

其中，P₁为当前步骤下的接触压力合力，p_ij ^new为修正后的接触压力分布；

S19、迭代是否结束的判断：根据给定的精度值ε₀，若ε＜ε₀，则得到最终的接触压力分布情况；若ε≥ε₀则将步骤S17更新后的接触压力分布重新传入步骤S7，进行循环计算；其中，ε根据公式19计算：

S20、由车轮和钢轨材料的单调拉伸应力-应变曲线，取相同的应力水平进行分段，可以得到：

其中，σ₁＝0，σ₂＝min(σ_s,车轮,σ_s,钢轨)；σ_n＝min(σ_b,车轮,σ_b,钢轨)；n为总的应力个数；

和

分别为每个应力对应的车轮、钢轨材料的应变；

S21、根据步骤S19得到的车轮、钢轨材料应力应变数据点，计算车轮、钢轨材料的硬化模量：

其中，H_r,车轮和H_r,钢轨分别为车轮、钢轨的硬化模量点对；

S22、根据公式22计算轮轨材料的等效弹性模量和等效硬化模量：

其中，G_eq和

分别为等效弹性模量和等效硬化模量点对；v_p,车轮和v_p,钢轨分别为车轮、钢轨材料的塑性条件下的泊松比，一般取0.5；

S23、计算修正系数：

β_r＝λ_νλ_r,Hλ_r,σ,r＝1～n-2 (公式23)

其中，β_i为修正系数；λ_ν、λ_r,H和λ_r,σ可分别表示为：

其中，C_k＝2.24；

步骤24、寻找弹塑性接触区I_p：

I_p＝{(i,j)∈I_c:p_ij ^new≥C_kσ₂}(公式27)

S24、寻找弹塑性接触区I_p：

I_p＝{(i,j)∈I_c:p_ij ^new≥C_kσ₂}(公式27)

S25、对弹塑性接触区I_p内的所有接触压力进行修正，如果满足C_kσ_t≤p_ij ^new＜C_kσ_t,t＝2～n-2，则对该接触压力进行如下修正：

优选的，还包括：

S26、计算每个网格放大系数ξ：

S27、根据步骤S26计算的网格放大系数，对所有网格尺寸和节点坐标进行放缩：

与现有技术相比，本发明具有如下有益的技术效果：

本发明基于真实三维轮轨接触几何廓形、接触载荷及轮轨弹性属性，采用共轭梯度法及快速傅里叶变换获取了轮轨弹性接触压力的精确分布，进而根据真实轮轨材料应力-应变响应曲线对精确弹性接触压力分布进行了塑性修正，得到考虑曲率过渡下三维弹塑性轮轨法向载荷分布情况。因此，本发明相比现场直接检测或实验研究而言，成本小；而相比建立车轮和钢轨接触模型进行计算而言，耗时少；而相比现有接触理论而言，本发明进一步考虑了材料的塑性行为及轮轨接触区的曲率过渡情况，更加接近轮轨接触的真实情形。

附图说明

图1为轮轨廓形初始位置图；

图2为曲面离散示意图；

图3为弹性法向接触压力示意图；

图4为处理后的车轮、钢轨材料应力应变曲线示意图；

图5为当前接触条件下弹塑性接触压力示意图；

图6和图7为接触压力在X轴(曲率过渡方向)上的分布情况示意图。

具体实施方式

根据上面的考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法，首先，根据步骤S1，选取LMa型标准踏面廓形车轮，直径为840mm，选取60kg/m标准廓形钢轨，二维轮廓如图1所示。进而确定了二者的三维曲面廓形函数f_车轮(x,y,z)和f_钢轨(x,y,z)。

根据步骤S2，确定了车轮和钢轨真实材料属性，见表1。

表1：车轮和钢轨材料基本属性

<![CDATA[E<sub>车轮</sub>＝203.125GPa，v<sub>e,车轮</sub>＝0.3，v<sub>p,车轮</sub>＝0.5，σ<sub>s,车轮</sub>＝465.83MPa，σ<sub>b,车轮</sub>＝840.02MPa]]>
	<![CDATA[E<sub>钢轨</sub>＝202.535GPa，v<sub>e,钢轨</sub>＝0.3，v<sub>p,钢轨</sub>＝0.5，σ<sub>s,钢轨</sub>＝325MPa，σ<sub>b,钢轨</sub>＝1223.44MPa]]>

根据步骤S3，将轮轨接触区在OXZ平面内离散为2⁸×2⁸的均匀矩形网格，如图2所示。

根据步骤S4，取横移量X₀＝0，计算了初始接触点M₀(0.005369，-0.03000，0)和初始间隙函数h(x,y,z)。

根据步骤S5，计算了车轮和钢轨每个坐标处的影响系数K(x,z)，进而根据公式2计算进一步得到影响系数矩阵K_ij。

根据步骤S6，取法向接触合力P₀＝100kN，在接触区I_g内给出满足公式31的任意初始接触压力分布p_ij：

其中，根据步骤3确定的网格，Δx＝4.5313e-04，Δz＝5.8086e-04。

根据步骤S7，进行影响系数矩阵K_ij和公式4的卷积计算，计算每个节点处的挠度组成挠度矩阵u_ij。

根据步骤S8，根据初始轮轨间隙函数h(x,y,z)得到初始轮轨间隙矩阵h_ij，进而结合步骤7计算的挠度矩阵u_ij，得到当前情况下每个节点的间隙组成间隙矩阵g_ij(公式5)。

根据步骤S9，可以得到接触斑内所有节点的平均间隙。

根据步骤S10，计算平均间隙

对间隙矩阵进行更新(公式7)，进而可以根据公式8求得此时法向间隙的平方和G。

根据步骤S11，计算梯度折减矩阵t_ij。

根据步骤S12，将G_old的值进行更新(公式10)，并计算影响系数矩阵K_ij和梯度折减矩阵t_ij卷积(公式11)。

根据步骤S13，对修正矩阵r_ij进行更新(公式12)。

根据步骤S14，计算折减系数τ(公式13)。

根据步骤S15，对当前步骤下的接触压力进行更新(公式14)。

根据步骤S16，计算下一步即将产生接触压力的集合I_ol(公式15)。

根据步骤S17，当

时，则下一步计算过程中δ＝1，否则δ＝0。并且，对I_ol内的节点压力进行更新(公式16)。

根据步骤S18，通过接触压力合力与法向外载荷之比对总体接触压力沿折减方向进行控制(公式17)。

根据步骤S19，迭代是否结束的判断：根据给定的精度值ε₀，若ε＜ε₀，则得到最终的接触压力分布情况。若ε≥ε₀则将步骤17更新后的接触压力分布重新传入步骤S7，进行循环计算。最终可以得到弹性法向接触压力分布规律如图3所示。其中，ε可根据公式19计算。

根据步骤S20，由车轮和钢轨材料的单调拉伸应力-应变曲线，取相同的应力水平进行分段(如图4所示)，可以得到σ₁＝0，σ₂＝325MPa；σ_n＝840.02MPa。

根据步骤S21，根据步骤19得到的车轮、钢轨材料应力应变数据点，计算车轮、钢轨材料的硬化模量(公式21)。

根据步骤S22，根据公式22计算轮轨材料的等效弹性模量G_eq和等效硬化模量

根据步骤S23，取C_k＝2.24，计算修正系数β_i。

根据步骤S24，寻找弹塑性接触区I_p。

根据步骤S25，对弹塑性接触区I_p内的所有接触压力进行修正：如果满足C_kσ_t≤p_ij ^new＜C_kσ_t,t＝2～n-2，则按照公式28对该接触压力进行修正。

根据步骤S26，计算得到每个网格放大系数ξ＝1.0683。

根据步骤S27，根据步骤25计算的网格放大系数，对所有网格尺寸和节点坐标进行放缩。

表2：本方法计算精度

通过上述计算，最终得到当前接触条件下考虑材料塑性的接触压力分布情况，如图5所示。表2、图6和图7为本方法与有限元计算结果的对比，可以发现本方法与有限元计算结果基本吻合，即本方法的精度较高。

本发明基于真实三维轮轨接触几何廓形、接触载荷及轮轨弹性属性，采用共轭梯度法及快速傅里叶变换获取了轮轨弹性接触压力的精确分布，进而根据真实轮轨材料应力-应变响应曲线对精确弹性接触压力分布进行了塑性修正，得到考虑曲率过渡下三维弹塑性轮轨法向载荷分布情况。因此，本发明相比现场直接检测或实验研究而言，成本小；而相比建立车轮和钢轨接触模型进行计算而言，耗时少；而相比现有接触理论而言，本发明进一步考虑了材料的塑性行为及轮轨接触区的曲率过渡情况，更加接近轮轨接触的真实情形。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于此，在所属技术领域的技术人员所具备的知识范围内，在不脱离本发明宗旨的前提下还可以作出各种变化。

Claims (2)

1.考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法，其特征在于，包括如下步骤：

S1、根据车轮和钢轨真实几何形貌，分别建立接触区附近车轮和钢轨三维曲面廓形函数f_车轮(x,y,z)和f_钢轨(x,y,z)，其中X轴为钢轨横向，Y轴垂直向上，Z轴沿列车行进方向，坐标原点为车轮标准轮径检测点；

S2、根据车轮和钢轨真实材料属性，获取其单调拉伸应力-应变曲线及其弹性模量E_车轮和E_钢轨、弹性泊松比v_e,车轮和v_e,钢轨、塑性泊松比v_p,车轮和v_p,钢轨、屈服强度σ_s,车轮和σ_s,钢轨、抗拉强度σ_b,车轮和σ_b,钢轨基本材料属性；

S3、根据步骤S1建立的车轮和钢轨的曲面函数f_车轮(x,y,z)和f_钢轨(x,y,z)对接触区I_g内进行离散，车轮和钢轨接触区均在OXZ平面内离散为2ⁿ×2ⁿ的均匀矩形网格，n为正整数，即离散后车轮和钢轨上下对应节点具有相同的X、Z坐标，仅Y方向坐标不同；

S4、考虑在横移量为X₀情况下的轮轨真实接触情形，根据步骤S1建立的车轮和钢轨的曲面函数f_车轮(x,y,z)和f_钢轨(x,y,z)计算初始接触点位置M₀(x₀,y₀,z₀)及初始轮轨间隙函数h(x,y,z)；

S5、根据公式1计算车轮和钢轨每个坐标处的影响系数K(x,z)，进而根据公式2计算进一步得到影响系数矩阵K_ij：

K_ij＝∫∫K(x_i-x',z_i-z')dx'dz',(i,j)∈I_g(公式2)

S6、在接触区I_g内给出满足公式3的任意初始接触压力分布p_ij：

其中，P₀为所考虑的真实法向外载荷，Δx为当前离散后在X方向的节点距离，Δz为当前离散后在Z方向的节点距离；

S7、根据步骤S5得到的影响系数矩阵K_ij和公式4的卷积计算过程计算每个节点处的挠度组成挠度矩阵u_ij：

S8、根据步骤S4得到的初始轮轨间隙函数h(x,y,z)得到初始轮轨间隙矩阵h_ij，进而结合步骤S7计算的挠度矩阵u_ij，得到当前情况下每个节点的间隙组成间隙矩阵g_ij：

g_ij＝u_ij+h_ij,(i,j)∈I_g (公式5)

S9、根据步骤S8，进而可以得到接触斑内所有节点的平均间隙：

其中，I_c为接触压力p_ij非负的接触斑内区域，I_c∈I_g；N_c为接触斑内接触压力p_ij非负的节点个数；

S10、根据步骤S9得到的平均间隙

通过公式7对间隙矩阵进行更新，进而根据公式8求得此时法向间隙的平方和G：

S11、计算梯度折减矩阵t_ij：

其中，在初始条件下，t_ij ^old＝0，G_old＝1，δ＝0，其值在迭代过程中会发生变化；

S12、根据公式10将G_old的值进行更新，并根据公式11计算影响系数矩阵K_ij 和梯度折减矩阵t_ij卷积：

G_old＝G(公式10)

S13、根据步骤S12，通过公式12对修正矩阵r_ij进行更新：

其中，

为在接触斑内修正矩阵r_ij的平均值；

S14、更新后的修正矩阵r_ij ^*用于在梯度折减矩阵t_ij方向上的折减系数τ：

S15、对当前步骤下的接触压力进行更新：

p_ij ^*＝p_ij-τt_ij,(i,j)∈I_c(公式14)

S16、计算下一步即将产生接触压力的集合I_ol：

I_ol＝{(i,j)∈I_g:p_ij＝0,g_ij ^*＜0}(公式15)

S17、当

时，则下一步计算过程中δ＝1，否则δ＝0，并且，通过公式16对I_ol内的节点压力进行更新：

p_ij ^#＝p_ij ^*-τg_ij ^*,(i,j)∈I_ol(公式16)

S18、通过接触压力合力与法向外载荷之比对总体接触压力沿折减方向进行控制：

其中，P₁为当前步骤下的接触压力合力，p_ij ^new为修正后的接触压力分布；

S19、迭代是否结束的判断：根据给定的精度值ε₀，若ε＜ε₀，则得到最终的接触压力分布情况；若ε≥ε₀则将步骤S17更新后的接触压力分布重新传入步骤S7，进行循环计算；其中，ε根据公式19计算：

S20、由车轮和钢轨材料的单调拉伸应力-应变曲线，取相同的应力水平进行分段，可以得到：

其中，σ₁＝0，σ₂＝min(σ_s,车轮,σ_s,钢轨)；σ_n＝min(σ_b,车轮,σ_b,钢轨)；n为总的应力个数；

和

分别为每个应力对应的车轮、钢轨材料的应变；

S21、根据步骤S19得到的车轮、钢轨材料应力应变数据点，计算车轮、钢轨材料的硬化模量：

其中，H_r,车轮和H_r,钢轨分别为车轮、钢轨的硬化模量点对；

S22、根据公式22计算轮轨材料的等效弹性模量和等效硬化模量：

其中，G_eq和H_e ^r _q分别为等效弹性模量和等效硬化模量点对；v_p,车轮和v_p,钢轨分别为车轮、钢轨材料的塑性条件下的泊松比，一般取0.5；

S23、计算修正系数：

β_r＝λ_νλ_r,Hλ_r,σ,r＝1～n-2(公式23)

其中，β_i为修正系数；λ_ν、λ_r,H和λ_r,σ可分别表示为：

其中，C_k＝2.24；

步骤24、寻找弹塑性接触区I_p：

I_p＝{(i,j)∈I_c:p_ij ^new≥C_kσ₂}(公式27)

S24、寻找弹塑性接触区I_p：

I_p＝{(i,j)∈I_c:p_ij ^new≥C_kσ₂}(公式27)

S25、对弹塑性接触区I_p内的所有接触压力进行修正，如果满足C_kσ_t≤p_ij ^new＜C_kσ_t,t＝2～n-2，则对该接触压力进行如下修正：

2.根据权利要求1所述的考虑曲率过渡的三维弹塑性轮轨法向载荷快速计算方法，其特征在于，还包括：

S26、计算每个网格放大系数ξ：

S27、根据步骤S26计算的网格放大系数，对所有网格尺寸和节点坐标进行放缩：

Images