汽车二自由度振动模型计算

1、kl fk2(x+12*e )建立系统的动力学方程建立系统的动力学方程的方法：1、牛顿力学：牛顿第二定律；2、分析力学：拉格朗日方程。x、。为两个变量建立二自由度系统动力学方程；3、影响系数法-张量算法1、根据牛顿第二定律,mx + kl<x-ll)+k2*(x + 12 * 0) = 02、6 - kl ll*(x-H*e)+k2*I2 » (x + 12 0) = 0拉格朗日方程，d .:L;LL 一 = Q i i =1,2,. n dt :qi ;qi拉格朗日第二类方程L =T -V ,称为拉氏函数，泛函;势能函数:=V qi ,q2 ,. qn动能函数:= T(qi,

2、q2,. qn )；广义坐标qi= 1,2,. n )对应的非保守力：Qi (i =1,2,.nd;:L;:L.-=0 i =1,2,. ndt m:4保守系统的拉氏方程利用上诉拉氏方程求解T1T-=q M q ; V21 T=-q2KqT12. 2=(mx : u J )2V1. 2=kjx -)2k 2 (x.2成)将L=T-V代入拉氏方程可解的：+k2* (x + 12 * 0) = 0(8-kl* lL*'：x-ll -0i +k2*12 * (x + 12 * 9) = 0可见：与牛顿第二定律求得的结果一致。M q +Kq =Q;对保守系统 Q=0;n寸 g q j 吟 jq

3、 J =Qj i =1,2,. nj 4m0M=()0jk1'k2k22 -k1K=(22 )k 2 l 2-k1 l 1kf - k2l2xq=(.)enM q +Kq =Q和£ (m q j +kj j q j) = Q i (i =1,2,. n )有明确的物理意义：弹性恢复力 j 文-Kq、惯性力-M q与保守力Q平衡。如张量理论，m、匕j可认为是张量的坐标， 匕j表 示：使系统仅产生沿 q j坐标的单位位移时，沿 qi坐标必须施加的外力 Qi，或者说Q的 i分量在q的j分量上的影响量的投影 。3、影响系数法-张量算法!下面由张量分量投影计算理论直接求解 !！设&qu

4、ot;mHm12 '*x”、+如1k12 '*xm21m22j 0 -+L k21k22L 9k11:仅当x动1单位时，x向的作用力为 k1*1+k2*1=k1+k2;k12:仅当0动1单位时, k21:仅当x动1单位时, k22:仅当0动1单位时,x0向的作用力为0向的作用力为向的作用力为-k1*(l1*1)+k2*(l2*1)=-k1*l1+k2*l2;-(k1*1)*l1+(k2*1)*l2=-k1*l1+k2*l2;k1*(l1*1)*l1+k2*(l2*1)*l2=k1*l1A2+k2*l2A2;m11:仅当x"动1单位时，x"向的作用力为 m*1

5、;k12:仅当0 ”动1单位时，x”向的作用力为0; %仅绕质心转动时不影响x ”向惯性力k21:仅当x”动1单位时，。”向的作用力为0; %仅平动时不影响。”向惯性 力？ ？不过质心时该怎么计算？若旋转中心偏离质心a,则变为ma°旋转中心'C质心绕质心的转动惯量为J,则绕0为J+m*a"2此时，k12:仅当。”动1单位时，x”向的作用力为 m*(a*1);k21:仅当x'动1单位时，。”向的作用力为 m*1*a;k22:仅当0 ”动1单位时，。”向的作用力为 (J+m*aA2)*1;k22:仅当0 ”动1单位时，。”向的作用力为 J*1;可见与上述结果一致

6、。对建立的动力学方程更换坐标求偏频对上述系统建立前后轮纵向位移x1、x2的动力学方程 :x1=x-l1* 0 ; x2=x+l2* 0。使用matlab 的 solve ( 'x1=x-l1* 0 ; x2=x+l2*。' x1' ,'x2')可以直接解出：0=(x1-x2)/(I1+I2);x=(x1*I2+I1*x2)/(I1+I2),代入前面创建的方程组：+k" (x + 12 * 0) = 0J6 - kl»+k2*I2 « (x + 12 0) = 0消去x1和x2,可得到如下的方程：l22m 一广 ：、2 .x1

7、l12 -m l2x2-m l1'2 2 x1 -k2x2 = 0所以l22 m l2 一 ：、2mnr"2 l1l2 - m ，2_£l2乂 I-Jk10"!°k2_ lx2 _=0Xi 1x2-22* x -2 * x2 =0式中,M -P2|22 +p2联系系数，表示两坐标之间的联系2l1；-k1l2ml 2 + 32 偏频，表示前后悬挂独立振动时的振动频率，即 m (1 2-)x1=0时的振动频率是w2,x2=0时的振动频率是 w1,不同于系统的固有频率(2自由度独立时才相等)k2l2m汽车绕质心轴的回转半径在汽车设计中，希望行车时一个悬

8、挂的振动不传到另一个悬挂上，为此，应使车身质量 分布和前后轮位置满足：质量分配系数l1l2=1，这时-=七=022 m2l对于一般质量分配系数 &尹1的耦合情况，可以用模态分析法求固有频率及其通解:X1 = An sin( pn ) - % sin( p?t 七) x2 = A21 sin( p1,:i 1) A22 sin( p2t , *2) -1 An sin( pj Q ：22 sin( p?t七)Xi11A11 sin（ Pi-:m）11xpi即虹底小12 sin（ P2t + 中2）_ B B2Xp2可见，特征向量阵（模态矩阵）巾组成坐标变换矩阵（由老基到新的主坐标基的坐标

9、变换矩阵），Xp=4）Tx为主坐标，主振动的坐标，在该坐标系上，各自由度独立振动（ 解耦）。 1 11 一"-f3分别是新的主坐标基的两个基矢量在老基下的投影坐标。在新的主坐标基下，Mp=4TM巾为主质量阵（主质量组成的对角阵）；Kp=（i）TK巾为主刚度阵（主刚度组成的对角阵）。这种矩阵变换的本质是张量的坐标变换。Mp x+Kp x =0主坐标方程组为解耦方程组。利用特征值分解找到系统的主坐标基，通过坐标变换进行解耦、简化计算、然后再变换 回去，这是坐标变换的意义所在。用matlab特征值分解法求平等与转动主模态（振型）%SH760小轿车空载主要参数m=1340;a=1.54;b=

10、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*aA2+k2*bA2;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)./m

11、ode_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<w2。求平动与平动主模态（振型）与解析法仿真计算%SH760小轿车空载主要参数clear;m=1340;a=1.54;b=1.29;l=a+b;Ic=2395; %绕质心的转动惯量 rou=sqrt(Ic/m);k1=

12、40*1000;k2=44*1000;M=m*(bA2+rouA2)/lA2,m*(a*b-rouA2)/lA2;m*(a*b-rouA2)/lA2,m*(aA2+rouA2)/lA2;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_

13、vec(:,2)=mode_vec(:,2)./mode_vec(1,2);w1=sqrt(k1*lA2)/(m*(bA2+rouA2);w2=sqrt(k2*lA2)/(m*(aA2+rouA2);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(:,

14、2)title(strcat('w_2=',num2str(omeg(2),'w_2pian=',num2str(w2_pian);%完全用 matlab sym解析法运算求解 syms A11 A12 phi1 phi2 tx1=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,&#

15、39;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

16、)+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

17、; % 时间向量A=zeros(2,2),eye(2);-MK,zeros(2,2);% 四阶参数矩阵 Y'=AY->Y=expm(A*t)*Y0Y=x1;x2;x1'x2'%expm(A)的意义是将坐标先变换到主坐标系，对对角值进行exp运算后再变换到原坐标系，如同张量坐标变换 help expmy0=x0;xd0; %四元变量的初始条件for i=1:round(tf/dt)+1%设定计算点，作循环计算tj(i)=dt*(i-1);y(:,i)=expm(A*tj(i)*y0;% 循环计算矩阵指数endsubplot(2,2,4),plot(tj,y(1,:)

18、',y(2,:)'),grid可见，坐标的选取对固有频率没有影响，但对振型有影响。W1_pianpin和w2_pianpin是前后的偏频（假设质量分配系数为1计算）。用matlab ode45()直接进行仿真计算%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*(bA2+rouA2)/lA2,m*(a*b-rouA2)/lA2;m*(a*b-rouA2)/lA2,m*(aA2+rouA2)/lA2;K=k1,0;0,

19、k2;%用 matlab ode45 数值解A=zeros(2,2),eye(2);-MK,zeros(2,2);% 四阶参数 4X4 矩阵 X'=AX->X=expm(A*t)*X0X=x1;x2;x1'x2'syms x1 x2 dx1 dx2df_sym=A*x1;x2;dx1;dx2;df_sym=subs(df_sym,'xT,'x2','dxT,'dx2','x(1)','x(2)','x(3)','x(4)');n=length(df_sy

20、m);i=1;ss=''%先定义好很重要，否则再循环体中定义时，每一循环ss不累加。while i<nss=strcat(ss,char(df_sym(i),'');i=i+1;endss=strcat(ss,char(df_sym(i);ss=strcat(ss,'');f=inline(ss,'t','x');t,x=ode45(f,0 10,1,0,0,0);% 初始 y=0,y'=1%subplot(2,2,1)plot(t,x(:,1),x(:,2) % 时间状态系列用s-function进

21、行仿真计算%sh760.m function sys,x0,str,ts=s_function(t,x,u,flag) switch flag, case 0, sys,x0,str,ts=mdlInitializeSizes; case 1, sys=mdlDerivatives(t,x,u); case 3, sys=mdlOutputs(t,x,u); case 2, 4, 9 sys =; otherwiseerror( 'Unhandled flag = ',num2str(flag);end function sys,x0,str,ts=mdlInitializeSizes sizes = simsizes; sizes.NumContStates =