简介：一套 MATLAB 程序包，专门用于轮轨关系计算，面向车辆工程、轨道工程及轮轨动力学方向的科研人员、研究生和工程师。程序以单轮对和多轮对为对象，能够求解轮轨接触斑的尺寸与形状、接触角随轮对横移的变化规律，并进一步计算等效锥度这一关键稳定性指标，适合用于踏面匹配分析、车辆动力学仿真前期参数准备以及轮轨磨耗研究。压缩包共 12 个文件，包含 11 个 .m 源程序和 1 个说明文本，总大小仅 15KB，代码轻量、入口清晰。核心模块涵盖单点与双点接触判断、法向接触几何、蠕滑力计算、轮对运动方程及滚动半径求解等，配套单轮对与转向架级仿真脚本，方便按需调用和二次开发。目前已有 2210 人学习下载，对于正在开展轮轨接触机理研究或进行车辆动力学仿真的读者具有直接参考价值。

1. 轮轨关系计算入门：为什么先做接触几何再做蠕滑力

把轮轨关系当黑匣子的项目，基本都翻在第一步。我拆过的这套用于计算轮轨关系的 MATLAB 程序，不追求上万个节点的有限元接触模型，而是把赫兹法向接触、接触几何扫描和 FASTSIM 切向蠕滑力计算串成一条完整可复现的计算链，几百行代码就能输出接触斑尺寸、接触位置、轮轨力和蠕滑力矩。对做车辆动力学、钢轨廓形打磨和磨耗评估的工程师来说，这套程序的真正价值是把计算顺序理顺：先算接触几何，再算接触力学，最后才谈得上力与滑移的关系。适合两类人：一类是刚接触轮轨关系、想搞懂参数怎么影响结果的新手，另一类是已经会用商业软件、但需要批量扫参或做数据前处理的一线工程师。

2. 数值求解的“地基”：Hertz法向接触与FASTSIM切向蠕滑

轮轨关系的难点不在单个公式，而在法向问题和切向问题必须耦合着看。法向力决定接触斑的大小、形状和压力分布，切向力则取决于接触斑内部的局部滑移状态。工程上很少直接调用完整边界元接触解，因为单工况计算时间太长，批量扫横移量、摇头角时根本不现实。常用的做法是把问题拆成两层：第一层用赫兹理论求法向接触，第二层用 Kalker 简化理论（FASTSIM）求切向蠕滑力。这套 MATLAB 程序走的也是这个路线。

2.1 接触力学选型：为什么工程上不用全量 CONTACT，先用 Hertz 加 FASTSIM

全量接触解在学术上最严谨，它对接触斑内部每个网格点做弹性半空间积分，能给出非常细腻的应力分布和微滑移状态，但代价是慢。一个接触工况要算几十秒甚至几分钟，整车动力学仿真一跑几千个时间步，这种开销在实际项目中很难接受。FASTSIM 是基于 Winkler 弹性基础假设的简化方法，把接触斑离散成若干条带，沿滚动方向递推切向应力，速度比全量解快两三个数量级，精度在工程误差范围内。对于磨耗预测、稳定性校核、轮轨廓形优选这类任务，FASTSIM 的误差完全可以接受，而且参数透明，能看清楚每个调节项到底影响了什么。

FASTSIM 的输入主要是三样东西：法向接触参数（接触斑长短半轴、最大法向压力）、蠕滑率（纵向、横向、自旋）、摩擦系数与柔度系数。其中蠕滑率的定义必须前后一致，否则算出来的轮轨力方向会整体反转，后面第 4 章会专门说这个坑。

2.2 接触几何求解：用踏面和轨面离散点找接触点

接触力学算得再准，接触点找错了也白搭。这套程序里接触几何的核心思路是离散化扫描：把轮对踏面轮廓和钢轨轮廓分别用一组二维坐标点表示，给定轮对横移量和摇头角，在横向坐标上逐点搜索踏面与轨面之间的最小法向间隙。间隙最小处就是当前工况下的轮轨接触点。

