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

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

S2、根据车轮和钢轨真实材料属性，获取其单调拉伸应力‑应变曲线及其弹性模量E车轮和E钢轨、弹性泊松比ve,车轮和ve,钢轨、塑性泊松比vp,车轮和vp,钢轨、屈服强度σs,车轮和σs,钢轨、抗拉强度σb,车轮和σb,钢轨基本材料属性；

S3、根据步骤S1建立的车轮和钢轨的曲面函数f车轮(x,y,z)和f钢轨(x,y,z)对接触区Ig内n n

进行离散，车轮和钢轨接触区均在OXZ平面内离散为2×2的均匀矩形网格，n为正整数，即离散后车轮和钢轨上下对应节点具有相同的X、Z坐标，仅Y方向坐标不同；

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

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

Kij＝∫∫K(xi‑x',zi‑z')dx'dz',(i,j)∈Ig(公式2)S6、在接触区Ig内给出满足公式3的任意初始接触压力分布pij：其中，P0为所考虑的真实法向外载荷，Δx为当前离散后在X方向的节点距离，Δz为当前离散后在Z方向的节点距离；

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

S8、根据步骤S4得到的初始轮轨间隙函数h(x,y,z)得到初始轮轨间隙矩阵hij，进而结合步骤S7计算的挠度矩阵uij，得到当前情况下每个节点的间隙组成间隙矩阵gij：gij＝uij+hij,(i,j)∈Ig (公式5)S9、根据步骤S8，进而可以得到接触斑内所有节点的平均间隙：其中，Ic为接触压力pij非负的接触斑内区域，Ic∈Ig；Nc为接触斑内接触压力pij非负的节点个数；

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

S11、计算梯度折减矩阵tij：old

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

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

Gold＝G(公式10)

S13、根据步骤S12，通过公式12对修正矩阵rij进行更新：其中， 为在接触斑内修正矩阵rij的平均值；

*

S14、更新后的修正矩阵rij用于在梯度折减矩阵tij方向上的折减系数τ：S15、对当前步骤下的接触压力进行更新：*

pij＝pij‑τtij,(i,j)∈Ic(公式14)S16、计算下一步即将产生接触压力的集合Iol：*

Iol＝{(i,j)∈Ig:pij＝0,gij＜0}(公式15)S17、当 时，则下一步计算过程中δ＝1，否则δ＝0，并且，通过公式16对Iol内的节点压力进行更新：

# * *

pij＝pij‑τgij ,(i,j)∈Iol(公式16)S18、通过接触压力合力与法向外载荷之比对总体接触压力沿折减方向进行控制：new

其中，P1为当前步骤下的接触压力合力，pij 为修正后的接触压力分布；

S19、迭代是否结束的判断：根据给定的精度值ε0，若ε＜ε0，则得到最终的接触压力分布情况；若ε≥ε0则将步骤S17更新后的接触压力分布重新传入步骤S7，进行循环计算；其中，ε根据公式19计算：S20、由车轮和钢轨材料的单调拉伸应力‑应变曲线，取相同的应力水平进行分段，可以得到：

其中，σ1＝0，σ2＝min(σs,车轮,σs,钢轨)；σn＝min(σb,车轮,σb,钢轨)；n为总的应力个数； 和分别为每个应力对应的车轮、钢轨材料的应变；

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

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

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

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

S23、计算修正系数：

βr＝λνλr,Hλr,σ,r＝1～n‑2(公式23)其中，βi为修正系数；λν、λr,H和λr,σ可分别表示为：其中，Ck＝2.24；

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

Ip＝{(i,j)∈Ic:pij ≥Ckσ2}(公式27)S24、寻找弹塑性接触区Ip：new

Ip＝{(i,j)∈Ic:pij ≥Ckσ2}(公式27)new

S25、对弹塑性接触区Ip内的所有接触压力进行修正，如果满足Ckσt≤pij ＜Ckσt,t＝2～n‑2，则对该接触压力进行如下修正：

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

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