简介:面向车辆工程学习者与悬架动力学入门研究者的MATLAB计算材料,围绕二自由度悬架模型,用传递函数刻画路面激励到车身振动的频域关系,并支持多阻尼参数对比。压缩包仅含1个.m脚本文件,体量约1KB,核心代码聚焦建模、传递函数求解与时域响应计算,便于在MATLAB中直接运行、修改参数并观察曲线变化。资源页显示已有1598人学习,适合作为课程作业、教学演示、论文复现或悬架性能对比实验的轻量参考。读者可据此核对车轮质量、车身质量、弹簧刚度与阻尼系数等参数设置,理解不同阻尼对共振峰、振幅衰减和振动抑制效果的影响,并以此为模板扩展频响曲线与时域响应的对比分析,形成可复用的悬架动态特性评估思路。
1. 二自由度悬架传递函数:多阻尼对比到底在比哪条曲线
车辆悬架的二自由度模型里,传递函数不是一张装饰图,它决定了你换阻尼系数后曲线到底怎么动。这份资源围绕“二自由度悬架传递函数-多阻尼对比”展开,把簧上/簧下质量、悬架刚度、轮胎刚度、阻尼系数放进同一套方程,输出车身位移、车身加速度、悬架动挠度和轮胎动载荷的频响,再对多个阻尼系数做批量扫描。它适合做车辆动力学课程设计、悬架参数初选、频响图复现的人,也适合已经有一版模型但不确定符号和单位有没有写错的从业者。多阻尼对比的价值不在多画几条线,而在于把“阻尼变大后哪个频段被压下去、哪个指标反而变差”摊开看,避免只盯一个峰值就下结论。
2. 从动力学方程到H(s):二自由度悬架传递函数的代码落地
2.1 模型假设与参数口径先对齐
这套资源用的是四分之一车辆二自由度模型:簧上质量代表车身,簧下质量代表车轮和部分悬架部件,悬架弹簧和阻尼连接两者,轮胎简化成径向弹簧接地。常见假设是轮胎阻尼忽略、悬架刚度线性、路面输入为位移激励、只研究垂向运动。这个模型不处理俯仰和侧倾,也不适合直接评价整车舒适性,但用来比较不同阻尼系数下的频响形态足够。
参数口径必须先定死。下面这组量级只用于跑通代码,不代表任何具体车型:
| 符号 | 含义 | 常用量级 | 单位 |
|---|---|---|---|
| m1 | 簧下质量 | 30~60 | kg |
| m2 | 簧上质量 | 250~450 | kg |
| k1 | 轮胎径向刚度 | 150000~250000 | N/m |
| k2 | 悬架刚度 | 20000~40000 | N/m |
| c | 悬架阻尼系数 | 500~3500 | N·s/m |
最容易翻车的是把 k1 和 k2 写反,或者把 c 的单位写成 N·s/mm。前者会让车轮共振频率跑到不合理的位置,后者会让阻尼效果差 1000 倍。我一般先算两个粗估频率:车身频率约 sqrt(k2/m2)/(2π),车轮频率约 sqrt((k1+k2)/m1)/(2π)。如果脚本跑出来的峰值离这两个量级很远,先别怀疑算法,回头查单位和符号。
2.2 传递函数推导与可复用实现
设车身位移为 z2,车轮位移为 z1,路面位移为 zr,悬架阻尼系数为 c,悬架刚度为 k2,轮胎刚度为 k1。忽略轮胎阻尼后,方程写成:
m2 z2'' + c(z2' - z1') + k2(z2 - z1) = 0
m1 z1'' + c(z1' - z2') + k2(z1 - z2) + k1(z1 - zr) = 0
整理成矩阵形式,令 s 为拉普拉斯变量:
A11 = m2 s^2 + c s + k2
A12 = -(c s + k2)
A22 = m1 s^2 + c s + k2 + k1
车身位移对路面位移的传递函数为:
H_z2(s) = k1(c s + k2) / [A11 A22 - A12^2]
车轮位移传递函数为:
H_z1(s) = k1 A11 / [A11 A22 - A12^2]
悬架动挠度等于 z2 - z1,传递函数为:
H_def(s) = -k1 m2 s^2 / [A11 A22 - A12^2]
下面这段代码把三个传递函数封装成可复用的函数,分子分母都用多项式系数表示,方便后面直接丢给 scipy 做频响:
import numpy as np def quarter_car_tf(m1, m2, k1, k2, c): """ 二自由度四分之一车辆模型。 返回车身位移、车轮位移、悬架动挠度对路面位移的传递函数。 多项式系数按 s 的高阶到低阶排列。 """ # A11 = m2*s^2 + c*s + k2 A11 = np.array([m2, c, k2], dtype=float) # A22 = m1*s^2 + c*s + k2 + k1 A22 = np.array([m1, c, k2 + k1], dtype=float) # A12 = -(c*s + k2) A12 = np.array([-c, -k2], dtype=float) # detA = A11*A22 - A12^2 det = np.polysub(np.polymul(A11, A22), np.polymul(A12, A12)) # H_z2 = k1*(c*s + k2) / detA num_z2 = k1 * np.array([c, k2], dtype=float) # H_z1 = k1*A11 / detA num_z1 = k1 * A11 # H_def = H_z2 - H_z1 = -k1*m2*s^2 / detA num_def = -k1 * m2 * np.array([1.0, 0.0, 0.0]) return { "det": det, "z2": (num_z2, det), "z1": (num_z1, det), "deflection": (num_def, det), }这段代码的关键不是把公式抄一遍,而是把 A12 的负号单独放进多项式。A12 实际是 -(c s + k2),平方后变成正号,如果手写时漏掉负号,行列式会变成 A11 A22 + A12^2,低频增益立刻不对。判断方法很直接:路面位移在 0 Hz 时是缓慢抬升,车身和车轮位移的传递函数在 s=0 处应该趋近 1,悬架动挠度应该趋近 0。跑完函数先做这三个检查,比看曲线快。
参数上,m1、m2 用 kg,k1、k2 用 N/m,c 用 N·s/m。函数里没有做单位换算,这是故意的:单位换算放在参数表入口,模型内部只认 SI。如果从文档里拿到 N/mm,先乘 1000 再传进来。我一般会在调用前加一行断言,把明显不合理的值拦住。
2.3 频响计算与阻尼初扫
传递函数是 s 域表达式,工程上更常看频率响应:令 s = jω,ω = 2πf,扫 0.1~20 Hz。下面的代码用 scipy.signal.freqs 计算车身加速度、悬架动挠度和轮胎动载荷的复频响。注意车身加速度不是直接对 H_z2 取峰值,而是乘上 s^2,即 -(2πf)^2 H_z2。这个细节在避坑章还会展开。
import numpy as np from scipy.signal import freqs def response_curves(params, f_hz): """ 计算单位路面位移激励下的频响。 返回车身加速度、悬架动挠度、轮胎动载荷的复频响。 """ tf = quarter_car_tf(**params) w = 2.0 * np.pi * f_hz # 车身位移频响 _, h_z2 = freqs(*tf["z2"], worN=w) # 车轮位移频响 _, h_z1 = freqs(*tf["z1"], worN=w) # 悬架动挠度频响 _, h_def = freqs(*tf["deflection"], worN=w) # 加速度 = s^2 * 位移,s = jw,所以乘 -(w^2) h_acc = -(w ** 2) * h_z2 # 轮胎动载荷 = k1*(z1 - zr),单位路面位移下就是 k1*(H_z1 - 1) h_load = params["k1"] * (h_z1 - 1.0) return { "acc": h_acc, "deflection": h_def, "load": h_load, "z2": h_z2, "z1": h_z1, }这里的参数说明要分清“输入量”和“观测量”。f_hz 是频率轴,常用 0.1~20 Hz 对数分布,800 点到 2000 点足够看峰值;点太少会把窄峰削平。h_acc 的单位不是米,而是“单位路面位移下的加速度增益”,做相对比较没问题,不要直接当成 m/s^2 的绝对舒适性指标。h_load 同样是在单位路面位移下的轮胎动载荷增益,乘以实际路面谱后才更接近工程评价。多阻尼对比阶段看的是趋势和峰值位置,不是拿单条曲线下绝对结论。
初扫时先取三档阻尼,比如 800、1500、3000 N·s/m,跑同一组频率。把三条车身加速度曲线画在同一张对数频率轴上,先看 1~2 Hz 附近的车身共振峰有没有被压下去,再看 10~13 Hz 附近的车轮共振有没有被抬起来。如果曲线形态和预期相反,先查 H_z2(0) 是否等于 1,再查加速度有没有乘 s^2。这个顺序能省很多时间。
3. 多阻尼对比实验:阻尼系数扫描、指标与曲线解读
3.1 阻尼系数与阻尼比的换算边界
多阻尼对比最容易写成“c 越大越好”或者“c 越小越舒服”,这两种说法都不稳。悬架阻尼系数 c 影响的是车身和车轮两个模态,单自由度阻尼比公式 zeta = c/(2 sqrt(k2 m2)) 只能当粗估。它把车轮侧当成固定基础,忽略了车轮模态和轮胎刚度,所以低阻尼时还能看,阻尼上来后预测的峰值位置会有偏差。
更稳的做法是把方程写成状态空间,求系统矩阵特征值,从复特征值里取模态频率和模态阻尼比。下面这段代码不依赖传递函数,直接对 M、C、K 矩阵做特征值分解:
import numpy as np def modal_damping(m1, m2, k1, k2, c): """ 通过状态空间特征值求两个垂向模态的频率和阻尼比。 返回按频率升序排列的 fn(Hz) 和 zeta。 """ M = np.array([[m2, 0.0], [0.0, m1]]) C = np.array([[c, -c], [-c, c]]) K = np.array([[k2, -k2], [-k2, k2 + k1]]) Minv = np.linalg.inv(M) A = np.block([ [np.zeros((2, 2)), np.eye(2)], [-Minv @ K, -Minv @ C] ]) eig = np.linalg.eigvals(A) # 取正虚部的共轭根,避免重复 pos = eig[eig.imag > 0] fn = np.abs(pos) / (2.0 * np.pi) zeta = -pos.real / np.abs(pos) idx = np.argsort(fn) return fn[idx], zeta[idx]这段代码给出的模态阻尼比比单自由度公式更接近二自由度模型的真实响应。参数仍然用 SI,M 矩阵里 m2 在前、m1 在后,和状态变量顺序一致。若把顺序写反,特征值频率不会变,但输出顺序和传递函数对不上,画图时会误判哪个峰是车身模态。我的习惯是把模态频率打印出来,和 2.1 节的粗估频率对一眼,再继续扫描。
3.2 批量扫描脚本与三类评价指标
多阻尼对比不要手工改一个参数跑一次。下面这段扫描脚本把阻尼列表、频率轴和三类指标放在一个循环里,输出可以直接转成表格。三类指标分别是车身加速度峰值、悬架动挠度均方根、轮胎动载荷均方根。这里的均方根按等带宽白噪声位移输入做相对比较,不是法规认证值,适合看趋势。
import numpy as np from scipy.signal import freqs def sweep_damping(base_params, c_list, f_hz): rows = [] for c in c_list: p = dict(base_params) p["c"] = float(c) tf = quarter_car_tf(**p) w = 2.0 * np.pi * f_hz _, h_z2 = freqs(*tf["z2"], worN=w) _, h_z1 = freqs(*tf["z1"], worN=w) _, h_def = freqs(*tf["deflection"], worN=w) h_acc = -(w ** 2) * h_z2 h_load = p["k1"] * (h_z1 - 1.0) # 等带宽相对指标:对功率谱做积分再除以带宽 df = f_hz[-1] - f_hz[0] acc_abs = np.abs(h_acc) def_rms = np.sqrt(np.trapz(np.abs(h_def) ** 2, f_hz) / df) load_rms = np.sqrt(np.trapz(np.abs(h_load) ** 2, f_hz) / df) rows.append({ "c": p["c"], "zeta_susp": p["c"] / (2.0 * np.sqrt(p["k2"] * p["m2"])), "acc_peak": float(np.max(acc_abs)), "acc_peak_f": float(f_hz[np.argmax(acc_abs)]), "def_rms": float(def_rms), "load_rms": float(load_rms), }) return rows调用时给 base_params 传 m1、m2、k1、k2、c 的字典,c 会被覆盖。f_hz 建议用 np.logspace(-1, np.log10(20), 1200)。acc_peak 看的是车身加速度最大幅值,acc_peak_f 告诉你峰值落在哪个频率;def_rms 反映悬架动挠度整体能量,太小可能意味着阻尼过硬;load_rms 反映轮胎动载荷整体能量,和接地性相关。三个指标经常互相拉扯:阻尼加大,车身共振峰值下降,但车轮高频段的动载荷可能上升。多阻尼对比的意义就是把这个拉扯关系摆出来。
3.3 曲线形态背后的物理含义
把扫描结果画成图,通常会看到两个明显区域。低频区 1~2 Hz 对应车身模态,阻尼增大时车身加速度峰值下降,悬架动挠度峰值也下降,这是最直观的收益。高频区 10~13 Hz 对应车轮模态,阻尼增大不一定继续变好,有时轮胎动载荷会升高,因为阻尼把更多运动传给了车轮侧。只看车身加速度会得出“阻尼越大越舒服”的片面结论,把轮胎动载荷和动挠度一起看,才会发现边界。
还有一条经验:阻尼系数扫描不要等间隔取 500、1000、1500 这种线性序列,最好在低阻尼段密一点,高阻尼段疏一点,或者直接按阻尼比取点。因为峰值变化在低阻尼段更敏感,到了高阻尼段曲线已经压平,线性等间隔会浪费很多点。我一般先按 zeta 从 0.1 到 0.7 取 12 个点,反算 c,再补几个边界值。这样做出来的峰值-阻尼曲线更容易看出拐点。
图表上,频率轴用对数,纵轴用 dB 或对数幅值都可以,但同一篇文章里不要混用。车身加速度、动挠度、轮胎动载荷分三个子图,共用频率轴,标注清楚单位是“单位路面位移下的增益”。如果只放一张归一化曲线,读者看不出量级,也不方便复核。曲线旁边最好附一张扫描表,把 c、zeta_susp、acc_peak、acc_peak_f、def_rms、load_rms 列出来,图看形态,表看数值。
4. 多阻尼对比排查:五个常见翻车点及修复
4.1 单位、符号与模型口径的坑
现象:车身位移传递函数在 0 Hz 不趋近 1,低频曲线直接漂了。
原因多半是 A12 的负号漏了,或者行列式写成了 A11 A22 + A12^2。二自由度模型的耦合项带负号,手推公式时很容易在平方那一步丢掉。解决办法是跑完 quarter_car_tf 后先算 s=0 的值:H_z2(0) 应接近 1,H_z1(0) 应接近 1,H_def(0) 应接近 0。三个检查不过,不要继续画图。
现象:车轮共振峰跑到 3 Hz 或 30 Hz,和粗估频率对不上。
原因常见于 k1 和 k2 写反,或者轮胎刚度用了 N/mm 却没乘 1000。轮胎刚度通常比悬架刚度大一个量级,写反后车轮模态频率会明显偏低。解决方式是先算 sqrt(k2/m2)/(2π) 和 sqrt((k1+k2)/m1)/(2π),把两个数写在脚本注释里,每次换参数先对量级。单位检查函数放在参数入口,别放在模型内部,否则后面批量扫描还会绕过。
现象:阻尼从 800 改到 3000,曲线几乎不动。
先查 c 的单位。N·s/m 和 N·s/mm 差 1000 倍,如果参数表里写的是 N·s/mm,直接传进模型会得到过阻尼或者几乎无阻尼的假象。另一个可能是扫描时只改了字典副本,实际调用仍用旧参数。解决方式是打印每次循环的 c 和 zeta_susp,确认扫描值真的进了 quarter_car_tf。
4.2 数值计算与图表表达的坑
现象:车身加速度峰值频率和车身位移峰值频率对不上,怀疑模型错了。
这不一定错。加速度是位移乘 s^2,幅值上多了一个 ω^2 权重,峰值会往高频偏移。若把 H_z2 的峰值直接当加速度峰值,就会得到偏低的频率。解决方式是单独封装 h_acc = -(w**2)*h_z2,画图时明确标注“加速度传递函数”,不要和位移传递函数共用一条曲线。
现象:多阻尼曲线看起来交叉混乱,结论反复。
常见原因是频率范围截断。只扫 0.5~8 Hz 会漏掉车轮模态,只扫 5~20 Hz 又会漏掉车身模态,两个频段分开看容易得出相反结论。解决方式是把频率轴固定为 0.1~20 Hz 对数分布,三个观测量的子图共用同一范围,并在峰值附近加密取点。还有一个小坑是纵轴用线性幅值又用 dB 混着画,视觉上峰值被压缩,误以为阻尼影响不大。同一批对比图统一口径,图注写清楚。
现象:扫描表里 acc_peak 很小,但图上看着峰很高。
检查是不是对复频响取了实部或者忘了取绝对值。np.max(h_acc) 对复数没有意义,应该用 np.max(np.abs(h_acc))。另一个可能是频率点太稀,峰值落在两个点之间被削平。把 f_hz 从 800 点加到 2000 点,或者用峰值附近局部加密,数值会稳定下来。这个坑很隐蔽,因为图看起来“差不多”,但表格里的峰值会偏差十几个百分点。
5. 传递函数结果复核:解析解、数值解与报告化输出
5.1 解析传递函数与状态空间频响交叉验证
一套多阻尼对比结果要能复核,最好用两条独立路径算同一个频响。第一条是第 2 章的传递函数多项式,第二条是状态空间模型。两条路径的方程来源相同,但数值实现不同,交叉验证能抓出符号和系数错误。下面这段代码构造状态空间,输出车身位移和车轮位移,再用 freqresp 和解析传递函数对比。
import numpy as np from scipy.signal import StateSpace, freqresp, freqs def state_space_model(m1, m2, k1, k2, c): """ 状态变量顺序:[z2, z1, v2, v1]。 输出:[z2, z1]。 """ M = np.diag([m2, m1]) C = np.array([[c, -c], [-c, c]], dtype=float) K = np.array([[k2, -k2], [-k2, k2 + k1]], dtype=float) Minv = np.linalg.inv(M) A = np.block([ [np.zeros((2, 2)), np.eye(2)], [-Minv @ K, -Minv @ C] ]) B = np.vstack([ np.zeros((2, 1)), Minv @ np.array([[0.0], [k1]]) ]) Cmat = np.array([ [1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0] ]) D = np.zeros((2, 1)) return StateSpace(A, B, Cmat, D) def check_tf_vs_ss(params, f_hz): tf = quarter_car_tf(**params) w = 2.0 * np.pi * f_hz _, h_z2_tf = freqs(*tf["z2"], worN=w) _, h_z1_tf = freqs(*tf["z1"], worN=w) sys = state_space_model(**params) _, h_ss = freqresp(sys, w) h_z2_ss = h_ss[0, 0, :] h_z1_ss = h_ss[1, 0, :] err_z2 = np.max(np.abs(h_z2_tf - h_z2_ss)) err_z1 = np.max(np.abs(h_z1_tf - h_z1_ss)) return err_z2, err_z1这段代码里 B 矩阵的输入项是 [0; k1],因为路面位移只直接作用在车轮方程上。状态变量顺序必须和 Cmat 对应,z2 在前、z1 在后,速度项不输出。err_z2 和 err_z1 在数值精度内应该很小,常见量级是 1e-10 到 1e-8。如果误差到了 1e-2,先查 Cmat 是否把 z2、z1 顺序写反,再查 B 矩阵是不是把 k1 放到了车身方程。交叉验证不需要画图,打印两个残差就够。
5.2 参数敏感性与批量导出
多阻尼对比做完后,最好再做一轮参数敏感性:阻尼固定几档,悬架刚度按比例缩放,观察峰值频率和峰值幅值怎么变。这样在报告里能说清楚“结论依赖哪几个参数”,而不是只给一组曲线。下面这段代码把扫描结果导出 CSV,字段和 3.2 节保持一致,方便后续画图或贴表。
import csv def export_sweep(rows, path="sweep_result.csv"): if not rows: raise ValueError("rows 为空,先跑 sweep_damping") fieldnames = list(rows[0].keys()) with open(path, "w", newline="", encoding="utf-8-sig") as f: writer = csv.DictWriter(f, fieldnames=fieldnames) writer.writeheader() writer.writerows(rows)导出前把浮点数保留 5~6 位有效数字即可,太多位数会让表格难读,也不代表精度。路径用相对路径,和脚本放同一目录,避免换机器后找不到文件。如果要批量跑多组 k2,建议文件名带参数标识,比如 sweep_k2_0p8.csv,别用 result1、result2 这种名字,过两天自己都分不清。
5.3 报告里该放哪些图表和指标
多阻尼对比的报告不需要堆图,四张图基本够用:
| 图表 | 内容 | 检查点 |
|---|---|---|
| 图1 | 车身加速度幅频对比 | 1~2 Hz 峰值是否随阻尼下降 |
| 图2 | 悬架动挠度幅频对比 | 低频峰值是否下降,高频是否抬升 |
| 图3 | 轮胎动载荷幅频对比 | 10~13 Hz 是否出现反向趋势 |
| 表1 | c、zeta、峰值、RMS 汇总 | 数值和图形态是否一致 |
表里至少要保留 c、zeta_susp、acc_peak、acc_peak_f、def_rms、load_rms 六列。若用了模态阻尼比,再加一列 zeta_modal_1 和 zeta_modal_2,并注明是通过特征值算的。图注写清楚“单位路面位移激励”“频率范围 0.1~20 Hz”“纵轴为增益”,不要把不同归一化方式的曲线放在同一张图里。复核时先看 0 Hz 渐近值,再看峰值趋势,最后看表格里的 RMS,顺序不要乱。
6. 把多阻尼对比做成模板:三个可复用技巧
6.1 参数表驱动与单位检查函数
我习惯把参数表单独放一个字典,所有脚本只读这个字典,不直接写数字。入口加一个 check_params,把单位、量级和符号拦住:
def check_params(p): assert p["m1"] > 0 and p["m2"] > 0, "质量必须为正" assert p["k1"] > p["k2"] > 0, "轮胎刚度通常大于悬架刚度" assert 0 < p["c"] < 3 * 2 * (p["k2"] * p["m2"]) ** 0.5, "阻尼系数超出粗估范围" assert 0.5 < p["m2"] / p["m1"] < 20, "质量比异常,先查单位"这段检查不追求覆盖所有工况,只拦住最常见的单位错和符号错。参数表里每个键后面写上单位注释,换人接手也能看懂。
6.2 图表固定口径与一键复跑
把频率轴、图幅、线型、图例位置写进一个 plot_style 字典,所有对比图共用。每次改阻尼列表只改 c_list,不改画图函数。这样出的图口径一致,报告里不会出现这张图 0.1~20 Hz、那张图 0.5~10 Hz 的混乱。我的习惯是跑完扫描先打印前三行和最后一行,确认 c、zeta、峰值没有明显跳变,再导出 CSV 和图片。
从那以后我每次换阻尼参数前,都先跑 0 Hz 渐近检查和状态空间残差,再画多阻尼曲线。这个习惯救过好几次“图看起来合理、表里峰值全错”的翻车。资源包拿到后,建议先跑通主脚本,再替换自己的参数表,别一上来就改传递函数。希望帮到你。
本文还有配套的精品资源,点击获取