汽车系统动力学Matlab.doc

汽车系统动力学Matlab

汽车系统动力学Matlab
作业报告

小组成员：

一、组内任务分配

二、Matlab程序及图形
1、不同转向特性车辆在不同车速下的系统特征根
m=1000;I=1500;a1=1.15;b1=1.35;Caf=53000;Car=53000;
i=1;R=[];
for uc=10:5:100;
D=(I*(Caf+Car)+m*(a1^2*Caf+b1^2*Car))/(m*I*uc);
S=(a1+b1)^2*Caf*Car/(m*I*uc^2)+(b1*Car-a1*Caf)/I;
P=[1 D S];
r=roots(P);
R(i,1)=r(1,1);R(i,2)=r(2,1);i=i+1;
end
plot(real(R(:,1)),imag(R(:,1)),'bo');
hold
a2=1.25;
b2=1.25;
t=1;
S=[];
for uc=10:5:100
P=[m 0;0 I];
Q=[(Caf+Car)/uc,m*uc+(a2*Caf-b2*Car)/uc;(a2*Caf-b2*Car)/uc,(a2^2*Caf+b2^2*Car)/uc];
R=[Caf;a2*Caf];
A=-P^(-1)*Q;
d=eig(A);
i=imag(d);
r=real(d);
S(t,1)=r(1);
S(t,2)=i(1);
t=t+1;
end
plot(S(:,1),S(:,2),'*')
a3=1.35;
b3=1.15;
for uc=10:5:100
P=[m 0;0 I];
Q=[(Caf+Car)/uc,m*uc+(a3*Caf-b3*Car)/uc;
(a3*Caf-b3*Car)/uc,(a3^2*Caf+b3^2*Car)/uc];
R=[Caf;a3*Caf];
A=-P^(-1)*Q;
d=eig(A);
i=imag(d);
r=real(d);
S(t,1)=r(1);
S(t,2)=i(1);
t=t+1;
end
grid on
plot(S(:,1),S(:,2),'d');
axis([-14 2 0 3]);
xlabel('实轴(Re)');
ylabel('虚轴(Im)');
text(-8,2.8,'不足转向');
text(0,0.2,'过多转向');
text(-3,0.2,'中性转向')
set(gca,'FontName','Helvetica','FontSize',10)
title(['不同转向特性车辆在不同车速下的系统特征根'],'FontSize',12);

