1曲柄连杆结构参数识别.docx
《1曲柄连杆结构参数识别.docx》由会员分享，可在线阅读，更多相关《1曲柄连杆结构参数识别.docx（16页珍藏版）》请在冰豆网上搜索。
1曲柄连杆结构参数识别
1曲柄连杆结构参数识别
所要识别的参数：
R（曲柄半径）
1R=52.5e-3;%曲柄半径
2L=184e-3;%连杆长度
3Ap=0.0071;%活塞面积
4J=0.3;%曲轴转动惯量
5m=0.4;%单个缸往复质量
1.1简化的转速动力学模型公式推导（把连杆简化成两部分，一部分随活塞上下运动，令一部分随曲轴旋转）
-------------------------------------------------------------
（1）
根据
其中
为连杆转角，
为曲柄夹角。
。
根据三角公式
则
（1）变为
.
得出
，带入下式
（2）
其中
利用整个曲柄连杆机构瞬时动能相等的原则得：
（3）
为换算后的往复和旋转的质量；将
（2）带入
（2）得到瞬时转动惯量
（4）
然后（3）对
求导得到整个曲柄连杆机构的惯性扭矩T（包括往复惯性扭矩
和旋转惯性扭矩
）；
（5）
分为旋转和往复惯性扭矩两部分；往复惯性扭矩
旋转惯性扭矩
单缸发动机气力扭矩如下
对于N缸发动机按照发火相位的总和气力扭矩
惯性扭矩为
其中
，
为转动部分惯量。
根据扭矩平衡关系得到
（6）
其中摩擦扭损失压力
从而得出摩擦扭矩为
将上面几个式子带入（6）
得到二阶转速动力学模型为
（8）式是一个关于t的非线性二阶微分方程直接求解较为困难，将其转化为关于角度的函数
设
，角域内的
和时域内的
相等，则
（9）；
（10），
令
把（10），（9）带入（7）得到一阶非齐次线性微分方程
（11）
其中
1.2用matlab求求微分方程表达式
symsoRJLApmTLp1p2p3p4y;
f=sin（o）+R*sin（2*o）/（2*（L^2-（R*sin（o））^2）^.5）;
fori=1:
4
a（i,:
）=（subs（f,o,o-（i-1）*pi））^2;
end;
fsum=sum（a）;
g=diff（f,o）;
fori=1:
4
fg（i,:
）=subs（f,o,o-（i-1）*pi）*subs（g,o,o-（i-1）*pi）;
end;
fgsum=sum（fg）;
f1=f;
f2=subs（f,o,o-pi）;
P=2*m*R^2*fgsum/（J+m*R^2*fsum）;
Q=2*（Ap*R*（（p1+p3）*f1+（p2+p4）*f2）-TL）/（J+m*R^2*fsum）;
dy1=Q-P*y;
dy1=（（R^2*sin（o）^2-L^2）*（2*TL+y*（（-8*m*cos（o）*L^2*R^2*sin（o）+32*m*cos（o）*R^4*sin（o）^3-8*m*cos（o）*R^4*sin（o）+2*J*cos（o）*R^2*sin（o））/（R^2*sin（o）^2-L^2）+（2*R^2*cos（o）*sin（o）*（4*m*L^2*R^2*sin（o）^2+J*L^2-8*m*R^4*sin（o）^4+4*m*R^4*sin（o）^2-J*R^2*sin（o）^2））/（R^2*sin（o）^2-L^2）^2）-2*Ap*R*（（p1+p3）*（sin（o）+（R*sin（2*o））/（2*（L^2-R^2*sin（o）^2）^（1/2）））-（p2+p4）*（sin（o）-（R*sin（2*o））/（2*（L^2-R^2*sin（o）^2）^（1/2））））））/（4*m*L^2*R^2*sin（o）^2+J*L^2-8*m*R^4*sin（o）^4+4*m*R^4*sin（o）^2-J*R^2*sin（o）^2）
然后令k
（1）=R;k
（2）=L;k（3）=Ap;k（4）=J;k（5）=m;k（6）=TL
functiondydo=sim_model（o,y,k,p1,p2,p3,p4）
dydo=-（2*k（6）-2*k（3）*k
（1）*（（p1+p3）*（sin（o）+（k
（1）*sin（2*o））/（2*（k
（2）^2-k
（1）^2*sin（o）^2）^（1/2）））-（p2+p4）*（sin（o）-（k
（1）*sin（2*o））/（2*（k
（2）^2-k
（1）^2*sin（o）^2）^（1/2）））））/（k（4）-（4*k
（1）^2*k（5）*sin（o）^2*（k
（2）^2+cos（2*o）*k
（1）^2））/（k
（1）^2*sin（o）^2-k
（2）^2））-（2*k
（1）^2*k（5）*y*（2*k
（2）^4*sin（2*o）+（5*k
（1）^4*sin（2*o））/4-k
（1）^4*sin（4*o）+（k
（1）^4*sin（6*o））/4-2*k
（2）^2*k
（1）^2*sin（2*o）+2*k
（2）^2*k
（1）^2*sin（4*o）））/（（k
（1）^2*sin（o）^2-k
（2）^2）^2*（k（4）-（4*k
（1）^2*k（5）*sin（o）^2*（k
（2）^2+cos（2*o）*k
（1）^2））/（k
（1）^2*sin（o）^2-k
（2）^2）））;
functionw=cal_speed_value（k,x）
globalp1p2p3p4;
y0=（76.9280）^2;
yy=y0;
tspan=x;
fori=1:
length（tspan）-1
[t,y]=ode45（@sim_model,[tspan（i）,tspan（i+1）],y0,[],k,p1（i）,p2（i）,p3（i）,p4（i））;
y0=y（end,:
）;
yy=[yy;y0];
end
w=yy.^0.5;
1.3用求解微分方程的方法识别曲柄连杆结构参数
真实值：
R=0.0485;%曲柄半径m
L=0.1410;%连杆长度m
Ap=0.0058;%活塞面积m^2
J=0.23;%转动惯量kg*m^2，带负载或者冷启动时，转动惯量稍大一些，达到0.38-0.4
m=0.33;%单个缸往复质量kg
TL=17.8;%随不同工况变动
1.3.1利用'cs#1.xls'数据
clear;
globalp1p2p3p4;
data=xlsread（'cs#1.xls'）;
pdata=data（:
5:
8）;
p11=pdata（203:
400,1）;p22=pdata（203:
400,4）;p33=pdata（203:
400,3）;p44=pdata（203:
400,2）;
p1=（p11+31.4346）*6.895e3;p2=（p22+31.0164）*6.895e3;p3=（p33+29.3030）*6.895e3;p4=（p44+21.1333）*6.895e3;
owdata=data（:
2:
3）;
xdata=（owdata（203:
400,1）-owdata（203））*pi/180;
ydata=owdata（203:
400,2）;
lb=[0.02,0.034,0.003,0.15,0.15,10];
ub=[0.1,0.35,0.02,1,2,100];
k0=[400,400,400,400,400,120]';
options=optimset（'TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',200,'Algorithm','trust-region-reflective','Display','iter'）;
kp=lsqcurvefit（@cal_speed_value,k0,xdata,ydata,lb,ub,options）
kp1=[0.0485;0.14;0.0058;0.25;0.15;74]
%kp=[0.03590.05570.00870.25091.024783.1993]’
w1=cal_speed_value（kp1,xdata）;
plot（xdata,ydata,'r',xdata,w1,'b'）;
1.3.2利用'ne_norm.xls'数据
clear;
globalp1p2p3p4;
data=xlsread（'ne_norm.xls'）;
pdata=data（:
5:
8）;
p11=pdata（363:
363+228-1,1）;p22=pdata（363:
363+228-1,2）;p33=pdata（363:
363+228-1,3）;p44=pdata（363:
363+228-1,4）;
p1=（p11+15.8812+14.5）*6.895e3;p2=（p22+13.8540+14.5）*6.895e3;p3=（p33+12.7006+14.5）*6.895e3;p4=（p44+5.7273+14.5）*6.895e3;
owdata=data（:
2:
4）;
xdata=（owdata（363:
363+228-1,1）-owdata（363,1））*pi/180;
ydata=owdata（363:
363+228-1,2）;
lb=[0.03,0.034,0.003,0.15,0.15,10];
ub=[0.1,0.35,0.02,1,2,40];
k0=[0.04,0.14,0.06,0.3,0.3,30];
options=optimset（'TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter'）;
kp=lsqcurvefit（@cal_speed_value,k0,xdata,ydata,lb,ub,options）
%kp=0.03960.06830.00940.62980.150324.0868
w1=cal_speed_value（kp,xdata）;
plot（xdata,ydata,'r',xdata,w1,'b'）;
1.3.3利用ne_850.xls
clear;
globalp1p2p3p4;
data=xlsread（'ne_850.xls'）;
towdata=data（:
1:
4）;
owdata=data（:
2:
4）;
xdata=（owdata（224:
1000,1）-owdata（224,1））*pi/180;
ydata=owdata（224:
1000,2）;
pdata=data（:
5:
8）;
p11=pdata（224:
1000,1）;p22=pdata（224:
1000,2）;p33=pdata（224:
1000,3）;p44=pdata（224:
1000,4）;
p1=（p11+16.4886）*6.895e3;p2=（p22+14.6024）*6.895e3;p3=（p33+13.2834）*6.895e3;p4=（p44+4.5723）*6.895e3;
lb=[0.03,0.034,0.003,0.15,0.15,10];
ub=[0.1,0.35,0.02,1,2,40];
k0=[0.04,0.14,0.06,0.3,0.3,30];
options=optimset（'TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter'）;
kp=lsqcurvefit（@cal_speed_value,k0,xdata,ydata,lb,ub,options）
%kp=0.03200.34680.01440.51810.150029.1136
w1=cal_speed_value（kp,xdata）;
plot（xdata,ydata,'r',xdata,w1,'b'）;
1.4利用自己改进的扭矩模型进行识别。
1.5finalmodel
Tf采用。
functiony=Tf_model_new211（k,o）
y=（k
（1）+k
（2）*sin（o）+k（3）*sin（2*o）+k（4）*sin（3*o）+k（5）*sin（4*o））+k（6）*sin（o）.^2+k（7）*sin（2*o）.^2+k（8）*sin（3*o）.^2+k（9）*sin（4*o）.^2;
functiondydo=sim_model_Tf1（o,y,k,p1,p2,p3,p4）
dydo=-（2*（（k（6）+k（7）*sin（o）+k（8）*sin（2*o）+k（9）*sin（3*o）+k（10）*sin（4*o））+k（11）*sin（o）.^2+k（12）*sin（2*o）.^2+k（13）*sin（3*o）.^2+k（14）*sin（4*o）.^2）-2*k（3）*k
（1）*（（p1+p3）*（sin（o）+（k
（1）*sin（2*o））/（2*（k
（2）^2-k
（1）^2*sin（o）^2）^（1/2）））-（p2+p4）*（sin（o）-（k
（1）*sin（2*o））/（2*（k
（2）^2-k
（1）^2*sin（o）^2）^（1/2）））））/（k（4）-（4*k
（1）^2*k（5）*sin（o）^2*（k
（2）^2+cos（2*o）*k
（1）^2））/（k
（1）^2*sin（o）^2-k
（2）^2））-（2*k
（1）^2*k（5）*y*（2*k
（2）^4*sin（2*o）+（5*k
（1）^4*sin（2*o））/4-k
（1）^4*sin（4*o）+（k
（1）^4*sin（6*o））/4-2*k
（2）^2*k
（1）^2*sin（2*o）+2*k
（2）^2*k
（1）^2*sin（4*o）））/（（k
（1）^2*sin（o）^2-k
（2）^2）^2*（k（4）-（4*k
（1）^2*k（5）*sin（o）^2*（k
（2）^2+cos（2*o）*k
（1）^2））/（k
（1）^2*sin（o）^2-k
（2）^2）））;
functionw=cal_speed_value_Tf1（k,x）
globalp1p2p3p4;
y0=（71.6349）^2;
yy=y0;
tspan=x;
fori=1:
length（tspan）-1
[t,y]=ode45（@sim_model_Tf,[tspan（i）,tspan（i+1）],y0,[],k,p1（i）,p2（i）,p3（i）,p4（i））;
y0=y（end,:
）;
yy=[yy;y0];
end
w=yy.^0.5;
clear;
globalp1p2p3p4;
data=xlsread（'ne_norm.xls'）;
pdata=data（:
5:
8）;
p11=pdata（363:
1000,1）;p22=pdata（363:
1000,2）;p33=pdata（363:
1000,3）;p44=pdata（363:
1000,4）;
p1=（p11+15.8812）*6.895e3;p2=（p22+13.8540）*6.895e3;p3=（p33+12.7006）*6.895e3;p4=（p44+5.7273）*6.895e3;
owdata=data（:
2:
4）;
xdata=（owdata（363:
1000,1）-owdata（363,1））*pi/180;
ydata=owdata（363:
1000,2）;
k0=[0.0485,0.14,0.0058,0.23,0.4,-10,0.1,0.1,0.1-1.4,1,-0.1,1,1];
lb=[0.02,0.034,0.003,0,0,-1e2,-1e3,-1e3,-100,-100,-100,-100,-100,-100];
ub=[0.1,0.35,0.026,1,2,1e2,1e3,1e3100,100,100,100,1e2,100];
options=optimset（'TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',400,'Algorithm','trust-region-reflective','Display','iter'）;
[kp,resnorm]=lsqcurvefit（@cal_speed_value_Tf1,k0,xdata,ydata,lb,ub,options）
kp=[0.04570.15130.00640.35500.10544.913319.45499.19262.73780.3356-0.8043-0.777720.96801.0000];
resnorm=15.6801;
w1=cal_speed_value_Tf1（kp,xdata）;
plot（xdata,ydata,'r',xdata,w1,'b'）;
2只识别J,m两个参数
2.1把Tf当做常数识别（效果不好）
functiondydo=two_para_model（o,y,k,p1,p2,p3,p4）
R=48.5e-3;%曲柄半径m
L=141e-3;%连杆长度m
Ap=0.0058;%活塞面积m^2
dydo=-（2*k（3）-2*Ap*R*（（p1+p3）*（sin（o）+（R*sin（2*o））/（2*（L^2-R^2*sin（o）^2）^（1/2）））-（p2+p4）*（sin（o）-（R*sin（2*o））/（2*（L^2-R^2*sin（o）^2）^（1/2）））））/（k
（1）-（4*R^2*k
（2）*sin（o）^2*（L^2+cos（2*o）*R^2））/（R^2*sin（o）^2-L^2））-（2*R^2*k
（2）*y*（2*L^4*sin（2*o）+（5*R^4*sin（2*o））/4-R^4*sin（4*o）+（R^4*sin（6*o））/4-2*L^2*R^2*sin（2*o）+2*L^2*R^2*sin（4*o）））/（（R^2*sin（o）^2-L^2）^2*（k
（1）-（4*R^2*k
（2）*sin（o）^2*（L^2+cos（2*o）*R^2））/（R^2*sin（o）^2-L^2）））;
functionw=cal_speed_value_two（k,x）
globalp1p2p3p4;
y0=（71.6349）^2;
yy=y0;
tspan=x;
fori=1:
length（tspan）-1
[t,y]=ode45（@sim_model,[tspan（i）,tspan（i+1）],y0,[],k,p1（i）,p2（i）,p3（i）,p4（i））;
y0=y（end,:
）;
yy=[yy;y0];
end
w=yy.^0.5;
clear;
globalp1p2p3p4;
data=xlsread（'ne_norm.xls'）;
pdata=data（:
5:
8）;
p11=pdata（363:
363+228-1,1）;p22=pdata（363:
363+228-1,2）;p33=pdata（363:
363+228-1,3）;p44=pdata（363:
363+228-1,4）;
p1=（p11+15.8812）*6.895e3;p2=（p22+13.8540）*6.895e3;p3=（p33+12.7006）*6.895e3;p4=（p44+5.7273）*6.895e3;
owdata=data（:
2:
4）;
xdata=（owdata（363:
363+228-1,1）-owdata（363,1））*pi/180;
ydata=owdata（363:
363+228-1,2）;
k0=[0.3,0.3,30];
lb=[-1,-2,-100];
ub=[1,2,100];
options=optimset（'TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflective','Display','iter'）;
kp=lsqcurvefit（@cal_speed_value_two,k0,xdata,ydata,lb,ub,options）
w1=cal_speed_value_two（kp,xdata）;
plot（xdata,ydata,'r',xdata,w1,'b'）;
方法2.2把Tf用自己所建立的模型（效果也不好）
Tf采用。
functiony=Tf_model_new211（k,o）
y=（（k（3）+k（4）*sin（o）+k（5）*sin（2*o）+k（6）*sin（3*o）+k（7）*sin（4*o））+k（8）*sin（o）^2+k（9）*sin（2*o）^2+k（10）*sin（3*o）^2+k（11）*sin（4*o）^2）;
functiondydo=two_para_model_Tf（o,y,k,p1,p2,p3,p4）
R=48.5e-3;%曲柄半径m
L=141e-3;%连杆长度m
Ap=0.0058;%活塞面积m^2
dydo=-（2*（（k（3）+k（4）*sin（o）+k（5）*sin（2*o）+k（6）*sin（3*o）+k（7）*sin（4*o））+k（8）*sin（o）^2+k（9）*sin（2*o）^2+k（10）*sin（3*o）^2+k（11）*sin（4*o）^2）-2*Ap*R*（（p1+p3）*（sin（o）+（R*sin（2*o））/（2*（L^2-R^2*sin（o）^2）^（1/2）））-（p2+p4）*（sin（o）-（R*sin（2*o））/（2*（L^2-R^2*sin（o）^2）^（1/2）））））/（k
（1）-（4*R^2*k
（2）*sin（o）^2*（L^2+cos（2*o）*R^2））/（R^2*sin（o）^2-L^2））-（2*R^2*k
（2）*y*（2*L^4*sin（2*o）+（5*R^4*sin（2*o））/4-R^4*sin（4*o）+（R^4*sin（6*o））/4-2*L^2*R^2*sin（2*o）+2*L^2*R^2*sin（4*o）））/（（R^2*sin（o）^2-L^2）^2*（k
（1）-（4*R^2*k
（2）*sin（o）^2*（L^2+cos（2*o）*R^2））/（R^2*sin（o）^2-L^2）））;
functionw=cal_speed_value_two_Tf（k,x）
globalp1p2p3p4;
y0=（71.6349）^2;
yy=y0;
tspan=x;
fori=1:
length（tspan）-1
[t,y]=ode45（@two_para_model_Tf,[tspan（i）,tspan（i+1）],y0,[],k,p1（i）,p2（i）,p3（i）,p4（i））;
y0=y（end,:
）;
yy=[yy;y0];
end
w=yy.^0
升级会员 