v1=xlsread('C:\Users\wtm22923\Desktop\MATLAB\中国典型城市工况数据CCBC循环');
m=2500;%满载质量（kg）
g=9.8;%重力加速度g （m/s2)
xl=0.9;%机械效率
r=0.336;%轮胎滚动半径 （m）
Cd=0.319;%风阻系数CD
Av=2.46;%迎风面积AV （m2）
rou=1.2258;%空气密度ρ （kg/m3）
f=0.0165;%最高车速下滚动阻力系数f1
seiguma=1.18;%旋转质量换算系数δm
i0=8.74;%主减速比
ig=[1 1 1 1];
V=v1(:,2);
t=v1(:,1);%t现在是s
v=V/3.6;
efficency=0.7938;
% plot(t,V)
% figure(2)
a=0;%a要求m/s^2
for e=2:length(V);
a1=(V(e)-V(e-1))/3.6/(t(e)-t(e-1));
a=[a;a1];
end
%plotyy(t,a,t,V)
%figure(3)
b=('rgbrkmcorw'); %颜色属性
%u=['<*oxv+pd*ox.sd^<>']; % 点标记属性
d={'<-',':','-.','--','-',':','-.','--'};%线型
for c=1:length(t)
if a(c)>=0&&V(c)>0;
Fe(c)=m*g*f+m*seiguma.*a(c)+(Av*Cd.*V(c).^2)/21.15;%滚阻加风阻加加速阻力
P(c)=Fe(c).*V(c)/3600/efficency;
n(c)=V(c).*i0/0.377/r;
T(c)=Fe(c).*r/i0;
else if a(c)>=0&&V(c)==0;
Fe(c)=m*g*f+m*seiguma.*a(c)+(Av*Cd.*V(c).^2)/21.15;%滚阻加风阻加加速阻力
P(c)=Fe(c).*V(c)/3600/efficency;
n(c)=V(c).*i0/0.377/r;
T(c)=0;
else a(c)<=0&&V(c)>=0;
Fe(c)=m*g*f+(Av*Cd.*V(c).^2)/21.15;
W1(c)=((0.5*m*((v(c))^2-(v(c-1))^2))-(t(c)-t(c-1))*Fe(c).*(v(c)+v(c-1))/2)*efficency;
P(c)=0;
n(c)=V(c).*i0/0.377/r;
T(c)=0;
end
end
end
plot(t,T,'--k')
hold on
plotyy(t,P,t,V)
hold on
legend('Torque（N.m）','power(kw)','speed(km/h)')
Recycle=sum(W1)/3600000
t1=t/3600;
Q=trapz(t1,P)
Qrecycle=Q+Recycle
S=trapz(t1,V)
C=1/Q*S
E=Q*S*100
