企业官网建设流程全解析

轮轨关系这四个字，写在纸上是一组教科书公式，落到代码里就是几百行的数值循环加一堆容易搞错的约定。很多做车辆动力学的人，日常工作是调用现成的多体软件，轮轨力那部分是个黑箱——输入轴重和蠕滑率，输出三个方向的力和力矩，曲线平滑、看着很舒服。可一旦需要自己写一个轮轨力模型接进自研求解器，或者只是想搞清楚那条蠕滑力曲线为什么在某个点突然拐弯，就会发现从 Hertz 接触理论到 FASTSIM 蠕滑力求解这一整条链路，中间隔着不少没人明说的细节。

这篇内容就是把这根链条完整拆一遍：法向接触怎么从两个曲面折算成一个等效问题、Hertz 的椭圆接触斑尺寸怎么用 MATLAB 算出来、Kalker 简化理论为什么要把接触斑切成网格沿滚动方向一条条扫、柔度系数和自旋项到底在算什么、最后蠕滑力和自旋力矩怎么积分出来。读者对象是学过材料力学和高等数学、能看懂微分方程但没真正手写过轮轨力代码的工程人员，也包括想复现论文里蠕滑力曲线验证自己算法的研究生。全部代码用 MATLAB，不依赖任何工具箱，复制粘贴就能跑。

我踩过的坑基本都在“公式都对，结果不对”这个区间里——单位、符号约定、归一化方式、系数表的版本，随便哪一处出问题，曲线形状看着还挺像，数值能差出两三倍。所以下面除了讲怎么算，更多篇幅会放在为什么这么算、哪里最容易错、怎么用已知解做反向校核上。

1. 一条从轴重到蠕滑力的计算链，卡在哪一环

1.1 商业软件里那个“黑箱”到底是什么

把轮轨力拆开看，它其实是两个物理上完全不同的子问题叠在一起。法向方向决定接触斑有多大、压力怎么分布，这是一个纯几何加弹性力学的问题，输入是轮轨型面、轴重、接触角，输出是接触椭圆半轴、峰值压力和法向接近量。切向方向决定轮对相对钢轨有微小相对运动时，接触斑里能传递多少纵向力、横向力和自旋力矩，输入是蠕滑率和摩擦系数，输出是三个力分量。

这两件事的数学性质差得很远。法向问题是椭圆积分，有闭式的近似解，算一次不到一毫秒；切向问题本质上是接触斑内的一个二维积分方程，严格求解需要解 Cauchy 型奇异积分方程，慢到没法放进实时仿真。所以工程上清一色走“法向用 Hertz 精确解 + 切向用简化理论”这条路线，FASTSIM 就是后者里最经典的一个实现。

理解了这一点，就能明白为什么很多自研模型在法向算得很准，切向却怎么调都对不上——法向只要几何和材料参数对了，结果基本不会错；切向的误差来源是模型假设本身加上离散化误差，需要单独处理。

1.2 Hertz 能解决什么，解决不了什么

Hertz 接触理论的前提假设有四条：两个接触体都是各向同性线弹性材料、接触区尺寸远小于物体特征尺寸、接触面光滑无摩擦、法向压力垂直于接触面。放到轮轨场景里逐条对照：材料假设基本成立，钢的线弹性范围足够覆盖轮轨的工作应力；尺寸假设在轮对踏面这种大曲率变化的地方会有点问题，但对常规接触斑尺寸（十毫米量级）仍然可用；光滑假设在打磨后的新轨上成立，磨损后期表面粗糙度上升会引入偏差；法向压力垂直这条，在轮对横移量较大、接触角不为零时需要用接触角做投影修正。

真正容易被忽略的是“无摩擦”这条。Hertz 解给出的压力分布 p(x,y) = p0·√(1 − x²/a² − y²/b²) 是在接触面完全光滑、没有切向力的情况下推出来的。实际有切向力的时候，压力分布会略微偏移，压力中心往滚动方向前移一点点。FASTSIM 这类简化理论直接沿用 Hertz 的压力分布，把它当作已知量喂给切向求解，这是工程上普遍接受的近似，但你得知道这个近似存在。

