企业官网建设流程全解析

简介：这是一套面向铁路车辆动力学与轮轨关系研究的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 是垂直距离阈值，单位与型面坐标一致，通常是毫米。阈值取得太大会把非接触点误判为接触，取得太小则可能漏掉轮缘贴靠时刻的接触点。一个快速自检方法是：当 yw=0 时，左右滚动半径应该非常接近且对称，如果左右接触点横向坐标之差超过一个采样间距，就要检查坐标原点约定。

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 如果只返回正值，调用端要把左侧接触角取反。验证方法是让轮对横移 y=0，左右接触角应接近相等但符号相反，否则就是坐标约定搞错了。另一个容易踩的坑是型面数据方向：有些文件把车轮型面的横向坐标从左到右定义，有些定义成从右到左，直接调用 gradient 会得到完全相反的斜率。

2.3 等效锥度：没有独立文件时的计算路径

等效锥度不是某个 .m 文件单独算出来的，而是通过 rolling_radius.m 得到滚动半径差 Δr = rl - rr，再按定义求平均斜率。UIC 519 线性化定义是在一段横移半幅值 A 内，取 Δr(y) 对 y 的线性拟合斜率。工程上常用 A=3mm，也就是 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 明显不是线性关系，拟合残差会很大，说明该横移范围内接触点已经跳到轮缘，单一等效锥度不再适用，应该改用非线性函数表示，否则动力学方程里的蛇行项会偏大。

下方是常见测试工况参数，可直接用于核对程序输出：

这个参数表也可以作为单元测试的基准。如果 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 系数的经验值：

这些是文献常引用的近似值，精确值应查 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 定律，轮缘接触点的位置由接触几何决定，法向力比例与两个接触点的侵入量有关。实际程序里通常用迭代求解，流程如下：

1. 根据轮对横移量和摇头角，用 rolling_radius.m 判断是否存在第二个接触点。
2. 若存在，初始化法向力 N1 = 0.4W、N2 = 0.6W。
3. 分别调用 onept_normal 计算两个接触斑尺寸。
4. 修正 N1、N2，使总法向力等于轮荷，同时满足垂向力平衡和几何约束。
5. 迭代到误差小于 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 往往会让接触几何模块在 y=0 附近来回振荡，建议先加一个 1 mm 的横移初值。运行后从位移时程里提取蛇行频率，再看不同速度下特征值实部的符号变化，实部由负变正的速度就是线性临界速度。

这一步能解释为什么在仿真曲线里高频振动消失得很快：悬挂阻尼把高频分量滤掉了，剩下的 1~3 Hz 分量是蛇行。若计算结果发散，先检查等效锥度是不是取到了轮缘贴靠以后的异常值，再检查蠕滑力方向是否与轮对运动方向一致。用 matlab 优化工具箱做参数扫描时，通常把临界速度作为目标函数，用一系纵向刚度作为设计变量，但要注意扫描步长不要越过非线性分岔点。

5. 把轮轨计算程序跑出可信结果的三个细节

细节一是接触点搜索前的型面对齐。rolling_radius.m对型面坐标原点极其敏感，钢轨型面通常以轨顶中心为原点，车轮型面以名义滚动圆为原点。从原始数据读入后，先画一条型面曲线，检查左右轮是否对称，钢轨轨顶是否在 y=0 处。常用做法是重采样到等弧长，采样间距取 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 程序从接触几何到蠕滑力的链路越可信。

本文还有配套的精品资源，点击获取