摇头角的影响不能忽略。轮对摇头后，踏面轮廓在纵向发生错动，接触点位置会沿踏面横向偏移。常见做法是用小角度近似，把轮对坐标系下的踏面点旋转到轨道坐标系再参与间隙计算。这样做的好处是计算量小，同时能覆盖直线、曲线和道岔区的基本工况。以下代码块是核心搜索函数，输入左右两侧的廓形离散点，输出接触点坐标和最小间隙。

function [contact] = find_contact_point(rail_y, rail_z, wheel_y, wheel_z, y_shift, dpsi) % 轮轨接触几何搜索：逐横向位置计算法向间隙 % rail_y/rail_z : 轨面廓形离散点 (mm) % wheel_y/wheel_z : 踏面廓形离散点，已放在轮对局部坐标系 (mm) % y_shift : 轮对相对轨道的横移量 (mm) % dpsi : 轮对摇头角 (rad)，小角度近似 r_y = rail_y(:) - y_shift; % 轨道坐标平移，统一到轮对中心 w_y = wheel_y(:); w_z = wheel_z(:); % 将踏面点按摇头角旋转到轨道坐标系 rot_z = [cos(dpsi), -sin(dpsi); sin(dpsi), cos(dpsi)]; pts_w = rot_z * [w_y'; w_z']; wy = pts_w(1, :)'; wz = pts_w(2, :)'; % 统一横向坐标，用线性插值求同一 y 位置下的 z 坐标 y_min = max(min(r_y), min(wy)); y_max = min(max(r_y), max(wy)); yq = linspace(y_min, y_max, 2000)'; z_rail = interp1(r_y, rail_z(:), yq, 'linear', NaN); z_wheel = interp1(wy, wz, yq, 'linear', NaN); gap = z_rail - z_wheel; % 轨面减踏面，最小正值即接触位置 gap(isnan(gap)) = inf; [gap_min, idx] = min(gap); contact.y = yq(idx); % 接触点横向坐标 contact.z = 0.5 * (z_rail(idx) + z_wheel(idx)); contact.gap = gap_min; end

这段代码的逻辑是：先把轨道廓形平移到轮对坐标系，再用摇头角旋转矩阵修正踏面点位置，最后在统一横向网格上插值并求间隙。参数方面需要注意两点：第一，yq 的密度直接决定接触点分辨率，2000 个点对常规 60 轨和磨耗踏面够用，但如果要做轮缘接触分析，建议加密到 4000 以上；第二，interp1 的边界要显式指定为 NaN，否则廓形端点外的假接触会被当成真实接触点。实际项目里我发现，很多接触几何算错都是因为左右轨廓形数据长度不一致，导致插值结果在端部失真。

3. MATLAB实现：从输入参数到轮轨力输出的完整代码

有了接触几何和接触力学两部分基础，现在可以拼装完整的计算流程了。这套程序的主体逻辑并不复杂：读入廓形和工况参数，计算接触点位置，然后由赫兹理论求出接触斑尺寸，最后用 FASTSIM 积分得到轮轨力。难的是把每一步的输入输出接口对整齐，尤其要注意单位、坐标系方向和数据结构的设计。

3.1 输入数据组织：踏面廓形、轨面廓形与运行工况

我一般建议把输入数据分成三个层次：廓形数据、工况参数、控制参数。廓形数据是固定不变的，通常从廓形仪或设计图纸导出后存成 CSV 文件；工况参数是每次仿真要扫的变量，比如横移量、摇头角、蠕滑率；控制参数则包括网格密度、摩擦系数、收敛阈值等。程序运行前先统一加载，再用结构体传给各个计算函数。

下面的代码展示了我常用的加载方式，兼容 CSV 和 MAT 两种格式，适合批处理场景。

