汽车系统动力学Matlab之欧阳化创编.docx
《汽车系统动力学Matlab之欧阳化创编.docx》由会员分享，可在线阅读，更多相关《汽车系统动力学Matlab之欧阳化创编.docx（17页珍藏版）》请在冰点文库上搜索。
汽车系统动力学Matlab之欧阳化创编
汽车系统动力学Matlab
时间：
2021.02.12
创作人：
欧阳化
作业报告
小组成员：
一、组内任务分配
二、Matlab程序与图形
1、不同转向特性车辆在不同车速下的系统特征根
m=1000;I=1500;a1=1.15;b1=1.35;Caf=53000;Car=53000;
i=1;R=[];
foruc=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=[1DS];
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=[];
foruc=10:
5:
100
P=[m0;0I];
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;
foruc=10:
5:
100
P=[m0;0I];
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
gridon
plot（S（:
1）,S（:
2）,'d'）;
axis（[-14203]）;
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;
forf=0.1:
0.1:
10
w=2*pi*f;
E=[10;01];
P=[m0;0I];
R1=[Caf;a1*Caf];
R2=[Caf;a2*Caf];
R3=[Caf;a3*Caf];
Q1=[（Caf+Car）/ucm*uc+（a1*Caf-b1*Car）/uc;（a1*Caf-b1*Car）/uc（a1^2*Caf+b1^2*Car）/uc];
Q2=[（Caf+Car）/ucm*uc+（a2*Caf-b2*Car）/uc;（a2*Caf-b2*Car）/uc（a2^2*Caf+b2^2*Car）/uc];
Q3=[（Caf+Car）/ucm*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））
gridon
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;
forf=0.1:
0.1:
10
w=2*pi*f;
E=[10;01];
P=[m0;0I];
R1=[Caf;a1*Caf];
R2=[Caf;a2*Caf];
R3=[Caf;a3*Caf];
Q1=[（Caf+Car）/ucm*uc+（a1*Caf-b1*Car）/uc;（a1*Caf-b1*Car）/uc（a1^2*Caf+b1^2*Car）/uc];
Q2=[（Caf+Car）/ucm*uc+（a2*Caf-b2*Car）/uc;（a2*Caf-b2*Car）/uc（a2^2*Caf+b2^2*Car）/uc];
Q3=[（Caf+Car）/ucm*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））
gridon
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））
gridon
set（gca,'Xtick',[0.1,0.3,1.5,3,10]）
set（gca,'FontName','Helvetica','FontSize',10）
legend（'不足转向','中性转向','过多转向'）
title（['具有不同转向特性车辆的横摆角速度相频响应'],'FontSize',12）;
xlabel（'频率/Hz'）
ylabel（'横摆角速度相位（°）'）
gridon
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;
forf=0.1:
0.1:
10
w=2*pi*f;
E=[10;01];
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）））
gridon
set（gca,'XTick',[0.10.31310]）
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））
gridon
legend（'不足转向','中性转向','过多转向'）
set（gca,'XTick',[0.10.31310]）
set（gca,'FontName','Helvetica','FontSize',10）
title（['具有不同转向特性车辆的侧向加速度幅频响应'],'FontSize',12）;
axis（[0.110-120120]）
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;
ifsb<=s0
u=uh/s0*sb;
elseu=（uh-ug*s0）/（1-s0）-rd*sb;
end
Fxb=u*m*g;
ifsb>0.22
Tb=Tb-kd*Ts;
elseifsb<0.18
Tb=Tb+ki*Ts;
elseTb=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）;
gridon
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]）;
gridon
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;
fori=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;
whilem<=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）
gridon
title（['单节点振型图（4.26Hz）'],'FontSize',12）
axis（[115-0.21.1]）
subplot（3,2,2）
plot（G（:
12））
xlabel（'质点号'）;
set（gca,'FontName','Helvetica','FontSize',10）
gridon
axis（[115-51.2]）
title（['双节点振型图（11.99Hz）'],'FontSize',12）
subplot（3,2,3）
plot（G（:
11））
xlabel（'质点号'）;
set（gca,'FontName','Helvetica','FontSize',10）
gridon
axis（[115-71.2]）
title（['三节点振型图（73.9Hz）'],'FontSize',12）
subplot（3,2,4）
plot（G（:
10））
xlabel（'质点号'）;
set（gca,'FontName','Helvetica','FontSize',10）
gridon
title（['四节点振型图（147.7Hz）'],'FontSize',12）
axis（[115-937]）
subplot（3,2,5）
plot（G（:
9））
xlabel（'质点号'）;
set（gca,'FontName','Helvetica','FontSize',10）
gridon
axis（[115-0.21.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=[];
forf=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（'悬架动行程增益'）
gridon
subplot（3,2,2）;
plot（f,B（2,:
））;
gridon
xlabel（'频率/Hz'）
ylabel（'不舒适性参数增益/[（m/s^2）/m]'）
subplot（3,2,3）;
plot（f,B（3,:
））;
gridon
xlabel（'频率/Hz'）
ylabel（'轮胎动载荷增益/[N/m]'）
subplot（3,2,4）
plot（f,B（6,:
））;
xlabel（'频率/Hz'）
ylabel（'悬架动行程功率谱密度/[m^2/Hz]'）
gridon
subplot（3,2,5）;
plot（f,B（5,:
））;
gridon
xlabel（'频率/Hz'）
ylabel（'不舒适性功率谱密度/[（m/s^2）^2/Hz]'）
subplot（3,2,6）;
plot（f,B（7,:
））;
gridon
xlabel（'频率/Hz'）
ylabel（'轮胎动载荷功率谱密度/[N^2/Hz]'）
6、魔术公式
xdata=[00.781.882.793.804.826.297.829.3111.80];
ydata=[07.9916.0520.6523.8226.0027.7828.8729.6530.09];
zdata=[01.651.610.77-0.21-1.03-1.90-2.59-2.61-2.68];
a0=[1110.81];
b0=[110.81];
a=lsqcurvefit（@fun1,a0,xdata,ydata）;
b=lsqcurvefit（@fun2,b0,xdata,zdata）;
yy=fun1（a,xdata）;
zz=fun2（b,xdata）;
subplot（2,1,1）
holdon
scatter（xdata,ydata）;
gridon
plot（xdata,yy）;
set（gca,'FontName','Helvetica','FontSize',10）
title（['MF公式拟合轮胎侧偏力-侧偏角曲线'],'FontSize',12）;
xlabel（'侧偏角α/（°）'）
ylabel（'侧偏力Fy/（N）'）
subplot（2,1,2）
holdon;
scatter（xdata,zdata）;
gridon
plot（xdata,zz）;
set（gca,'FontName','Helvetica','FontSize',10）
title（['MF公式拟合轮胎回正力矩-侧偏角曲线'],'FontSize',12）;
xlabel（'侧偏角α/（°）'）
ylabel（'回正力矩Mz/（N·m）'）
functiony=fun1（a,xdata）;
y=a
（1）*sin（a
（2）*ata
升级会员 