基于matlab的大学生方程式汽车动力性仿真【最高车速、百公里加速时间、75m加速时间、最优传动比】

分为两个程序，第一个给出发动机参数拟合曲线、驱动力-阻力曲线、最高车速、加速度倒数曲线、百公里加速时间、75m加速时间。第二个程序给出最优主减速比。主减速比一般取2.9-3.8，取值太离谱会报错。大多数报错都是因为拟合之后的曲线产生了错误的交点，可以通过改变拟合阶次来尝试解决。

% 程序一drivetrain for 21C
% 创建版本：matlab2018b
% 车辆参数
nmin=4000;      % 直线加速换挡过程中发动机转速最小值
nmax=8500;      % 直线加速换挡过程中发动机转速最大值
n=nmin:0.1:nmax;        % 发动机转速区间
ig1=5.475;      % 一档传动比
ig2=3.8325;     % 二档传动比
ig3=2.9199;     % 三档传动比
ig4=2.3985;     % 四档传动比
ig5=2.0947;     % 五档传动比
ig6=1.9044;     % 六档传动比
i0=3.5;     % 主减速比
nt=0.9;     % 传动效率
r=0.2286;       % 车轮半径
f=0.013;        % 滚阻系数
m=250;      % 质量（带车手）
g=9.8;      % 重力加速度
Cd=1.5374998;       % 空气阻力系数
A=1.040632;     % 迎风面积
CDA=Cd*A;       % 乘一起命名新变量
L=1.58;     % 轴距
a_1=0.8374;     % 前轴至质心的距离
hg=0.3;     % 质心高度
Iw1=0.42964;        % 两个前轮转动惯量和
Iw2=0.42964;        % 通过catia量轮胎、轮辋、轮辐的转动惯量加一起，是这个数值(两个后轮的)。后轮因为还要加上差速器半轴等会略大一些，不过影响甚小可以忽略不计。
If=0.009;       % 飞轮转动惯量
fd=1.2;     % 弹射起步较大滑移率时轮胎的纵向附着系数。建议取值之后在图4中和绿色线对比，二者应相近，此时fd的取值和实际值可能有出入，不过绿色线是HRT20C在练习时的实际数据
% （路面条件一般），所以建议以这个值为准。
ah20c=0.9*g;        % HRT20C的MOTEC测量值。测量取样4次，20年赛场路面保守加速度。
trantime=0.1;       % 换挡时间
n0=3500:500:9000;       % 标定数据给出的转速点
Tq0=[46.4 46.2 48.9 48.9 52.6 52.5 54 53.7 50.6 46 40.2 33.8];      % 标定数据给出的扭矩值
Pe0=[17.2 19.3 23.2 25.6 30.1 33 36.7 39.4 39.4 38.5 35.2 31.5];       % 标定数据给出的功率值% 弹射起步
uhmax=0.377*nmax*r/ig2/i0;      % 定义弹射起步曲线绘制区间上限
uh=0:0.01:uhmax;        % 定义弹射起步曲线绘制自变量uh
FN2=fd*m*g*a_1/(L-fd*hg);      % 弹射起步驱动力
Fwh=CDA/21.15*uh.^2;        % 空气阻力
ah=(FN2-m*g*f-Fwh)/m;       % 弹射起步加速度随速度的变化函数
adh=1./ah;      % 加速度取倒数% 拟合出扭矩和功率的曲线表达式，最高阶设置为3阶
tqq=polyfit(n0,Tq0,3);      % 拟合转矩随转速变化曲线
Tq=polyval(tqq,n);      % 根据polyfit函数得到的多项式系数向量求得不同n取值时的转矩值
pee=polyfit(n0,Pe0,3);      % 拟合功率随转速变化曲线
Pe=polyval(pee,n);      % 根据polyfit函数得到的多项式系数向量求得不同n取值时的功率值
subplot(2,3,1);     % 划分绘图窗格（两行三列，当前取第一个窗格）
plot(n0,Tq0,'bo',n0,Pe0,'mo')       % 绘制标定出的扭矩、功率点
axis([0,10000,0,60]);       % 坐标轴范围定义
grid on;        % 显示网格线
hold on;        % 在下面的绘制过程中保持已经画出的部分（扭矩、功率点）
plot(n,Tq,n,Pe)     % 绘制扭矩、功率随转速拟合之后的曲线
ylabel('Tq(Nm),Pe(kw)')     % y轴名称
xlabel('n(r/min)')      % x轴名称
title('发动机外特性曲线')       % 图像标题% 驱动力-行驶阻力平衡图,求出最高车速
u1=0.377*n*r/ig1/i0;        % 写出不同挡位车速随转速n的表达式
u2=0.377*n*r/ig2/i0;
u3=0.377*n*r/ig3/i0;
u4=0.377*n*r/ig4/i0;
u5=0.377*n*r/ig5/i0;
u6=0.377*n*r/ig6/i0;% 驱动力表达式
Ft1=Tq*ig1*i0*nt/r;     % 写出不同挡位驱动力随转速n的表达式
Ft2=Tq*ig2*i0*nt/r;
Ft3=Tq*ig3*i0*nt/r;
Ft4=Tq*ig4*i0*nt/r;
Ft5=Tq*ig5*i0*nt/r;
Ft6=Tq*ig6*i0*nt/r;Ff=f*m*g;       % 滚阻% 各档阻力表达式
Fw1=CDA/21.15*u1.^2+Ff;     % 写出不同挡位阻力随转速n的表达式
Fw2=CDA/21.15*u2.^2+Ff;
Fw3=CDA/21.15*u3.^2+Ff;
Fw4=CDA/21.15*u4.^2+Ff;
Fw5=CDA/21.15*u5.^2+Ff;
Fw6=CDA/21.15*u6.^2+Ff;% 绘制第二张图
subplot(2,3,2);
plot(u1,Ft1,u2,Ft2,u3,Ft3,u4,Ft4,u5,Ft5,u6,Ft6,u1,Fw1,'b',u2,Fw2,'b',u3,Fw3,'b',u4,Fw4,'b',u5,Fw5,'b',u6,Fw6,'b');
grid on;        % 在第二个窗格绘制驱动力、阻力随车速变化曲线
hold on;
title('驱动力-行驶阻力平衡图')
xlabel('ua(km/h)');
ylabel('Ft(N)');
legend('Ft1','Ft2','Ft3','Ft4','Ft5','Ft6','Ff+Fw');        % 图例
hold on;
[~,ux]=min(abs(Ft6-Fw6));       % 找出驱动力和阻力曲线的近似交点，返回交点对应的下标
umax=u6(ux);        % u6此时是一个数组的形式，此时对应的车速就是u6这个数组取ux下标的值
plot(u6(ux),Ft6(ux),'m*')       % 最高车速对应点
str=['umax=',num2str(umax),'km/h'];     % 图上注释内容
text(40,750,str);       % 定义注释位置
disp('umax=')       % 在命令行窗口显示最高车速
disp(umax)% 加速度倒数曲线，积分求解加速时间
nfin=umax*ig6*i0/0.377/r;       % 定义六档加速度倒数曲线绘制区间的速度上限（上面求出的最高速度）
n6=nmin:1:round(nfin-1);        % 定义六档加速度倒数曲线绘制区间的速度区间
u6=0.377*n6*r/ig6/i0;       % 重写六档速度随转速表达式
Ft6=Tq(n6)*ig6*i0*nt/r;     % 重写六档驱动力随转速表达式
Fw6=CDA/21.15*u6.^2+Ff;     % 重写六档阻力随转速表达式% 防止五档驱动力曲线和阻力曲线有交点导致加速度产生负值
if(min(abs(Ft5-Fw5))<1)     % 判定如果阻力曲线和五档驱动力曲线有交点[~,u5x]=min(abs(Ft5-Fw5));u5max=u5(u5x);      % 找出驱动力阻力相等时对应的车速n5fin=u5max*ig5*i0/0.377/r;     n5=nmin:1:round(n5fin);     % 重新定义五档加速度倒数曲线绘制区间的速度区间u5=0.377*n5*r/ig5/i0;       % 重写表达式Ft5=Tq(n5)*ig5*i0*nt/r;Fw5=CDA/21.15*u5.^2+Ff;
end% 旋转质量换算系数
delta1=1+(Iw1+Iw2+If*ig1.^2*i0.^2*nt)/m/r.^2;       % 旋转质量换算系数
a1=1/m/delta1*(Ft1-Fw1);        % 一档时的加速度表达式
delta2=1+(Iw1+Iw2+If*ig2.^2*i0.^2*nt)/m/r.^2;
a2=1/m/delta2*(Ft2-Fw2);
delta3=1+(Iw1+Iw2+If*ig3.^2*i0.^2*nt)/m/r.^2;
a3=1/m/delta3*(Ft3-Fw3);
delta4=1+(Iw1+Iw2+If*ig4.^2*i0.^2*nt)/m/r.^2;
a4=1/m/delta4*(Ft4-Fw4);
delta5=1+(Iw1+Iw2+If*ig5.^2*i0.^2*nt)/m/r.^2;
a5=1/m/delta5*(Ft5-Fw5);
delta6=1+(Iw1+Iw2+If*ig6.^2*i0.^2*nt)/m/r.^2;
a6=1/m/delta6*(Ft6-Fw6);% 加速度取倒数
ad1=1./a1;      % 加速度取倒数
ad2=1./a2;
ad3=1./a3;
ad4=1./a4;
ad5=1./a5;
ad6=1./a6;
subplot(2,3,3)
plot(uh,adh,u2,ad2,u3,ad3,u4,ad4,u5,ad5,u6,ad6);    % 加速度倒数随车速变化曲线
% axis([0,120,0,4]);
% plot(uh,adh,u2,ad2);
% axis([0,70,0,0.4]);
hold on;
% plot(u2(1),ad2(1),'r*')
grid on;
hold on;
title('加速度倒数曲线')
xlabel('ua(km/h)')
ylabel('1/a')% 更换自变量重新拟合
uu=0:0.01:120;      % 重新定义自变量
uu2=polyfit(u2,ad2,6);      % 对加速度倒数-速度曲线在更换自变量uu之后进行重新拟合，也就是复刻，拟合阶次过高或过低都会产生过大误差，一般在3-8阶选取，6阶经过验证比较合适。
ad2=polyval(uu2,uu);        % 此时，加速度倒数曲线就是以uu为自变量的了。由于是拟合出来的，因此和原曲线没有完全重合，这一步将产生一定误差，不过如果拟合阶次选取合适，换挡
% 车速的误差在1km/h以内，求出结果的误差在0.02秒以内，可以接受。
uu3=polyfit(u3,ad3,6);
ad3=polyval(uu3,uu);
uu4=polyfit(u4,ad4,6);
ad4=polyval(uu4,uu);
uu5=polyfit(u5,ad5,7);
ad5=polyval(uu5,uu);
uu6=polyfit(u6,ad6,6);
ad6=polyval(uu6,uu);% 如果各档加速度倒数曲线没有交点，用以下模块
% [uuu23,uu23]=min(abs(ad3-ad2));       % 找2、3档交点，返回误差值和对应下标
% u23=uu(uu23);     % 将u23换挡车速从uu的数组里提取出来
% plot(uu(uu23),ad2(uu23),'mo')     %绘制换挡点，和上面绘制出的加速度倒数曲线的交点对比，可以评估在拟合这一步产生的误差
% [uuu34,uu34]=min(abs(ad4-ad3));
% u34=uu(uu34);
% plot(uu(uu34),ad3(uu34),'go')
% [uuu45,uu45]=min(abs(ad5-ad4));
% u45=uu(uu45);
% % plot(u45,ad4(uu45(end)),'ro')
% plot(uu(uu45),ad4(uu45),'ro')
% [uuu56,uu56]=min(abs(ad6-ad5));
% u56=uu(uu56);
% plot(uu(uu56),ad5(uu56),'bo')
% legend('adh','ad2','ad3','ad4','ad5','ad6','23换挡点','34换挡点','45换挡点','56换挡点')% 如果各档加速度倒数曲线有交点，用以下模块
% 寻找换挡点，重新定义各档曲线绘制范围
uuuu23=find(abs(ad3-ad2)<0.001);       % 用find函数找到交点（两条曲线相差0.001以内的点，会有很多个。阈值过大会有过多错误交点，过小可能找不到交点（因为n的取值不是连续的））
uuu23=uu(uuuu23);       % 将下标对应的值提取出来
uu23=find(uuu23>45&uuu23<70);       % 筛选交点，因为拟合出的曲线会在其他不正确的范围里产生额外错误的交点，我们只取正确范围内的
u23=uuu23(uu23(end));       % 找正确范围内的最后一个符合find中条件的值，认为是交点（有很小很小的误差，忽略不计），也就是换挡点
uuuu34=find(abs(ad4-ad3)<0.001);
uuu34=uu(uuuu34);
uu34=find(uuu34>60&uuu34<85);
u34=uuu34(uu34(end));
uuuu45=find(abs(ad5-ad4)<0.002);
uuu45=uu(uuuu45);
uu45=find(uuu45>75&uuu45<105);
u45=uuu45(uu45(end));
uuuu56=find(abs(ad6-ad5)<0.003);
uuu56=uu(uuuu56);
uu56=find(uuu56>88&uuu56<110);
u56=uuu56(uu56(end));
legend('adh','ad2','ad3','ad4','ad5','ad6');
% legend('adh','ad2','转换点','ad4','ad5','ad6');% 舍去6档
if(i0<3.2)      % 如果主减速比小于3.2，5、6档的加速度倒数曲线离的太近，会产生很多难以分辨的错误交点。因为传动比较小时将在5档完成75m，故六档可以不画，不再求5、6换挡点u56=u5(end);u6=[0 0];       % u6被定义为1*2的全零数组，不再参与下面的计算
end% 重新划分各档速度区间，为后面的积分做准备
umin=0.377*nmin*r/ig2/i0;       % 定义2档最低车速
u2=umin:0.01:u23;       % 重新定义二档图像绘制区间
ad2=polyval(uu2,u2);        % 将ad2表达式以u2为自变量重写
u3=u23:0.1:u34;
ad3=polyval(uu3,u3);
u4=u34:0.1:u45;
ad4=polyval(uu4,u4);
u5=u45:0.1:u56;
ad5=polyval(uu5,u5);
u6=u56:0.1:110;
ad6=polyval(uu6,u6);% 弹射起步区间划分，表达式重写
uhmax=0.377*nmax*r/ig2/i0;
uh=umin:0.01:u23;
FN2=fd*m*g*a_1/(L-fd*hg);
Fwh=CDA/21.15*uh.^2;
ah=(FN2-m*g*f-Fwh)/m;
adh=1./ah;% 弹射起步工况判别模块
if(min(abs(ad2-adh))<0.001)     % 如果弹射起步和2档的加速度倒数曲线有交点[~,u2uh]=min(abs(ad2-adh));     uh2=u2(u2uh);       % 找出交点对应车速[uuu2mid,uu2mid]=min(abs(ad2));u2mid=u2(uu2mid);       % 找出2档加速度曲线最低点对应的车速if u2mid<uh2       % 如果交点在u2mid之后，将这个交点忽略，认为在2档最低速时弹射滑转过程结束u2=umin:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:umin;else if u2mid>=uh2         % 如果交点在u2mid之前，弹射起步过程至交点结束，2档曲线自交点开始u2=uh2:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:uh2;endend
else if((min(abs(ad2-adh))>=0.001))     % 如果弹射起步和2档的加速度倒数曲线没有交点if(adh(1)>ad2(1))       % 如果弹射起步滑转过程加速度倒数曲线在2档加速度倒数曲线上方，认为在2档最低速时弹射滑转过程结束u2=umin:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:umin;else if(adh(1)<ad2(1))      % 如果弹射起步滑转过程加速度倒数曲线在2档加速度倒数曲线下方，最复杂的工况t2=cumtrapz(u2/3.6,ad2);        % 使用cumtrapz函数，ad2对u2做积分得出时间。由du/dt=a,得dt=du*(1/a),两边同时积分，左边是时间t，右边是1/a（加速度倒数）%对速度的积分（积分前统一单位，故u应除以3.6）Fw2=CDA/21.15*u2.^2+Ff;     % 重写2档阻力表达式T2=(1./ad2*m*delta2+Fw2)*r;     % 2档驱动力表达式，就是用ad2倒推回去了alpha=(FN2*r-T2)/(Iw2+If*ig2.^2*i0.^2*nt);      % 滑转过程中的角加速度，（弹射起步驱动力（其实就是车轮在滑转过程中的阻力）-发动机驱动力）/转动惯量deltaw=cumtrapz(t2,alpha);      % 角速度的变化值就是角加速度对时间积分w=nmax/ig2/i0*pi/30-deltaw;     % 车轮转速随时间的表达式，以2档发动机最高转的车轮转速为初速度uw=w*r*3.6;[~,u2uh]=min(abs(u2-uw));       % 寻找交点，一方面是车速增加，一方面是车轮转速降低，二者速度相等时车轮不再有大的滑移率，弹射滑转过程结束s=(size(u2));       % 将u2变量的数组大小（就是u2这个数组里有多少变量）赋给sdisp(s(end));if(u2uh==s(end))        % 这一个if结构是为了防止上述方程在给定的区间里无解，既已经到了升档的车速车轮仍在滑转，这将导致u2无解报错u2=linspace(0,0,2);     % 把u2定义成一个1*2的全零数组，防止报错，此时相当于2档升3档时滑转过程结束ad2=polyval(uu2,u2);        uh=0:0.1:u23;else if(u2uh~=s(end))       % 如果滑转结束时的车速小于u23uh2=u2(u2uh);       % 滑转结束的速度可以用上面计算出的速度（使用u2数组里的第u2uh个数据）u2=uh2:0.1:u23;     % 重新划定区间ad2=polyval(uu2,u2);uh=0:0.1:uh2;endendendendend
end% 弹射起步
FN2=fd*m*g*a_1/(L-fd*hg);      %划定弹射滑转区间后，重写弹射过程表达式
Fwh=CDA/21.15*uh.^2;
ah=(FN2-m*g*f-Fwh)/m;
adh=1./ah;% 绘制第四张图
subplot(2,3,4);
plot(uh,adh,u2,ad2,u3,ad3,u4,ad4,u5,ad5,u6,ad6);        % 划分区间之后再次绘图，对比第3、4张图，可以检查以下区间划分的精度和正确性
grid on;
hold on;
plot([0,uhmax],[1/ah20c,1/ah20c],'g')       % 绘制HRT20C在练习过程中测得的加速度数据参考线
plot([0,120],[1/g,1/g],'--k');      % 绘制重力加速度参考线
title('加速度倒数积分曲线');
xlabel('ua(km/h)');
ylabel('1/a');
legend('adh','ad2','ad3','ad4','ad5','ad6','20c参考线','1/g');% 积分求解各档及弹射起步过程加速时间
th=cumtrapz(uh/3.6,adh);         % 使用cumtrapz函数，adh对uh做积分得出时间。由du/dt=a,得dt=du*(1/a),两边同时积分，左边是时间t，右边是1/a（加速度倒数）对速度的积分%（积分前统一单位，故u应除以3.6），求得弹射起步过程时间，不难想到，曲线对横轴围成的面积越小，用时越短，这也就是为什么要在曲线交点换挡。
t2=cumtrapz(u2/3.6,ad2);
t3=cumtrapz(u3/3.6,ad3);
t4=cumtrapz(u4/3.6,ad4);
t5=cumtrapz(u5/3.6,ad5);
t6=cumtrapz(u6/3.6,ad6);
disp(th(end));      % 将每个挡位的加速时间在命令行窗口显示出来
disp(t2(end));
disp(t3(end));
disp(t4(end));
disp(t5(end));% 绘制第5张图，由于cumtrupz函数会把积分过程每一步的结果计算出来，故可直接去plot每一个时间点和相应的速度值，后面几档的别看那么长，其实只是为了去平移曲线，在时间上加上前面每
% 一段时间的值，让曲线首尾相接
subplot(2,3,5);
plot(th,uh,'k');
hold on;
plot(th(end)+t2,u2,'k');
hold on;
plot(th(end)+t2(end)+t3+trantime,u3,'k');
hold on;
plot(th(end)+t2(end)+t3(end)+t4+2*trantime,u4,'k');
hold on;
plot(th(end)+t2(end)+t3(end)+t4(end)+t5+3*trantime,u5,'k');
hold on;
plot(th(end)+t2(end)+t3(end)+t4(end)+t5(end)+t6+4*trantime,u6,'k');
hold on;
plot([th(end)+t2(end),th(end)+t2(end)+trantime],[u2(end),u2(end)],'k');
hold on;
plot([th(end)+t2(end)+t3(end)+trantime,th(end)+t2(end)+t3(end)+2*trantime],[u3(end),u3(end)],'k');
hold on;
plot([th(end)+t2(end)+t3(end)+t4(end)+2*trantime,th(end)+t2(end)+t3(end)+t4(end)+3*trantime],[u4(end),u4(end)],'k');
hold on;
plot([th(end)+t2(end)+t3(end)+t4(end)+t5(end)+3*trantime,th(end)+t2(end)+t3(end)+t4(end)+t5(end)+4*trantime],[u5(end),u5(end)],'k');
hold on;
plot([0,7],[100,100],'--b');
hold on;
grid on;
title('速度-时间曲线');
xlabel('t(s)');
ylabel('ua(km/h)');% 求解百公里加速时间
if(min(abs(u5-100))<0.2)        % 由于并不知道在4、5、6档那个挡位破百，故需要分情况找交点，如果和5档有交点[~,vvb]=min(abs(u5-100));       % 求交点在数组中的下标Tb=th(end)+t2(end)+t3(end)+t4(end)+3*trantime+t5(vvb);      % 求出相应时间
else if(min(abs(u4-100))<0.2)       % 如果和4档有交点[~,vvb]=min(abs(u4-100));Tb=th(end)+t2(end)+t3(end)+2*trantime+t4(vvb);else if(min(abs(u6-100))<0.2)       % 如果和6档有交点[~,vvb]=min(abs(u6-100));Tb=th(end)+t2(end)+t3(end)+t4(end)+t5(end)+4*trantime+t6(vvb);endend
end
str=['百公里加速时间=',num2str(Tb),'s'];
text(3,50,str);if(u2(end)==0)      % 如果u2满足了弹射起步判别模块的倒数第二种工况（259行）u2(end)=uh(end);        % 把弹射起步过程的最后速度赋给u2的最后速度
end% 换挡时间内的区间划分，方便后面积分（其实是个矩形，因为换挡过程中认为速度不变，不过还要划分速度区间，为了绘制之后的位移曲线）
step=0.001;     % 步进量
ttran=0:step:trantime;      % 划分时间段
u2end=linspace(u2(end),u2(end),trantime/step+1);        % 划分换挡时间内的积分小区间（小矩形）
u3end=linspace(u3(end),u3(end),trantime/step+1);
u4end=linspace(u4(end),u4(end),trantime/step+1);
u5end=linspace(u5(end),u5(end),trantime/step+1);% 时间对速度进行积分，得出位移
xh=cumtrapz(th,uh/3.6);     % 弹射起步滑转过程的位移
x2=cumtrapz(t2,u2/3.6);     % 2档走过的位移
x23=cumsum(u2end.*(ttran(2)-ttran(1))/3.6);     % 2-3档换挡过程走过的位移
x3=cumtrapz(t3,u3/3.6);
x34=cumsum(u3end.*(ttran(2)-ttran(1))/3.6);
x4=cumtrapz(t4,u4/3.6);
x45=cumsum(u4end.*(ttran(2)-ttran(1))/3.6);
x5=cumtrapz(t5,u5/3.6);
x56=cumsum(u5end.*(ttran(2)-ttran(1))/3.6);
x6=cumtrapz(t6,u6/3.6);% 第6张图绘制，和第5张图一样，后面几档的曲线要平移，要不然他们都是从（0，0）点开始的，由于是位移曲线，在这里曲线在纵轴方向也要平移，加上之前各档走过的位移
subplot(2,3,6);
plot(th,xh,'k');
hold on;
plot(th(end)+t2,xh(end)+x2,'k');
hold on;
plot(th(end)+t2(end)+t3+trantime,xh(end)+x2(end)+x23(end)+x3,'k');
hold on;
plot(th(end)+t2(end)+t3(end)+t4+2*trantime,xh(end)+x2(end)+x3(end)+x23(end)+x34(end)+x4,'k');
hold on;
plot(th(end)+t2(end)+t3(end)+t4(end)+t5+3*trantime,xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45(end)+x5,'k');
hold on;
plot(th(end)+t2(end)+t3(end)+t4(end)+t5(end)+t6+4*trantime,xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56(end)+x6,'k');
hold on;
plot([th(end)+t2(end),th(end)+t2(end)+trantime],[xh(end)+x2(end),xh(end)+x2(end)+x23(end)],'k');
hold on;
plot([th(end)+t2(end)+t3(end)+trantime,th(end)+t2(end)+t3(end)+2*trantime],[xh(end)+x2(end)+x3(end)+x23(end),xh(end)+x2(end)+x3(end)+x23(end)+x34(end)],'k');
hold on;
plot([th(end)+t2(end)+t3(end)+t4(end)+2*trantime,th(end)+t2(end)+t3(end)+t4(end)+3*trantime],[xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end),xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45(end)],'k');
hold on;
plot([th(end)+t2(end)+t3(end)+t4(end)+t5(end)+3*trantime,th(end)+t2(end)+t3(end)+t4(end)+t5(end)+4*trantime],[xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end),xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56(end)],'k');
hold on;
plot([0,5.8],[75,75],'--b');
hold on;
grid on;
title('位移-时间曲线');
xlabel('t(s)');
ylabel('x(m)');% 求解75m加速时间，和求百公里加速时间一样的，这里分为5种情况：在4、5、6档和在4-5、5-6换挡点达到75m
if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45(end)+x5-75))<0.2)      % 5档达到75m[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45(end)+x5-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+3*trantime+t5(xxb);
else if(min(abs(xh(end)+x2(end)+x3(end)+x23(end)+x34(end)+x4-75))<0.2)      % 4档达到75m[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x23(end)+x34(end)+x4-75));Tx=th(end)+t2(end)+t3(end)+2*trantime+t4(xxb);else if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56(end)+x6-75))<0.2)        % 6档达到75m[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56(end)+x6-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+t5(end)+4*trantime+t6(xxb);else if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45-75))<0.2)     % 4-5换挡点达到75m[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+2*trantime+step*(xxb-1);else if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56-75))<0.2)        % 5-6换挡点达到75m[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+t5(end)+3*trantime+step*(xxb-1);endendendend
end
str2=['75m加速时间=',num2str(Tx),'s'];
text(4,35,str2);% 程序说明：主要误差还是在对弹射起步的处理上，一方面忽略了空套下压力随速度变化而对弹射的加速度产生影响（我认为这个值很小可以忽略掉），一方面对于滑转过程何时结束处理的比较
% 粗糙，比如弹射的加速度倒数曲线和2档加速度倒数曲线有交点的工况，选取交点作为转换时刻是最快的，也是可以通过tc或者稍微松油门达到的，就是难度比较大，另外，过了二者的交点我们
% 看到二档曲线是更低的，这意味着更大的加速度也就是需要更大的附着力，这段程序这里默认附着系数是给够的，最好当然是建立一个2档的最大加速度有没有达到极限附着系数的判据，我这里
% 就偷懒了，一方面影响不大，一方面2档能突破极限附着系数的时候（间）并不多。其次，误差也出现在拟合和用find求交点的时候，拟合出的曲线并不是完全和原曲线贴合的，所以会对积分
% 产生一个微小的影响，可以忽略不计，用find求交点（换挡车速）的时候，不难看出也会有误差，经过测量，如果151-161行用6、6、6、7、6的阶次拟合时，各换挡车速误差一般小于1km/h，
% 粗估换挡时间会有0.005s的误差，四次换挡的话就对应一个小于0.02s的误差，还是可以接受的。为什么要把第3、4张图都画出来，就是方便对比这些误差，有离谱的误差直接能看出来，此外
% 第4张图就是积分用的图像，可以放大去看看各曲线之间，你会发现他们是没有完美接上的，这也就是前面说的用find函数时的误差，如果误差过大，可以把182-198行find函数里的范围再缩小
% （过小就找不到交点了），目前的范围已经可以达到较高的精度。总的来说，这段程序对于大多数工况是能做到误差0.1s之内的，前提是参数给的准哈哈，还是可以接受的。如果程序报错，大概
% 率是拟合出的曲线产生了错误的交点，导致换挡车速计算错了，导致有的挡位根本没有对应的速度区间（就像5-6档换挡点早于4-5档这种错误），这时候解决方案就是尝试改改151-162那段程序
% 的拟合阶次，一般主要改4档之后的，一般是往低阶去改（极少数情况改高阶），最高不建议超过8阶，最低3阶。所有图画出来都有它的意义，记得好好看看，图6换挡过程中走过的位移曲线我并
% 没有定义成黑色，也是为了能方便看出来什么时候换的档、现在是几档。
% HRT-JYZ

% 程序二drivetrain for 21C
% 车辆参数
nmin=4000;
nmax=8500;
n=nmin:0.1:nmax;
ig1=5.475;
ig2=3.8325;
ig3=2.9199;
ig4=2.3985;
ig5=2.0947;
ig6=1.9044;
nt=0.9;
r=0.2286;
f=0.013;
m=250;
g=9.8;
Cd=1.5374998;
A=1.040632;
CDA=Cd*A;
L=1.58;
a_1=0.8374;
hg=0.3;
Iw1=0.42964;
Iw2=0.42964;
If=0.009;
fd=1.4;
ah20c=0.8*g;
trantime=0.1;
n0=3500:500:9000;
Tq0=[46.4 46.2 48.9 48.9 52.6 52.5 54 53.7 50.6 46 40.2 33.8];
Pe0=[17.2 19.3 23.2 25.6 30.1 33 36.7 39.4 39.4 38.5 35.2 31.5];
imin=2.9        % 定义传动比计算下限
imax=3.8        % 定义传动比计算上限
iz=imin:0.05:imax      % 定义所需计算的传动比（建议范围小于2.9-3.8，一方面考虑到耐久高避取不到，一方面程序拟合之后的曲线会有难以去除的多余交点，会报错）
i=1;        % 循环结构中的变量，用来换行for i0=iz      % for循环计算不同主减速比下的75m加速时间uhmax=0.377*nmax*r/ig2/i0;uh=0:0.01:uhmax;FN2=fd*m*g*a_1/(L-fd*hg);Fwh=CDA/21.15*uh.^2;ah=(FN2-m*g*f-Fwh)/m;adh=1./ah;% 拟合出扭矩和功率的曲线表达式，最高阶设置为3阶tqq=polyfit(n0,Tq0,3);Tq=polyval(tqq,n);pee=polyfit(n0,Pe0,3);Pe=polyval(pee,n);% 驱动力-行驶阻力平衡图,求出最高车速u1=0.377*n*r/ig1/i0;u2=0.377*n*r/ig2/i0;u3=0.377*n*r/ig3/i0;u4=0.377*n*r/ig4/i0;u5=0.377*n*r/ig5/i0;u6=0.377*n*r/ig6/i0;Ft1=Tq*ig1*i0*nt/r;Ft2=Tq*ig2*i0*nt/r;Ft3=Tq*ig3*i0*nt/r;Ft4=Tq*ig4*i0*nt/r;Ft5=Tq*ig5*i0*nt/r;Ft6=Tq*ig6*i0*nt/r;Ff=f*m*g;Fw1=CDA/21.15*u1.^2+Ff;Fw2=CDA/21.15*u2.^2+Ff;Fw3=CDA/21.15*u3.^2+Ff;Fw4=CDA/21.15*u4.^2+Ff;Fw5=CDA/21.15*u5.^2+Ff;Fw6=CDA/21.15*u6.^2+Ff;[~,ux]=min(abs(Ft6-Fw6));umax=u6(ux);% 加速度倒数曲线，积分求解加速时间nfin=umax*ig6*i0/0.377/r;n6=nmin:1:round(nfin-1);u6=0.377*n6*r/ig6/i0;Ft6=Tq(n6)*ig6*i0*nt/r;Fw6=CDA/21.15*u6.^2+Ff;if(min(abs(Ft5-Fw5))<1)[~,u5x]=min(abs(Ft5-Fw5));u5max=u5(u5x);n5fin=u5max*ig5*i0/0.377/r;n5=nmin:1:round(n5fin);u5=0.377*n5*r/ig5/i0;Ft5=Tq(n5)*ig5*i0*nt/r;Fw5=CDA/21.15*u5.^2+Ff;enddelta1=1+(Iw1+Iw2+If*ig1.^2*i0.^2*nt)/m/r.^2;a1=1/m/delta1*(Ft1-Fw1);delta2=1+(Iw1+Iw2+If*ig2.^2*i0.^2*nt)/m/r.^2;a2=1/m/delta2*(Ft2-Fw2);delta3=1+(Iw1+Iw2+If*ig3.^2*i0.^2*nt)/m/r.^2;a3=1/m/delta3*(Ft3-Fw3);delta4=1+(Iw1+Iw2+If*ig4.^2*i0.^2*nt)/m/r.^2;a4=1/m/delta4*(Ft4-Fw4);delta5=1+(Iw1+Iw2+If*ig5.^2*i0.^2*nt)/m/r.^2;a5=1/m/delta5*(Ft5-Fw5);delta6=1+(Iw1+Iw2+If*ig6.^2*i0.^2*nt)/m/r.^2;a6=1/m/delta6*(Ft6-Fw6);ad1=1./a1;ad2=1./a2;ad3=1./a3;ad4=1./a4;ad5=1./a5;ad6=1./a6;uu=0:0.01:120;uu2=polyfit(u2,ad2,6);ad2=polyval(uu2,uu);uu3=polyfit(u3,ad3,6);ad3=polyval(uu3,uu);uu4=polyfit(u4,ad4,6);ad4=polyval(uu4,uu);uu5=polyfit(u5,ad5,7);ad5=polyval(uu5,uu);uu6=polyfit(u6,ad6,6);ad6=polyval(uu6,uu);uuuu23=find(abs(ad3-ad2)<0.001);uuu23=uu(uuuu23);uu23=find(uuu23>45&uuu23<70);u23=uuu23(uu23(end));uuuu34=find(abs(ad4-ad3)<0.001);uuu34=uu(uuuu34);uu34=find(uuu34>60&uuu34<85);u34=uuu34(uu34(end));uuuu45=find(abs(ad5-ad4)<0.002);uuu45=uu(uuuu45);uu45=find(uuu45>75&uuu45<105);u45=uuu45(uu45(end));uuuu56=find(abs(ad6-ad5)<0.003);uuu56=uu(uuuu56);uu56=find(uuu56>88&uuu56<110);u56=uuu56(uu56(end));if(i0<3.2)u56=u5(end);u6=[0 0];endumin=0.377*nmin*r/ig2/i0;u2=umin:0.01:u23;ad2=polyval(uu2,u2);u3=u23:0.1:u34;ad3=polyval(uu3,u3);u4=u34:0.1:u45;ad4=polyval(uu4,u4);u5=u45:0.1:u56;ad5=polyval(uu5,u5);u6=u56:0.1:110;ad6=polyval(uu6,u6);uhmax=0.377*nmax*r/ig2/i0;uh=umin:0.01:u23;FN2=fd*m*g*a_1/(L-fd*hg);Fwh=CDA/21.15*uh.^2;ah=(FN2-m*g*f-Fwh)/m;adh=1./ah;if(min(abs(ad2-adh))<0.001)[~,u2uh]=min(abs(ad2-adh));uh2=u2(u2uh);[uuu2mid,uu2mid]=min(abs(ad2));u2mid=u2(uu2mid);if u2mid<uh2u2=umin:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:umin;else if u2mid>=uh2u2=uh2:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:uh2;endendelse if((min(abs(ad2-adh))>=0.001))if(adh(1)>ad2(1))u2=umin:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:umin;else if(adh(1)<ad2(1))t2=cumtrapz(u2/3.6,ad2);Fw2=CDA/21.15*u2.^2+Ff;T2=(1./ad2*m*delta2+Fw2)*r;alpha=(FN2*r-T2)/(Iw2+If*ig2.^2*i0.^2*nt);deltaw=cumtrapz(t2,alpha);w=nmax/ig2/i0*pi/30-deltaw;uw=w*r*3.6;[~,u2uh]=min(abs(u2-uw));s=(size(u2));disp(s(end));if(u2uh==s(end))u2=linspace(0,0,2);ad2=polyval(uu2,u2);uh=0:0.1:u23;else if(u2uh~=s(end))uh2=u2(u2uh);u2=uh2:0.1:u23;ad2=polyval(uu2,u2);uh=0:0.1:uh2;endendendendendendFN2=fd*m*g*a_1/(L-fd*hg);Fwh=CDA/21.15*uh.^2;ah=(FN2-m*g*f-Fwh)/m;adh=1./ah;th=cumtrapz(uh/3.6,adh);t2=cumtrapz(u2/3.6,ad2);t3=cumtrapz(u3/3.6,ad3);t4=cumtrapz(u4/3.6,ad4);t5=cumtrapz(u5/3.6,ad5);t6=cumtrapz(u6/3.6,ad6);if(min(abs(u5-100))<0.2)[~,vvb]=min(abs(u5-100));Tb=th(end)+t2(end)+t3(end)+t4(end)+3*trantime+t5(vvb);else if(min(abs(u4-100))<0.2)[~,vvb]=min(abs(u4-100));Tb=th(end)+t2(end)+t3(end)+2*trantime+t4(vvb);else if(min(abs(u6-100))<0.2)[~,vvb]=min(abs(u6-100));Tb=th(end)+t2(end)+t3(end)+t4(end)+t5(end)+4*trantime+t6(vvb);endendendif(u2(end)==0)u2(end)=uh(end);endstep=0.001;ttran=0:step:trantime;u2end=linspace(u2(end),u2(end),trantime/step+1);u3end=linspace(u3(end),u3(end),trantime/step+1);u4end=linspace(u4(end),u4(end),trantime/step+1);u5end=linspace(u5(end),u5(end),trantime/step+1);xh=cumtrapz(th,uh/3.6);x2=cumtrapz(t2,u2/3.6);x23=cumsum(u2end.*(ttran(2)-ttran(1))/3.6);x3=cumtrapz(t3,u3/3.6);x34=cumsum(u3end.*(ttran(2)-ttran(1))/3.6);x4=cumtrapz(t4,u4/3.6);x45=cumsum(u4end.*(ttran(2)-ttran(1))/3.6);x5=cumtrapz(t5,u5/3.6);x56=cumsum(u5end.*(ttran(2)-ttran(1))/3.6);x6=cumtrapz(t6,u6/3.6);if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45(end)+x5-75))<0.2)[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45(end)+x5-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+3*trantime+t5(xxb);else if(min(abs(xh(end)+x2(end)+x3(end)+x23(end)+x34(end)+x4-75))<0.2)[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x23(end)+x34(end)+x4-75));Tx=th(end)+t2(end)+t3(end)+2*trantime+t4(xxb);else if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56(end)+x6-75))<0.2)[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56(end)+x6-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+t5(end)+4*trantime+t6(xxb);else if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45-75))<0.2)[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x23(end)+x34(end)+x45-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+2*trantime+step*(xxb-1);else if(min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56-75))<0.2)[~,xxb]=min(abs(xh(end)+x2(end)+x3(end)+x4(end)+x5(end)+x23(end)+x34(end)+x45(end)+x56-75));Tx=th(end)+t2(end)+t3(end)+t4(end)+t5(end)+3*trantime+step*(xxb-1);endendendendendplot(i0,Tx,'ro');       % 绘制不同主减速比的时间，在图上以红色圆圈标注出来hold on;disp(th(end));      % 在命令行窗口显示各档的时间disp(t2(end));disp(t3(end));disp(t4(end));disp(t5(end));Tz(i,:)=[Tx,0];     % 将每一次循环过程中计算出的时间储存在Tz中形成一个数组i=i+1;      % i值加1，即在上述数组中换行以准备保存下一个时间
end         % for循环结束
T=Tz(:,1);      % 调用Tz数组的第一列，因为上述循环形成的是一个两列的数组（第二列全0）
ii=imin:0.01:imax;      % 定义变量ii
tt=polyfit(iz,T',3);        % 将T的转置和iz进行3阶拟合（因为T是一列数，要转过来形成行才能和iz形成一一对应关系）
TT=polyval(tt,ii);      % 以ii为自变量，输出TT，即对于每一个ii相应的时间
[~,ttmin]=min(TT);      % 找最小值下标
ibest=ii(ttmin);        % 把最优传动比表达出来
tmin=TT(ttmin);     % 把最短时间表达出来
plot(ii,TT);        % 绘制ii和TT的关系曲线（离散点的拟合曲线）
hold on;
plot(ibest,tmin,'r*');      % 图上标注最优点
ylabel('t(s)');
xlabel('主减速比');
title('加速时间-主减速比');
disp('最优主减速比');
disp(ibest);%个别传动比下因为各种判据的巧合会产生比较离谱的值，比如明显远离拟合曲线的值，此时该结果误差较大。

基于matlab的大学生方程式汽车动力性仿真【最高车速、百公里加速时间、75m加速时间、最优传动比】相关推荐

matlab故障识别,基于Matlab的电力系统故障分析与仿真（V2.1）最新版
<基于Matlab的电力系统故障分析与仿真.doc>由会员分享,可免费在线阅读全文,更多与<基于Matlab的电力系统故障分析与仿真(V2.1)>相关文档资源请在帮帮文库(ww ...
matlab中stms和taylor,基于Matlab的电力系统故障分析与仿真V2.1（手机版）
<基于Matlab的电力系统故障分析与仿真.doc>由会员分享,可免费在线阅读全文,更多与<基于Matlab的电力系统故障分析与仿真(V2.1)>相关文档资源请在帮帮文库(ww ...
matlab中stms和taylor,基于Matlab的电力系统故障分析与仿真V2.1（模版2）
<基于Matlab的电力系统故障分析与仿真.doc>由会员分享,可免费在线阅读全文,更多与<基于Matlab的电力系统故障分析与仿真(V2.1)>相关文档资源请在帮帮文库(ww ...
matlab中stms和taylor,基于Matlab的电力系统故障分析与仿真V2.1（网络分享版）
<基于Matlab的电力系统故障分析与仿真.doc>由会员分享,可免费在线阅读全文,更多与<基于Matlab的电力系统故障分析与仿真(V2.1)>相关文档资源请在帮帮文库(ww ...
matlab频分复用,基于MATLAB的频分复用系统的仿真_.doc
基于MATLAB的频分复用系统的仿真_ 基于MATLAB的频分复用系统的仿真_毕业论文(设计) Abstract With the development of communication techn ...
matlab非线性系统频域标识,基于MATLAB的最小二乘法系统辨识与仿真
第 29卷第 2期许昌学院学报 Vol. 29. No. 2 2010年 3月 JOURNAL OF XUCHANG UN IVERSITY Mar. 2010 收稿日期: 2008 - 10 - ...
tcsc工作原理matlab仿真,基于Matlab的TCSC建模与仿真研究.doc
基于Matlab的TCSC建模与仿真研究基于Matlab的TCSC建模与仿真研究第17卷第5期 2006年1O月巾原T学院 JOURNALOFZHONGYUANINSTIT[ITEOFTECHN ...
【基于MATLAB的火灾疏散模拟仿真】——安全隐患提前发现，疏散方案优化
[基于MATLAB的火灾疏散模拟仿真]--安全隐患提前发现,疏散方案优化随着城市化进程的不断加速,人口密度越来越大,特别是在高层建筑中.万一发生火灾,往往会对人的生命和财产造成严重损失.因此,火灾疏 ...
matlab模拟Fraunhofer衍射,基于Matlab的夫琅禾费衍射光学仿真.doc
基于Matlab的夫琅禾费衍射光学仿真基于Matlab的夫琅禾费衍射光学仿真摘要计算机仿真技术是以多种学科和理论为基础,以计算机及其相应的软件为工具,通过虚拟试验的方法来分析和解决问题的一门综合 ...

基于matlab的大学生方程式汽车动力性仿真【最高车速、百公里加速时间、75m加速时间、最优传动比】

基于matlab的大学生方程式汽车动力性仿真【最高车速、百公里加速时间、75m加速时间、最优传动比】相关推荐

最新文章

热门文章