% load_inputs.m % 统一加载廓形与工况，返回结构体 function [par] = load_inputs(base_dir, case_id) par.rail = readmatrix(fullfile(base_dir, 'rail_profile.csv')); % 第一列横向坐标，第二列垂向坐标，单位 mm par.wheel = readmatrix(fullfile(base_dir, 'wheel_profile.csv')); % 工况参数也可以从表格批量读入 case_tbl = readtable(fullfile(base_dir, 'cases.csv')); row = case_tbl(case_id, :); par.y_shift = row.y_shift; % 横移量，mm par.dpsi = row.dpsi; % 摇头角，rad par.xi = row.xi; % 纵向蠕滑率 par.eta = row.eta; % 横向蠕滑率 par.spin = row.spin; % 自旋蠕滑率，1/mm par.Fn = row.Fn; % 法向轮轨力，N par.Rw = row.Rw; % 车轮滚动圆半径，mm par.Rr = row.Rr; % 轨顶横向曲率半径，mm par.mu = row.mu; % 摩擦系数 par.E = 2.06e5; % 钢弹性模量，MPa par.v = 0.3; % 泊松比 end

这种加载方式的好处是领域工程师能直接看懂，不需要额外熟悉数据格式。参数方面有几个容易忽略的点：Rw 要取当前接触点处的滚动圆半径，而不是名义轮径，轮缘贴靠时滚动圆半径会明显变化；Rr 是钢轨横向曲率半径，轨顶打磨后实际曲率半径和标准廓形并不一样，最好从实测廓形反算。蠕滑率的单位也需要提前说清楚：纵向和横向蠕滑率无量纲，自旋蠕滑率应该带 1/mm 或 1/m，两个单位制都会出现在论文里，但程序内部必须统一成一种，我习惯用 1/mm，因为廓形坐标都是毫米。

3.2 Hertz接触参数计算：确定接触斑长短半轴并校验压力分布

赫兹接触在轮轨问题里的作用是给出一个等效接触椭圆。长半轴 a 沿滚动方向，短半轴 b 沿横向，最大法向压力 p0 按椭圆抛物面分布。工程实现上，最稳妥的是按曲率差查赫兹系数表，但表格式数据不利于批量编程。我用的是简化迭代法：先假设一个椭圆偏心率，计算辅助函数，再更新长半轴和短半轴，直到前后两次偏心率之差小于收敛阈值。

function [a, b, delta, p0] = hertz_contact(Fn, Rw, Rr, E, v) % 赫兹法向接触参数计算：轮轨正交圆柱接触简化模型 % Fn : 法向力 (N) % Rw : 接触点处滚动圆半径 (mm) % Rr : 轨顶横向曲率半径 (mm) % E : 弹性模量 (MPa) % v : 泊松比 Eprime = E / (1 - v^2); % 等效弹性模量，平面应变假设 SigmaRho = 1/Rw + 1/Rr; % 主曲率之和，1/mm % 两圆柱正交接触下偏心率随曲率差变化，这里采用迭代求解 e = 0.8; % 偏心率初始猜测 for iter = 1:50 [F_e, E_e] = ellipke(e^2); % 第一类、第二类完全椭圆积分 factor = F_e / E_e; a = (3 * Fn / (4 * Eprime * SigmaRho) * factor)^(1/3); b = a * sqrt(1 - e^2); e_new = sqrt(1 - (b/a)^2); if abs(e_new - e) < 1e-4 e = e_new; break; end e = e_new; end delta = Fn / (Eprime * (pi * a * b / 2)); % 法向接近量近似 p0 = 1.5 * Fn / (pi * a * b); % 最大法向压力，MPa end

这里要说明两点。第一，ellipke 在 MATLAB 中直接用数值积分结果，不需要你自己查表，但要注意参数是 e^2，不是偏心率本身，传错参数会让迭代发散。第二，delta 的表达式只适用于椭圆接触的简化估算，如果后续要算接触力做功或接触刚度，建议用更精细的数值积分；如果只是把 delta 作为输出参数记录，这个近似完全够用。实际使用中，我通常会在函数外部加一个校验逻辑：接触椭圆面积乘以平均压力应该近似等于 Fn，如果偏差超过 5%，说明曲率输入或单位制有误，程序应提前报警。

3.3 FASTSIM切向蠕滑力计算：条带网格上的剪应力积分

