前一阵子做极限工况仿真折腾了好几个礼拜最后发现最卡脖子的不是控制算法而是车辆模型本身。普通乘用车模型在常规工况下精度还行但只要把车扔进漂移状态模型就开始“胡来”——后轮侧偏角一超过20度线性轮胎模型给出的侧向力早就飞到天上去了仿真结果根本不可能是对的。这篇文章就记录一下我搭建的一套适用于漂移工况的非线性三自由度车辆模型从坐标定义、轮胎建模到Simulink子系统实现全部展开讲透顺便把我踩过的坑也一并交代清楚。这套模型适合谁如果你在做车辆动力学仿真、漂移控制算法开发、极限工况下的底盘策略验证或者纯粹想搞明白“漂移时车辆到底处于什么受力状态”那你应该会需要它。我会尽量用工程师之间聊天的口吻来讲不绕弯子该上公式就上公式该贴参数就贴参数。1. 为什么线性模型一到漂移工况就失效1.1 线性二自由度模型的边界在哪里经典线性二自由度模型也叫“自行车模型”是车辆动力学入门的第一课。它把左右轮合并成一个前轴和一个后轴只保留横向速度和横摆角速度两个自由度侧向力直接用“侧偏刚度 × 侧偏角”来算Fy -Cα · α这个式子成立的前提是轮胎工作在线性区也就是侧偏角通常不超过3到5度。日常驾驶、高速变道、稳态圆周行驶基本都落在这个范围内。此时侧偏刚度Cα可以看作常数整车横摆响应也能用传递函数很漂亮地表达出来。问题是漂移工况完全不满足这个前提。漂移时后轮的侧偏角经常到15度甚至40度轮胎早就过了线性区进入饱和区甚至下降区。你再用一个常数刚度去近似算出来的后轴侧向力可能比真实值大出两三倍横摆力矩完全失真。1.2 漂移工况带来的三个核心挑战除了轮胎非线性漂移工况还有几个线性模型从根本上就没法处理的问题。第一个是载荷转移。漂移时侧向加速度通常超过1g四个车轮的垂直载荷变化幅度非常大。举个例子一台1400kg的车1g侧向加速度下前后轴之间的载荷转移量可以达到几百牛到上千牛同轴左右轮的垂直载荷可能相差一倍以上。轮胎侧向力本质上是“摩擦系数 × 垂直载荷”载荷一变抓地力就变线性模型里固定不变的等效侧偏刚度自然失效。第二个是纵向力与侧向力的耦合。漂移不是纯转向状态它通常需要后轮大油门甚至拉手刹来维持后轴滑移。这时后轮同时承受大纵向驱动力和大侧偏力轮胎工作点落在“摩擦椭圆”的边界附近侧向力会被纵向力挤占。线性模型把纵向和侧向完全解耦在这类工况下给出的侧向力会严重偏高。第三个是运动学非线性项。线性模型通常会做小角度假设把速度变化率方程里的一些交叉项直接省略。但漂移时横摆角速度可能到1 rad/s以上纵向和横向速度的耦合项比如 r·vx、r·vy非常大不能忽略。所以结论很明确要分析漂移机理至少需要纵向、横向、横摆这三个自由度的非线性模型并且轮胎、载荷转移、纵滑侧偏联合工况这三个环节一个都不能省。1.3 为什么三自由度而不是七自由度或十四自由度做仿真的人看到这里可能会问现在商业软件里动不动就十四自由度甚至更高为什么这里选三自由度我的理由是漂移的核心机理主要体现在车身平面运动和后轮轮胎力饱和上更高自由度的悬架运动、俯仰侧倾动力学对解释“为什么要反打方向盘、为什么后轮会抱死滑移”这个问题贡献有限。三自由度模型保留了纵向速度、横向速度、横摆角速度这三个最关键的状态量同时通过准静态载荷转移公式把侧倾、俯仰的影响折算到四个车轮的垂直载荷上。这样做的优势很明显参数少、计算量小、Simulink里搭起来一天就能跑通而且因为模型简单异常行为容易被定位。更重要的是后续如果要做漂移平衡点分析或者闭环控制器设计三自由度的模型维度和状态量都足够支撑线性化处理和反馈控制设计。2. 三自由度模型的坐标系与动力学方程定义2.1 坐标系选取和状态量定义我自己习惯用ISO坐标系车辆前进方向为x轴正方向左侧为y轴正方向z轴垂直向上。横摆角速度r以逆时针为正俯视车辆逆时针旋转为正。这个坐标系在多数车辆动力学文献里通用性最好和轮胎模型参数的符号约定也不容易打架。模型保留三个自由度对应的三个状态量是vx车体坐标系下的纵向速度单位 m/svy车体坐标系下的横向速度单位 m/sr绕z轴的横摆角速度单位 rad/s注意这里vx和vy都是车体坐标系的不是大地坐标系。如果你习惯用大地坐标系下的 GPS 速度来换算需要额外做坐标变换但在这个模型内部直接用车体坐标系是最方便的。2.2 车身平面运动方程推导三个自由度的运动方程在惯性系下写出来是这样的m·(dvx/dt - r·vy) ΣFx m·(dvy/dt r·vx) ΣFy Iz·dr/dt ΣMz这三条式子对应的物理意义要讲清楚第一条是纵向力平衡加离心效应修正第二条是横向力平衡第三条是横摆力矩平衡。左边括号里的-r·vy和r·vx就是上面提到的运动学非线性项漂移工况下它们对受力平衡的影响非常显著一定不能省略。这四个车轮的纵向力和侧向力需要先合成到车体坐标系再代进上面的三个方程里。以四轮独立建模为例前轮转角为δ后轮不转向合成后的车体纵向力、横向力和横摆力矩可以写成ΣFx (Fx_fl Fx_fr)cosδ - (Fy_fl Fy_fr)sinδ Fx_rl Fx_rr ΣFy (Fx_fl Fx_fr)sinδ (Fy_fl Fy_fr)cosδ Fy_rl Fy_rr ΣMz a·[(Fy_fl Fy_fr)cosδ (Fx_fl Fx_fr)sinδ] - b·(Fy_rl Fy_rr) (tw_f/2)·[(Fx_fr - Fx_fl)cosδ (Fy_fr - Fy_fl)sinδ] (tw_r/2)·(Fx_rr - Fx_rl)这里a是质心到前轴的距离b是质心到后轴的距离tw_f和tw_r分别是前后轮距。横摆力矩表达式里那个车轮距相关的部分就是左右轮纵向力差异产生的直接横摆力矩——漂移工况里通过后轮左右驱动力矩差来“扭屁股”正是靠这一项。实际搭建Simulink模型时我建议把这三个合成公式单独封装成一个“力合成”子系统方便检查和调试后面章节会细说。2.3 轮胎侧偏角与滑转率的计算四个轮胎各自的侧偏角和纵向滑转率是连接车身状态和轮胎模型的桥梁。侧偏角定义为车轮实际速度方向与车轮朝向之间的夹角在ISO坐标系下前轮 α_fl δ - atan2(vy a·r, vx - tw_f/2·r) α_fr δ - atan2(vy a·r, vx tw_f/2·r)后轮 α_rl -atan2(vy - b·r, vx - tw_r/2·r) α_rr -atan2(vy - b·r, vx tw_r/2·r)注意前轮转角δ本身对侧偏角的影响是减去的驾驶员的转向角反打方向时δ为负侧偏角会向反方向偏移这正是漂移中“前轮反打”的数学体现。纵向滑转率定义的符号约定各家常不一样我自己习惯定义为驱动时κ (ω·R - vx) / (ω·R)制动时κ (vx - ω·R) / vx这样处理可以保证滑转率始终在0到1附近区间方便调用魔术公式。漂移仿真中后轮的驱动力矩通常很大ω·R会明显大于vxκ会比较接近1轮胎工作点在纵向力饱和区需要轮胎模型有足够的数值稳定性。3. 非线性轮胎模型与载荷转移的计算3.1 魔术公式的基本形式和参数含义提到非线性轮胎模型绕不开魔术公式。它的基本形式是y D·sin(C·arctan(B·x - E·(B·x - arctan(B·x))))式子里的x是输入侧偏角或者滑转率y是输出力四个系数含义如下B刚度因子决定曲线初始斜率C曲线形状因子决定曲线整体形态侧向力曲线通常在1.3左右D峰值因子决定曲线峰值大小通常等于 μ·FzE曲率因子决定峰值附近的弯曲程度和下降段形态对于漂移仿真来说峰值过后的下降段特别重要。漂移时后轮侧偏角很大很可能越过峰值点进入“微降段”。如果C值、E值选得不合适后轴侧向力随时间的变化趋势就会失真。我一般在拟合轮胎数据时会在20度到40度区间多给些权重确保大侧偏角下的力变化趋势是对的。这组参数从哪里来最优解是轮胎厂的Pacejka参数文件但多数项目拿不到。退而求其次用同尺寸轮胎的公开测试数据做最小二乘拟合也能得到能用的参数。实在不行就先从经验值起步再拿整车的稳态转向试验去校核这个调参流程后面会专门写。下面给一组我在某紧凑型后驱车仿真里用过的参数仅供参考参数前轮后轮说明C1.301.25后轮峰值过后下降更明显D1.10·Fz1.15·Fz后轮附着系数略高B0.110.09单位与弧度保持一致E0.850.90后轮更“尖”一些实际上这四个参数也跟轮胎的垂直载荷相关严格来说应该做成Fz的函数。更精细的做法是每个Fz点一组参数再用插值表。但对于漂移机理分析固定的B/C/D/E已经能复现绝大多数关键现象了。3.2 四轮垂直载荷的实时计算前面提到漂移工况下载荷转移是影响轮胎力的关键变量。我用的是准静态载荷转移模型不引入悬架自由度而是直接根据纵向加速度ax和横向加速度ay来估算四个轮的垂直载荷。先算静态载荷分配Fz_f_static m·g·b / L Fz_r_static m·g·a / L其中L a b。然后考虑纵向载荷转移。加速时ax 0后轴加载、前轴减载ΔFz_long m·ax·h / L再考虑横向载荷转移。这里我用一个简化假设前后轴的侧倾力矩分配等于其静态轴荷分配比例即前轴承担的比例是b/L后轴是a/L。于是ΔFz_front_cross m·ay·h / L · b / L ΔFz_rear_cross m·ay·h / L · a / L有了这些量四个轮的垂直载荷就能写出来了。以向右转向ay为正为例左侧车轮加载右侧减载Fz_fl (Fz_f_static ΔFz_long) / 2 ΔFz_front_cross / 2 Fz_fr (Fz_f_static ΔFz_long) / 2 - ΔFz_front_cross / 2 Fz_rl (Fz_r_static - ΔFz_long) / 2 ΔFz_rear_cross / 2 Fz_rr (Fz_r_static - ΔFz_long) / 2 - ΔFz_rear_cross / 2注意这里ay本身就是车身状态的函数ay dvy/dt r·vx。模型里会出现“载荷转移依赖轮胎力轮胎力依赖垂直载荷垂直载荷依赖加速度加速度依赖轮胎力”的反馈环。在Simulink里直接搭这种代数环容易卡死后面第五章会有处理办法我建议在这里用上一时刻的ay做近似或者干脆把ay直接从上一仿真步的数值里取。另外还要加一道限制垂直载荷不能为负值。极端漂移情况下内侧车轮可能完全离地Fz为负说明物理上就已经失效了模型这时候基本也没法继续跑。我在代码里会加一句max(Fz, 0)防止数值崩溃。3.3 纵向力与侧向力的摩擦椭圆约束单独由魔术公式算出来的纯侧向力和纯纵向力并不能直接叠加使用。因为轮胎与地面的总附着力是有限的存在一个“摩擦椭圆”边界。漂移时后轮大油门驱动纵向力已经占据了相当一部分摩擦裕度侧向力必须打折。工程上最常用的做法是对纯侧向力乘以一个修正系数Fy_combined Fy_pure · sqrt(1 - (Fx / (μ·Fz))²)这个公式的几何意义就是椭圆约束纵向力越大可用的侧向力越小。μ·Fz就是当前轮载下的最大总附着力。实际实现时需要注意两点。第一Fx本身也是由魔术公式算出来的因此要先算纯纵向力再加椭圆修正顺序不要反。第二如果Fx超过μ·Fz括号里面变成负数sqrt就出问题了。我在代码里会对这个比值做限幅确保不超过0.99。漂移工况下后轮大油门时Fx接近甚至达到附着极限Fy会被压得很低这就解释了为什么后轴那么容易失去横向抓地力——不是“轮胎抓地力突然没了”而是纵向力把侧向裕度吃光了。这个机理在纯侧偏模型里完全体现不出来所以在搭建模型时千万不能省掉椭圆约束这一步。4. Simulink建模架构与子系统实现4.1 顶层架构与信号流设计Simulink里搭这个模型我强烈建议按照信号流的方向来组织子系统而不是把一堆公式塞进一个MATLAB Function里。我的顶层结构是这样分的输入层驾驶员转向角δ、后轴驱动转矩Tq也可以拆成四个轮各自的控制转矩车辆运动学子系统由状态量计算侧偏角、滑转率、加速度等中间量轮胎力子系统根据侧偏角、滑转率、垂直载荷计算四个轮各自的Fx和Fy载荷转移子系统计算四个轮的垂直载荷Fz力合成子系统把轮胎力合成为车体ΣFx、ΣFy、ΣMz状态积分层由运动方程和合成力计算出dvx/dt、dvy/dt、dr/dt再积分得到vx、vy、r用这种模块化结构的好处是后期调试时可以单独给某个子系统喂一组测试信号快速定位问题。比如说怀疑轮胎模型参数不对就不需要跑完整个整车仿真直接把轮胎子系统拎出来输入一组常值侧偏角看一眼输出曲线就行了。信号类型方面建议全部用double所有角度统一用弧度制。单位制不一致导致的仿真结果异常是我见过频率最高的低级错误。4.2 车身动力学子系统搭建车身动力学部分我不用MATLAB Function直接用积分器加增益模块来搭。三个运动方程写成积分形式dvx/dt (ΣFx m·r·vy) / m dvy/dt (ΣFy - m·r·vx) / m dr/dt ΣMz / IzSimulink里就是用三个Integrator模块分别积分vx、vy、r反馈回路里有加减法和增益模块。这种加法器加积分器的搭法虽然看起来原始但信号流直观而且积分器的初始值可以直接设置成漂移工况的初始状态用起来很方便。有一点需要注意积分器输出必须限制幅值尤其是vy。如果不加限制模型在发散时vy可能跑到几十米每秒远远超出物理可能费很大劲找原因才发现是数值积分飘了。我给vy加的正负限幅一般取±20 m/s足够覆盖真实漂移工况。4.3 轮胎子系统与载荷转移的实现轮胎和载荷转移这两个子系统相互依赖我一般把它们放在同一个MATLAB Function里实现输入是四个轮的侧偏角、滑转率和垂直载荷输出是四个轮的Fx和Fy。用MATLAB Function的好处是魔术公式的复杂表达式写起来更舒服调试时也方便加断点查中间量。魔术公式的实现代码大概长这样以纯侧向力为例function Fy magic_formula(alpha, B, C, D, E) x tan(alpha); Bx B * x; Fy D * sin(C * atan(Bx - E * (Bx - atan(Bx)))); end注意x这里用tan(alpha)而不是直接alpha是为了在侧偏角比较大的时候让曲线更贴试验数据这是很多文献里的做法我实际对比过小侧偏角下差异不大大侧偏角下更合理。载荷转移子系统我放在一个独立模块里输入是当前状态vx、vy、r输出是四个Fz。ay用上一时刻的vx和r计算我在函数里把ay定义成ay vdoty r * vx;这里vdoty是上一仿真步的dvy/dt我是通过在Simulink里给这个信号加一个Memory模块或者Unit Delay来避免代数环的。虽然严格来说有一步延迟但实际仿真结果显示在1ms步长下这个延迟带来的误差可以忽略。4.4 参数初始化与基本仿真配置我习惯用matlab脚本统一管理参数而不是把数值写到模型里。这样每次改参数只需要编辑脚本模型文件本身保持干净。下面是一个最小参数集的初始化脚本示例% 车辆基本参数 m 1400; % 整车质量 kg Iz 2350; % 横摆转动惯量 kg·m^2 a 1.10; % 质心到前轴距离 m b 1.50; % 质心到后轴距离 m L a b; twf 1.62; % 前轮距 m twr 1.58; % 后轮距 m h 0.55; % 质心高度 m g 9.81; % 轮胎模型参数前/后轮 tire.Cf 1.30; tire.Bf 0.11; tire.Df 1.10; tire.Ef 0.85; tire.Cr 1.25; tire.Br 0.09; tire.Dr 1.15; tire.Er 0.90; % 初始状态低速直线 x0 [15.0; 0.0; 0.0]; % [vx; vy; r]仿真配置方面我用的求解器是固定步长ode4四阶Runge-Kutta步长1ms。变步长求解器理论上更快但漂移工况轮胎力变化很快变步长在步长调整时容易引入数值振荡反而更难调试。固定步长1ms在普通笔记本上跑几秒仿真只要几十秒完全够用。如果你要跑控制器在环建议把模型输出保存为MATLAB timeseries或者直接接Simulink Control Design工具箱做线性化。我后面做漂移平衡点分析时就直接用这个模型做了线性化处理。5. 漂移工况仿真的调试经验与常见问题5.1 仿真发散的五种典型原因第一个原因魔术公式参数的单位不对。B值如果按“度”的量级带入而你的侧偏角是弧度曲线形态完全错掉侧向力可能瞬间变成上千牛的反向力整车直接原地起飞。我的检查方式就是跑一个纯侧偏角的开环测试固定vx和vy给定一个从0到40度变化的侧偏角看轮胎力曲线是否平滑、是否单调或者按要求出现峰值下降段。第二个原因初始状态给得离谱。漂移工况的初始状态不是一个猜测值它应该大致满足漂移平衡点附近的速度和横摆角速度关系。你如果上来就给vx30, vy8, r2.5系统很可能直接积分发散。建议先用较小的漂移强度起步比如vx15, vy3, r0.8逐渐往深漂状态调。第三个原因垂直载荷变成负数。内侧轮离地时Fz小于零魔术公式里D μ·Fz变成负的力的方向就反了整车模型瞬间崩溃。前面说的Fz限幅一定不能省另外要留意仿真过程中内侧后轮的Fz是否长时间接近0如果接近0实际上说明你给的工况早就不在模型适用范围内了。第四个原因摩擦椭圆约束里开根号出问题。Fx超过μ·Fz时sqrt计算错误Simulink会提示NaN或者是复数计算错误。我处理的办法是把sqrt内部值先clip到0到1之间。第五个原因纵向驱动力矩给的过大导致后轮完全打滑空转轮胎模型在极高滑转率下数值不稳。这时候可以限制驱动转矩的爬升速率避免阶跃给矩。5.2 代数环的处理方法前面提到了由加速度驱动的载荷转移会引入代数环。Simulink在检测到代数环时通常会给出提示而且求解速度明显变慢。我处理的方式有三种按推荐程度排序办法一给ay信号加单位延迟Unit Delay模块引入一个仿真步长的滞后。这个最简单对结果的影响在1ms步长下可以接受。办法二用MATLAB Function里persistent变量缓存上一时刻的ay值等效于做一步延迟好处是不需要在Simulink图里加额外模块。办法三如果模型必须完全没有延迟那就需要用迭代求解器来解这个隐式方程Simulink里可以对代数环设置迭代求解但我实际跑下来发现迭代收敛在漂移工况下并不稳定尤其是轮胎力不光滑的区间。所以我的实践建议就一条接受一步延迟它带来的误差远远小于模型本身参数不确定性的误差。5.3 一套可复用的调参流程我总结了一套从零开始调通这个模型的操作流程照着走基本不会卡壳。第一步先把轮胎模型单独拿出来测。给固定Fz遍历不同侧偏角画出前后轮的侧向力曲线确认弧度和参数匹配、曲线形态正确、极限值在μ·Fz附近。第二步整车上电不带转向验证直线加速和制动工况。给一个较小的阶跃驱动力矩看vx是否合理变化、vy和r是否保持在接近0的数值。如果直线工况都飘回头查力合成公式的符号问题。第三步给一个中等幅值的阶跃转向角比如20度的阶跃输入观察车辆是否出现合理的转向响应。此时如果vy和r的响应有明显的高频振荡大概率是代数环或者步长问题。第四步开始做漂移激励。一种很实用的开环激励方式在1秒内给后轮一个很大的驱动力矩比如2000Nm同时把前轮转角打到-15度反打方向看模型能否自然进入持续滑动状态。如果整车直接掉头说明驱动力矩太大或者初始速度不够需要降低力矩或者提高初速。第五步把仿真数据和实车实验或者公开数据集里的极限工况数据对比调轮胎参数。漂移工况最敏感的参数是后轮的C值和E值它们决定了大侧偏角下的侧向力下降速率对漂移能否稳定维持影响最大。5.4 一个值得记录的仿真坑最后分享一个我实际遇到的、比较隐蔽的坑。做漂移维持的闭环仿真时发现横摆角速度一直有大约0.2 Hz的振荡怎么调控制器都不管用。后来排查了大半天发现是轮胎模型里的摩擦椭圆修正系数在作怪——我写的修正系数在前一步的Fx和Fy之间绑定得太紧导致两个力之间形成了一个自激振荡的小闭环。解决方式是在摩擦椭圆修正函数里做一阶低通滤波把高频增益降下来振荡就消失了。这个问题的启示是轮胎模型的输出不仅要“准”还要“平滑”。漂移工况下轮胎力处于饱和边界非线性很强可能在很小的输入变化下产生很大的输出跳变。模型输出的高频毛刺会和车身动力学形成耦合振荡所以建议在轮胎力输出的环节加一阶低通滤波器时间常数在5到20ms之间。这不是为了掩盖问题而是模拟真实轮胎的松弛特性——轮胎本身就是一个低通滤波器这一步在机理上也是讲得通的。我个人现在搭极限工况模型时已经习惯性地把轮胎松弛模型也一并考虑进去最简单的做法就是给轮胎力输出的每个力加一步低通滤波时间常数按前后轮分别取10ms和15ms左右效果不错。如果你按这个思路搭完了模型下一步就可以开始做漂移平衡点分析、线性化和闭环控制器设计了。那一部分牵扯到的东西更多留着下一篇再写。