还有一条：Hertz 只解决单点接触。轮对在轨道上横移时可能同时出现踏面接触和轮缘接触，也就是所谓两点接触。这种情况不能简单地把两个 Hertz 解加起来，因为两个接触区的法向力分配取决于轮对的平衡位置，需要联立求解。刚开始做的时候建议先把工况限制在单点接触、直线轨道、无冲角，把这条链路跑通再说。

1.3 从法向到切向：为什么必须引入简化理论

假设你现在已经拿到了接触椭圆半轴 a、b 和峰值压力 p0，也知道轮对相对钢轨的蠕滑率。接下来要算切向力，最直接的想法是：接触斑里每个点都满足库仑摩擦，粘着区不滑、滑移区滑到底。这个思路叫“刚性接触模型”，它的问题是——在小蠕滑时给出的切向力是跳变的，要么全粘着给不出力，要么一滑就跳到摩擦极限，中间那段连续的过渡完全丢失。

实际轮轨的蠕滑力曲线在小蠕滑区是近似线性的，斜率可观，然后在 1% 到 2% 的蠕滑率附近逐渐饱和到 μN。这条曲线的形状决定了车辆系统的稳定性：斜率太陡，转向架容易蛇行失稳；饱和太快，曲线通过时导向力不够。所以必须有一个能描述“粘着区逐渐缩小、滑移区逐渐扩大”这个连续过程的模型。

Kalker 的贡献在于用一组变分不等式严格描述了这个问题，然后 FASTSIM 用一种工程上极其讨巧的方式做了离散化：接触斑切成网格，沿着滚动方向一条条往前走，每走一格就把切向应力的增量累加一次，同时检查有没有超过摩擦极限。这个算法把所有非线性都收敛到了一行判断语句里，跑一次二维网格循环就能出结果，速度比严格解快两三个数量级，精度对工程应用来说足够。

2. Hertz 接触的 MATLAB 落地：曲率折算与椭圆系数

2.1 把车轮和钢轨折算成一个等效弹性体

两个弹性体接触，直接处理需要同时跟踪两个物体的变形，很麻烦。经典处理是把其中一个当作刚体，把两个物体的弹性常数和曲率都折算到另一个上，变成一个等效弹性体压在刚性平面上。折算规则是：

曲率方面，纵向和横向各定义一个等效曲率：

• A = ½·(1/R₁ₓ + 1/R₂ₓ)
• B = ½·(1/R₁ᵧ + 1/R₂ᵧ)

材料方面，等效弹性模量：

• 1/E* = (1 − ν₁²)/E₁ + (1 − ν₂²)/E₂

轮轨都是钢，E = 206 GPa，ν = 0.3，所以 1/E* = 2×0.91/206e9，得到 E* ≈ 113.2 GPa。这里有个特别容易错的地方：很多人直接用单体的 E 或者用 E/(1−ν²)，那是“刚体压弹性半空间”的公式，用在轮轨这种两个弹性体接触的场景会大出去一倍。

曲率的符号约定同样容易翻车。车轮踏面横断面是凸的，钢轨轨顶也是凸的，两个凸面接触时曲率半径取正值，等效曲率是两个正值相加的一半。如果遇到轮缘根部与钢轨侧面接触这类“凸对凹”的情形，凹面的曲率半径要取负值。直线轨道上钢轨的纵向曲率半径是无穷大，1/R = 0，这一项直接不参与计算。

举个具体算例：车轮滚动圆半径 0.42 m，踏面横断面半径取 0.5 m，钢轨轨顶半径 0.3 m，单轮法向载荷 80 kN。代入得 A = ½×2.381 = 1.190，B = ½×(2.0 + 3.333) = 2.667，A/B = 0.446。这个比值决定了后面椭圆系数怎么查。

2.2 椭圆积分系数 m、n 的求解与查表的取舍

Hertz 解的最终形式里有两个无量纲系数 m 和 n，它们是 A/B 的函数：

• a = m·[3P/(2E*(A+B))]^(1/3)
• b = n·[3P/(2E*(A+B))]^(1/3)

这两个系数的严格求法涉及第一类和第二类完全椭圆积分，需要数值迭代。工程上最常用的做法是查表加插值，下面这张表是我用了很多年的一组近似值，来源是经典弹性力学手册，用来做工程计算足够：

