做轮胎模型这件事,工程圈子里绕不开PAC。哪怕现在各家都有自己标定好的高精度模型,聊起车辆动力学仿真里的轮胎,大家第一反应还是那句“魔术公式”。这篇就基于我实际用Matlab/Simulink 2021a搭建PAC轮胎模型的经验,把公式理解、Simulink实现、参数拟合、踩坑记录一次讲透。适合刚接触整车动力学仿真的工程师、车辆工程专业的学生,以及准备把轮胎模型接进ABS/ESC/自动驾驶控制算法里的开发者。
1. 为什么做轮胎模型,首选还是PAC
1.1 不是所有轮胎模型都叫“魔术”
先聊一个基础问题:轮胎模型那么多,为什么PAC能成为业界事实标准?
轮胎模型本质上是轮胎输入输出关系的数学映射。输入是滑移率、侧偏角、垂直载荷、路面附着系数,输出是纵向力、横向力、回正力矩。线性轮胎模型只能描述小侧偏角范围内的近似关系,适合自行车模型做操纵稳定性理论推导,但侧偏角一大就完全失真。Dugoff模型是纯解析模型,公式简洁、参数少,适合做实时控制算法嵌入,但精度有限,尤其是联合工况(同时纵向滑移和侧偏)下误差明显。国内也有几种半经验模型方案,工程上能跑通,但文献资料、标定经验、工具链支持都不如PAC成熟。
PAC的优势在于三点:一是数学结构是半经验的,也就是说公式形式来自大量轮胎实验数据的规律总结,每个参数都有明确的物理意义,不像纯黑箱模型那样不可解释;二是参数数量适中,全工况参数也就几十个,配合实验台架数据识别起来不算难;三是生态成熟,从论文、书籍到Simulink工具箱都有直接支持,网上能找到大量标定好的参数示例。
我见过不少项目初期尝试自己从零搭一个“更精确”的轮胎模型,最后都被参数辨识和数值稳定性折磨到放弃,老老实实切回PAC。不是说别的不行,而是PAC的精度和工程效率比最平衡。
1.2 常见版本:从PAC94到PAC2002,怎么选
PAC有几个版本,工程中最常遇到的是PAC94和PAC2002,以及衍生出来的魔术公式简化版。
PAC94是最经典的版本。纵向力公式只依赖滑移率,横向力公式只依赖侧偏角,垂直载荷作为修正项挂在参数外面。公式结构简单,单工况拟合精度高,做纯纵滑或纯侧偏的研究完全够用。但PAC94有一个明显短板:纵向力和侧向力是分开计算的,没有显式的联合工况耦合处理,在同时存在制动/驱动和转向的工况下,计算出来的轮胎力可能超出摩擦圆约束。
PAC2002是荷兰TNO改进并推广的版本,最大的变化是引入了联合工况修正,用无量纲滑移率把纵向滑移和侧向滑移统一处理,同时增加了大滑移区的松弛行为修正。做整车操稳仿真、极限工况、电子稳定控制策略验证,建议直接用PAC2002。Simulink的Vehicle Dynamics Blockset里内置的Magic Formula Tire帮助文档就基于PAC2002,参数结构也是按这个版本组织的。
我做整车控制器仿真时习惯直接选PAC2002,虽然初始参数多了一列,但公式里自带防“不合理外推”的保护逻辑,整车的仿真稳定性好很多。如果你只是做单轮台架数据拟合练习,PAC94上手更简单,后面再升级也不难,两个版本的核心变量描述是相通的。
2. 核心公式不是背出来的,是“理解”出来的
很多人拿到PAC参数表一上来就背公式,这是最容易走弯路的地方。PAC公式看着复杂,但只要你理解了每个字母在曲线形状上“画”的是什么,拟合调参的时候就能少调白发。
2.1 纵向力与滑移率:D、C、B、E各干一件事
PAC纵向力公式的基本形式是:
F_x = D * sin(C * atan(B * x - E * (B * x - atan(B * x)))) + S_v这里的x通常是滑移率加上一个水平偏移项。初看那个多层嵌套的atan很唬人,其实拆开就清楚了:
- D是峰值因子,代表曲线的最大纵向力值。直观理解就是轮胎能提供的最大抓地力,主要由垂直载荷和路面附着决定。
- C是形状因子,控制输出曲线是“正弦型”还是“反正切型”,纵向力的C通常在1.6附近,这个值决定了曲线头部是否平缓、尾部是否下压。
- B是刚度因子,它和C、D一起决定曲线原点附近的斜率,也就是轮胎的纵向滑移刚度。原点斜率等于B * C * D,所以三者是相乘关系,不是独立的。这个约束在做参数辨识时特别重要,初值给得不好,拟合结果经常在刚度项上飘。
- E是曲率因子,用来修正曲线的整体弯曲程度。E为正时曲线峰值附近的过渡更平滑,E为负时曲线中段容易出现一个下凹拐弯。
实际影响直观感受一下:如果你在Simulink里用同一组实验数据拟合出来的纵向力曲线在滑移率0.1附近总是低于实测点,优先考虑调E,而不是去动B。因为B动一点点,原点的斜率就全变了,整体形态会歪掉。
2.2 横向力与侧偏角:重点理解侧偏刚度
横向力公式形式和纵向力几乎一样,只是x换成了侧偏角,同时多了水平偏移S_h和垂直偏移S_v,用来处理由于帘布层转向效应造成的曲线不对称。
F_y = D * sin(C * atan(B * x - E * (B * x - atan(B * x)))) + S_v侧偏工况下D就是峰值侧向力,B * C * D对应的是侧偏刚度。这个侧偏刚度在小侧偏角区(通常5度以内)几乎就是线性段,整车操稳分析里用的就是这段斜率。你拿PAC参数表算出来的侧偏刚度如果和实验报告的数值差异很大,大概率是B因子的归一化基准没搞对。
这里有一个工程上很常见的误解:以为垂直载荷只是一个乘进去的系数。其实PAC里F_z的影响是通过各系数随F_z的“二次函数拟合”实现的,不是简单线性缩放。也就是说,同一套B、C、D参数,在不同载荷下会各自变成不同的值。这就是为什么PAC参数表不是一个常数,而是一张“参数随载荷变化”的表格。Simulink实现时,你要么暴力查表,要么在线计算这些系数。
2.3 回正力矩:大滑移区会“翻转”,千万别忽略
回正力矩M_z的公式结构更复杂,但工程上最重要的是记住一个现象:侧偏角增大到一定程度后,回正力矩会从正变负,也就是所谓“翻转区”。这个现象描述的是轮胎自对准扭矩在大侧偏角下开始反向,直接影响驾驶员手力感和EPS标定逻辑。
如果只是做车辆轨迹仿真,回正力矩往往可以先用简化模型顶住。但如果你在做转向手感模拟或者EPS力矩控制,这个M_z的翻转特性必须用PAC完整描述,否则转向回正力矩曲线在高侧偏角区会和实车完全对不上号。
3. 在Simulink 2021a里搭PAC:三条路各有什么坑
Matlab/Simulink 2021a这个版本下做PAC轮胎模型,实操路径大体分三条。我把每条路的实现方式、适用场景、典型坑都过一遍。
3.1 路径一:Vehicle Dynamics Blockset自带的魔术公式模块
2021a的Vehicle Dynamics Blockset里已经内置了Magic Formula Tire模块,双击进去可以看到它按PAC2002组织的参数选项。这条路最简单,适合整车架构已经用Simulink搭好的场景。
实际操作步骤:
- 打开Simulink库浏览器,在Vehicle Dynamics Blockset的Tires分类下找到Magic Formula Tire模块。
- 把模块拖到模型中,连接到车轮的纵向滑移率、侧偏角、垂直载荷输入端口。
- 双击模块,在参数页选择Tire model version为PAC2002,然后加载一个预置的轮胎参数.mat文件或者手动输入参数向量。
- 连接输出端的F_x、F_y、M_z信号到车身动力学模块。
坑点有两个。第一,内置模块对输入信号的单位有严格要求,滑移率用的是无量纲,侧偏角要用弧度而不是角度。很多初学者直接把角度信号连进去,仿真出来的力直接爆炸,还以为是参数问题。第二,内置模块坐标系约定和自定义模型可能不一致,F_x正方向是车轮前进方向还是胎体坐标系下的方向,一定要对着文档确认,不然转弯工况整车横摆方向会反。
3.2 路径二:MATLAB Function手写公式:全流程可控
如果不想被工具箱的固定框架限制,自己手写一个PAC轮胎模型也就几十行代码。我日常做算法验证时反而更喜欢这条路,因为参数、状态、接口完全自己说了算,配合Simulink的MATLAB Function模块能随时改逻辑。
下面是一份基于PAC2002简化的实现示例,只做纵向力部分,展示核心逻辑:
function Fx = PAC_longitudinal(kappa, Fz, params) % kappa: 纵向滑移率 (无量纲) % Fz: 垂直载荷 (N) % 从params结构体取PAC2002纵向力系数 D = params.pD1 * Fz.^2 + params.pD2 * Fz; % 峰值因子随载荷变化 BCD = params.pB1 * Fz + params.pB2; % 纵向刚度随载荷变化 C = params.pC1; % 形状因子近似常数 B = BCD / (C * D + eps); % 刚度因子反解 E = params.pE1 * Fz + params.pE2; % 曲率因子随载荷变化 % 魔术公式计算 % 这里已经包含水平偏移Sv和垂直偏移Sh的简化(取0) x = kappa; Fx = D * sin(C * atan(B*x - E*(B*x - atan(B*x)))); end这个函数放在MATLAB Function模块里,外部通过结构体传入标定好的参数即可。Simulink调用时注意数据类型和维度,MATLAB Function默认输出double,如果下游模块是single,中间需要加Data Type Conversion。
自己手写的好处是不受模块版本限制,2021a和2023a都能跑。坏处是细节容易错。我踩过的坑包括:atan还是atan2用混、角度弧度忘记换算、B被C和D除时没加eps导致零载荷除零。零载荷除零这个事特别隐蔽,整车仿真在车轮离地瞬间Fz变零,模型直接NaN,后面所有积分全乱。我建议在函数入口加一个Fz下限钳位,比如Fz = max(Fz, 10),既保证数值稳定又不影响正常工况。
3.3 路径三:查表法替代:数据够多时更快
如果手里有台架实测的轮胎力数据表,而不是一套现成的PAC参数,可以直接用Lookup Table(查表模块)来做。这种做法在2021a的Simulink里最容易实现,本质是把“公式拟合”这一步省略,让模型直接用插值输出力。
三维查表的典型维度是:滑移率、侧偏角、垂直载荷。输出是纵向力Fx、侧向力Fy。Simulink里的n-D Lookup Table可以设置断点维度和插值方法。建议插值方法用线性,外推用“Clip”模式,避免外推发散。
查表法最大的坑在于数据稀疏区。台架实验不可能把所有载荷、所有滑移率组合都采一遍,如果整车仿真工况跑到了数据边界之外,查表会直接输出边界值,整车响应看起来“卡死”了一样。所以我的习惯是:查表法只用来做快速原型验证,最终交付的模型还是回归到PAC公式。
4. 参数拟合:没有实验数据,PAC就是一堆字母
PAC模型精度上限由实验数据决定。参数拟合这条路基本绕不开,我拆成三个环节讲。
4.1 实验数据清洗与坐标规范化
很多台架导出的轮胎数据文件里,纵向力和滑移率符号约定每个台架都不一样。第一步永远是统一坐标系。建议统一成车辆动力学标准:纵向力正方向为前进方向,滑移率全滑移为1,纯滚动接近0,制动滑移率为负。
数据清洗至少要处理三件事:剔除明显异常尖峰(比如传感器打滑瞬间的跳变)、对力信号做低通滤波(截止频率建议在20-50Hz,看采样率)、把重复工况的几组数据做平均处理。别小看这一步,数据不干净,后续拟合出来的参数再漂亮也是废的,因为代进去仿真根本收敛不了。
4.2 用lsqcurvefit做初值拟合:别直接扔给工具箱
Matlab里做PAC参数拟合最常用的就是lsqcurvefit,但直接扔给优化器十个参数乱搜,结果几乎必是局部最优解。我的做法分两步走。
第一步,固定物理约束,缩小搜索范围。比如已知这台轮胎的峰值附着系数大约是0.9,那么这个载荷下的D初值就按0.9 * Fz给;侧偏刚度按实验曲线小角度段的斜率直接量出来,再反算B。这样一来,D、B都有物理初值,C取1.3到1.6之间的典型值,E取0,S_v和S_h取0。整套初值已经能画出一条和实验数据大差不差的曲线。
第二步,设置边界约束。比如B不能小于0,C在1到2之间,E一般在-1到2之间。这些边界是我从几百组标定数据里总结的常见取值范围,设好之后优化器就不会跑到荒唐区域去了。
拟合代码如下:
% kappa_data: 滑移率实验数据 % Fx_data: 纵向力实验数据 % Fz_nom: 实验载荷 % 初值结构 [D, B, C, E] x0 = [0.9 * Fz_nom, 15, 1.6, 0.2]; lb = [0, 0, 1, -1]; ub = [1.2 * Fz_nom, 50, 2, 2]; % 定义拟合函数 fun = @(x, kappa) x(1) * sin(x(3) * atan(x(2)*kappa - x(4)*(x(2)*kappa - atan(x(2)*kappa)))); % 拟合 [x_opt, resnorm] = lsqcurvefit(fun, x0, kappa_data, Fx_data, lb, ub);fminsearch别用,就是个玩具;lsqnonlin也可以,但lsqcurvefit直接传数据向量更顺手。
4.3 摩擦圆与联合工况约束
单工况拟合完,纵向力和横向力参数可能各自看起来都挺好,但联合工况下轮胎合力会突破摩擦圆,也就是sqrt(Fx^2 + Fy^2) > mu * Fz。这说明两套参数的峰值因子D不一致导致的。
解决思路是:把纵向力峰值和横向力峰值用同一个摩擦系数约束起来。具体做法是在联合工况下引入摩擦圆修正系数,让纵向力形函数和横向力形函数通过当前位置的滑移状态分配附着能力。不要用简单乘法把Fx和Fy各自乘一个折扣系数,那样低速小滑移区会失真。PAC2002自带的无量纲滑移率组合公式就是干这个的,所以前面推荐直接上PAC2002。
5. 我踩过的坑:仿真发散、标定乱跳、结果不自洽
最后整理一份排查速查表,都是我在实际项目中反复遇到的高频问题。
| 现象 | 常见原因 | 排查方向 |
|---|---|---|
| 轮胎力瞬间NaN/Inf | Fz出现0或负值 | 在模型入口加Fz下限钳位,检查悬架模型是否输出负载荷 |
| 低速大滑移时Fx剧烈振荡 | 求解器步长过大,或滑移率计算分母接近0 | 使用变步长求解器并限制最大步长,滑移率计算分母加eps |
| 拟合收敛但曲线形态明显错误 | 初值里B和D比例失调 | 按物理刚度反算B初值,不要随机给 |
| 转角输入下横摆方向反向 | 侧偏角正负号约定与模型不一致 | 核对坐标系定义,侧偏角正方向统一 |
| 查表法在外推区整车响应“锁死” | Lookup Table外推设置成了Clip | 改用线性外推,或者确认工况不超数据边界 |
| 联合工况下同时制动和转向,车辆过于“推头” | 纵向和横向D没做摩擦圆一致约束 | 统一峰值因子,检查摩擦圆约束条件 |
| 仿真速度慢到跑不动 | 模型用了连续变步长+高精度公式轮胎力,且没做离散化 | 用Tustin或零阶保持将轮胎力计算离散化,步长设到1ms以上看能否接受 |
其中数值振荡这个事我想多说几句。Simulink默认变步长求解器在轮胎模型进入高刚度区时,很容易发生步长骤降,仿真速度变成爬行。我在做ABS仿真时遇到过2秒仿真跑了20分钟的情况。后来把轮胎模型单独抽出来做成离散模块,或者把最大步长限制到1ms并选用ode23t,速度立刻提升了几个量级。轮胎模型本质上是代数方程而不是状态方程,没必要参与连续积分器的刚性求解,离散化完全能满足工程精度。
还有一个标定“乱跳”的问题,本质是使用了过拟合的拟合算法。轮胎实验数据本身有噪声,如果你用很高阶的多项式去拟合Fx曲线,在数据点之间会画蛇添足地扭来扭去,仿真时车辆响应也跟着抖。PAC之所以用固定的sin(atan())结构,就是防止过度自由拟合。遇到标定结果“乱跳”,回去审视数据质量而不是加自由度。
6. 从单胎模型到整车仿真:怎么把PAC接进去
6.1 接口设计
PAC模型本身只负责输出轮胎力,把它接到整车动力学模型前,你需要定义好一套稳定接口。我惯用的接口包括输入:纵向滑移率kappa、侧偏角alpha、垂直载荷Fz、路面附着系数mu(可选);输出:纵向力Fx、侧向力Fy、回正力矩Mz。这样设计的好处是无论后面换查表模型还是换成别的高精度轮胎模型,整车框架不用动。
6.2 给ABS/ESC留接口
做ABS/ESC策略仿真时,轮胎模型的输出往往还要加一个“路面附着估计”信号。虽然PAC模型内部用的是真实mu,但控制器算法里估计的mu应该走另一个通道,不能把仿真里的真值直接喂给控制器,否则你验证的是理想环境而不是算法性能。
我常在胎模型外层包一个使能和故障注入逻辑,模拟传感器丢帧、轮胎力延迟输出等故障场景。这个小成本的做法对后面做控制器鲁棒性验证帮助很大。
6.3 想再进一步,还能加点什么
基础PAC跑通之后,你可以考虑扩展三个方向:一是加轮胎热模型,用胎温修正峰值附着系数,模拟连续激烈驾驶下轮胎性能衰减;二是加磨损模型,里程增加导致侧偏刚度下降;三是加动态松弛长度,用一个一阶惯性环节把轮胎力延迟一下,模拟瞬态响应。前两个方向在真实赛车仿真、耐久测试里有用,第三个方向对底盘控制器的相位裕度分析有价值。
在实际使用中我的体会是,PAC模型本身的Simulink实现并不难,难的是参数和实验数据之间的闭环迭代。每次拿到一批新的台架数据,都要重新走一遍“数据清洗 -> 初值估算 -> 拟合 -> 整车验证 -> 反馈修正”这个循环。模型搭好只是开始,花时间和数据较劲才是日常。给新手一个建议:不要迷信看到的任何一组“通用参数”,轮胎参数和轮胎型号、路面、胎压都绑定,换一个场景就必须重新标定一轮。这听起来麻烦,却是走向靠谱仿真必须的过程。
最后分享一个小技巧。在Simulink模型里做一个PAC模型版本切换开关,用Manual Switch在“查表快速原型”和“PAC公式精算”之间切换。前期调整车控制策略的时候用查表跑速度,后期做精度验证时切到PAC公式,两个场景都不用改整车框架。这个习惯帮我省了大量重复调试时间,第一次搭模型就顺手把开关做好,后面会非常舒服。