FASTSIM 的离散实现是整套程序里最容易写错的一环。核心思想是把接触斑沿滚动方向切成若干条带，从接触入口开始，逐条带累积弹性位移，再由位移计算局部剪应力。剪应力超过摩擦极限的点进入滑动区，剪应力被截断到摩擦锥边界；未超过的点属于粘着区。整个过程可以理解成一种逐行递推的“积分器”。

function [Fx, Fy, Mz] = fastsim(a, b, p0, xi, eta, spin, mu, Lx, Ly, nx, ny) % FASTSIM 简化蠕滑力计算 % a/b : 接触斑长半轴(滚动方向)与短半轴(横向), mm % p0 : Hertz最大法向压力, MPa % xi/eta/spin : 纵向、横向、自旋蠕滑量 % Lx/Ly : 柔度系数，控制切向力-滑移关系的刚度 % nx/ny : 滚动方向与横向网格数 % 条带划分 dx = 2*a / nx; dy = 2*b / ny; Fx = 0; Fy = 0; Mz = 0; for ky = 1:ny yv = -b + (ky - 0.5) * dy; % 条带中心横向坐标 % 每个条带内的纵向递推 sx = 0; sy = 0; % 弹性位移累积量，初始为0 for kx = 1:nx xv = -a + (kx - 0.5) * dx; % 条带中心纵向坐标，入口到出口 % 法向压力椭圆分布 pz = p0 * sqrt(1 - (xv/a)^2) * sqrt(1 - (yv/b)^2); pz = max(pz, 0); % 自旋项对局部蠕滑的贡献，符号按车辆前进方向右手系 wxi = xi + spin * yv; weta = eta - spin * xv; % 弹性位移增量，柔度系数把应变转换成应力 sx = sx + wxi * dx / Lx; sy = sy + weta * dx / Ly; px = -sx; py = -sy; pt = sqrt(px^2 + py^2); % 摩擦极限截断 limit = mu * pz; if pt > limit && pt > 0 scale = limit / pt; px = px * scale; py = py * scale; % 滑动区内重置位移，避免误差往下游传递 sx = -px; sy = -py; end Fx = Fx + px * dx * dy; Fy = Fy + py * dx * dy; Mz = Mz + (xv * py - yv * px) * dx * dy; end end

这段代码的关键在递推顺序：纵向从接触入口（xv 为负）向出口（xv 为正）遍历，弹性位移逐步累积，这样粘着区到滑动区的过渡是自然发生的。柔度系数 Lx、Ly 实际控制着切向力-滑移曲线的初始斜率和饱和点，取值过大会让蠕滑力过早饱和，过小则让粘着区范围异常偏大。规范做法是从半空间弹性模型推导一组初值，再用纯纵向蠕滑工况标定一次，修正系数在 0.8 到 1.2 之间。网格数 nx 和 ny 取 40 和 20 就能得到稳定结果，继续加密对总力影响很小，只会增加耗时。自旋项对高速转向架尤其重要，忽略自旋会让摇头频率附近的稳定性判断明显偏乐观。

3.4 主流程封装：一次横移量扫描的批处理入口

有了以上三个核心函数，最后还需要一个批量扫参的入口。最典型的工程场景是固定摇头角和蠕滑率，扫一组横移量，观察接触点位置和轮轨横向力变化趋势。主流程里先预分配结果数组，再对每个横移量依次调用接触几何、赫兹和 FASTSIM。

y_shift_list = -12:0.5:12; % 横移量扫描范围，单位 mm n = length(y_shift_list); result = struct('y_shift', [], 'Fx', [], 'Fy', [], 'Mz', [], 'a', [], 'b', []); for k = 1:n c = find_contact_point(rail_y, rail_z, wheel_y, wheel_z, y_shift_list(k), 0); [a, b, ~, p0] = hertz_contact(Fn, Rw, Rr, E, v); [Fx, Fy, Mz] = fastsim(a, b, p0, xi, eta, spin, mu, Lx, Ly, 40, 20); result(k).y_shift = y_shift_list(k); result(k).Fx = Fx; result(k).Fy = Fy; result(k).Mz = Mz; result(k).a = a; result(k).b = b; end

