简介这是一套面向铁路车辆动力学与轮轨关系研究的MATLAB代码包对需要分析轮轨接触特性、预测蛇行稳定性、评估曲线通过能力的工程师和研究生尤为适用。压缩包内共12个文件以11个m脚本为主体辅以1个txt说明文件整体仅15KB。代码从接触斑、接触角、等效锥度等基础几何量切入逐步覆盖单轮对/双轮对模型的蠕滑力与法向力计算方程脚本进一步给出轮轨动力学的核心数学描述。所有脚本围绕轮轨接触问题紧凑组织模块边界清楚既适合按步骤深入学习也便于直接改造嵌入自己的仿真流程帮助研究者快速得到接触应力分布、轮轨作用力等关键数据为车轮磨耗和脱轨安全性分析提供支撑。已有965人学习使用对相关课程设计、科研实验或工程优化都是轻量而实用的参考资料。1. 轮轨关系的数值计算为什么绕不开接触几何轮轨关系里最容易被低估的是接触几何。很多人把用于计算轮轨关系的 MATLAB 程序当成纯函数库读入型面文件、调用 contact_angle.m 出个角度就算完成了一次计算。但真正上线路工况时接触角变化率在轮缘贴靠附近很容易出现量级跳变导致蛇行失稳临界速度的仿真结果差出 20 km/h。这里的差别不是算法原理不够而是程序包里 rolling_radius.m、onept_normal.m、twopt_creep.m 等文件各自承担的物理环节没有被拆开理解。这套以轮轨关系为核心的 matlab 程序正好覆盖了接触几何、法向接触、蠕滑力到轮对动力学方程的完整链路。适合正在做车辆动力学仿真、轮轨磨耗评估的工程师也适合刚从各类资源站把 matlab 程序下载下来、却不知道从哪个文件开始读的初学者。读完你至少能回答三个问题每个 .m 文件该喂什么参数、输出是什么物理量、跑出来的数值能不能直接拿去算临界速度。2. 从 rolling_radius.m 到接触角接触几何怎么算2.1 接触点搜索rolling_radius.m 的核心假设接触几何计算要区分“型面几何”和“接触力学”。rolling_radius.m 处理的是几何问题给定轮对横移量 y左右车轮和两根钢轨各自的空间位置找出一对接触点然后读取该点的滚动半径和横向位置。常见做法是把轮轨型面离散成点列用最小距离搜索替代解析求交因为实际型面是实测点云不存在闭合的曲线方程。function [rl, rr, xl, xr] rolling_radius(yw, wheel_z, wheel_y, rail_z, rail_y) % 轮对横移 yw 下的左右滚动半径和接触点位置 % wheel_y, wheel_z : 车轮型面点列相对名义滚动圆 % rail_y, rail_z : 钢轨型面点列轨顶坐标系 dy wheel_y(2) - wheel_y(1); % 型面采样间距 [~, idx] min(abs(rail_y - yw)); % 粗略定位钢轨上与轮对中点对应的点 rl 1e9; rr 1e9; xl NaN; xr NaN; for k 1:numel(wheel_y) d sqrt((wheel_y(k) - yw - (rail_y - idx*dy)).^2 ... (wheel_z(k) - rail_z).^2); [dmin, ~] min(d); % 对每个车轮点找最近钢轨点 if dmin 1e-3 wheel_y(k) 0 % 左侧接触阈值可调 if dmin rl, rl dmin; xl k; end elseif dmin 1e-3 wheel_y(k) 0 % 右侧接触 if dmin rr, rr dmin; xr k; end end end % 返回的仍是距离真正滚动半径需要叠加轮对自身的滚动圆半径 end这段代码用双重循环完成最小距离搜索优点是结构直观方便核对几何关系缺点是计算量随型面点数平方上升。实际工程程序里会把型面重采样到等弧长再用最近邻搜索或二分法提速。输入轮对横移量 yw 必须是标量如果要算等效锥度就要在外层对 yw 扫掠滚动半径差 Δr rl - rr 是后面所有非线性接触参数的基础。代码里的 1e-3 是垂直距离阈值单位与型面坐标一致通常是毫米。阈值取得太大会把非接触点误判为接触取得太小则可能漏掉轮缘贴靠时刻的接触点。一个快速自检方法是当 yw0 时左右滚动半径应该非常接近且对称如果左右接触点横向坐标之差超过一个采样间距就要检查坐标原点约定。2.2 接触角contact_angle.m 为什么用数值微分接触角是指轮轨接触点处公切面与水平面之间的夹角更准确地说是接触点处轮缘切线方向与竖直方向的夹角。这个角度影响轮轨法向力的方向也决定了两点接触出现的时机。由于型面是离散点列程序里通常不保留解析解而是用型面坐标的梯度近似。function alpha contact_angle(wheel_y, wheel_z, point_idx) % 从离散型面计算接触角 dy wheel_y(2) - wheel_y(1); dz gradient(wheel_z) ./ gradient(wheel_y); % 踏面斜率 slope dz(point_idx); alpha atan(slope); % 弧度 end这段代码只用了一行 gradient但实际使用时要注意三点。第一原始型面往往有测量噪声直接 gradient 会产生毛刺我一般先用 sgolayfilt 做局域多项式平滑阶数取 3窗口长度取 15 到 31 个点。第二轮缘根部斜率接近 90 度atan 的数值误差会被放大建议改用两侧差分再查表。第三接触角随横移量的导数在轮缘贴靠点附近不连续这个不连续性正是两点接触模型的触发条件后续 twopt_normal.m 会用到这个判据。2.2.1 左右接触角的符号约定左轮和右轮的接触角符号相反因为坐标系方向不同。程序包里的 contact_angle.m 如果只返回正值调用端要把左侧接触角取反。验证方法是让轮对横移 y0左右接触角应接近相等但符号相反否则就是坐标约定搞错了。另一个容易踩的坑是型面数据方向有些文件把车轮型面的横向坐标从左到右定义有些定义成从右到左直接调用 gradient 会得到完全相反的斜率。2.3 等效锥度没有独立文件时的计算路径等效锥度不是某个 .m 文件单独算出来的而是通过 rolling_radius.m 得到滚动半径差 Δr rl - rr再按定义求平均斜率。UIC 519 线性化定义是在一段横移半幅值 A 内取 Δr(y) 对 y 的线性拟合斜率。工程上常用 A3mm也就是 y ∈ [-3, 3] mm 范围因为在这个范围内接触点多数还处于踏面区域等效锥度能反映蛇行稳定性。function lambda equivalent_conicity(y_range, rl_func, rr_func, dy) % 在给定横移范围内拟合等效锥度 y -y_range:dy:y_range; dr arrayfun((yy) rl_func(yy) - rr_func(yy), y); p polyfit(y, dr, 1); % 一次线性拟合 lambda p(1); % 斜率即等效锥度 end参数说明y_range 是轮对横移半幅值单位毫米dy 是扫掠步长一般取 0.1 mm 或 0.2 mm。若 dr 与 y 明显不是线性关系拟合残差会很大说明该横移范围内接触点已经跳到轮缘单一等效锥度不再适用应该改用非线性函数表示否则动力学方程里的蛇行项会偏大。下方是常见测试工况参数可直接用于核对程序输出参数典型值说明型面采样间距0.1~0.5 mm越小接触点越稳定但计算量增大横移扫掠范围±6 mm覆盖踏面接触到轮缘贴靠步长0.1~0.2 mm等效锥度拟合至少需要 30~60 个点接触角平滑窗口15~31 点sgolayfilt 平滑窗口长度必须是奇数最小距离阈值0.1~1e-3 mm阈值过大会把非接触点误判为接触这个参数表也可以作为单元测试的基准。如果 rolling_radius.m 在 ±6mm 扫掠范围内输出的滚动半径差曲线没有单调性多半是型面基准点没对齐而不是程序算法有问题。对齐型面时通常把钢轨轨顶中心或名义滚动圆作为坐标原点左右对称的型面数据要各自翻转一次。3. 蠕滑力与法向力onept 与 twopt 系列程序在算什么3.1 法向接触问题onept_normal.m 的 Hertz 背景接触几何给出接触点位置后下一步是计算接触斑上的法向力分布。onept_normal.m 处理的是单个接触点的法向问题理论基础是 Hertz 接触理论。Hertz 理论假设接触体为弹性半空间、接触面光滑、接触区域为椭圆输入是轮轨接触点处的主曲率、法向力和材料参数输出是接触椭圆的长半轴 a、短半轴 b 和最大接触压力 p0。function [a, b, p0] onept_normal(N, Rw, Rr, E, nu) % N 为法向力Rw/Rr 为车轮和钢轨在接触点处的曲率半径 % E 为弹性模量nu 为泊松比 E_star E / (1 - nu^2); rhos 1/Rw 1/Rr; % 曲率和 delta (3*N/(4*E_star))^(1/3) * rhos^(2/3); % 法向接近量 % 实际 a/b 需要由曲率差查表得到这里给出极限简化 a sqrt(4*N*Rr/(pi*E_star*delta)); % 长半轴近似 b a * 0.8; % 短半轴按接触椭圆率近似 p0 3*N / (2*pi*a*b); end这段代码是极端简化实际 onept_normal.m 里会数值求解椭圆积分并根据接触点处主曲率差计算椭圆率。很多初学者直接把轴重除以轮对数填进 N算出的接触斑会小于实测因为忽略了轨道不平顺引起的动轮载。我一般会先用动力学仿真输出轮轨垂向力再把这个时变力传给法向接触模块。若不想引入动力学仿真至少要把静轮重乘以 1.2~1.5 的动载系数。3.2 切向蠕滑onept_creep.m 的线性理论切向问题处理的是车轮和钢轨之间“既滚动又滑动”的蠕滑现象。蠕滑率分为纵向蠕滑 vx、横向蠕滑 vy 和自旋蠕滑 w。Kalker 线性理论给出了蠕滑力与蠕滑率在小编程范围内的线性关系。onept_creep.m 的文件名表明它实现的是单接触点蠕滑力计算输入为接触椭圆尺寸、材料剪切模量和蠕滑率输出为纵向力 Fx、横向力 Fy 和自旋力矩 Mz。function [Fx, Fy, Mz] onept_creep(vx, vy, w, a, b, G, C) % vx 纵向蠕滑率, vy 横向蠕滑率, w 自旋蠕滑率 % a, b 接触椭圆长短半轴, G 剪切模量, C [C11 C22 C23] Fx -G*a*b*C(1)*vx; Fy -G*a*b*(C(2)*vy C(3)*b*w); Mz -G*a*b*(C(2)*b*w) * 0.01; % 自旋项很小通常忽略 end这里的 C11、C22、C23 是 Kalker 系数与泊松比和接触椭圆率 a/b 有关。程序里如果只给了一组固定常数线性范围之外的结果会明显偏大只能在蠕滑率小于 1% 时使用。实际应用中蠕滑率会超过 1%此时要引入饱和修正常见的做法是沈志云-Hedrick-Elkins 方法用合成蠕滑力把线性值压缩到库仑极限内。twopt_creep.m 的输出如果与 onept_creep.m 的差异超过 30%多半是两点接触时切向力非线性叠加没有处理好。为方便核对系数范围下面给出一组 Kalker 系数的经验值泊松比a/bC11C22C230.251.04.054.050.8980.252.05.203.841.550.281.04.124.120.8910.282.05.303.911.52这些是文献常引用的近似值精确值应查 Kalker 数据表。使用时要留意 onept_creep.m 里参数的顺序有些程序把蠕滑率按 [vx, vy, w] 排列有些按 [w, vy, vx] 排列顺序错了会产生一个很隐蔽的符号错误导致横向力方向反相。3.3 两点接触twopt_normal.m 和 twopt_creep.m 的分工在曲线通过或轮对横移较大时轮缘和钢轨侧面同时接触出现两个接触点。twopt_normal.m 负责把总法向力分配到两个接触点twopt_creep.m 负责分别计算两个接触点上的蠕滑力再合成为轮对受到的力和力矩。分配法向力的常见方法是假设两个接触点各自遵循 Hertz 定律轮缘接触点的位置由接触几何决定法向力比例与两个接触点的侵入量有关。实际程序里通常用迭代求解流程如下根据轮对横移量和摇头角用 rolling_radius.m 判断是否存在第二个接触点。若存在初始化法向力 N1 0.4W、N2 0.6W。分别调用 onept_normal 计算两个接触斑尺寸。修正 N1、N2使总法向力等于轮荷同时满足垂向力平衡和几何约束。迭代到误差小于 0.1%再调用 twopt_creep.m。这个流程里最容易被忽略的是第 4 步的约束两个接触点的垂向合力必须等于轮对垂向力横向力的分配则与轮缘角有关。twopt_creep.m 的输出会包含轮缘力曲线通过仿真里轮缘力出现台阶式上升就是两点接触被激活的标志。如果运行后发现力不守恒先检查迭代初值再看 onept_normal.m 输出的接触斑尺寸是否在一开始就发生了跳变。4. 从单轮对到转向架wheelset.m 与 wheelset_suspension.m 里的动力学方程4.1 轮对运动方程equations.m 在组装什么接触几何和蠕滑力只是轮轨接触局部的结果要评估车辆稳定性需要把它们放进轮对的动力学方程。wheelset.m 和 equations.m 的作用是把接触参数转变成方程系数。一个单轮对有横移 y、摇头 ψ、侧滚 φ、垂向 z 和纵向 x 自由度简化模型可以只考虑横移和摇头。运动方程的一般形式为M·q̈ C·q̇ K(q) F_contact(q, q̇)其中刚度项来自轮轨接触几何的等效刚度和悬挂刚度F_contact 来自蠕滑力。在 equations.m 中矩阵 M 是质量矩阵K 往往是随横移量变化的非线性矩阵因为等效锥度和接触角不是常数。function [M, C, K] equations(mw, Iwx, Iwz, lambda, g, e, Kc, Cc) % 单轮对横移/摇头模型的线性化方程矩阵 % mw 轮对质量, Iwx 侧滚惯量, Iwz 摇头惯量, lambda 等效锥度 % e 轮对滚动圆横向跨距之半, g 重力加速度, Kc/Cc 一系悬挂刚度和阻尼 M diag([mw, Iwz]); C [Cc, 0; 0, 0]; % 只给横向加阻尼 K [mwg*lambda/e, -2*Kc; ... % 重力刚度项 悬挂刚度 Kc, 2*Kc*e]; % 摇头自由度对角项 end这是线性化示例实际 wheelset.m 里还会有侧滚自由度与摇头的耦合项。注意这里的重力刚度 λ/e 来自轮对中心升高与横移的关系而不是传统意义上的结构刚度。等效锥度越大这个项越大轮对越容易发生蛇行失稳。若直接修改程序的参数矩阵会发现临界速度随等效锥度增大而下降这与实测规律方向一致。4.2 wheelset_suspension.m一系悬挂怎么影响接触力反馈wheelset_suspension.m 把轮对和转向架构架之间的弹簧阻尼连接加了进去对应一系悬挂的纵、横、垂向刚度和阻尼。悬挂参数不仅提供恢复力也改变了蠕滑力的反馈路径。简单模型的悬挂力可以表示为function F wheelset_suspension(q, dq, Ks, Cs, track_irregularity) % q,dq 轮对位移和速度; Ks,Cs 悬挂刚度阻尼矩阵 % track_irregularity 是轨道不平顺激励向量 F -Ks*(q - track_irregularity) - Cs*(dq - track_irregularity); end实际工程中一系悬挂纵横向刚度在 5~50 MN/m 范围内阻尼在 5~50 kNs/m。程序注释里如果写了单位提示不要忽略把 MN/m 当成 N/m 使用会导致特征值数量级完全错误。wheelset_suspension.m 里还可能包含抗蛇行减振器参数那是更高频的稳定性控制元件单独建模时要额外加一个串联刚度。4.3 用 ode45 求解和提取临界速度拿到 wheelset.m 和 wheelset_suspension.m 之后还要组合成状态方程才能求解。我的习惯是把二阶方程改写成一阶状态空间然后用 ode45 求解再对速度参数扫掠计算特征值。function dstate wheelset_state(~, state, params) % state [y; psi; dy; dpsi] y state(1); psi state(2); dy state(3); dpsi state(4); [M, C, K] equations(params.mw, params.Iwx, params.Iwz, ... params.lambda, params.g, params.e, ... params.Kc, params.Cc); rhs -C*[dy;dpsi] - K*[y;psi]; acc M \ rhs; % 解算加速度 dstate [dy; dpsi; acc(1); acc(2)]; end调用代码[t, s] ode45((t,s) wheelset_state(t,s,params), tspan, [0; 0.001; 0; 0]);参数说明tspan 要覆盖至少 10 个蛇行周期初速设为 0 往往会让接触几何模块在 y0 附近来回振荡建议先加一个 1 mm 的横移初值。运行后从位移时程里提取蛇行频率再看不同速度下特征值实部的符号变化实部由负变正的速度就是线性临界速度。这一步能解释为什么在仿真曲线里高频振动消失得很快悬挂阻尼把高频分量滤掉了剩下的 1~3 Hz 分量是蛇行。若计算结果发散先检查等效锥度是不是取到了轮缘贴靠以后的异常值再检查蠕滑力方向是否与轮对运动方向一致。用 matlab 优化工具箱做参数扫描时通常把临界速度作为目标函数用一系纵向刚度作为设计变量但要注意扫描步长不要越过非线性分岔点。5. 把轮轨计算程序跑出可信结果的三个细节细节一是接触点搜索前的型面对齐。rolling_radius.m对型面坐标原点极其敏感钢轨型面通常以轨顶中心为原点车轮型面以名义滚动圆为原点。从原始数据读入后先画一条型面曲线检查左右轮是否对称钢轨轨顶是否在 y0 处。常用做法是重采样到等弧长采样间距取 0.1 mm这样既保证接触点搜索稳定又不会让滚动半径差曲线出现锯齿。细节二是接触角平滑的窗口选择。用sgolayfilt平滑时窗口长度影响接触角导数的峰值。窗口太短轮缘贴靠点的接触角跳变仍然存在窗口太长会把真实的轮缘几何圆角抹掉。我一般先用 15 点窗口计算一次看接触角曲线在横移 3~5 mm 附近是否有平滑过渡如果还有毛刺逐步加到 31 点。注意接触角导数曲线是判断两点接触触发的关键最好单独输出一张图检查。细节三是用已知型面的参考值校验。以 LMA 踏面和 CN60 钢轨为例名义工况下 3 mm 等效锥度通常在 0.05~0.15 之间轮重 140 kN 时接触椭圆长半轴一般不超过 10 mm。把 computed_conicity 的结果和这些参考值对比如果差距超过一倍多半是接触几何模块的坐标或符号错误。更严格的验证是把横向力-横移曲线与多体软件 SIMPACK 或 UM 的输出做对比两者之间的差异应小于 5%。把两条曲线画在同一个图里重合度越高说明这套轮轨关系 matlab 程序从接触几何到蠕滑力的链路越可信。本文还有配套的精品资源点击获取