注意：不同手册给出的 m、n 定义并不完全一致，有的把系数定义在 (3PR/4E*)^(1/3) 这个形式里，有的定义在 (3P/2E*Σρ)^(1/3) 里。混用两套表会得到差一个 2^(1/3) 倍的结果，看起来还在合理范围内，非常隐蔽。建议选定一套之后，用球-球接触的解析解做一次校核，确认系数体系自洽。

自洽性检验很简单：取 A = B（球对球，或者两个相同半径的圆柱正交），此时 A/B = 1，表里 m = n = 1。代入得 a = b = [3P/(2E*·2/R)]^(1/3) = [3PR/(4E*)]^(1/3)，正好是标准 Hertz 球接触公式。这一步对上了，后面就不会有大问题。

查表插值的精度问题：A/B 在 0.2 到 1.0 之间时，用线性插值误差在 1% 到 2%，用三次样条能压到 0.3% 以内。如果 A/B 小于 0.1，也就是接触斑非常狭长，这张表就不够用了，得换成椭圆积分直接数值求解。轮轨常规工况下 A/B 一般在 0.2 到 0.8，表够用。

2.3 一个能直接跑的 hertz_solve 函数

把上面的逻辑写成一个函数，输入载荷、四个曲率半径、两个材料的弹性常数，输出半轴、峰值压力和法向接近量：

function out = hertz_solve(P, R1x, R1y, R2x, R2y, E1, nu1, E2, nu2) % Hertz 椭圆接触求解 % P 为法向载荷(N)，R 为曲率半径(m)，R 取 inf 表示该方向无曲率 % E 单位 Pa % 等效曲率 kx = 0; if isfinite(R1x), kx = kx + 1/R1x; end if isfinite(R2x), kx = kx + 1/R2x; end ky = 0; if isfinite(R1y), ky = ky + 1/R1y; end if isfinite(R2y), ky = ky + 1/R2y; end A = 0.5 * kx; B = 0.5 * ky; % 等效弹性模量 Estar = 1 / ((1 - nu1^2)/E1 + (1 - nu2^2)/E2); % 椭圆系数查表插值 ratio = A / B; ratioTab = [0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.0]; mTab = [1.884 1.546 1.361 1.245 1.166 1.109 1.038 1.000]; nTab = [0.319 0.451 0.545 0.622 0.690 0.753 0.872 1.000]; if ratio < 0.1 || ratio > 1.0 error('A/B = %.3f 超出查表范围，需改用椭圆积分数值求解', ratio); end m = interp1(ratioTab, mTab, ratio, 'pchip'); n = interp1(ratioTab, nTab, ratio, 'pchip'); % 接触半轴 c = (3 * P / (2 * Estar * (A + B)))^(1/3); a = m * c; b = n * c; % 峰值压力与接近量 p0 = 3 * P / (2 * pi * a * b); out.a = a; out.b = b; out.p0 = p0; out.A = A; out.B = B; out.Estar = Estar; out.m = m; out.n = n; end

用前面那组数据跑一下：P = 80000，R1x = 0.42，R1y = 0.5，R2x = inf，R2y = 0.3，E = 206e9，ν = 0.3。得到 A = 1.190，B = 2.667，A/B = 0.446，插值出 m ≈ 1.25、n ≈ 0.66，c ≈ 6.50 mm，最终 a ≈ 8.13 mm、b ≈ 4.29 mm、p0 ≈ 1.10 GPa。

这组数跟文献里的轮轨接触实测结果量级一致：接触斑长度十毫米以内、宽度四到五毫米、峰值压力一个 GPa 上下。如果算出来 a 是 80 mm 或者 p0 是 100 MPa，基本可以确定是某一项的单位搞错了。

2.4 用球形接触做反向校核

写完函数别急着往轮轨上套，先用一个能解析验证的算例跑一下。两个半径 20 mm 的钢球互相压，载荷 1000 N，按解析式 a = b = [3PR/(4E*)]^(1/3)：

3×1000×0.02 = 60，4×113.2e9 = 4.528e11，比值 1.325e-10，开三次方得 5.10e-4 m，也就是 0.51 mm。峰值压力 p0 = 3P/(2πa²) = 3000/(2π×2.6e-7) = 1.84e9 Pa。

把参数喂进 hertz_solve，R1x = R1y = R2x = R2y = 0.02，得到的 a、b 应该都是 0.51 mm，p0 一致。差超过 1% 就回头查等效模量和曲率的正负号。