这里要注意接触几何函数里的 yq 密度要和扫描步长匹配，如果扫描步长是 0.5 mm，而接触点搜索网格稀疏到 5 mm 一档，接触点位置会呈现阶梯状跳变。工程上 0.5 mm 步长配合 2000 个插值点是比较稳的组合，既能看出轮缘贴靠的突变，又不至于让计算量失控。

4. 轮轨计算高频翻车点：坐标、单位与网格条带的避坑排查

代码能跑通和结果能用是两码事。我自己在轮轨关系项目里踩过不少坑，有些问题看起来是数值误差，实际是坐标约定或单位制错位。这一章把高频翻车点按“现象、原因、解决”写清楚，每条都是可以直接对照排查的经验。

4.1 蠕滑力与位移方向相反：坐标系约定不一致

现象：纯纵向蠕滑工况下，算出来的纵向力 Fx 和蠕滑率 xi 符号正好相反，而且横向力 Fy 也不为零，像是程序内部发生了“串扰”。

原因：轮对前进方向、纵向坐标轴和 FASTSIM 内部的 x 轴约定不一致。有的文献把接触斑前进方向定义为正 x，有的则定义车辆前进方向为负 x；如果几何函数里用摇头角旋转，而 FASTSIM 里又直接用原始蠕滑率，符号就会错位。

解决：在程序入口处统一约定：纵向 x 正向取车辆前进方向，垂向 z 向上，横向 y 按右手定则确定。所有输入蠕滑率、接触点坐标和自旋项符号都基于这个基准。加入一个自检函数，输入纯纵向蠕滑率，输出 Fx 应该与 xi 同号，否则立即报错。

提示：这个自检我固定在每个工程项目的入口脚本里，运行任何工况前都先跑一遍，几十毫秒的成本换一次方向翻车，非常划算。

4.2 接触椭圆参数反复跳变：曲率半径正负号没处理

现象：横移量扫描时，赫兹接触的 a、b 值在某几个位置突然增大或减小，出现明显跳变，甚至 a 比 b 还大，明显不符合轮轨接触常识。

原因：轮轨接触中，车轮踏面在接触点处是凸面，轨顶在局部也是凸面，两者曲率方向相同；但一旦进入轮缘贴靠区，接触点移到轮缘圆弧与轨距角附近，曲率半径的正负号会变化。如果程序里把所有曲率半径都按正值处理，接触几何就会失真。

解决：为曲率半径定义符号协议：凸面接触取正，凹面接触取负。轮轨接触几何中，接触点处的等效曲率要和接触力方向对应。建议在计算 Hertz 参数前，先检查接触点处踏面弧段的凹凸性，用一个前后点差分判断曲率符号，而不是直接采用名义半径。

4.3 FASTSIM 摩擦锥截断后总切向力越界

现象：程序输出的 Fx 和 Fy 合成后大于 μFn，从物理上看已经突破库仑摩擦极限，说明滑动区处理有漏洞。

原因：常见原因是接触斑网格划分太粗，法向压力 pz 在条带中心取值后，边缘处 pz 很小，但切向应力累积值偏大，截断后被重置的弹性位移没有和网格尺寸匹配，导致错误累积到相邻条带。

解决：把 nx 从 20 提高到 50 以上，同时把滑动区重置逻辑改成“重新计算当前条带内的等效位移”，而不是简单截断应力。另一个实用检查是在程序结果里附上总接触面积、粘着区域面积和滑动区域面积，三者之和应等于理论接触斑面积；如果相差超过 2%，就说明网格步长太大或边界处理有问题。

4.4 接触几何搜索无解：廓形坐标原点和长度不一致

现象：find_contact_point 函数返回的 gap 全部是 NaN，或者接触点横坐标跑到廓形端部之外。