2.1、具有不同转向特性车辆的横摆角速度幅频和相频响应
m=1000;I=1500;
a1=1.15;a2=1.25;a3=1.35;
b1=1.35;b2=1.25;b3=1.15;
Caf=53000;Car=53000;
D=[];C=[];M=[];
uc=20;
i=1;
for f=0.1:0.1:10
w=2*pi*f;
E=[1 0;0 1];
P=[m 0;0 I];
R1=[Caf;a1*Caf];
R2=[Caf;a2*Caf];
R3=[Caf;a3*Caf];
Q1=[(Caf+Car)/uc m*uc+(a1*Caf-b1*Car)/uc;(a1*Caf-b1*Car)/uc (a1^2*Caf+b1^2*Car)/uc];
Q2=[(Caf+Car)/uc m*uc+(a2*Caf-b2*Car)/uc;(a2*Caf-b2*Car)/uc (a2^2*Caf+b2^2*Car)/uc];
Q3=[(Caf+Car)/uc m*uc+(a3*Caf-b3*Car)/uc;(a3*Caf-b3*Car)/uc (a3^2*Caf+b3^2*Car)/uc];
A1=-inv(P)*Q1;
A2=-inv(P)*Q2;
A3=-inv(P)*Q3;
B1=inv(P)*R1;
B2=inv(P)*R2;
B3=inv(P)*R3;
Hw1=-inv(A1-1i*w*E)*B1;
Hw2=-inv(A2-1i*w*E)*B2;
Hw3=-inv(A3-1i*w*E)*B3;
D(i,1)=angle(Hw1(2))*180/pi;D(i,2)=abs(Hw1(2));
C(i,1)=angle(Hw2(2))*180/pi;C(i,2)=abs(Hw2(2));
M(i,1)=angle(Hw3(2))*180/pi;M(i,2)=abs(Hw3(2));
i=i+1;
end
subplot(2,1,1)
f=0.1:0.1:10;
semilogx(f,D(:,2),'k',f,C(:,2),'r-.',f,M(:,2))
grid on
set(gca,'Xtick',[0.1,0.3,1.5,3,10])
set(gca,'FontName','Helvetica','FontSize',10)
legend('不足转向','中性转向','过多转向')
title(['具有不同转向特性车辆的横摆角速度幅频响应'],'FontSize',12);
xlabel('频率/Hz')
ylabel('横摆角速度增益/{(°/s)/(°)}')
subplot(2,1,2)
f=0.1:0.1:10;
semilogx(f,D(:,1),'k',f,C(:,1),'r-.',f,M(:,1))m=1000;I=1500;
a1=1.15;a2=1.25;a3=1.35;
b1=1.35;b2=1.25;b3=1.15;
Caf=53000;Car=53000;
D=[];C=[];M=[];
uc=20;
i=1;
for f=0.1:0.1:10
w=2*pi*f;
E=[1 0;0 1];
P=[m 0;0 I];
R1=[Caf;a1*Caf];
R2=[Caf;a2*Caf];
R3=[Caf;a3*Caf];
Q1=[(Caf+Car)/uc m*uc+(a1*Caf-b1*Car)/uc;(a1*Caf-b1*Car)/uc (a1^2*Caf+b1^2*Car)/uc];
Q2=[(Caf+Car)/uc m*uc+(a2*Caf-b2*Car)/uc;(a2*Caf-b2*Car)/uc (a2^2*Caf+b2^2*Car)/uc];
Q3=[(Caf+Car)/uc m*uc+(a3*Caf-b3*Car)/uc;(a3*Caf-b3*Car)/uc (a3^2*Caf+b3^2*Car)/uc];
A1=-inv(P)*Q1;
A2=-inv(P)*Q2;
A3=-inv(P)*Q3;
B1=inv(P)*R1;
B2=inv(P)*R2;
B3=inv(P)*R3;
Hw1=-inv(A1-1i*w*E)*B1;
Hw2=-inv(A2-1i*w*E)*B2;
Hw3=-inv(A3-1i*w*E)*B3;
D(i,1)=angle(Hw1(2))*180/pi;D(i,2)=abs(Hw1(2));
C(i,1)=angle(Hw2(2))*180/pi;C(i,2)=abs(Hw2(2));
M(i,1)=angle(Hw3(2))*180/pi;M(i,2)=abs(Hw3(2));
i=i+1;
end
subplot(2,1,1)
f=0.1:0.1:10;
semilogx(f,D(:,2),'k',f,C(:,2),'r-.',f,M(:,2))
grid on
set(gca,'Xtick',[0.1,0.3,1.5,3,10])
set(gca,'FontName','Helvetica','FontSize',10)
legend('不足转向','中性转向','过多转向')
title(['具有不同转向特性车辆的横摆角速度幅频响应'],'FontSize',12);
xlabel('频率/Hz')
ylabel('横摆角速度增益/{(°/s)/(°)}')
subplot(2,1,2)
f=0.1:0.1:10;
semilogx(f,D(:,1),'k',f,C(:,1),'r-.',f,M(:,1))
grid on
set(gca,'Xtick',[0.1,0.3,1.5,3,10])
set(gca,'FontName','Helvetica','FontSize',10)
legend('不足转向','中性转向','过多转向')
title(['具有不同转向特性车辆的横摆角速度相频响应'],'FontSize',12);
xlabel('频率/Hz')
ylabel('横摆角速度相位（°）')

grid on
set(gca,'Xtick',[0.1,0.3,1.5,3,10])
set(gca,'FontName','Helvetica','FontSize',10)
legend('不足转向','中性转向','过多转向')
title(['具有不同转向特性车辆的横摆角速度相频响应'],'FontSize',12);
xlabel('频率/Hz')
ylabel('横摆角速度相位（°）')

