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

建立系统的动力学方程

建立系统的动力学方程的方法：1、牛顿力学：牛顿第二定律；2、分析力学：拉格朗日方程。 以x 、θ为两个变量建立二自由度系统动力学方程；3、影响系数法-张量算法 1、 根据牛顿第二定律，

2、 拉格朗日方程，

()n ,.....2,1Q i ==??-??i q L q

L dt d i

i -------------拉格朗日第二类方程

V T L -=，称为拉氏函数，泛函；

势能函数：()n q q q V V ,.......,21=；

动能函数：()n q q q

,.......,T T 21=； 广义坐标()n i i ,.....2,1q =对应的非保守力：()n i ,.....2,1Q i =

()n i q L q

L dt d i

i ,.....2,10==??-?? -------------保守系统的拉氏方程

-------------------利用上诉拉氏方程求解-------------------------------------

q q 21T T M = ；q q 2

1V T

K =

)J (mx

21T 2

2

θ +=

])l (x k )l -(x [k 2

1V 2

22211θθ++=

将L=T-V 代入拉氏方程可解的：

可见：与牛顿第二定律求得的结果一致。

M q

+K q =Q;对保守系统Q=0； ()n i q k q

m

n

j j i j ij

,.....2,1Q )(i 1

j ==+∑= M =(

J

0m )

K=(

2

2

22

111

122112221l k l k l k l k l k l k k k +--+)

q =(θx

)

M q

+K q =Q 和()n i q k q m n

j j i j ij ,.....2,1Q )(i 1j ==+∑= 有明确的物理意义：弹性恢复力-Kq 、惯性力-M q

与保守力Q 平衡。如张量理论，j i j i k m 、可认为是张量的坐标，j i k 表示：使系统仅产生沿j q 坐标的单位位移时，沿i q 坐标必须施加的外力i Q ，或者说Q 的i 分量在q 的j 分量上的影响量的投影。

3、影响系数法-张量算法

！！！！！！！下面由张量分量投影计算理论直接求解-----------！！！！！！！！ 设

m11 m12 * x” + k11 k12 * x

m21 m22 θ” + k21 k22 θ

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,可得到如下的方程：

01122

2

21122

2

2=+-++x k x

l

l l m x

l

l m ρρ

02212

2

2122

2

2

1=+-++x k x

l

l l m x

l

l m

ρρ

所以0002121212

2

2122

2122

212

222=??

?

??????????+????????????

??????

?+--+x x k k x x l l m l l l m l l l m l

l m

ρρ

ρρ 0x *x

*x 12

1211=++ωη 0x *x

*x 22

2122=++ωη

式中，2222

211ρρ

η+-=l l l 联系系数，表示两坐标之间的联系

2

2

12212ρ

ρη+-=

l l l

)

(2

2

22

11ρω+=

l m l

k 偏频，表示前后悬挂独立振动时的振动频率，即x1=0时的振动频

率是w2,x2=0时的振动频率是w1,不同于系统的固有频率（2自由度独立时才相等）。 )(2

212

22ρω+=

l m l

k

m J

=

ρ 汽车绕质心轴的回转半径

在汽车设计中，希望行车时一个悬挂的振动不传到另一个悬挂上，为此，应使车身质量分布和前后轮位置满足：

质量分配系数1212

==l l ρ

ε，这时021==ηη

11

211m k ml l k =

=

ω 22

1

22m k ml l k =

=

ω

2

2

22

112

l m l m m J +==ρ

对于一般质量分配系数1

≠ε的耦合情况，可以用模态分析法求固有频率及其通解：

)

sin()sin()sin()sin()sin()sin(2212211111222211212221211111?β?β????+++=+++=+++=t p A t p A t p A t p A x t p A t p A x

即??

?

?

????????=??????++??????=??????212122121111212111)sin()sin(11p p x x t p A t p A x x ββ??ββ

可见，特征向量阵(模态矩阵)ф组成坐标变换矩阵（由老基到新的主坐标基的坐标变换矩阵），x p =фT

x 为主坐标，主振动的坐标，在该坐标系上，各自由度独立振动（解耦）。

??

?

?????????2111ββ分别是新的主坐标基的两个基矢量在老基下的投影坐标。

在新的主坐标基下，M p =фT

M ф为主质量阵（主质量组成的对角阵）；K p =фT

K ф为主刚度阵（主刚度组成的对角阵）。这种矩阵变换的本质是张量的坐标变换。

M p x

+K p x =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)));

%完全用matlab sym解析法运算求解------------------------------------------------

syms A11 A12 phi1 phi2 t

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)*Y0 Y=[x1;x2;x1';x2']

%expm(A)的意义是将坐标先变换到主坐标系，对对角值进行exp运算后再变换到原坐标系，如同张量坐标变换help expm

y0=[x0;xd0]; %四元变量的初始条件

for i=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计算）。

用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*(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 ode45数值解------------------------------------------------

A=[zeros(2,2),eye(2);-M\K,zeros(2,2)]; %四阶参数4X4矩阵X'=AX-->X=expm(A*t)*X0 X=[x1;x2;x1';x2']

syms x1 x2 dx1 dx2

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不累加。

while i

ss=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,[0 10],[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)

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 = [];

otherwise

error(['Unhandled flag = ',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=[1 0 0 0];

str=[];

ts=[];

function sys=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;

function sys=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))