★ FEATURED ARTICLE 轮轨关系这四个字写在纸上是一组教科书公式落到代码里就是几百行的数值循环加一堆容易搞错的约定。很多做车辆动力学的人日常工作是调用现成的多体软件轮轨力那部分是个黑箱——输入轴重和蠕滑率输出三个方向的力和力矩曲线平滑、看着很舒服。可一旦需要自己写一个轮轨力模型接进自研求解器或者只是想搞清楚那条蠕滑力曲线为什么在某个点突然拐弯就会发现从 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.190B ½×(2.0 3.333) 2.667A/B 0.446。这个比值决定了后面椭圆系数怎么查。2.2 椭圆积分系数 m、n 的求解与查表的取舍Hertz 解的最终形式里有两个无量纲系数 m 和 n它们是 A/B 的函数a m·[3P/(2E*(AB))]^(1/3)b n·[3P/(2E*(AB))]^(1/3)这两个系数的严格求法涉及第一类和第二类完全椭圆积分需要数值迭代。工程上最常用的做法是查表加插值下面这张表是我用了很多年的一组近似值来源是经典弹性力学手册用来做工程计算足够A/B1.00.80.60.50.40.30.20.1m1.0001.0381.1091.1661.2451.3611.5461.884n1.0000.8720.7530.6900.6220.5450.4510.319注意不同手册给出的 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 80000R1x 0.42R1y 0.5R2x infR2y 0.3E 206e9ν 0.3。得到 A 1.190B 2.667A/B 0.446插值出 m ≈ 1.25、n ≈ 0.66c ≈ 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 604×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_solveR1x R1y R2x R2y 0.02得到的 a、b 应该都是 0.51 mmp0 一致。差超过 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/∂xw_y νy φ·x ∂u_y/∂x其中 νx、νy 是纵向和横向蠕滑率φ 是自旋蠕滑率单位 rad/mx、y 是单元中心坐标。自旋项 φ·y 和 φ·x 反映的是刚体转动在接触斑内不同位置造成的相对速度差异——轮对绕垂直轴转一下接触斑前后两端的滑动方向是相反的。把位移梯度项用柔度系数替换掉u_x L1·q_x就得到递推式q_x(i1) q_x(i) − (νx − φ·y)·Δx/L1q_y(i1) 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-10.5); yc -b dy*(0.5:Ny-10.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
阅读完成 · 觉得有帮助？
