主振型与固有频率;偏频(假设分配系数=1)
模态分析:可以清楚地看到模态和主模态之间的关系!!!!!!
当取初值x
10=1时,仿真的结果为:
x1=A11sin(ω1t+φ1)+A12sin(ω2t+φ2)
=0.8275 sin(8.4293t+π/2)+0.1725sin(7.8622t+π/2)
x2=β1A11sin(ω1t+φ1)+ β2A12sin(ω2t+φ2)
=-0.4353*0.8275 sin(8.4293t+π/2)+2.0883*0.1725sin(7.8622t+π/2)
=0.4353*0.8275 sin(8.4293t-π/2)+2.0883*0.1725sin(7.8622t+π/2)
主振动(不排序):系统以某一阶固有频率按其相应的主振型进行振动。
主振动在原坐标基上的投影为:
X
11=0.8275 sin(8.4293t+π/2)
x
2
1=0.4353*0.8275 sin(8.4293t-π/2)
在ω
1=8.4293的固有频率下,
主振动相位相反,振幅比为0.4353,X
11振动为主,在振型图上,
节点(振幅为0的点)在X11 、x21之间;
X
12=0.1725sin(7.8622t+π/2)
x
2
2=2.0883*0.1725sin(7.8622t+π/2)
在ω
1=7.8622的固有频率下,
主振动相位相同,振幅比为2.0883,X
22振动为主;
用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<w2
。
求平动与平动主模态(振型)与解析法仿真计算
%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<n
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);
%%%%%%%%%%--------------------------------------------------------------
附上模态做图代码:
subplot(2,2,1)
plot(tj,[y(1,:)',y(2,:)']),grid
xlabel('t');ylabel('x_1,x_2');
subplot(2,2,2)
plot3(y(1,:)',y(2,:)',tj),grid
xlabel('x_1');ylabel('x_2');zlabel('t');
subplot(2,2,3)
plot(y(1,:)',y(3,:)'),grid
xlabel('x_1');ylabel('x_1''');
subplot(2,2,4)
plot(y(2,:)',y(4,:)'),grid
xlabel('x_2');ylabel('x_2''');
能量传递分析:
初值x
10=1时,输入能量0.5*k1*1^2=0.5*40000*1^2=20000, x
2振幅最大为0.7189,势能为0.5*44000*0.7189^2=11370;传到x
2的能量为56.85%
初值x
20=1时,输入能量0.5*k2*1^2=0.5*44000*1^2=22000, x
1振幅最大为0.7908,势能为0.5*40000*7908^2=12507;传到x
1的能量为56.85%。
故通过系统传递的能量比例不变。
但是当前后都输入能量时,
能量传递有方向性?????
当前后位移均输入1单位位移时,看看系统振动与能量传递的情况,再增加后偏频(使k2由44000增加到54000)使振型改变,看看有什么影响。
%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(3,3,1)
plot([1;2],mode_vec(:,1))
title(strcat('w_1=',num2str(omeg(1)),';w_1pian=',num2str(w1_pian)));
subplot(3,3,2)
plot([1;2],mode_vec(:,2))
title(strcat('w_2=',num2str(omeg(2)),';w_2pian=',num2str(w2_pian)));
x0=[1;1];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(3,3,3)
plot([0;1],[0,0;mode_vec(2),mode_vec(4)]),hold
plot(y(1,1:400),y(2,1:400))
xlabel('x_1');ylabel('x_2')
subplot(3,3,4),plot(tj,[y(1,:)',y(2,:)']),grid
xlabel('t');ylabel('x_1,x_2')
subplot(3,3,5)
plot3(y(1,:)',y(2,:)',tj),grid
xlabel('x_1');ylabel('x_2');zlabel('t');
subplot(3,3,7)
plot(y(1,:)',y(3,:)'),grid
xlabel('x_1');ylabel('x_1''');
subplot(3,3,8)
plot(y(2,:)',y(4,:)'),grid
xlabel('x_2');ylabel('x_2''');
可见当x
10=1、x
20=1时,x
2(与低偏频对应)振幅会增加;
当x
10=1降到0.48,x
20=1不变,x
1x
2振幅均不变,以低频振动,如下图;
当x
10=0.48继续下降到0,x
20=1不变,x
1振幅增大,x
2振幅减小,如最上面的图;
当x
10=0继续下降到-2.3,x
20=1不变,x
1x
2振幅均不变,以高频振动,如下图;
当x
10=-2.3继续下降,x
20=1不变,x
2振幅增大,x
1振幅减小;
总之,x20=1不变,x10在0.48以上时,x2(与低偏频对应)振幅增大,
x10在0.48到-2.3之间,x2(与低偏频对应)振幅减小,
x10在-2.3以下时,x2(与低偏频对应)振幅增大。
其原因从振型图上可以看出来!!!1/0.48=2.0833约等于β2;1/(-2.3)=-0.4348约等于β1
当增加后刚度k2到54000以增加后偏频改变振型后的结果为:
可见:
能量传递与初始值有关??????初始值在平衡点附近倾向于接受能量,而远离平衡点则倾向于发送能量,从而使系统稳定并在平衡点附近振动。