简介面向铁路轮轨接触力学研究者的Python简化模型专门处理超出经典赫兹理论适用范围的轮轨接触问题覆盖非线性变形、滑动接触、黏着特性与滚动接触疲劳等非赫兹因素。模型基于Piotrowski-Kik理论实现可模拟轮轨动态响应预测接触压力分布、磨损区域及疲劳损伤演化为列车设计、运行条件优化与轨道维护策略提供量化依据。压缩包共13个文件包括5个Python核心模块、配置文件与文档以及车轮钢轨廓形文件整体仅49KB轻量易用。Python脚本负责建模与数值计算Markdown与reStructuredText文档提供使用说明和理论背景便于快速上手。目前已有263人学习下载适合具备力学基础与Python编程经验的科研人员、轨道工程师及车辆工程专业学生既可用于教学演示也能作为二次开发与参数校核的起点。1. 非赫兹轮轨接触问题简化模型在解决什么轮轨接触分析里Hertz 解总是最先被想起但整车动力学、磨耗预测真正面对的是磨耗车轮和轮缘贴靠下的接触接触斑远离椭圆间隙函数无法用二次曲面近似Hertz 解会同时算错接触斑面积与压力峰值。非赫兹问题的轮轨接触力学简化模型用条带式半赫兹假设替代精确数值接触解在损失少量精度的前提下把单次接触计算压到毫秒级。这套以 Python 实现、zip 分发的工程适合两类人做动力学仿真的工程师要快速稳定的接触力做代理模型的研究者要可批量调用的函数接口。2. 非赫兹接触简化模型的数学原理虚拟贯入与条带降维要正确使用并修改这类工程代码先要理解它把完整接触问题压缩成了哪几步。核心就是两条用虚拟贯入量代替求解真实接触面积用条带分解把二维接触问题拆成一串一维线接触问题。2.1 Hertz 假设为什么在轮轨接触里撑不住Hertz 接触理论建立在几个前提下接触体在接触域内可近似为二次曲面、接触斑尺寸远小于曲率半径、法向与切向解耦、材料均匀且表面无摩擦。轮轨接触对这五条几乎逐条破坏。磨耗后的车轮踏面在接触域内高阶项明显等效曲率半径法会把非椭圆项直接抹掉轮缘接触时同时存在踏面和轮缘两个接触块小半径曲线和大横向曲率让接触斑宽度达到曲率半径的百分之几不再是小尺寸问题。最典型的现象是马鞍形或带缺口的接触斑中间贴合、前后缘翘起Hertz 椭圆解会把根本不接触的区域也包进去。用实测廓形算出的间隙函数往往呈多谷底结构直接套 Hertz 公式得到的接触斑和压力椭球只能是看着像。和 CONTACT 精确解相比斑面积与压力峰值误差超过 20% 是常态。这不是 Hertz 公式的错而是前提失效工程上需要一条介于完全数值解和Hertz 近似之间的路径。2.2 虚拟贯入法用穿透量确定接触斑边界Kik-Piotrowski 的虚拟贯入法是这类简化模型最常见的地基。思路是让两个接触体沿法向虚拟重叠一个微小量 Δ把几何间隙小于该穿透量的区域当作候选接触域。设车轮滚动方向为 x、横向为 y未变形间隙可写成g(x, y) g₀(y) x² / (2R_w)其中 g₀(y) 是横向位置 y 处的最小间隙R_w 是车轮在该处的滚动半径钢轨沿纵向视为等截面所以纵向曲率只由车轮贡献。施加虚拟贯入量 Δ 后接触条件 g(x,y) ≤ Δ 给出每个条带的纵向半长 a(y) sqrt(2 R_w (Δ - g₀(y)))当 Δ g₀(y)。代码里就是一条向量化计算a np.sqrt(np.maximum(2.0 * R_w * (penetration - gap), 0.0))说明np.maximum把负值截成 0对应间隙大于贯入量的条带无接触R_w取该横向位置的车轮滚动半径若轮对存在摇头角纵向曲率还要叠加轮轨相对曲率。注意Δ 不是真实压入量而是数值参数。Δ 太小会让接触斑对廓形噪声敏感太大又会把实际不接触的谷底卷进来常见做法是取 0.020.1 mm并配合最小条带宽度阈值过滤孤立点。2.3 条带法把二维接触拆成一串一维线接触拿到每个条带的半长 a(y) 后法向压力近似为沿 x 的半椭圆分布 p(x,y) C(y) sqrt(a(y)² - x²)。比例系数 C(y) 由材料参数和该条带局部等效曲率决定工程实现里更常见的是先按统一刚度假设计算压力形态再整体缩放使总法向力等于给定轮重。两种做法对接触斑边界形状没有影响区别只在压力幅值分布。方法接触斑形状压力分布单次耗时量级适用场景Hertz 解析解椭圆半椭球微秒级新廓形、快速筛选、初值条带虚拟贯入法任意形状条带拼合条带半椭圆 整体缩放毫秒级磨耗廓形、大样本迭代、代理模型CONTACT 边界元任意形状二维离散满足边界积分方程秒级以上基准校核、疲劳机理研究这张表的边界要清楚Hertz 用来做退化极限验证CONTACT 用来当裁判真正跑批量计算的是中间的简化模型。条带分解的代价是压力分布不严格满足弹性半空间的 Boussinesq 约束表现为斑边缘压力偏大、峰值偏低如果分析目标是疲劳裂纹萌生这类对压力细节敏感的问题应在简化模型后面再接一个局部精算步骤。3. Python 实现从轮轨廓形到非赫兹接触斑与压力分布下面给出最小可跑的实现路径。工程包里通常会拆成廓形处理、条带扫描、压力组装三个模块按这个顺序把关键函数写出来。3.1 廓形读取与预处理插值、平滑、统一网格轮轨廓形文件一般是两列文本横向坐标 y 和垂向坐标 z单位都是 mm。原始测点间距通常 0.51 mm直接拿来算接触斑会把条带边界切成台阶所以先重采样import numpy as np from scipy.interpolate import CubicSpline def load_profile(path): 读取两列廓形横向坐标 y(mm)、垂向坐标 z(mm) return np.loadtxt(path, unpackTrue) def resample_profile(y, z, y_start, y_end, dy0.1): cs CubicSpline(y, z) # 默认边界条件 y_new np.arange(y_start, y_end dy, dy) return y_new, cs(y_new)说明CubicSpline对离散点做三次样条插值返回可调用的插值函数重采样间距 dy 取 0.1 mm 时接触斑横向能分到几十到上百个条带。磨耗廓形里有高频测量噪声时三次样条会产生小幅振荡建议先做 35 点滑动平均或改用scipy.signal.savgol_filter否则噪声会被虚拟贯入法放大成虚假条带。然后是间隙计算。车轮和钢轨各自重采样到同一组横向坐标后间隙取垂向坐标之差gap z_wheel - z_rail。注意符号约定要统一正值表示分离负值表示初始过盈实际程序里按廓形文件的正负方向核对。3.2 条带扫描接触斑边界求解核心函数是按横向位置扫描出每个条带的纵向半长def contact_strips(y_grid, gap, penetration, R_w): 按条带扫描非赫兹接触斑返回 [(y, a), ...] strips [] for y, g in zip(y_grid, gap): if g penetration: # 虚拟贯入后存在重合 a np.sqrt(2.0 * R_w * (penetration - g)) if a 0.1: # 过滤数值噪声条带 strips.append((y, a)) return strips说明g penetration与公式 a(y) sqrt(2R_w(Δ - g₀(y))) 完全对应0.1 mm 的半长阈值是工程经验值用来过滤浮点误差和廓形毛刺产生的孤立条带。如果接触斑本身很小薄轮缘贴靠阈值要下调到 0.01 mm否则会把真实接触滤掉。3.3 压力组装与法向力缩放每个条带按半椭圆压力分布组装再整体缩放满足目标轮重def assemble_pressure(strips, dy, N_target, C_ref1.0): nx 201 rows, volume [], 0.0 for y, a in strips: x np.linspace(-a, a, nx) p np.sqrt(np.maximum(a * a - x * x, 0.0)) # 半椭圆形态 volume np.trapz(p, x) * dy # 压力体积 rows.append((y, a, x, p)) scale N_target / volume if volume 0 else 0.0 return [(y, a, x, p * scale) for y, a, x, p in rows]说明np.trapz是梯形积分NumPy 2.0 之后官方改名np.trapezoid两个名字均可。C_ref1.0表示先算单位比例系数下的压力形态scale保证总法向力严格等于 N_target这个缩放不改变压力形态是简化模型的固定近似。如果做位移控制给定贯入量求轮重跳过缩放直接把volume当轮重输出此时 C_ref 要换成由材料参数算出的等效刚度。数据项常见单位说明y / zmm廓形坐标接触点附近 z 一般取 0R_wmm滚动半径与廓形坐标单位一致N_targetN目标轮重缩放基准dymm条带宽度通常 0.10.2提示整体缩放不改变压力形态这是简化模型的主要近似误差来源。若目标是疲劳寿命这类对压力峰值敏感的指标应在条带模型后再加一步局部精算。3.4 载荷控制下的贯入量二分迭代法向力随贯入量单调增加所以给定轮重反求贯入量用二分法最稳def integrate_normal_force(strips, dy): return dy * sum(0.5 * np.pi * a**2 for _, a in strips) def solve_penetration(y_grid, gap, R_w, N_target, dy): lo, hi 0.0, 0.2 # 贯入量区间单位 mm for _ in range(40): mid 0.5 * (lo hi) N_now integrate_normal_force( contact_strips(y_grid, gap, mid, R_w), dy) if N_now N_target: lo mid else: hi mid return 0.5 * (lo hi)说明条带半椭圆压力对 x 的积分是 πa²/2乘上条带宽度 dy 就是该条带对总法向力的贡献这是integrate_normal_force的解析依据。二分法的收敛条件是 N(Δ) 单调这在本模型里天然成立40 次迭代把区间压到 2⁻⁴⁰ 量级工程上可以在区间宽度小于 1e-4 mm 时提前退出。4. 切向力计算与 FASTSIM 落地参数标定和 Python 工程目录4.1 法向解之后为什么必须接切向迭代法向解得到的是接触斑形状和压力分布但动力学需要的是蠕滑力磨耗和黏着分析依赖黏滑边界不接切向解简化模型只完成了一半。非赫兹条带斑上做切向解最常见的是 Kalker 简化理论里的 FASTSIM 算法从接触斑入口向前积分每一步用库仑界限 μp 截断剪应力。它的优势和条带式法向模型天然同构每个条带独立求解计算量可以忽略。4.2 FASTSIM 条带适配柔度系数与黏滑边界FASTSIM 原本假设椭圆接触斑和一个统一柔度系数非赫兹适配的做法是每个条带用各自的半长 a(y) 和局部压力def fastsim_strip(x_grid, p, xi, mu, L1): 纵向蠕滑的条带 FASTSIM x_grid: 从 -a 到 a 的纵向网格 p: 该条带的法向压力 xi: 纵向蠕滑率 L1: 纵向柔度常见量级 1e-9 ~ 1e-8 m^3/N q np.zeros_like(x_grid) q_acc 0.0 for i in range(1, len(x_grid)): dx x_grid[i] - x_grid[i-1] q_acc - xi / L1 * dx # 柔度关系的增量形式 q_acc np.clip(q_acc, -mu * p[i], mu * p[i]) # 库仑界限 q[i] q_acc return q说明入口处剪应力为零每一格按 dτ -ξ/L1·dx 累加超过 ±μp 就截断截断点即黏滑边界负蠕滑对应负剪应力符号由蠕滑率方向决定。L1 越大剪应力增长越慢斑内更容易全黏着L1 偏小则过早饱和蠕滑力被低估。纵向柔度按 Kalker 系数表 L 8a/(3GC₁₁) 标定横向柔度类似但用 C₂₂。4.3 必调参数表与三个常见误用参数常见取值主要影响调参方向横向条带间距 dy0.10.2 mm接触斑边界锯齿磨耗廓形取 0.1虚拟贯入量 Δ0.020.1 mm斑面积、压力峰值偏大高估面积偏小易振荡纵向柔度 L1/L2按 C₁₁/C₂₂ 标定黏滑边界、蠕滑力与实验或 CONTACT 对比定摩擦系数 μ0.30.6饱和区大小雨天、污染轨面下调条带纵向点数 nx≥201黏滑边界锯齿接触斑细长时加大常见误用有三个一是把 Hertz 等效椭圆直接当条带半长分布这等于没做非赫兹修正二是忽略自旋蠕滑曲线通过工况下自旋对横向蠕滑力影响显著FASTSIM 里要补自旋项不能只算纵向三是网格太粗黏滑边界出现锯齿状跳变纵向点数少于 100 时尤其明显。4.4 zip 工程包的目录组织和运行入口这类 zip 工程常见的组织方式是把廓形放数据目录、算法按模块拆分根目录留一个命令行入口。解压和运行的最小命令序列是unzip 非赫兹问题的轮轨接触力学简化模型_Python_下载.zip -d wheel_rail_model cd wheel_rail_model python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate pip install -r requirements.txt python main.py --wheel profiles/wheel.txt --rail profiles/rail.txt --load 80000说明先解压再建虚拟环境避免依赖污染系统 Python。如果提示python was not found; run without arguments to install是解释器没进 PATH要么重装时勾选 Add to PATH要么用py -3调用 Windows 的 Python 启动器VSCode 里跑之前按 CtrlShiftP 选解释器否则终端用的可能是另一个环境。Linux 下没有 unzip 时先sudo apt install unzip也可以用7z x后者对损坏 zip 的容忍度更高。--load 80000是目标轮重N按实际轴重换算。5. 验证方法与 zip 解压排错把模型的边界圈出来验证简化模型我一般分三层做先跑 Hertz 退化自检再与 CONTACT 对比最后用实测蠕滑力标定柔度。这套顺序能把离散误差、模型近似误差和参数标定误差分开排查。5.1 Hertz 退化自检对二次曲面间隙条带模型会自动退化为 Hertz 椭球压力因为半椭圆条带拼合出来的正是椭球面。构造球面车轮加平面钢轨的算例python main.py --wheel profiles/sphere_r460.txt --rail profiles/plane.txt --check-hertz以面积误差 ±5%、压力峰值 ±8% 为验收线超差说明插值、条带离散或缩放环节有系统性错误先查廓形单位是否混用mm 与 m。5.2 与 CONTACT 对比时看什么与 CONTACT 对比时不要只盯压力峰值。先按横向坐标对比接触斑宽度分布再对比压力沿横剖面的积分曲线最后看黏滑边界。简化模型系统性高估接触面积 10%20% 属于模型固有行为因为条带半椭圆压力不满足弹性半空间的边界积分条件误差主要来自斑边缘的孤立条带按压力贡献剔除面积占比小于 1% 的条带能明显改善。5.3 zip 伪加密与 EOCD 报错的快速处理解压报invalid zip archive: could not find EOCD时多数情况是文件截断或格式伪装Linux 下先file model.zip看真实格式再重新下载。提示要密码但文件是伪加密时内容实际没加密只是通用位标志第 0 位被置位清掉即可import zipfile def clear_fake_encryption(src, dst): with zipfile.ZipFile(src) as zin, \ zipfile.ZipFile(dst, w, zipfile.ZIP_DEFLATED) as zout: for info in zin.infolist(): data zin.read(info.filename) info.flag_bits ~0x0001 # 清除加密标志位 zout.writestr(info, data) clear_fake_encryption(model.zip, model_fixed.zip)注意zipfile 读取时同样会被伪加密标志拦住所以要先按原样读出再写入新文件如果连读都报错那就是真加密或文件损坏优先重新下载不要尝试爆破。本文还有配套的精品资源点击获取