企业官网建设流程全解析

轮轨关系计算入门这件事，说穿了就两个核心：一个是 Hertz 接触给出的法向压力分布，另一个是 FASTSIM 算出来的切向蠕滑力。前者决定接触斑有多大、压得多狠，后者决定轮对在钢轨上滚动时到底被"拽"了多少、会不会打滑。这两块拼起来，才是一套能用的轮轨关系计算流程。我最早接触这块内容是在做车辆动力学仿真的时候，当时对着公式推导看了一个礼拜，真正动手用 MATLAB 写出来才发现，纸上看着顺的公式，一进代码全是符号和边界条件的坑。这篇文章就把这条从 Hertz 到 FASTSIM 的完整链路拆开讲，包含曲率半径怎么取、椭圆积分怎么解、条带怎么推进、黏滑怎么判，最后给一套能直接跑通的 MATLAB 代码。不管你是刚进这个方向的学生，还是做了几年仿真想回头补基本功的工程师，应该都能从里面捞出点能用的东西。

1. 先想清楚：轮轨关系计算到底在算什么

1.1 一个接触斑，两个性质完全不同的问题

车轮压在钢轨上，接触区域小得可怜，典型情况下长轴十几毫米、短轴六七毫米，面积也就几十平方毫米。但就是这几十平方毫米，要传递整辆车的重量，还要传递牵引、制动和导向产生的切向力。这里有个新手特别容易混淆的点：法向和切向虽然发生在同一个接触斑上，但它们的求解逻辑完全是两套东西。

法向问题相对"干净"。轮轨材料都是弹性体，接触区尺寸远小于构件曲率半径，满足 Hertz 假设的前提，所以可以用 Hertz 接触理论解析地算出接触斑的椭圆形状、半轴长度和压力分布。它本质上一个几何加弹性力学的封闭问题，输入是曲率半径和法向力，输出是压力场，中间不需要迭代求解历史状态。

切向问题就麻烦得多。轮轨之间存在相对滑动趋势，也就是蠕滑。蠕滑率很小的时候，接触斑里一部分区域黏着不动、一部分区域发生微滑，黏着区和滑动区的边界还依赖于加载历史。这种"有记忆"的特性，让切向问题必须沿着滚动方向一步步推进求解，这也是 FASTSIM 这类算法的由来。我在实际项目里见过有人拿 Hertz 压力乘个摩擦系数就当切向力用，短蠕滑工况下误差还能忍，一旦进入饱和区，误差能到百分之几十。

所以这一整套计算的标准顺序是：先用 Hertz 把法向边界条件定下来（接触斑尺寸、压力分布），再把这套几何喂给 FASTSIM 算切向力。顺序不能反，因为 FASTSIM 的每一条计算条带都需要知道当地的极限摩擦力，而极限摩擦力等于摩擦系数乘以当地法向压力。

1.2 为什么我不建议一上来就调商业软件

商业多体动力学软件里，轮轨接触模块通常是黑箱，点几下鼠标就出蠕滑力曲线。用是能用，但一旦结果不对劲，你根本不知道从哪儿查。我踩过最典型的一次坑：某次仿真里横向力明显偏大，查了两天才发现是踏面横曲率半径填成了直径值，接触斑短轴直接缩了一半，接触压力翻倍，切向力自然跟着虚高。

自己用 MATLAB 把这套流程写一遍，最大的收益不是代码本身，而是你会清楚每一个输入参数从哪来、影响到哪一步、量级大概是多少。写完一次之后再看任何软件的接触设置界面，心里都有底。而且 MATLAB 的调试体验在矩阵运算类问题上确实顺手，中间变量随便看，画个图验证一下分布形态也很方便。

当然，自研代码有边界。它适合做参数敏感性分析、算法理解、教学演示，真要算整车曲线通过性能，还是得靠成熟的动力学平台。我的建议是：先用 MATLAB 把单点接触算明白，再回到软件里用同样的参数核对一遍，两边对得上，你才算真的掌握了这个模块。

1.3 MATLAB 在这套计算里的定位

具体到工具层面，这套计算对 MATLAB 的依赖其实很轻：一个椭圆积分的内置函数、一个一维求根、再加上几个矩阵循环就够了。不需要任何工具箱，基础版本就能跑。椭圆积分用ellipke，求根用fzero，这两个是核心。

真正需要花心思的是怎么组织代码结构。我的习惯是拆成三个文件：一个主脚本负责参数设定和结果绘图，一个hertz_contact函数专门处理法向，一个fastsim函数处理切向。这样你在调切向参数的时候，法向那块不用重复计算，也方便单独验证每一部分的正确性。后面第 4 节会把这几个文件完整给出来。

还有一个实操层面的建议：所有长度单位统一用米，力用牛顿，弹性模量用帕斯卡。轮轨这块参数跨度很大，滚动圆半径是零点几米，接触斑半轴是零点零几米，接触压力是几百兆帕。单位一旦混用，比如曲率半径填毫米、力填千牛，结果会错得离谱而且很难从数量级上看出来。我自己的习惯是在每个函数开头写一行注释，标明输入输出的单位，看着啰嗦，救过我好几次。