2.2、具有不同转向特性车辆的侧向加速度和相频响应
m=1000;I=1500;
a1=1.15;a2=1.25;a3=1.35;
b1=1.35;b2=1.25;b3=1.15;
Caf=53000;Car=53000;
D=[];C=[];M=[];
uc=20;L=a1+b1;
i=1;
for f=0.1:0.1:10
w=2*pi*f;
E=[1 0;0 1];
Vi=w*I*Caf;
Rr=L*Caf*Car/uc;
Vr1=(L*b1*Caf*Car/uc-m*a1*Caf*uc);
Vr2=(L*b2*Caf*Car/uc-m*a2*Caf*uc);
Vr3=(L*b3*Caf*Car/uc-m*a3*Caf*uc);
Ri1=w*m*a1*Caf;
Ri2=w*m*a2*Caf;
Ri3=w*m*a3*Caf;
Dr1=-w^2*m*I+L^2*Caf*Car/(uc^2)+m*(b1*Car-a1*Caf);
Dr2=-w^2*m*I+L^2*Caf*Car/(uc^2)+m*(b2*Car-a2*Caf);
Dr3=-w^2*m*I+L^2*Caf*Car/(uc^2)+m*(b3*Car-a3*Caf);
Di1=w*(I*(Caf+Car)+m*(a1^2*Caf+b1^2*Car))/uc;
Di2=w*(I*(Caf+Car)+m*(a2^2*Caf+b2^2*Car))/uc;
Di3=w*(I*(Caf+Car)+m*(a3^2*Caf+b3^2*Car))/uc;
Hay1=1i*w*(Vr1+1i*Vi)/(Dr1+1i*Di1)+uc*(Rr+1i*Ri1)/(Dr1+1i*Di1);
Hay2=1i*w*(Vr2+1i*Vi)/(Dr2+1i*Di2)+uc*(Rr+1i*Ri2)/(Dr2+1i*Di2);
Hay3=1i*w*(Vr3+1i*Vi)/(Dr3+1i*Di3)+uc*(Rr+1i*Ri3)/(Dr3+1i*Di3);
D(i,1)=angle(Hay1)*180/pi;D(i,2)=abs(Hay1);
C(i,1)=angle(Hay2)*180/pi;C(i,2)=abs(Hay2);
M(i,1)=angle(Hay3)*180/pi;M(i,2)=abs(Hay3);
i=i+1;
end
subplot(2,1,1)
f=0.1:0.1:10;
semilogx(f,D(:,2)*(pi/(180*9.8)),'k',f,C(:,2)*(pi/(180*9.8)),'r-.',f,M(:,2)*(pi/(180*9.8)))
grid on
set(gca,'XTick',[0.1 0.3 1 3 10])
set(gca,'FontName','Helvetica','FontSize',10)
legend('不足转向','中性转向','过多转向')
title(['具有不同转向特性车辆的侧向加速度幅频响应'],'FontSize',12);
xlabel('频率/Hz')
ylabel('侧向加速度增益/[g/(°)]')
subplot(2,1,2)
f=0.1:0.1:10;
semilogx(f,D(:,1),'k',f,C(:,1),'r-.',f,M(:,1))
grid on
legend('不足转向','中性转向','过多转向')
set(gca,'XTick',[0.1 0.3 1 3 10])
set(gca,'FontName','Helvetica','FontSize',10)
title(['具有不同转向特性车辆的侧向加速度幅频响应'],'FontSize',12);
axis([0.1 10 -120 120])
xlabel('频率/Hz')
ylabel('侧向加速度相位（°）')

