ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

半车悬架模型Simulink仿真:从微分方程到工程应用

半车悬架模型Simulink仿真:从微分方程到工程应用 做车辆动力学仿真的人一般最早接触的都是四分之一悬架模型——一个簧上质量、一个簧下质量加一组弹簧阻尼简单归简单但真拿它去分析俯仰、制动点头、双轴路面输入就明显不够用了。所以当工程师把研究对象放在二分之一车辆悬架半车模型上时用Simulink搭建仿真平台几乎是最顺手的路径。半车模型在Simulink里既能保留足够的动力学细节又比整车模型清爽很多调试起来不费劲特别适合做悬架参数匹配、主动控制算法验证和俯仰特性研究。这篇文章我打算把半车悬架模型的Simulink仿真拆开讲透从为什么选半车模型、运动方程怎么推导到模块怎么搭、路面激励怎么做、结果怎么看、常见坑怎么避最后聊一下怎么往CarSim联合仿真和代码生成方向扩展。如果你手头正好要交一个悬架仿真课题或者刚开始接触半车模型这篇可以直接照着抄作业。1. 从四分之一到半车建模维度的一次关键升级1.1 半车模型的物理结构与自由度半车模型又叫二分之一车辆模型本质上是从整车纵向对称面上“切一刀”把前轴和后轴两侧的悬架等效到同一个平面里。简化后的结构包含一个车身俯仰刚体、两个簧下质量前轴和后轴各一个以及对应两组悬架弹簧、减振器和轮胎等效刚度。自由度方面常见的半车模型有四个车身垂向位移、车身俯仰角、前悬架簧下质量垂向位移、后悬架簧下质量垂向位移。如果你还要考虑纵向前后方向的耦合可以再加纵向自由度但绝大多数悬架性能研究不需要加了反而干扰主要结论。四个自由度对应的就是四阶微分方程组状态变量不多用Simulink搭起来层级清楚后期如果要改成六自由度整车模型结构上也能平滑过渡。这里有一个很多人会忽略的点半车模型的几何参数比如轴距、质心到前后轴的距离直接影响俯仰惯量的分配。你用的俯仰转动惯量如果是从整车参数里反推的一定要结合质心位置修正不能直接拿经验值硬套否则仿真出来的俯仰角响应会产生明显的相位偏差。1.2 为什么要用半车而不是整车讨论这个问题之前先看各自能干什么。四分之一模型只能研究单轮垂向跳动它对悬架刚度、阻尼匹配的评估是有用的但完全看不出俯仰动态也无法处理前后轴关联的输入激励。整车模型倒是信息全但所需的参数数量非常多轮胎、转向、横向稳定杆、衬套刚度、质心高度、转动惯量张量——每一个参数的误差都可能让仿真结果偏离真实车调试周期很长。半车模型正好落在中间它保留了前后轴的耦合关系能够反映制动或加速时的俯仰运动、不同车速下经过凸起路面时前后轮依次响应的过程也能评估悬架动行程和轮胎动载荷的前后分配情况。同时参数数量可控前悬架刚度、后悬架刚度、减振器阻尼系数这些核心参数基本都能从悬架设计手册或者台架试验数据里拿到。我在实际做课题的时候习惯把半车模型当作“金字塔的中间层”先用四分之一模型做初步参数扫描确定刚度阻尼的大致区间然后升级到半车模型验证俯仰表现最后如果项目需要再扩展成整车模型做精细化调校。这样比直接跳到整车模型效率高得多排查问题也容易定位是前悬还是后悬的问题。2. 数学模型与Simulink模块搭建2.1 运动微分方程推导建模先列方程。半车模型常用的坐标定义是以车辆静止平衡位置为原点车身质心垂向位移记为 $z_c$俯仰角记为 $\theta$前轴簧下质量垂向位移记为 $z_1$后轴记为 $z_2$。前轴到质心的距离是 $a$后轴到质心的距离是 $b$轴距 $Lab$。先看车身。车身受到前、后悬架力的作用垂向运动方程是$$m_c \ddot{z}c F{s_1} F_{s_2}$$其中 $F_{s_1}$ 和 $F_{s_2}$ 是前、后悬架传递给车身的力。这里方向要注意按照通常的符号约定压缩方向为正我习惯在Simulink里通过增益模块和加法器直接整理成标准形式避免符号混乱。俯仰运动方程则要考虑前后悬架力对质心的力矩$$I_y \ddot{\theta} -a F_{s_1} b F_{s_2}$$$I_y$ 是车身俯仰转动惯量。注意力矩方向前悬架在质心前方后悬架在质心后方所以两个力矩符号相反。这一点在搭建时非常容易出错我在初学时曾经把符号写反结果仿真出来的俯仰角完全反相查了一个多小时才发现是加减号的问题。前、后悬架力本身由弹簧和减振器产生。设前悬架弹簧刚度为 $k_{s1}$减振器阻尼系数为 $c_{s1}$悬架的相对位移由车身前轴位置与簧下质量位移决定。车身前轴处的垂向位移可以表示为$$z_{c1} z_c - a \sin\theta \approx z_c - a\theta$$小角度假设下 $\sin\theta \approx \theta$线性化处理在常规悬架分析中足够精确。于是前悬架力$$F_{s_1} k_{s1}(z_{c1} - z_1) c_{s1}(\dot{z}_{c1} - \dot{z}_1)$$后悬架同理$$F_{s_2} k_{s2}(z_{c2} - z_2) c_{s2}(\dot{z}_{c2} - \dot{z}_2)$$其中 $z_{c2} z_c b\theta$。最后是前、后簧下质量 $m_1$ 和 $m_2$ 的垂向运动方程。它们一方面通过悬架与车身相连另一方面通过轮胎等效刚度 $k_{t1}$、$k_{t2}$ 与路面接触$$m_1 \ddot{z}1 -F{s_1} k_{t1}(z_{r1} - z_1)$$$$m_2 \ddot{z}2 -F{s_2} k_{t2}(z_{r2} - z_2)$$$z_{r1}$、$z_{r2}$ 分别是前、后轮的路面激励位移。这四个方程合在一起构成一个典型的二输入前、后路面激励、多输出车身垂向加速度、俯仰角加速度、悬架动行程、轮胎动载荷的线性系统。建模到这一步Simulink里的事情就变成“把微分方程翻译成积分器和加减运算”。2.2 模块化搭建思路Simulink里搭建半车模型有两种路线一是完全用积分器、增益、加法器逐项搭建优点是可以直观对照微分方程适合教学和初学阶段二是用状态空间模块直接把方程整理成 $\dot{x}AxBu$ 的形式一个状态空间模块搞定优点是结构极其简洁适合参数化扫描。我最开始做课题时用的是第一种因为要反复检查方程有没有写错积分器连线和微分方程一一对应出问题好查。后期做优化迭代才改成了状态空间形式。下面分别说一下每条路线的要点。用积分器搭建时推荐的做法是把状态变量设为 $z_c$、$\theta$、$z_1$、$z_2$以及它们的导数然后在每个积分器前汇集对应的加速度表达式。以车身垂向加速度为例在加法器里把 $F_{s1}$ 和 $F_{s2}$ 相加后除以 $m_c$再送进积分器得到 $\dot{z}_c$再积一次得到 $z_c$。悬架力的计算则用增益模块乘以相对位移和相对速度可以用Bus Creator把相关信号打包层次更清晰。用状态空间模块时状态向量取 $x[z_c, \theta, z_1, z_2, \dot{z}_c, \dot{\theta}, \dot{z}1, \dot{z}2]^T$控制输入 $u[z{r1}, z{r2}]^T$。系统矩阵 $A$ 和输入矩阵 $B$ 可以直接从方程中整理出来。我建议用MATLAB脚本写一个初始化脚本把参数赋值和矩阵组装全部放到脚本里每次改参数只需要改脚本头部的变量不用动模型文件。2.3 参数选取与初始化脚本参数是仿真能否靠谱的基础。以一辆典型中型轿车为例我常用的参数如下簧上质量车身等效质量$m_c$ 950 kg俯仰转动惯量 $I_y$ 1200 kg·m²前悬架刚度 $k_{s1}$ 28000 N/m后悬架刚度 $k_{s2}$ 30000 N/m前减振器阻尼 $c_{s1}$ 2200 N·s/m后减振器阻尼 $c_{s2}$ 2100 N·s/m前簧下质量 $m_1$ 45 kg后簧下质量 $m_2$ 40 kg前轮胎刚度 $k_{t1}$ 260000 N/m后轮胎刚度 $k_{t2}$ 240000 N/m前轴到质心距离 $a$ 1.2 m后轴到质心距离 $b$ 1.4 m这里有个经验值得分享悬架刚度并不是前后随便选的前后悬架的刚度比要和质心位置匹配也就是要让前后悬架的偏频接近。前悬架偏频 $f_1 \frac{1}{2\pi}\sqrt{\frac{k_{s1}}{m_{c1}}}$后悬架偏频 $f_2 \frac{1}{2\pi}\sqrt{\frac{k_{s2}}{m_{c2}}}$其中 $m_{c1}$ 和 $m_{c2}$ 是根据质心位置分配到前、后轴上的车辆质量。如果前后偏频相差太多车辆路过障碍后会出现明显的俯仰振荡体感很差这也是半车模型能暴露而四分之一模型看不见的问题。初始化脚本我会写成下面这样的结构放到模型的InitFcn回调里或者直接放在脚本中先于sim()执行% 半车悬架模型参数初始化 mc 950; % 簧上质量, kg Iy 1200; % 俯仰转动惯量, kg*m^2 a 1.2; b 1.4; % 质心到前/后轴距离, m L a b; % 轴距, m ks1 28000; ks2 30000; % 前后悬架刚度, N/m cs1 2200; cs2 2100; % 前后减振器阻尼, N*s/m m1 45; m2 40; % 前后簧下质量, kg kt1 260000; kt2 240000; % 前后轮胎刚度, N/m % 状态空间矩阵组装 % 状态: [zc theta z1 z2 zc_dot theta_dot z1_dot z2_dot] A zeros(8,8); A(1,5) 1; A(2,6) 1; A(3,7) 1; A(4,8) 1; % 具体元素按方程填入 % ... B zeros(8,2); % 路面输入映射到轮胎力项 % ...模型里到前、后轮的路面激励通常有时间延迟关系同一路面轮廓后轮要比前轮晚 $L/v$ 到达。在Simulink里可以用Transport Delay模块实现输入设为 $z_{r1}$延迟时间设为 $L/v$得到 $z_{r2}$。这个细节特别重要如果你前后轮用同一个路面输入模拟出来的俯仰运动会完全失真因为前后轮同时压过障碍和先后压过障碍明明是两种完全不同的工况。3. 路面激励与仿真工况设计3.1 随机路面生成悬架仿真的输入不外乎两种随机路面和确定性障碍。随机路面用来评估平顺性确定性障碍用来看特定工况的瞬态响应比如过减速带、过坑。随机路面常用滤波白噪声法生成也就是对白噪声进行整形滤波逼近标准路面功率谱密度。国标GB/T 7031或ISO 8608把路面等级按谱密度大小分成A到H级B、C级是常见公路。这里的关键是生成一条时域的路面不平度序列再用空间采样转成时间序列。假设车速 $v 20 \mathrm{m/s}$空间频率 $n_00.1 \mathrm{m^{-1}}$路面不平度系数 $G_q(n_0)64\times10^{-6} \mathrm{m^3}$C级路面。滤波白噪声的时域模型可以写成$$\dot{z}_r(t) -2\pi n_0 v \cdot z_r(t) 2\pi n_0 \sqrt{G_q(n_0) v} \cdot w(t)$$其中 $w(t)$ 是白噪声。在Simulink里实现这个公式的常用办法是搭一个闭环子系统白噪声源接到增益再经过一阶惯性环节。也可以用MATLAB预先生成好时间序列然后用Repeating Sequence Interpolate模块或者From Workspace模块导入。我的习惯是后者因为能看到生成的路面曲线心里有底。生成随机路面的脚本参考v 20; % 车速, m/s Gq0 64e-6; % C级路面不平度系数 n0 0.1; % 参考空间频率 dt 0.001; % 采样步长 T_end 20; % 仿真时长 N T_end / dt; rng(2024); w randn(N, 1) / sqrt(dt); % 离散白噪声 zr zeros(N, 1); for k 2:N zr(k) zr(k-1) * (1 - 2*pi*n0*v*dt) ... 2*pi*n0*sqrt(Gq0*v) * w(k) * dt; end t (0:N-1) * dt;注意白噪声在离散域的方差要和连续域匹配所以除以 $\sqrt{dt}$否则仿真结果会随步长变化。3.2 减速带凸块输入确定性障碍我常用两种正弦半波凸块和矩形凹坑。正弦半波减速带适合做频率特性的直观验证其表达式为$$z_r(t) \begin{cases} h \sin\left(\frac{2\pi v}{L_p} t\right) 0 \le t \le \frac{L_p}{v} \ 0 \text{otherwise} \end{cases}$$$h$ 是减速带高度$L_p$ 是减速带宽度。比如 $h0.05\mathrm{m}$$L_p0.3\mathrm{m}$车速 $v10\mathrm{m/s}$则前轮经过的时间只有0.03秒后轮在此基础上延迟 $L/v0.26$ 秒。这个时间差造成的俯仰激励是半车模型特别典型的工况。在Simulink里可以用Clock模块、比较器和Switch组合出这个波形也可以直接用Signal Builder或Signal Editor画出来。我建议用Signal Editor因为它能直接编辑分段波形还能导出改起来方便。3.3 求解器配置求解器这块容易被忽视但实际影响非常大。半车模型是一个刚性程度适中的动力学系统轮胎刚度比悬架刚度大一个数量级所以高频分量不小。我的经验是做随机路面平顺性分析时用ode45配合最大步长限制在0.001秒以内精度够速度也还好做减速带这类瞬态冲击工况时改用ode15s因为短时间内的冲击容易让ode45自适应步长变得很激进反而拉长计算时间如果只是验证模型是否稳定可以先用Fixed-step-discrete步长0.001秒跑20秒看看曲线趋势。还有一点Transport Delay模块是连续时间模块在固定步长离散求解器下可能会产生延迟量的量化误差。如果后轮激励延迟时间不是步长的整数倍仿真波形会出现微小的台阶。解决方法是把延迟放到路面生成阶段统一处理先算出前轮时间序列再按样本数精确偏移 $N_d \mathrm{round}(L/v/dt)$ 得出后轮序列然后用From Workspace导入。这样做出来的前后轮激励天然同步也省掉了Transport Delay带来的数值麻烦。4. 结果分析垂向、俯仰与悬架行程的联合解读4.1 时域响应曲线判读Simulink模型跑完第一件事不是左看右看那些花花绿绿的Scope输出而是先把关键信号导出到工作区用MATLAB脚本统一处理。我的标准流程是把车身垂向加速度、俯仰角加速度、前后悬架动行程、前后轮胎动载荷这六路信号全部导出存成结构体或时间表然后写脚本集中画图和分析。看曲线的时候优先关注三个点第一车身垂向加速度的稳态幅值。随机路面激励下垂向加速度均方根值如果超过 $2.5 \mathrm{m/s^2}$体感已经开始不舒服了。这个值可以作为悬架参数好坏的初步判断依据。第二俯仰角的变化幅度和振荡收敛速度。过减速带工况下俯仰角峰值如果超过3度说明前后悬架阻尼协调性不够好如果振荡经过四五次衰减还没稳定下来可能是后悬架阻尼偏小。第三悬架动行程有没有撞限位。悬架动行程是悬架相对位移如果它的峰值超过了设计行程限位值说明悬架刚度偏低或者阻尼匹配不当实车上就会出现悬架击穿。半车模型虽然不直接模拟限位块但你可以根据动行程峰值提前判断会不会接近限位。我手动处理时习惯把信号分成两组画图一组是车身状态垂向加速度、俯仰角一组是悬架和轮胎状态动行程、轮胎动载荷。上下对称排版一眼就能看出参数调整对舒适性和操稳性的影响。4.2 悬架评价指标计算仿真只是过程最后还是要落到指标上。悬架性能评价有几个常用指标车身垂向加速度均方根值反映乘坐舒适性。俯仰角加速度均方根值反映俯仰抑制能力。悬架动行程均方根值和峰值反映悬架行程利用率。轮胎动载荷均方根值和峰值反映车轮接地性动载荷峰值超过静载就意味着车轮可能离地。这些指标可以在MATLAB脚本里直接计算az_rms sqrt(mean(az.^2)); theta_rms sqrt(mean(theta.^2)); susp_rms_front sqrt(mean(susp_f.^2)); DTL_front_max max(dtl_f);做参数优化时把这四个指标当作多目标来权衡。舒适性和操稳性通常是矛盾的阻尼调大车身加速度和俯仰抑制会变好但悬架动行程和轮胎动载荷会恶化。半车模型的优势恰恰在于能看到这种矛盾在前后轴之间的分配关系。这里我想说一个很实用的技巧在分析时把前后悬架动载荷和动行程的相位差也看一眼。前后轴响应之间的相位关系决定车辆整体的俯仰姿态变化模式。如果后悬架动行程峰值明显晚于前悬架而且幅值较大说明后悬架的相对阻尼可能要重新调整。5. 仿真过程中的典型问题与排查5.1 代数环问题代数环是Simulink里做物理建模时最常见的问题之一。半车模型中如果你把悬架力直接建模为相对位移和相对速度的代数函数同时又把悬架力反馈到加速度计算那么在不落地建模的情况下很容易出现一个“输出依赖输入、输入又依赖输出”的闭环。出现代数环时Simulink会在诊断窗口里提示而且仿真速度会变慢甚至出现数值震荡。消除代数环的办法有几个一是用积分器把环路打破。物理上加速度先积分得到速度速度再积分得到位移力由位移和速度计算后反馈给加速度。这个“积分延迟”天然打破代数环。所以用积分器组搭建时代数环基本不会出现反而是你用状态空间模块直接组装代数方程时一定要确保没有把增益模块首尾相接。二是如果必须保留代数结构可以在模型中插入Unit Delay或者Memory块但这个方法会引入一步延迟误差对精度要求高的课题不建议用。三是检查模型里有没有把同一个信号既作为输出又作为输入比如把悬架力同时送到受力端和测量端时误连了信号线。这种低级错误也会触发代数环提示。我用状态空间模块搭建时从来不直接写矩阵A里的代数关系而是把 $\dot{z}$ 变量显式放在状态里A矩阵中包含位移和速度的耦合这样可以保证状态矩阵是严格的常微分方程形式从根上避开代数环。5.2 仿真发散原因与对策仿真发散是另一个高频问题。现象是曲线一开始很正常跑了零点几秒后突然指数式飞掉或者NaN、Inf直接出现。最常见的原因是参数数量级差异过大。轮胎刚度是 $2\times10^5$ 级别悬架刚度是 $2\times10^4$ 级别车身质量是 $10^2$ 级别这些数字混在一起如果积分步长太大收敛性会出问题。解决办法是把步长调到0.001秒以下同时检查是否存在刚性过大的环节。第二个原因是初始条件设得不合适。模型刚搭建时微分方程里的积分器初始值如果设成0而对应的平衡位置又不是0系统会先经历一个很大的瞬态调整如果激励叠加在这个瞬态上很容易超出数值范围。解决办法是先让模型空跑一次不加路面激励等平衡稳定后再接入激励。更好的做法是直接计算出平衡状态下的弹簧压缩量把它作为积分器初始值。第三个原因是阻尼被设成负值。减振器阻尼系数虽然在物理上不允许为负但参数扫描时如果扫描区间不小心包含负值仿真必然发散。排查这类问题我习惯在初始化脚本里加assert检查比如assert(cs1 0 cs2 0, 阻尼系数必须为正); assert(all([ks1 ks2 kt1 kt2] 0), 刚度参数必须为正);第四路面激励的幅值或者频率超出模型的线性假设范围。滤波白噪声生成的随机路面如果车速过高路面输入的高频分量会被放大线性悬架模型弹簧刚度不变、阻尼特性线性可能给出不合理的大幅响应。这种情况下要意识到不是模型错了而是线性模型适用范围到了。你可以先降低车速或者采用小一些的路面等级或者干脆把模型升级为非线性悬架模型比如加入变阻尼特性。5.3 单位与符号约定混乱单位问题看似简单实际统计下来大部分仿真结果对不上根子都在单位换算上。最常见的是把毫米当成米或者把千牛当成牛。轮胎刚度如果写成260而不是260000模型就会变成“轮胎是软的”悬架静平衡状态完全不对。符号约定方面我见过团队内部因为正方向定义不同导致最后联合调试时完全对不上数据。做半车模型之前最好先约定一个统一的符号规范并用注释写在模型里。我的习惯是垂向位移向上为正俯仰角抬头为正弹簧压缩力为正。所有力和位移的符号在方程推导阶段就固定下来再到Simulink的每个增益模块里检查一遍。检查符号是否正确有一个快速方法用重力阶跃输入测试车身加速度方向。垂直系统在中性位置自然下落如果车身加速度方向和重力方向一致说明符号是对的如果方向反了车会像被吸上去一样飞离路面这个错误在Scope里非常明显。6. 从Simulink到工程应用扩展与进阶6.1 悬架参数优化与DoE分析模型搭建完成且验证可靠后最常见的工作就是把悬架刚度和阻尼参数优化一遍。半车模型状态变量少、仿真速度快非常适合做基于仿真的参数优化。我自己常用的流程是确定设计变量比如前悬架阻尼系数 $c_{s1}$ 和后悬架阻尼系数 $c_{s2}$给定变化范围然后写一个循环脚本批量调用sim()函数采集每一组参数下的性能指标最后画Pareto前沿图综合权衡舒适性和操稳性后选出折中参数。这里有个小经验批量仿真前先把Simulink模型的Rapid Accelerator模式打开配合并行计算工具可以大幅缩短仿真时间。半车模型本身很轻量100组参数仿真通常也就几分钟但如果你把整车模型拿来做同样的事时间成本会上一个量级。6.2 主动悬架与半主动控制接入半车模型是控制算法验证的好载体。你可以保留机械悬架模型把前、后悬架力替换成可主动调节的力作动器表达式比如$$F_s k_s(z_{c1}-z_1) c_s(\dot{z}_{c1}-\dot{z}1) F{act}$$半主动悬架则是把阻尼系数做成实时可变的控制量常见的天棚阻尼Skyhook控制策略落在Simulink里就是% Skyhook控制伪代码 if (zdot_c * (zdot_c1 - zdot_1) 0) c_sky c_max; else c_sky c_min; end用半车模型做控制验证比用四分之一模型更能看出控制策略对俯仰动态的影响。天棚阻尼在四分之一模型上表现良好但放在半车模型上有时会导致前后悬架控制力矩相互干扰产生新的俯仰模态。这种跨轴耦合效应是半车模型独有的价值。6.3 CarSim联合仿真与代码生成最后聊一下扩展方向。半车模型稳定之后如果你还想让路面模型、轮胎模型更真实可以走CarSim与Simulink联合仿真的路线。CarSim提供高精度的整车动力学环境但它内部的悬架详细参数未必方便修改Simulink这边则适合做控制算法快速开发。典型的做法是在CarSim里搭建整车模型把六个悬架力或控制作动器接口暴露给Simulink这样既能保留Simulink的控制算法开发便利性又获得了CarSim的整车级精度。另一个方向是Simulink模型生成C代码。半车模型本身是一个线性多输入多输出系统生成嵌入式C代码后可以部署到硬件在环HIL仿真台架或者快速原型控制器上。用Simulink的嵌入式编码器生成代码再编写简单的I/O驱动程序就能在原车控制单元类似的硬件环境中测试悬架控制策略的实时性。这一步在院校课题和企业预研阶段都很加分。我记得有一次在硬件在环平台上调试半主动悬架控制策略代码生成后跑出来和离线仿真结果有一点偏差排查半天发现是浮点类型转换问题——模型默认双精度嵌入式硬件默认单精度关键信号被截断后控制效果就变了。所以如果你要走代码生成路线早期就要把模型里的数据类型显式设置为单精度或定点数并在Simulink里启用数据字典统一管理。我个人在实际操作中的体会是半车模型这个中间层做透之后的复用价值极高。模型架构、驱动脚本、后处理指标计算这些代码几乎可以原样迁移到后续的整车模型和联合仿真项目里。你唯一需要改的只是把矩阵维数从八乘八扩出来加上横向自由度本质方法没变。所以别觉得半车模型“简单”就不值得花时间反而是把这个层级的细节抠清楚后面的路会顺很多。
返回列表