原因：踏面廓形和轨面廓形来自不同设备或不同导出批次，横向坐标零点没有对齐；或者廓形点长度差得太多，统一插值后有效区间很小，导致搜索区间内没有有效点。

解决：在加载廓形数据时做强制归一化处理：把轨面顶部顶点和踏面基准点对齐到同一横向坐标零点，然后对两组廓形在公共区间内重新插值到等距网格。我一般会把对齐逻辑写成一个单独的函数，前端加载一步到位，后续几何计算不再碰原始坐标。

4.5 批处理内存爆炸：中间应力场保存过多

现象：扫 500 个横移量工况时，MATLAB 内存涨到几个 GB，运行越来越慢，最后强制关进程。

原因：循环里把 FASTSIM 每个工况的中间剪应力场全部存进了结构体，而实际工程分析只需要总力和总力矩，应力场仅用于少数调试场景。

解决：默认只保存 Fx、Fy、Mz 和 a、b、接触点坐标；调试时单独设置一个 debug 开关，只有打开时才保存应力场矩阵。循环前用 zeros 预分配所有结果数组，尽量避免结构体动态扩张。如果工况数上万，建议把 parfor 打开，并把廓形数据作为切片变量传入，减少并行进程间的数据复制。

5. 把计算结果用起来：等效锥度提取与多体仿真数据对接

轮轨力算出来只是中间产物，工程上更关心的往往是两个输出：一是接触几何随横移量的变化规律，二是可供整车动力学模型直接调用的轮轨力谱。前者通常用等效锥度来表征，后者需要把计算结果整理成规范的时序数据。最后一章说清楚这两件事怎么做。

等效锥度是衡量轮对自导向能力最直观的参数，定义上近似等于滚动圆半径差随横移量的变化率除以名义滚动圆半径。实际计算时，不需要直接测量滚动圆半径，只需要在接触几何求解里记录每个横移量下左右车轮接触点的滚动圆半径差。

% 从扫掠结果计算等效锥度 % result(i).y_shift 是横移量，rr(i)/rl(i) 分别记录左右接触滚动圆半径 y_shift = [result.y_shift]; delta_r = [result.rr] - [result.rl]; % 滚动圆半径差 idx = abs(y_shift) < 3; % 通常在±3mm内取线性段 p = polyfit(y_shift(idx), delta_r(idx), 1); lambda = p(1) / (2 * Rw_nom); % 等效锥度，无量纲

这段代码的思路是线性拟合线性段斜率，再除以两倍名义滚动圆半径。为什么是两倍？因为左右轮半径差的斜率和轮对横移量之间有一个几何放大关系，定义式中分母取 2R。注意 polyfit 得到的是一次项系数，单位是 mm/mm，和 Rw_nom 单位必须一致；如果 Rw_nom 用米而 y_shift 用毫米，锥度结果会差 1000 倍。这个坑我见不止一次，强烈建议在代码里强制统一。

轮轨力谱导出到多体动力学软件时，数据结构比数据本身更重要。整车模型通常按时间步或位置步调用轮轨力，程序里应把扫描工况整理成一张长表，每列是一个可识别的物理量，并带上单位信息。

% 导出为 CSV 长表格式 T = table(y_shift', Fx_list', Fy_list', Mz_list', ... 'VariableNames', {'y_shift_mm', 'Fx_N', 'Fy_N', 'Mz_Nm'}); writetable(T, 'wheel_rail_force_table.csv');

导出文件后，我一般会在多体软件里做一个单轮对静态校核：把程序输出的 Fy 随横移量变化曲线和已知轮缘接触试验曲线对比，趋势一致再往下走。第一次做轮轨关系项目时，我就是因为接触几何搜索函数里插值网格太稀疏，导致轮缘贴靠点偏移了 2 mm，整个横向稳定性判据失真，整车模型白跑了两周。从那以后，我每次跑轮轨程序都强制走一遍坐标单位自检、纯纵向蠕滑符号自检和接触点曲线平滑度检查，三步缺一不可。这套流程虽然多花十分钟，但能挡掉大部分低级错误，希望帮到你。

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