2.3、 ABS控制器设计
w0=120;v0=30;Tb=600;ki=4500;kd=5000;Ts=0.05;Iw=12;
rd=0.25;uh=0.8;ug=0.6;s0=0.2;m=300;g=9.8;
w=[];v=[];s=[];
w(1)=120;v(1)=30;s(1)=0;
k=0;i=2;
while(v0>0)
k=k+1;
sb=(v0-rd*w0)/v0;
s(i)=sb;
if sb<=s0
u=uh/s0*sb;
else u=(uh-ug*s0)/(1-s0)-rd*sb;
end
Fxb=u*m*g;
if sb>0.22
Tb=Tb-kd*Ts;
elseif sb<0.18
Tb=Tb+ki*Ts;
else Tb=Tb;
end
w0=w0+(Fxb*rd-Tb)/Iw*Ts;
w(i)=w0;
v0=v0-Fxb/m*Ts;
v(i)=v0;
i=i+1;
end
t=[0:Ts:Ts*k];
subplot(2,1,1)
plot(t,v,t,w*rd);
grid on
legend('车轮前进速度','车轮线速度');
set(gca,'FontName','Helvetica','FontSize',10)
title(['车轮前进速度及车轮线速度关系曲线'],'FontSize',12);
xlabel('时间/（s）');
ylabel('速度/（m/s）');
axis([0,Ts*k,0,32]);
subplot(2,1,2)
plot(t,s);
axis([0,Ts*k,0,1]);
grid on
set(gca,'FontName','Helvetica','FontSize',10)
title(['ABS控制的滑移率时域仿真结果'],'FontSize',12);
xlabel('时间/（s）');
ylabel('滑移率');
set(gca,'Ytick',0:0.2:1)

3、扭振系统振型图
J1=[1.986*10^(-3);1.910*10^(-3);1.931*10^(-3);1.931*10^(-3);1.910*10^(-3);1.924*10^(-3);7.8426*10^(-2);2.258*10^(-3);2.641*10^(-2);2.91*10^(-3);
2.51*10^(-3);1.77*10^(-3);7.836*10^(-2);3.238];
K1=[7.95*10^4;7.95*10^4;6.95*10^4;7.95*10^4;7.95*10^4;6.90*10^4; 8.93*10^3;1.41*10^4;1.02*10^4;4.4*10^3;1.38*10^4;1.616*10^2;
2.51*10^2];
J=diag(J1);
K(1,1)=K1(1);K(1,2)=-K1(1);
i=2;
for i=2:13;
K(i,i-1)=-K1(i-1);
K(i,i)=K1(i-1)+K1(i);
K(i,i+1)=-K1(i);
end
K(14,13)=-K1(13);K(14,14)=K1(13);
A=inv(J)*K;
[G,D]=eig(A);
f=sqrt(D)/(2*pi);
m=1;
while m<=14
G(:,m)=G(:,m)/G(1,m)
m=m+1
end
subplot(3,2,1)
plot((G(:,13)))
xlabel('质点号');
set(gca,'FontName','Helvetica','FontSize',10)
grid on
title(['单节点振型图（4.26Hz）'],'FontSize',12)
axis([1 15 -0.2 1.1])
subplot(3,2,2)
plot(G(:,12))
xlabel('质点号');
set(gca,'FontName','Helvetica','FontSize',10)
grid on
axis([1 15 -5 1.2])
title(['双节点振型图（11.99Hz）'],'FontSize',12)
subplot(3,2,3)
plot(G(:,11))
xlabel('质点号');
set(gca,'FontName','Helvetica','FontSize',10)
grid on
axis([1 15 -7 1.2])
title(['三节点振型图 (73.9Hz)'],'FontSize',12)
subplot(3,2,4)
plot(G(:,10))
xlabel('质点号');
set(gca,'FontName','Helvetica','FontSize',10)
grid on
title(['四节点振型图 (147.7Hz)'],'FontSize',12)
axis([1 15 -9 37])
subplot(3,2,5)
plot(G(:,9))
xlabel('质点号');
set(gca,'FontName','Helvetica','FontSize',10)
grid on
axis([1 15 -0.2 1.1])
title(['五节点振型图 (252Hz)'],'FontSize',12)

