SH760模态分析-多种解析与数字计算方法

 

 

详细分析参考：https://wenku.baidu.com/view/6915c300f08583d049649b6648d7c1c708a10b62

 

 

 

主振型与固有频率；偏频（假设分配系数=1）

 

模态分析：可以清楚地看到模态和主模态之间的关系！！！！！！

 

当取初值x ₁₀=1时，仿真的结果为：

x₁=A₁₁sin(ω₁t+φ₁)+A₁₂sin(ω₂t+φ₂)

 =0.8275 sin(8.4293t+π/2)+0.1725sin(7.8622t+π/2)

 

 

x₂=β₁A₁₁sin(ω₁t+φ₁)+ β₂A₁₂sin(ω₂t+φ₂)

=-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 ₁₁=0.8275 sin(8.4293t+π/2)

x _{2 1}=0.4353*0.8275 sin(8.4293t-π/2)

在ω ₁=8.4293的固有频率下， 主振动相位相反，振幅比为0.4353，X ₁₁振动为主，在振型图上， 节点（振幅为0的点）在X₁₁ 、x₂₁之间；

 

X ₁₂=0.1725sin(7.8622t+π/2)

x _{2 2}=2.0883*0.1725sin(7.8622t+π/2)

在ω ₁=7.8622的固有频率下， 主振动相位相同，振幅比为2.0883，X ₂₂振动为主；

 

 

 

 

用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 ₁₀=1时，输入能量0.5*k1*1^2=0.5*40000*1^2=20000, x ₂振幅最大为0.7189，势能为0.5*44000*0.7189^2=11370；传到x ₂的能量为56.85%

 

初值x ₂₀=1时，输入能量0.5*k2*1^2=0.5*44000*1^2=22000, x ₁振幅最大为0.7908，势能为0.5*40000*7908^2=12507；传到x ₁的能量为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 ₁₀=1、x ₂₀=1时，x ₂(与低偏频对应）振幅会增加；

当x ₁₀=1降到0.48，x ₂₀=1不变，x ₁x ₂振幅均不变，以低频振动，如下图；

 

 

 

当x ₁₀=0.48继续下降到0，x ₂₀=1不变，x ₁振幅增大，x ₂振幅减小，如最上面的图；

当x ₁₀=0继续下降到-2.3，x ₂₀=1不变，x ₁x ₂振幅均不变，以高频振动，如下图；

 

 

当x ₁₀=-2.3继续下降，x ₂₀=1不变，x ₂振幅增大，x ₁振幅减小；

 

总之，x₂₀=1不变，x₁₀在0.48以上时，x₂(与低偏频对应）振幅增大，

x₁₀在0.48到-2.3之间，x₂(与低偏频对应）振幅减小，

x₁₀在-2.3以下时，x₂(与低偏频对应）振幅增大。

其原因从振型图上可以看出来！！！1/0.48=2.0833约等于β_2；1/（-2.3）=-0.4348约等于β₁

 

 

 

 

当增加后刚度k2到54000以增加后偏频改变振型后的结果为：

 

 

可见： 能量传递与初始值有关？？？？？？初始值在平衡点附近倾向于接受能量，而远离平衡点则倾向于发送能量，从而使系统稳定并在平衡点附近振动。