1、1 行驶动力学计算机建模，仿真及主动悬架控制器设计1.1 建立 M 文件%车辆模型参数mb=320; %簧载质量(kg)mw=40; %非簧载质量(kg)ks=20000; %悬架刚度（N/m）kt=200000; %轮胎刚度（N/m）swsc=+-100; %悬架工作空间（mm ）%仿真路面输入参数g0=5*10-6; %路面不平度系数（m3/cycle）u=20; %车速（m/s）f0=0.1; %下截止频率（Hz ）%性能指标加权系数q1=80000; %轮胎动位移q2=5; %悬架动行程q3=1; %车身加速度%系统输入x0=wgn(10001,1,20); %生成高斯白噪声（采样点数
2、 列数 功率）t=0:0.005:50; %仿真时间Tw=t x0; %系统输入序列仿真时间序列 白噪声序列%系统状态方程矩阵A=0 0 -ks/mb ks/mb 0;0 0 ks/mw -(ks+kt)/mw kt/mw;1 0 0 0 0;0 1 0 0 0;0 0 0 0 -2*pi*f0;B=1/mb;-1/mw;0;0;0;F=0;0;0;0;2*pi*(g0*u)0.5；%二次最优控制器设计Q=0 0 0 0 0;0 0 0 0 0;0 0 q2+ks2/mb2 -q2-ks2/mb2 0;0 0 -q2-ks2/mb2 q1+q2+ks2/mb2 -q1;0 0 0 -q1 q1
3、;R=1/mb2;N=0;0;-ks/mb2;ks/mb2;0;K,S,E=LQR(A,B,Q,R,N); %K 为最优控制反馈增益矩阵以上程序执行后可得到的结果为：A=0 0 -62.5 62.5 0;0 0 500 -5500 5000;1 0 0 0 0;0 1 0 0 0;0 0 0 0 -0.6;B=0.0031;-0.025;0;0;0;F=0;0;0;0;0.0628;K =711.88 -1241.5 -19284 -2038.5 20864;Tw时间序列 白噪声幅值序列1.2 建立 simulink 模型设置输入模块：（Tw 序列）设置状态空间方程参数：（其中 B 矩阵为 M
4、 执行结果中 B 与 F 矩阵的和矩阵维数：42）设置增益模块：（K 矩阵）1.3 执行结果根据以上设置执行仿真，可得如下结果：车身加速度（m/s2）-时间（s）悬架动行程(m)-时间（s ）轮胎动位移(m)-时间（s ）路面位移输入(m)-时间（s）1.4 结果分析根据上图可分别导出车身加速度 BA，悬架动行程 SWS，轮胎动位移 DTD 的数据。三者都是 61362 的矩阵，第一列为幅值对应的时间，第二列为三者的幅值，对其第二列分别求均方根值，程序为：Y=sqrt(sum(X(:,2).2)/6136)执行程序，可得主动悬架性能指标的均方根值如表 1。在相同仿真条件下，可对一个被动系统进行
5、分析，取悬架刚度 Ks=22000N/m，阻尼系数 Cs=1000N*s/m。其他参数均与主动系统相同，计算可得被动系统性能指标如表 1。表 1 主动悬架与被动悬架性能指标均方根值的比较均方根值性能指标 单位主动悬架 被动悬架车身加速度 m/s2 1.4288 1.7533悬架动行程 mm 44.3 17.6轮胎动位移 mm 5.9 6.0由表可见，在轮胎动位移基本相同的情况下，所设计的最优主动悬架显著地降低了车身的垂向振动加速度，与被动悬架系统相比，其均方根值减少了 18.5%。而且，主动悬架系统的悬架动行程被很好的控制在设计的范围（100mm）内，且比被动悬架的行程大了许多，意味着许用的悬
6、架工作空间得到了更充分的利用。2 操纵动力学计算机建模及分析2.1 操纵动力学模型建立%模型参数m=2045 1008; %质量Izz=5428 1031; %横摆转动惯量a=1.488 1.234; %前轴到质心距离b=1.712 1.022; %后轴到质心距离Caf=77850 117440; %前轮侧偏刚度Car=76510 144930; %后轮侧偏刚度uc=40; %速度%建立状态方程矩阵%矩阵 AA11=-(Caf+Car)./(m*uc);A12=-(a.*Caf-b.*Car)./(m*uc)-uc;A21=-(a.*Caf-b.*Car)./(Izz*uc);A22=-(a.
7、2.*Caf+b.2.*Car)./(Izz*uc);A1=A11(1) A12(1);A21(1) A22(1);A2=A11(2) A12(2);A21(2) A22(2);%矩阵 BB11=Caf./m;B22=a.*Caf./Izz;B1=B11(1);B22(1);B2=B11(2);B22(2);%矩阵 CDC=0 1;D=0;执行结果为：A1 =-1.8870 -39.8149;0.06981 -1.8267A2 =-6.5072 -39.9207;0.0775 -8.0070B1=38.0685;21.3413B2=116.5079;140.5635C=0 1;D=0;2.2
8、时域分析建立 Simulink 模型阶跃输入参数设置：初始时刻（t=0s）输入为 0rad，t0s 时输入为 0.0058rad。别克模型状态方程参数设置：(矩阵 A1 B1 C D)法拉利模型状态方程参数设置：（矩阵 A2 B2 C D）运行结果：时间/s角阶跃输入下的横摆角速度时域响应( 紫色为法拉利，黄色为别克)由图可知，在同样的转向盘转角输入下，法拉利跑车的瞬态响应比别克轿车的要好，主要体现在较短的响应时间，较小的超调量以及更好的阻尼特性等。横摆角速度/rad/s2.3 频域分析执行程序如下：figure()sys1=ss(A1,B1,C,D);sys2=ss(A2,B2,C,D);b
9、ode(sys1,sys2);结果为：由图可知，对于幅频响应在同一行驶车速下，法拉利跑车的响应带宽大于别克轿车的响应带宽，从而也说明了前者具有更好的频率响应特性，对于相频特性，法拉利跑车的系统响应之后要比别克车的少，系统延迟小。2.3 稳定性分析求解车辆在车速 uc=1560m/s 的特征值执行程序如下：figure()u=15:5:60;for j=1:length(u);A=A11(1)*uc/u(j) (A12(1)+uc)*uc/u(j)-u(j);A21(1)*uc/u(j) A22(1)*uc/u(j);plot(real(eig(A),imag(eig(A),go);hold o
10、n;endfor j=1:length(u);A=A11(2)*uc/u(j) (A12(2)+uc)*uc/u(j)-u(j);A21(2)*uc/u(j) A22(2)*uc/u(j);plot(real(eig(A),imag(eig(A),bx);hold on;end结果为：由图可知，随着车速的增加，系统的特征值由复平面的左侧不断向虚轴靠近，因此系统越来越趋向于不稳定，对两种车型的比较可见，他们的特征跟随车速变化而变化的趋势也不同，由于法拉利跑车具有相对较大的稳定裕度，因而其特征根位置与别克轿车相比更远离虚轴，而别克轿车的转向特性的稳定裕量较小且容易失稳。3 半车模型的推导分析3.1
11、 模型建立%参数设置mhb=690; %二分之一车身质量（kg）Ihp=1222; %二分之一转动惯量(kg*m 2)mwf=40.55; %前轮非簧载质量(kg)mwr=45.4; %后轮非簧载质量(kg)Ktf=192000; %前轮胎刚度(kN/m)Ktr=192000; %后轮胎刚度(kN/m)Ksf=17000; %前悬架刚度(kN/m)Ksr=22000; %后悬架刚度(kN/m)Csf=1500; %前悬架阻尼系数(kN*s/m)Csr=1500; %后悬架阻尼系数(kN*s/m)a=1.25; %车身质心至前轴距离（ m）b=1.51; %车身质心至后轴距离（ m）L=2.76
12、; %轴距（m）i=(-1)0.5; %虚数u=20; %车速(m/s)t=L/u; %前后轮输入时间间隔（s）%建立动力学模型w=sym(2*pi*f)c1=1/mhb+a2/Ihp;c2=1/mhb-a*b/Ihp;c3=1/mhb+b2/Ihp; A1=Csf*i*w+Ksf+Ktf-mwf*w2;c1*(Ksf+Csf*i*w);0;c2*(Ksf+Csf*i*w);A2=-Csf*i*w-Ksf;w2-c1*(Ksf+Csf*i*w);0;-c2*(Ksf+Csf*i*w);A3=0;c2*(Ksr+Csr*i*w);Csr*i*w+Ksr+Ktr-mwr*w2;c3*(Ksr+Cs
13、r*i*w);A4=0;-c2*(Ksr+Csr*i*w);-Csr*i*w-Ksr;w2-c3*(Ksr+Csr*i*w);B=Ktf;0;Ktr/exp(i*w*t);0; G1=det(B A2 A3 A4)/det(A1 A2 A3 A4);G2=det(A1 B A3 A4)/det(A1 A2 A3 A4);G3=det(A1 A2 B A4)/det(A1 A2 A3 A4);G4=det(A1 A2 A3 B)/det(A1 A2 A3 A4);3.2 仿真程序及结果%求悬架动行程增益的频响函数Gswsf=G1-G2;Gswsr=G3-G4;fsf=inline(Gswsf);
14、fsr=inline(Gswsr);f=0:0.01:15;figure();plot(f,abs(fsf(f),b,f,abs(fsr(f),g);结果：%求垂向加速度增益的频响函数Gvaf=-w2*G2;Gvar=-w2*G4;fvf=inline(Gvaf);fvr=inline(Gvar);f=0:0.01:15;figure();plot(f,abs(fvf(f),b,f,abs(fvr(f),g);结果：%求轮胎动载荷增益的频响函数Gdtlf=Ktf*(G1-1);Gdtlr=Ktr*(G3-1);fdf=inline(Gdtlf);fdr=inline(Gdtlr);f=0:0.
15、01:15;figure();plot(f,abs(fdf(f),b,f,abs(fdr(f),g);结果：%求质心处垂向加速度Gvac=( Gvar*a+ Gvaf*b)/(a+b)fvc= inline(Gvac);f=0:0.01:15;figure();plot(f,abs(fvc(f);结果： %求俯仰角增益Gp=(Gvar-Gvaf)/(a+b)/-w2;frl= inline(Gp);f=0:0.01:15;figure();plot(f,abs(frl(f);结果：结果分析：与单轮车辆模型相比，半车模型不同的地方在于，后轮的输入将会比前轮延迟时间t=L/u。故前后两轮输入存在相
16、位差，前后轮输入的共同作用决定了质心垂向加速度和俯仰角的频率响应特性，且影响规律并不是简单的线性叠加。4 整车模型的推导及分析4.1 模型建立%设置参数mb=1380; %整车质量（kg）Ip=2440; %俯仰转动惯量(kg*m 2)Ir=380; %侧倾转动惯量(kg*m 2)mwa=40.5; %A 处的簧下质量质量（kg）mwb=40.5; %B 处的簧下质量质量（ kg）mwc=45.4; %C 处的簧下质量质量（ kg）mwd=45.4; %D 处的簧下质量质量（kg）kta=192000; %A 处的轮胎刚度（N/m）ktb=192000; %B 处的轮胎刚度（N/m）ktc=1
17、92000; %C 处的轮胎刚度（N/m）ktd=192000; %D 处的轮胎刚度（N/m）ksa=17000; %A 处的悬架刚度（N/m）ksb=17000; %B 处的悬架刚度（N/m）ksc=22000; %C 处的悬架刚度（N/m）ksd=22000; %D 处的悬架刚度（N/m）Csa=1500; %A 处的阻尼系数（N*s/m）Csb=1500; %B 处的阻尼系数（N*s/m）Csc=1500; %C 处的阻尼系数（N*s/m）Csd=1500; %D 处的阻尼系数（N*s/m）a=1.25; %车身质心至前轴距离（ m）b=1.51; %车身质心至后轴距离（ m）L=2.7
18、6; %轴距（m）i=(-1)0.5; %虚数u=20; %车速(m/s)t=L/u; %前后轮输入时间间隔（s）Bf=1.48; %前轴轮距Br=1.48; %后轴轮距%建立动力学模型ea=Csa*i*w+ksa;eb=Csb*i*w+ksb;ec=Csc*i*w+ksc;ed=Csd*i*w+ksd;A1=ea+kta-mwa*w2;0;0;0;ea;-a*ea;Bf/2*ea;A2=0;eb+ktb-mwb*w2;0;0;eb;-a*eb;-Bf/2*eb;A3=0;0;ec+ktc-mwc*w2;0;ec;b*ec;Br/2*ec;A4=0;0;0;ed+ktd-mwd*w2;ed;b
19、*ed;-Br/2*ed;A5=-ea;-eb;-ec;-ed;mb*w2-(ea+eb+ec+ed);a*(ea+eb)-b*(ec+ed);Bf/2*(-ea+eb)+Br/2*(-ec+ed);A6=a*ea;a*eb;-b*ec;-b*ed;a*(ea+eb)-b*(ec+ed);Ip*w2-a2*(ea+eb)+b2*(ec+ed);a*Bf/2*(ea-eb)-b*Br/2*(ec-ed);A7=-Bf/2*ea;Bf/2*eb;-Br/2*ec;Br/2*ed;Bf/2*(-ea+eb)+Br/2*(-ec+ed);a*Bf/2*(ea-eb)+b*Br/2*(-ec+ed);-
20、Bf2/4*(ea+eb)-Br2/4*(ec+ed); B=Kta;Ktb/exp(i*w*2*t);Ktc/exp(i*w*t);Ktd/exp(i*w*3*t);0;0;0;G1=det(B A2 A3 A4 A5 A6 A7)/det(A1 A2 A3 A4 A5 A6 A7);G2=det(A1 B A3 A4 A5 A6 A7)/det(A1 A2 A3 A4 A5 A6 A7);G3=det(A1 A2 B A4 A5 A6 A7)/det(A1 A2 A3 A4 A5 A6 A7);G4=det(A1 A2 A3 B A5 A6 A7)/det(A1 A2 A3 A4 A5 A
21、6 A7);G5=det(A1 A2 A3 A4 B A6 A7)/det(A1 A2 A3 A4 A5 A6 A7);G6=det(A1 A2 A3 A4 A5 B A7)/det(A1 A2 A3 A4 A5 A6 A7);G7=det(A1 A2 A3 A4 A5 A6 B)/det(A1 A2 A3 A4 A5 A6 A7);4.2 仿真程序及结果%侧倾角增益Groll=G7;frl= inline(Groll);f=0:0.01:15;figure();plot(f,abs(frl(f);结果：%俯仰角增益Gp=G6;fp= inline(Gp);f=0:0.01:15;figure();plot(f,abs(fp(f);结果：%质心加速度增益Gvc=-w2*G5;fvc= inline(Gvc);f=0:0.01:15;figure();plot(f,abs(fvc(f);结果：由图可知，对于整车模型而言，当四轮输入均不相同时，车身某点的响应可以认为是四者同时作用的主模态各自频率响应的叠加。故其频响分析较为复杂，但当前轮输入相同，后轮输入相同，且前后轮间输入存在时间延迟时，即转化为半车模型，通过设置相应的输入，可以通过频响函数进行验证。