2. Hertz 接触：法向参数的计算细节

2.1 曲率半径怎么取，符号怎么定

Hertz 接触的输入是接触点处两个物体的主曲率半径。轮轨接触这个场景下，通常简化成四个曲率：车轮的纵向曲率、车轮的横向曲率、钢轨的纵向曲率、钢轨的横向曲率。

车轮纵向曲率就是滚动圆半径，新轮 0.42 米左右，磨耗到限大概 0.38 米，这个值好取。车轮横向曲率指的是踏面在接触点附近的横向圆弧半径，这个就麻烦了：踏面是变曲率的，锥形踏面名义上是个直线段，实际运行后磨耗成凹形，接触点位置不同，局部曲率半径差很多。工程上一般取接触点附近拟合出来的等效半径，量级在 0.15 到 0.5 米之间。

钢轨纵向曲率通常直接忽略，因为钢轨沿线路方向的曲率半径是几百米量级，对接触斑的影响可以忽略不计。钢轨横向曲率取轨顶圆弧半径，新轨 300 毫米，磨耗后可能降到 80 到 150 毫米，这个变化对接触斑形状影响非常大，是现场磨耗分析里最敏感的参数之一。

符号约定是最容易出错的地方。我用的是这套规则：车轮的两个主曲率都取正，因为踏面是凸的；钢轨横向曲率取负，因为钢轨轨顶也是凸的，但从车轮的角度看过去，它的弯曲方向和车轮相反。这个符号一旦搞反，A 和 B 会直接变号，后面求出来的接触斑形状就完全不对。

提示：如果你的计算结果里接触斑长轴跑到了横向，而实际应该是纵向，八成就是钢轨横向曲率的符号写反了。先改符号，再怀疑别的。

2.2 A、B 两个曲率和怎么合成

有了四个曲率，按下面的方式合成两个中间量。设车轮纵向曲率为 rho11、横向为 rho12，钢轨纵向为 rho21、横向为 rho22，则

• A = (rho11 + rho12 + rho21 + rho22) / 2
• B = 0.5 × |(rho11 - rho12) + (rho21 - rho22)|

这里假设两个物体的主曲率平面重合，也就是夹角为零。轮轨在轨顶中心附近接触时这个假设基本成立，横移量大的时候会有接触角，严格算需要引入夹角项，但入门阶段先按零处理完全够用。

算出来的 A 和 B 要保证 A 大于等于 B，如果反了就交换一下。这一步看着是技术细节，实际上是给后面的椭圆积分定方向：A 对应接触斑短轴方向，B 对应长轴方向。交换之后，你算出来的 m 和 n 才是配套的。

拿一组典型参数试一下：车轮滚动圆半径 0.42 米，踏面横曲率半径 0.20 米，轨顶横曲率半径 0.30 米。四个曲率是 2.381、5.000、0、-3.333，单位是每米。A 等于 2.024，B 等于 0.357。A 比 B 大约 5.7 倍，说明接触斑会是一个明显的长椭圆，长轴沿纵向。这个比值直接决定了后面 m、n 系数的大小，也决定了接触斑有多扁。

2.3 椭圆积分求 m、n 的数值做法

Hertz 接触的半轴公式是

• a = m × (3P / (2E* × Σρ))^(1/3)
• b = n × (3P / (2E* × Σρ))^(1/3)

其中 Σρ 等于 A 加 B，E* 是等效弹性模量，两个物体材料相同时等于 E / (2(1-ν²))。m 和 n 是椭圆积分的函数，只跟 A/B 这个比值有关，跟载荷和材料都无关。

求 m、n 的标准路径是先解一个椭圆离心率 e。定义

• ratio(e) = [E(e)/(1-e²) - K(e)] / [K(e) - E(e)]

这个 ratio 随 e 从 0 增到 1 时，从 1 单调增到无穷。所以给定 A/B 之后，用一维求根把 e 反解出来就行。MATLAB 里 K 和 E 都用ellipke，注意它的入参是模数的平方，也就是 e²，不是 e，这一点特别容易搞错，我第一次写的时候直接把 e 传进去，结果接触斑尺寸偏了百分之十几。

拿到 e 之后：

• m = [2E(e) / (π(1-e²))]^(1/3)
• n = [2E(e)√(1-e²) / π]^(1/3)

用这组公式验证一个特例：A/B 等于 1 时，e 等于 0，K 和 E 都等于 π/2，代进去 m 和 n 都等于 1。这时候 a 等于 b，接触斑退化成圆，公式回到经典的球-平面接触解，对得上，说明公式没问题。这种"拿已知特例反查"的习惯建议你保留，比单纯背公式靠谱得多。

求根的时候有个数值细节要注意：e 太接近 1 时ellipke会掉精度，所以求根区间上界建议卡在 0.9999，不要给 1。另外 e 接近 0 时 ratio 是 0 除 0，函数里要