简介一套 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单工况耗时秒级到分钟级毫秒级接触斑内应力细节完整简化抛物线分布与小蠕滑试验对比误差基准约 10%20%代码可读性低数百行矩阵积分高一个循环即可实现适合场景研究性分析、验证基准批量扫参、实时仿真、工程回归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整个横向稳定性判据失真整车模型白跑了两周。从那以后我每次跑轮轨程序都强制走一遍坐标单位自检、纯纵向蠕滑符号自检和接触点曲线平滑度检查三步缺一不可。这套流程虽然多花十分钟但能挡掉大部分低级错误希望帮到你。本文还有配套的精品资源点击获取