这个校核花不了五分钟，但能省掉后面几个小时的排查。我见过太多次因为 E* 用了单体的 206 GPa 而不是 113 GPa，接触斑尺寸算小了四分之一，然后切向力整体偏小，最后花了一整天才定位到问题在法向这一层。

3. FASTSIM 的核心：沿着滚动方向逐条扫描

3.1 Kalker 简化理论究竟简化掉了什么

严格的三维弹性接触问题是这样的：接触斑内每一点的切向位移 u(x,y) 都与整个接触区内的切向应力分布 q(x',y') 通过一个积分核关联，这是一个二维积分方程。同时每个点还要满足库仑摩擦约束：粘着区 |q| < μp 且局部滑动为零，滑移区 |q| = μp 且滑动方向与切向力方向相反，两种区域的边界事先未知，需要作为解的一部分求出来。

Kalker 的简化理论砍掉了积分核的长程效应，假设接触斑内某一点的切向位移只由同一条“流线”上它前面那些点的应力决定。物理上的解释是：车轮滚过去的时候，接触斑里的材料是连续进入和离开的，材料点的变形历史沿着滚动方向累积，横向相邻点之间的耦合相对次要。这个假设让二维积分方程退化成一维递推，计算量从 O(N⁴) 降到 O(N²)。

这个简化对长宽比不算极端的接触椭圆精度很好，纵向力和横向力的误差通常在 5% 到 10% 以内，自旋力矩的误差会大一些。对于车辆动力学仿真来说这个精度完全够，因为我们本来就在跟一大堆参数不确定性打交道。

3.2 柔度系数 L1、L2、L3 的物理意义

简化理论里引入了三个柔度系数，物理意义是“单位切向应力引起的切向位移”：

• L1 对应纵向（滚动方向）
• L2 对应横向
• L3 对应自旋

它们的表达式形式是 L1 = 8a/(3G·c11)，其中 a 是纵向半轴，G 是剪切模量，c11 是 Kalker 系数，取决于接触椭圆的长宽比和泊松比。c11、c22、c23 这三个系数 Kalker 给出了表格，工程实现里通常做成查表插值。

注意：这三个系数表在文献里有好几套版本，有的按 a/b 排、有的按 b/a 排，有的把泊松比固定成 0.3 只给一列。混用会直接导致柔度系数差一个量级，最终算出来的蠕滑力斜率完全不对。建议在代码里把这三个系数单独封装成一个函数，注释写清楚用的是哪一套表、泊松比取多少，方便以后核对。

这里我的建议是：初期先用一组示意值跑通流程（比如 a/b ≈ 1.9 时取 c11 ≈ 0.55、c22 ≈ 0.45、c23 ≈ 0.35），把整条链路打通之后，再花时间去抠系数表的精度。反过来先抠系数表，流程没通，根本判断不出误差是系数带来的还是代码逻辑带来的。

三个柔度的量纲都是 长度/(力/面积) = 面积/力。物理含义是：施加 1 Pa 的切向应力，材料点位移 4.9e-13 m。数值很小是因为弹性模量很大，这很正常。

3.3 网格、步进与粘着-滑移判定的递推

FASTSIM 的核心循环逻辑可以用一句话概括：把接触椭圆切成 Nx × Ny 个矩形单元，对每一行（固定 y），沿着 x 从入口 −a 走到出口 +a，每一步根据局部滑动量累加切向应力，累加完立刻检查是否超过摩擦极限，超了就截断。

局部滑动量的表达式是：

• w_x = νx − φ·y + ∂u_x/∂x
• w_y = νy + φ·x + ∂u_y/∂x

其中 νx、νy 是纵向和横向蠕滑率，φ 是自旋蠕滑率（单位 rad/m），x、y 是单元中心坐标。自旋项 φ·y 和 φ·x 反映的是刚体转动在接触斑内不同位置造成的相对速度差异——轮对绕垂直轴转一下，接触斑前后两端的滑动方向是相反的。

把位移梯度项用柔度系数替换掉（u_x = L1·q_x），就得到递推式：

• q_x(i+1) = q_x(i) − (νx − φ·y)·Δx/L1
• q_y(i+1) = q_y(i) − (νy + φ·x)·Δx/L2

