汽车二自由度振动模型计算.docx
《汽车二自由度振动模型计算.docx》由会员分享，可在线阅读，更多相关《汽车二自由度振动模型计算.docx（16页珍藏版）》请在冰豆网上搜索。
汽车二自由度振动模型计算
汽车二自由度振动模型计算
建立系统的动力学方程
建立系统的动力学方程的方法：
1、牛顿力学：
牛顿第二定律；2、分析力学：
拉格朗日方程。
以x、θ为两个变量建立二自由度系统动力学方程；3、影响系数法-张量算法
1、根据牛顿第二定律，
2、拉格朗日方程，
-------------拉格朗日第二类方程
k11:
仅当x动1单位时，x向的作用力为k1*1+k2*1=k1+k2;
k12:
仅当θ动1单位时，x向的作用力为-k1*（l1*1）+k2*（l2*1）=-k1*l1+k2*l2;
k21:
仅当x动1单位时，θ向的作用力为-（k1*1）*l1+（k2*1）*l2=-k1*l1+k2*l2;
k22:
仅当θ动1单位时，θ向的作用力为k1*（l1*1）*l1+k2*（l2*1）*l2=k1*l1^2+k2*l2^2;
m11:
仅当x”动1单位时，x”向的作用力为m*1;
k12:
仅当θ”动1单位时，x”向的作用力为0;%仅绕质心转动时不影响x”向惯性力
k21:
仅当x”动1单位时，θ”向的作用力为0;%仅平动时不影响θ”向惯性力？
？
？
？
？
？
不过质心时该怎么计算？
？
若旋转中心偏离质心a,则变为ma
此时，k12:
仅当θ”动1单位时，x”向的作用力为m*（a*1）;
k21:
仅当x”动1单位时，θ”向的作用力为m*1*a;
k22:
仅当θ”动1单位时，θ”向的作用力为（J+m*a^2）*1;
k22:
仅当θ”动1单位时，θ”向的作用力为J*1;
可见与上述结果一致。
对建立的动力学方程更换坐标求偏频
对上述系统建立前后轮纵向位移x1、x2的动力学方程:
x1=x-l1*θ；x2=x+l2*θ。
使用matlab的solve（‘x1=x-l1*θ；x2=x+l2*θ’，‘x1’,‘x2’）可以直接解出:
θ=（x1-x2）/（l1+l2）;x=（x1*l2+l1*x2）/（l1+l2）,代入前面创建的方程组：
消去x1和x2,可得到如下的方程：
所以
式中，
联系系数，表示两坐标之间的联系
偏频，表示前后悬挂独立振动时的振动频率，即x1=0时的振动频率是w2,x2=0时的振动频率是w1,不同于系统的固有频率（2自由度独立时才相等）。
汽车绕质心轴的回转半径
在汽车设计中，希望行车时一个悬挂的振动不传到另一个悬挂上，为此，应使车身质量分布和前后轮位置满足：
质量分配系数
，这时
对于一般质量分配系数
的耦合情况，可以用模态分析法求固有频率及其通解：
即
可见，特征向量阵（模态矩阵）ф组成坐标变换矩阵（由老基到新的主坐标基的坐标变换矩阵），xp=фTx为主坐标，主振动的坐标，在该坐标系上，各自由度独立振动（解耦）。
分别是新的主坐标基的两个基矢量在老基下的投影坐标。
在新的主坐标基下，Mp=фTMф为主质量阵（主质量组成的对角阵）；Kp=фTKф为主刚度阵（主刚度组成的对角阵）。
这种矩阵变换的本质是张量的坐标变换。
Mp
+Kp
=0主坐标方程组为解耦方程组。
利用特征值分解找到系统的主坐标基，通过坐标变换进行解耦、简化计算、然后再变换回去，这是坐标变换的意义所在。
用matlab特征值分解法求平等与转动主模态（振型）
%SH760小轿车空载主要参数
m=1340;
a=1.54;
b=1.29;
Ic=2395;%绕质心的转动惯量
k1=40000;
k2=44000;
M=[m,0;0,Ic];
K=[k1+k2,-（k1*a-k2*b）;-（k1*a-k2*b）,k1*a^2+k2*b^2];
[eig_vec,eig_val]=eig（inv（M）*K）;
[omeg,w_order]=sort（sqrt（diag（eig_val）））;%频率
mode_vec=eig_vec（:
w_order）;%振型
T=2.*pi./omeg;%周期
mode_vec（:
1）=mode_vec（:
1）./mode_vec（1,1）;
mode_vec（:
2）=mode_vec（:
2）./mode_vec（1,2）;
subplot（2,1,1）
plot（[1;2],mode_vec（:
1））
title（strcat（'w1=',num2str（omeg
（1））））;
subplot（2,1,2）
plot（[1;2],mode_vec（:
2））
title（strcat（'w2=',num2str（omeg
（2））））;
因为对特征值进行了排序，所以w1
求平动与平动主模态（振型）与解析法仿真计算
%SH760小轿车空载主要参数
clear;
m=1340;
a=1.54;
b=1.29;
l=a+b;
Ic=2395;%绕质心的转动惯量
rou=sqrt（Ic/m）;
k1=40*1000;
k2=44*1000;
M=[m*（b^2+rou^2）/l^2,m*（a*b-rou^2）/l^2;m*（a*b-rou^2）/l^2,m*（a^2+rou^2）/l^2];
K=[k1,0;0,k2];
%用matlab特征值分解法求主振型------------------------------------------------
[eig_vec,eig_val]=eig（inv（M）*K）;
[omeg]=（sqrt（diag（eig_val）））;%频率不用sort排序
mode_vec=eig_vec;%（:
w_order）;%振型
T=2.*pi./omeg;%周期
mode_vec（:
1）=mode_vec（:
1）./mode_vec（1,1）;
mode_vec（:
2）=mode_vec（:
2）./mode_vec（1,2）;
w1=sqrt（（k1*l^2）/（m*（b^2+rou^2）））;
w2=sqrt（（k2*l^2）/（m*（a^2+rou^2）））;
w1_pian=sqrt（（k1*l）/（m*b））;
w2_pian=sqrt（（k2*l）/（m*a））;
subplot（2,2,1）
plot（[1;2],mode_vec（:
1））
title（strcat（'w_1=',num2str（omeg
（1））,';w_1pian=',num2str（w1_pian）））;
subplot（2,2,2）
plot（[1;2],mode_vec（:
2））
title（strcat（'w_2=',num2str（omeg
（2））,';w_2pian=',num2str（w2_pian）））;
%完全用matlabsym解析法运算求解------------------------------------------------
symsA11A12phi1phi2t
x1=A11*sin（omeg
（1）*t+phi1）+A12*sin（omeg
（2）*t+phi2）;
x2=mode_vec（2,1）*A11*sin（omeg
（1）*t+phi1）+mode_vec（2,2）*A12*sin（omeg
（2）*t+phi2）;
dx1=diff（x1）;
dx2=diff（x2）;
x1_0=subs（x1,'t',0）;
x2_0=subs（x2,'t',0）;
dx1_0=subs（dx1,'t',0）;
dx2_0=subs（dx2,'t',0）;
eq=[sym（strcat（char（x1_0）,'=1'））;sym（strcat（char（dx1_0）,'=0'））;
sym（strcat（char（x2_0）,'=0'））;sym（strcat（char（dx2_0）,'=0'））];
s=solve_sym（eq）;
x1=s.A11
（1）*sin（omeg
（1）*t+s.phi1
（1））+s.A12
（1）*sin（omeg
（2）*t+s.phi2
（1））;
x2=mode_vec（2,1）*s.A11
（1）*sin（omeg
（1）*t+s.phi1
（1））+mode_vec（2,2）*s.A12
（1）*sin（omeg
（2）*t+s.phi2
（1））;
ti=0:
0.02:
10;
x1i=subs（x1,'t',ti）;
x2i=subs（x2,'t',ti）;
subplot（2,2,3）
plot（ti',[x1i',x2i']）
%用matlab指数运算求解----------------------------------------------------------
x0=[1;0];xd0=[0;0];%初始条件
tf=10;dt=0.02;%时间向量
A=[zeros（2,2）,eye
（2）;-M\K,zeros（2,2）];%四阶参数矩阵Y'=AY-->Y=expm（A*t）*Y0Y=[x1;x2;x1';x2']
%expm（A）的意义是将坐标先变换到主坐标系，对对角值进行exp运算后再变换到原坐标系，如同张量坐标变换helpexpm
y0=[x0;xd0];%四元变量的初始条件
fori=1:
round（tf/dt）+1%设定计算点，作循环计算
tj（i）=dt*（i-1）;
y（:
i）=expm（A*tj（i））*y0;%循环计算矩阵指数
end
subplot（2,2,4）,plot（tj,[y（1,:
）',y（2,:
）']）,grid
可见，坐标的选取对固有频率没有影响，但对振型有影响。
W1_pianpin和w2_pianpin是前后的偏频（假设质量分配系数为1计算）。
用matlabode45（）直接进行仿真计算
%SH760小轿车空载主要参数
clear;
m=1340;
a=1.54;
b=1.29;
l=a+b;
Ic=2395;%绕质心的转动惯量
rou=sqrt（Ic/m）;
k1=40*1000;
k2=44*1000;
M=[m*（b^2+rou^2）/l^2,m*（a*b-rou^2）/l^2;m*（a*b-rou^2）/l^2,m*（a^2+rou^2）/l^2];
K=[k1,0;0,k2];
%用matlabode45数值解------------------------------------------------
A=[zeros（2,2）,eye
（2）;-M\K,zeros（2,2）];%四阶参数4X4矩阵X'=AX-->X=expm（A*t）*X0X=[x1;x2;x1';x2']
symsx1x2dx1dx2
df_sym=A*[x1;x2;dx1;dx2];
df_sym=subs（df_sym,{'x1','x2','dx1','dx2'},{'x
（1）','x
（2）','x（3）','x（4）'}）;
n=length（df_sym）;
i=1;
ss='[';%先定义好很重要，否则再循环体中定义时，每一循环ss不累加。
whileiss=strcat（ss,char（df_sym（i））,';'）;
i=i+1;
end
ss=strcat（ss,char（df_sym（i）））;
ss=strcat（ss,']'）;
f=inline（ss,'t','x'）;
[t,x]=ode45（f,[010],[1,0,0,0]）;%初始y=0,y'=1
%subplot（2,2,1）
plot（t,[x（:
1）,x（:
2）]）%时间状态系列
用s-function进行仿真计算
%sh760.m
function[sys,x0,str,ts]=s_function（t,x,u,flag）
switchflag,
case0,
[sys,x0,str,ts]=mdlInitializeSizes;
case1,
sys=mdlDerivatives（t,x,u）;
case3,
sys=mdlOutputs（t,x,u）;
case{2,4,9}
sys=[];
otherwise
error（['Unhandledflag=',num2str（flag）]）;
end
function[sys,x0,str,ts]=mdlInitializeSizes
sizes=simsizes;
sizes.NumContStates=4;
sizes.NumDiscStates=0;
sizes.NumOutputs=2;
sizes.NumInputs=1;
sizes.DirFeedthrough=0;
sizes.NumSampleTimes=0;
sys=simsizes（sizes）;
x0=[1000];
str=[];
ts=[];
functionsys=mdlDerivatives（t,x,u）
m=1340;
a=1.54;
b=1.29;
l=a+b;
Ic=2395;%绕质心的转动惯量
rou=sqrt（Ic/m）;
k1=40*1000;
k2=44*1000;
M=[m*（b^2+rou^2）/l^2,m*（a*b-rou^2）/l^2;m*（a*b-rou^2）/l^2,m*（a^2+rou^2）/l^2];
K=[k1,0;0,k2];
A=[zeros（2,2）,eye
（2）;-M\K,zeros（2,2）];%四阶参数4X4矩阵X'=AX
sys=A*x;
functionsys=mdlOutputs（t,x,u）
sys
（1）=x
（1）;
sys
（2）=x
（2）;
坐标x1、x2与主坐标x1p、x2p的关系
plot（[0;1],[0,0;mode_vec
（2）,mode_vec（4）]）
plot（y（1,1:
40）,y（2,1:
40））
升级会员 