5、福特Granada轿车后悬架单轮模型频率响应函数和系统相应输入功率谱密度
mb=317.5;mw=45.4;ks=22000;kt=192000;cs=1500;Go=5*10^(-6);Uc=20;
n=1;B=[];D=[];
for f=0:0.01:15
w=2*pi*f;
Sf=4.47*10^(-4)*power(f,-2.5);
A=[i*cs*w+(ks+kt-w^2*mw),-i*cs*w-ks;-i*cs*w-ks,i*cs*w+(ks-mb*w^2)];
C=[kt;0];
D=-A^(-1)*C;
B(1,n)=abs(D(2)-D(1));
B(2,n)=abs(-D(2)*w^2);
B(3,n)=abs(kt*(D(1)-1));
B(4,n)=Sf;
B(5,n)=B(2,n)^2*Sf;
B(6,n)=B(1,n)^2*Sf;
B(7,n)=B(3,n)^2*Sf;
n=n+1;
end
f=0:0.01:15
subplot(3,2,1)
plot(f,B(1,:));
xlabel('频率/Hz')
ylabel('悬架动行程增益')
grid on
subplot(3,2,2);
plot(f,B(2,:));
grid on
xlabel('频率/Hz')
ylabel('不舒适性参数增益/[(m/s^2)/m]')
subplot(3,2,3);
plot(f,B(3,:));
grid on
xlabel('频率/Hz')
ylabel('轮胎动载荷增益/[N/m]')
subplot(3,2,4)
plot(f,B(6,:));
xlabel('频率/Hz')
ylabel('悬架动行程功率谱密度/[m^2/Hz]')
grid on
subplot(3,2,5);
plot(f,B(5,:));
grid on
xlabel('频率/Hz')
ylabel('不舒适性功率谱密度/[(m/s^2)^2/Hz]')
subplot(3,2,6);
plot(f,B(7,:));
grid on
xlabel('频率/Hz')
ylabel('轮胎动载荷功率谱密度/[N^2/Hz]')

6、魔术公式
xdata = [0 0.78 1.88 2.79 3.80 4.82 6.29 7.82 9.31 11.80];
ydata = [0 7.99 16.05 20.65 23.82 26.00 27.78 28.87 29.65 30.09];
zdata = [0 1.65 1.61 0.77 -0.21 -1.03 -1.90 -2.59 -2.61 -2.68];
a0 = [11 1 0.8 1];
b0 = [1 1 0.8 1];
a = lsqcurvefit(@fun1, a0, xdata, ydata);
b = lsqcurvefit(@fun2, b0, xdata, zdata);
yy = fun1(a,xdata);
zz = fun2(b,xdata);
subplot(2,1,1)
hold on
scatter(xdata, ydata);
grid on
plot(xdata, yy);
set(gca,'FontName','Helvetica','FontSize',10)
title(['MF公式拟合轮胎侧偏力-侧偏角曲线'],'FontSize',12);
xlabel('侧偏角 α/（°）')
ylabel('侧偏力 Fy/(N)')
subplot(2,1,2)
hold on;
scatter(xdata,zdata);
grid on
plot(xdata, zz);
set(gca,'FontName','Helvetica','FontSize',10)
title(['MF公式拟合轮胎回正力矩-侧偏角曲线'],'FontSize',12);
xlabel('侧偏角 α/（°）')
ylabel('回正力矩 Mz/(N·m)')

function y = fun1(a,xdata);
y=a(1)*sin(a(2)*atan(a(3)*xdata-a(4)*(a(3)*xdata-atan(a(3)*xdata))))

function z = fun2(b,xdata);
z=b(1)*sin(b(2)*atan(b(3)*xdata-b(4)*(b(3)*xdata-atan(b(3)*xdata))))

三、结论
通过实际Matlab编程，不仅巩固了课堂知识，并且学以致用，更加深入地了解到了应该如何运用，解决问题。通过编程，也掌握了Matlab的一些基本语言，为以后研究阶段的软件应用打下基础。

33 / 33