这个递推是在“假定该点粘着”的前提下做的，因为粘着意味着局部滑动为零，位移梯度完全由刚体蠕滑提供。算完之后算出该点的合成切向应力幅值 |q| = √(q_x² + q_y²)，如果超过 μ·p(x,y)，就把它按比例缩放到 μ·p，并把这个单元标记为滑移单元。

缩放操作就是整个算法里唯一的非线性环节。别小看这两行代码，粘着区和滑移区的分界线就是靠它自然浮现出来的，不需要事先判断，这正是 FASTSIM 最漂亮的地方。

步进方向有讲究。必须从接触斑入口（x = −a）往出口（x = +a）扫，不能反过来，因为材料是从入口进入接触区的，它的变形历史从入口开始累积。反着扫出来的结果物理上说不通，数值上可能还挺像，但饱和段会明显偏大。

3.4 完整主循环与力/力矩的积分

把上面的逻辑写成函数，输入接触参数和蠕滑率，输出三个力分量：

function out = fastsim(a, b, p0, G, mu, nu_x, nu_y, phi, Nx, Ny, c11, c22) % FASTSIM 简化切向接触求解 % a, b 接触半轴(m)，p0 峰值压力(Pa)，G 剪切模量(Pa) % nu_x, nu_y 纵/横向蠕滑率，phi 自旋蠕滑率(rad/m) dx = 2*a / Nx; dy = 2*b / Ny; % 柔度系数 L1 = 8*a / (3*G*c11); L2 = 8*a / (3*G*c22); % 单元中心坐标 xc = -a + dx*(0.5:Nx-1+0.5); yc = -b + dy*(0.5:Ny-1+0.5); qx = zeros(Nx, Ny); qy = zeros(Nx, Ny); for j = 1:Ny y = yc(j); for i = 1:Nx x = xc(i); % 参考点压力 r2 = (x/a)^2 + (y/b)^2; if r2 >= 1 p = 0; else p = p0 * sqrt(1 - r2); end if i == 1 qx(i,j) = 0; qy(i,j) = 0; else % 粘着假设下的增量 qx(i,j) = qx(i-1,j) - (nu_x - phi*y) * dx / L1; qy(i,j) = qy(i-1,j) - (nu_y + phi*x) * dx / L2; end % 摩擦极限截断 qmag = sqrt(qx(i,j)^2 + qy(i,j)^2); qlim = mu * p; if qmag > qlim && qmag > 0 scale = qlim / qmag; qx(i,j) = qx(i,j) * scale; qy(i,j) = qy(i,j) * scale; end end end % 积分得到合力 dA = dx * dy; Fx = sum(qx(:)) * dA; Fy = sum(qy(:)) * dA; Mz = 0; for j = 1:Ny for i = 1:Nx Mz = Mz + (xc(i)*qy(i,j) - yc(j)*qx(i,j)) * dA; end end out.Fx = Fx; out.Fy = Fy; out.Mz = Mz; out.qx = qx; out.qy = qy; end

用前面算出的参数跑一次：a = 8.13e-3、b = 4.29e-3、p0 = 1.10e9、G = 8.1e10、μ = 0.3、νx = 1e-3、νy = 0、φ = 0、Nx = Ny = 20、c11 = 0.55、c22 = 0.45。算出来的 Fx 应该在 −1500 N 到 −2000 N 之间（负号取决于蠕滑率的符号约定），跟 Kalker 线性理论的预测值 Fx = −G·c11·a·b·νx ≈ 1554 N 在同一个量级。

如果跑出来是 −1.5e6 N，多半是柔度系数的量纲错了；如果跑出来是 −15 N，多半是 Δx 用了归一化坐标但柔度系数没跟着改。这两个方向差六个数量级，很好判断。

4. 复现最容易翻车的几处参数与约定

4.1 蠕滑率的定义和千分号

蠕滑率的定义是 ν = (V_轮 − V_轨)/V_名义，本质是一个无量纲的速度比。实际轮轨的蠕滑率量级很小，纵向蠕滑率在直线巡航时接近零，制动和牵引时会到千分之几，曲线通过时横向蠕滑率能到百分之一。

问题在于工程界习惯用千分号（‰）甚至百分号（%）来表达，而代码里必须是无量纲小数。表格里写“蠕滑率 1.5‰”，代码里要填 0