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

% 程序一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=uh2:0.1:u23;
            ad2=polyval(uu2,u2);
            uh=0:0.1:uh2;
        end
    end
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)

% 程序二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;
    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;
    
    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];
    end
    
    umin=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=uh2
                u2=uh2:0.1:u23;
                ad2=polyval(uu2,u2);
                uh=0:0.1:uh2;
            end
        end
    else 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)