简介这份基于Python的非赫兹轮轨接触力学简化模型面向铁路接触力学研究人员、轨道工程及车辆工程相关技术人员专门应对高速、重载或非线性变形下传统赫兹理论难以准确预测的接触行为。模型围绕Piotrowski-Kik思路实现覆盖非线性变形、滑动接触、黏着力学、滚动接触疲劳等关键因素借助NumPy、SciPy等科学计算库完成数值求解与可视化代码结构清晰便于扩展和二次开发。压缩包内共13个文件以5个Python脚本为主辅以2个RST说明、2个Markdown文档以及轮轨型面数据、配置和许可证文件整体仅49KB。通过模型代码与配套文档用户可学习非赫兹轮轨接触建模思路掌握接触压力、磨损及疲劳寿命的预测方法为列车设计、运行优化和轨道维护提供参考。已有263人学习下载适合具备一定力学与Python基础的读者深入研究。1. 非赫兹轮轨接触Hertz假设失效的典型工况某地铁线路车轮踏面磨耗后在横移量2mm时轮轨接触斑从单一椭圆变成“哑铃”形Hertz椭圆近似算出的最大接触压力比线路实测高出22%。这不是个例。曲线通过、长期磨耗、型面不匹配场景下接触区主曲率不再是常数接触斑可为不规则形状甚至两个独立小斑。Hertz公式要求接触体表面连续、曲率恒定并忽略摩擦切向影响这些条件在实际轮轨匹配中很难同时满足。这个模型包是一份Python实现的轮轨非赫兹接触简化模型覆盖型面离散、法向压力迭代和FASTSIM风格切向求解目标是把一次接触计算控制在数十毫秒量级同时保留压力峰和蠕滑力的工程精度。适合轨道车辆动力学方向的研究生、线路轮轨匹配工程师也是数值接触力学的入门样例。解压zip后直接运行main.py即可复现下文所有结果。2. 轮轨接触几何离散化型面坐标与曲率参数的Python表达2.1 型面数据读取与样条重采样轮轨型面数据通常来自CAD或激光轮对测量仪最常见的格式是两列数值相对基准点的横向坐标x和竖向坐标z。直接拿原始测点做接触求解会带来两个问题。一是测点密度不均匀在轨距角附近加密而在轨底区域稀疏二是离散点求二阶导数做曲率时噪声会被放大。所以一般先用样条插值把型面重采样到固定步距常见步距为0.05mm到0.1mm。这样做还能统一左侧和右侧型面的采样坐标系方便后续做左右轨同时求解。import numpy as np from scipy.interpolate import CubicSpline def load_and_resample(path, step0.05): # 输入文件两列x(mm) z(mm)#开头为注释 raw np.loadtxt(path, comments#) x, z raw[:, 0], raw[:, 1] # 去除重复横坐标否则CubicSpline会抛异常 keep np.r_[True, np.diff(x) 1e-6] x, z x[keep], z[keep] # 自然边界样条二阶导在端点平滑 cs CubicSpline(x, z, bc_typenatural) # 在原始范围上均匀重采样 x_new np.arange(x[0], x[-1], step) z_new cs(x_new) return x_new, z_new这段代码把散点型面变成均匀网格。参数step是重采样步距轮轨接触计算里常用0.010.1mm。步距太小会使后续接触网格节点数爆炸步距太大又丢失轨距角处的局部曲率突变。bc_typenatural表示两端二阶导为零对于型面这种开曲线比默认的not-a-knot更适合。若型面有垂直断点比如辙叉区CubicSpline可能出现过冲此时改用scipy.interpolate.Akima1DInterpolator更稳。2.2 接触点搜索与曲率计算给定轮对横向位置接触点由轮轨几何最小间隙决定。假设钢轨型面在全局坐标中固定轮对踏面相对钢轨有一个横向偏移dy那么最小间隙对应的横向坐标就是滚动圆附近的名义接触点。搜索范围要覆盖轮缘根部到轨距角通常从-15mm到15mm步长0.05mm。接触点处的曲率由型面样条的二阶导决定直接逐点计算即可。def find_contact_y(wheel_cs, rail_cs, dy, y_search): # y_search: 横向扫描网格mm # 轮对横移dy后踏面横向坐标比轨面横向坐标小dy dist rail_cs(y_search) - wheel_cs(y_search - dy) return y_search[np.argmin(dist)] def curv_radius(cs, y): d1 cs(y, 1) # 一阶导数 d2 cs(y, 2) # 二阶导数 k abs(d2) / (1 d1**2) ** 1.5 if k 1e-12: return 1e6 return 1.0 / kfind_contact_y把踏面曲线沿横移方向平移再与轨面曲线做垂直差最小值对应名义接触点。curv_radius返回曲率半径单位与坐标一致都是mm。当接触角小于8°时用竖直方向近似法向误差可以忽略如果接触角更大需要旋转型面建立局部接触坐标系代码量会增加不少。工程简化模型通常接受这个误差因为轮缘贴靠时主要关注的是接触斑位置而不是精确法向压力。下表是某磨耗踏面与60N轨在横向坐标上的曲率半径变化。这里的“相对接触点横向位置”以名义滚动圆为原点负值指向轮缘侧。相对接触点横向位置mm踏面横向曲率半径mm轨面横向曲率半径mm等效横向曲率半径mm-1.5320.513.0212.510.0385.213.0112.841.5410.880.4067.15等效曲率半径从12.5跳到67说明接触区几何与Hertz经典假设“曲率半径近似常数”矛盾这是非赫兹问题的最直接来源。如果只用名义接触点的曲率去套Hertz误差就会被放大到不可接受。2.3 接触候选网格与间隙矩阵在名义接触点附近建立一个矩形网格x方向为滚动方向y方向为横向。网格范围只要覆盖可能接触区域即可无需覆盖整个接触斑迭代中压力为零的单元自然退出。单元矩形的尺寸受计算量限制一般取0.20.5mm。间隙矩阵表示两个弹性体在无载荷状态下沿法线方向的距离。将轮轨型面在接触点处展开为二次曲面间隙可表达为g(y) 0.5 * (1 / R_eq_y) * y^2 0.5 * (1 / R_eq_x) * x^2其中R_eq_x由车轮名义滚动半径和轨面纵向半径合成R_eq_y由型面横向曲率合成。这个二次近似在接触斑尺寸远小于曲率半径时成立对非赫兹问题的接触斑宽度一般小于15mm而曲率半径最低也在12mm左右误差会稍大但作为简化模型的初始条件足够。def build_grid_and_gap(cx, cy, R_eq_x, R_eq_y, w, l, dx, dy): # cx, cy: 接触中心坐标w/l为半宽半长 xs np.arange(cx - l, cx l dx, dx) ys np.arange(cy - w, cy w dy, dy) XX, YY np.meshgrid(xs, ys) gap 0.5 * (1.0 / R_eq_x) * XX**2 0.5 * (1.0 / R_eq_y) * YY**2 return XX, YY, gap该函数的gap是初始刚性间隙单位是mm。注意R_eq_x和R_eq_y单位必须与坐标一致如果坐标是mm曲率半径也必须是mm否则计算出的gap会差三个数量级。我经常看到有人把半径单位混成米结果压力迭代出来是负数。网格步距看计算目标动态仿真里跑几百个接触点网格用0.5mm比较划算只研究单一静态接触时0.2mm可以把压力峰算得很准。完成几何离散后输出为四个量网格X/Y坐标、间隙矩阵gap、等效曲率R_eq_x/R_eq_y。它们作为第3章法向迭代的输入后续所有计算都在这套网格上进行。3. 法向接触非Hertz解影响系数矩阵与迭代求解3.1 离散影响系数矩阵对于弹性半空间法向表面位移由接触压力积分得到。对矩形单元上均布压力离散后每个单元j对单元i中心的位移影响系数为C_ij。严格做法是求矩形载荷的Love解析解但工程简化模型通常把远距离单元近似为集中力。这个模型源码采用点源近似并将自影响系数用等面积圆盘公式替换。这样既避免复杂的解析积分公式又保证对角线元素不因奇异性发散。def build_influence_matrix(xs, ys, dx, dy, E, nu): # xs, ys: 单元中心坐标(mm)E(MPa)nu泊松比 A dx * dy # 单元面积 (mm^2) a_eq np.sqrt(A / np.pi) # 等面积圆半径 n len(xs) coord_x xs.reshape(-1, 1) - xs.reshape(1, -1) coord_y ys.reshape(-1, 1) - ys.reshape(1, -1) r np.hypot(coord_x, coord_y) # 单元中心距离 # 点源近似u F * (1-nu^2) / (pi*E*r) C (1 - nu**2) / (np.pi * E) * A / np.where(r 1e-12, a_eq, r) # 自影响均布压力圆盘中心 u 2*(1-nu^2)*p*a_eq / E idx np.diag_indices(n) C[idx] 2.0 * (1 - nu**2) * a_eq / E return C逻辑非对角线元素把单元内总力视为作用在中心的集中力位移公式来自Boussinesq点载荷解。对角线元素由于不能把单元看成零面积点用等面积圆盘受均布压力时中心位移公式。等效替换后自影响与近邻单元的误差随网格细化而消失。参数E为弹性模量钢轨和轮对都取210000MPanu取0.3。注意单位统一坐标为mm压力用MPa时面积mm^2位移mm量级一致。如果压力用Pa而坐标用mm计算结果会差几个数量级。3.2 压力迭代与载荷平衡法向接触条件是一个线性互补问题接触区内“间隙 弹性变形 - 接近量”等于零非接触区压力为零。这个模型源码采用“接触集仿射求解法”。若接触集已知则在该集合上方程为C_sub p_sub δ * 1 - gap_sub因此p_sub是δ的仿射函数。先求解两次线性方程组得到两个基准解再根据总载荷约束∑p W计算δ紧接着更新接触集。该算法比逐次置换法收敛更快因为压力分布随δ变化是线性仿射只要接触集判别正确一次外迭代就能得到精确解。def solve_normal(C, gap, W, max_iter100): n len(gap) p np.zeros(n) contact gap 0.0 for it in range(max_iter): Ci C[np.ix_(contact, contact)] gi gap[contact] # 解 C_sub * v 1, C_sub * w g_sub v np.linalg.lstsq(Ci, np.ones(len(gi)), rcondNone)[0] w np.linalg.lstsq(Ci, gi, rcondNone)[0] # 载荷条件sum(p) W delta (W np.sum(w)) / max(np.sum(v), 1e-30) p_i delta * v - w p_new np.zeros(n) p_new[contact] p_i # 计算间隙 弹性变形 u C p_new h gap u - delta new_contact (h 0.0) | (p_new 0.0) if np.array_equal(new_contact, contact) and np.all(p_i -1e-12): break contact new_contact return p, delta, it逻辑第4行和第5行分别解出v和wv表示单位接近量下的压力分配w表示初始间隙引起的压力贡献。delta从总载荷条件解出这样任意迭代步压力总和都严格等于W。h是最终间隙小于等于零表示两个表面互相嵌入应加入接触p_new中为负的单元实际上应脱离。循环直到接触集和压力稳定。这里用lstsq代替solve因为接触矩阵在多个连通域时可能接近奇异最小二乘比直接求逆稳定得多。该算法每次迭代需要求解一次接触子矩阵子矩阵大小随接触面积变化通常在500×500到4000×4000之间。对于Python直接解在2000×2000时约1秒作为静态分析可接受若要嵌入动力学循环需要用Cholesky分解加预处理器加速。见第5章。3.3 压力云图与接触斑输出迭代结束后把一维压力reshape回二维并画出分布。也可以提取接触斑面积和最大压力用于与Hertz解比对。绘图时注意把X轴设为滚动方向Y轴设为横向这样能直观看到压力斑是否偏斜。import matplotlib.pyplot as plt def plot_contact(XX, YY, p_2d): plt.figure(figsize(6, 4)) cf plt.contourf(XX, YY, p_2d, levels25, cmapviridis) plt.colorbar(cf, labelContact pressure (MPa)) plt.xlabel(x (mm) 滚动方向) plt.ylabel(y (mm) 横向) plt.axis(equal) plt.show()参数说明XX和YY来自第2.3节的build_grid_and_gapp_2d形状必须和网格形状一致。建议把最大压力、接触斑面积、压力积分值打印出来自检。压力积分和输入W之间的相对误差应小于0.1%如果误差大说明迭代没收敛到正确接触集先检查接触集判定条件。下面是一个缩样本的收敛比较展示不同网格密度下法向求解的数值表现。网格尺寸mm单元数最大接触压力MPa迭代次数载荷相对误差%0.53721103870.090.310404109290.060.2228011120120.04可以看到压力峰对网格尺寸敏感0.5mm网格低估约7%所以研究最大接触应力时必须做网格收敛性验证。迭代次数受接触集变化影响接触斑越复杂比如多个不连通区域收敛越慢。4. 切向滚动接触简化模型FASTSIM思路在Python中的落地4.1 为什么用FASTSIM简化轮轨切向完全精确解需要Kalker的CONTACT程序该程序基于边界元求非Hertz切向解计算代价高且不适合多步动力学仿真。工程上常用FASTSIM它把接触斑按滚动方向分成条带用一维流动模型逼近切向应力分布。FASTSIM要求法向压力分布已知且对每个条带假设切向柔度系数固定。虽然单个单元上的误差存在但总蠕滑力误差一般控制在5%10%以内已满足绝大多数车辆动力学仿真。这里的代码是FASTSIM风格的离散版本突出递推和饱和两个核心逻辑。柔度系数的选取原始FASTSIM用接触椭圆半轴和材料参数组合出柔度系数L。为简化这个实现把切向应变增量直接映射到剪应力上用等效倒柔度乘以步长。不同文献对L的取值差20%不会影响最终蠕滑力饱和值但对蠕滑-力线性段有明显影响。所以在低蠕滑率工况下若结果偏差过大优先检查L的数值。4.2 纵向/横向蠕滑递推实现输入蠕滑率xi是纵向蠕滑率eta是横向蠕滑率phi是自旋蠕滑率。单元坐标xs沿滚动方向递增。从入口开始每一步按蠕滑引起的位移增量更新切向应力候选值若超过库伦摩擦极限μp则按比例拉回摩擦圆。这里的代码按一维条带递推每个单元保留当前累计应力。def fastsim_tangential(xs, ys, p, xi, eta, phi, mu, G, nu, dx): # xs, ys: 单元中心坐标(mm)p: 法向压力(MPa) n len(xs) tx np.zeros(n) ty np.zeros(n) # 简化倒柔度单位经过dx吸收让算法对步距不敏感 K_x G / (2.0 * (2.0 - nu)) * 0.01 K_y K_x * (2.0 - nu) / (2.0 * (1.0 - nu)) * 0.01 order np.argsort(xs) # 按滚动方向从入口到出口 qx 0.0 qy 0.0 prev_x None for i in order: if prev_x is not None: ds xs[i] - prev_x # 蠕滑累积位移纵向蠕滑自旋横向蠕滑-自旋 sx (xi - phi * ys[i]) * ds sy (eta phi * xs[i]) * ds qx K_x * sx qy K_y * sy prev_x xs[i] # 库伦摩擦极限 limit mu * p[i] q np.hypot(qx, qy) if q limit: scale limit / q qx * scale qy * scale tx[i] qx ty[i] qy return tx, ty代码逻辑先排序保证从接触斑入口向出口递推。每个单元先按蠕滑率累加剪应力再把总剪应力限制在摩擦圆内。这样的结果是压力大、蠕滑小的区域处于粘着态压力小或蠕滑大的区域被截断形成滑动区。参数dx没有直接参与计算但传入是为了提醒使用者如果修改坐标单位柔度系数中的步距换算也要同步调整。求和得到总蠕滑力单元面积是dx*dy单位为mm^2dx XX[0, 1] - XX[0, 0] dy YY[1, 0] - YY[0, 0] Fx np.sum(tx) * dx * dy * 0.001 # MPa*mm^2 mN转成N Fy np.sum(ty) * dx * dy * 0.001tx和ty的单位是MPa乘面积得到mN0.001再转成N。结果可直接用于车辆动力学求解。4.3 饱和效应与误差对照切向应力一旦达到库伦极限就按比例缩减等价于摩擦圆约束。实际接触斑内压力分布不对称饱和边界也不是一条直线可能在压力谷值处先滑动。这是简化模型最有价值的地方它能反映非赫兹压力斑对切向饱和区的直接影响。下面用一个单接触斑算例对照数据来自同参数下的简化模型与CONTACT结果网格尺寸0.3mm。纵向蠕滑率‰简化模型FxkNCONTACT FxkNFx偏差%17.98.46.0512.112.64.01513.213.41.5小蠕滑率时误差偏大因为柔度系数近似导致线性段刚度偏低大蠕滑率时大部分区域饱和简化模型由摩擦极限主导偏差明显变小。所以FASTSIM风格更适合蠕滑率较大的牵引和制动工况而不是直线稳态蠕滑极小的问题。5. 参数敏感性分析与加速技巧5.1 网格收敛性门槛第3章压力云图显示最大压力对网格尺寸高度敏感。做法是取同一工况从0.8mm依次加密到0.15mm观察最大压力变化率。当两次相邻尺寸的压力峰值变化小于1%认为该网格合格。若在不同横移量下做收敛曲线所需网格也不同大横移时接触斑更扁平需要更细网格。def convergence_check(W, sizes[0.8, 0.5, 0.3, 0.2, 0.15]): p_peaks [] for dx in sizes: # 复用2.3节网格构建步骤 XX, YY, gap generate_grid(dx) # 需要自己实现 C build_influence_matrix(XX.flatten(), YY.flatten(), dx, dx, 210000, 0.3) p, delta, it solve_normal(C, gap.flatten(), W) p_peaks.append(p.max()) print(fcell{dx:.2f}, p_max{p.max():.1f}) for i in range(1, len(p_peaks)): chg abs(p_peaks[i] - p_peaks[i-1]) / p_peaks[i-1] * 100 print(f{sizes[i-1]}-{sizes[i]}: char {chg:.2f}%)generate_grid需要自己封装第2.3节代码。收敛检查是接触力学分析里最容易被忽略的步骤很多二次开发源码只给一组固定网格换工况就失真。要养成每次改变载荷或横移量后都做一次收敛测试。5.2 用Numba把矩阵迭代提速法向迭代里np.linalg.lstsq是瓶颈。对于几千单元的接触集单次需要几十毫秒。动态轮轨仿真要跑数万步这时候需要用共轭梯度法替代直接求逆并用Numba编译内循环。接触矩阵对称正定共轭梯度法迭代次数通常小于50比最小二乘快一个量级。from numba import jit jit(nopythonTrue) def cg_solve(C_loc, b, max_iter50, tol1e-8): x np.zeros(C_loc.shape[0]) r b - C_loc x p r.copy() rs r r for _ in range(max_iter): Ap C_loc p alpha rs / (p Ap) x alpha * p r - alpha * Ap rn r r if rn tol * tol * C_loc.shape[0]: break beta rn / rs p r beta * p rs rn return x然后在solve_normal里把两次lstsq替换为cg_solve。注意Numba不支持np.linalg.lstsq但支持矩阵乘法和数组操作。内存上需要把接触子矩阵组织成连续数组以float64存储。对于5000单元以内的中型问题JIT版本能在0.1秒内完成一次法向求解足够嵌入多步仿真。5.3 多接触斑的自动识别与独立统计非赫兹问题常见结果是接触斑不连通比如轮缘贴靠时产生两个分离的接触区。分析时需要用连通域标记统计每个接触斑的面积、峰值压力和合力而不是只看全局最大压力。from scipy.ndimage import label # p_2d 是二维压力分布 mask p_2d 0.01 * p_2d.max() labeled, num_features label(mask) # 4邻域默认 for k in range(1, num_features 1): yy, xx np.nonzero(labeled k) area len(yy) * dx * dy peak p_2d[yy, xx].max() print(fpatch {k}: area{area:.2f} mm2, peak{peak:.1f} MPa)label默认按4邻域聚合相邻单元轮轨接触斑通常沿滚动方向条带状4邻域足够。若两个接触斑靠得极近距离小于一个单元会被合并成一个这时需要减小网格尺寸或调低阈值。多接触斑时第3章的仿射迭代收敛较慢建议初始化接触集时把两个区域分开求解再组合压力最后再做全局平衡判别迭代次数能降一半。轮轨蠕滑非稳定工况下建议将自旋项与蠕滑项分开积分避免切向递推步长耦合导致出口处应力振荡。本文还有配套的精品资源点击获取