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

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

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

x₂₂=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以增加后偏频改变振型后的结果为：

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

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

实对称矩阵的特征值求法_机械振动理论(3)-解析实模态分析
模态分析是一种研究系统振动特性的分析方法,可以分为:解析模态分析和试验模态分析. 解析法,在事先知道结构的几何形状.边界条件和材料特性的前提下,将结构的质量分布.刚度分布和阻尼分布分别用质量矩阵.刚度 ...
什么是模态分析？什么是振型？
模态和振型是两个比较难懂的概念,涉及的理论比较多,我想通过一句话引出,然后通过逐步解释的方法去阐释这两个概念. 以一根梁为例,通过理论计算寻找其固有频率.阻尼比.振型的过程就是解析模态分析,通过实验得 ...
基于空间分析的铁路征地数量计算方法
摘要:为减少铁路设计项目中土地外业调查工作量,提高征地数量统计准确性,提出利用土地利用数据库和基本农田数据库,基于空间相交分析和空间擦除分析计算征地数量的方法,实现快捷.准确的铁路设计项目用地数量的计 ...
一阶电路误差分析_电动涡旋压缩机转子的模态分析及试验研究
摘要: 对制冷压缩机的曲轴结构进行模态分析,可以使压缩机的额定工作频率避开共振频率区域,从而降低压缩机的振动与噪声.为了确定几何排量为28 mL的电动汽车空调涡旋压缩机的额定转速,文章通过UG软件建立 ...
JS 调试分析 + 字体解析（汽车之家）
JS 调试分析 + 字体解析(汽车之家) 当你看到这篇文章,讲一堆理论和基础,你一定会很烦..直接开始,上图!!(需要使用一个工具:FontCreator..如何下载,自己搜!) 惊不惊 ...
python数字规律分析_【小白学爬虫】用Python分析福彩3D|发现数字的秘密
2).我们用chrome浏览器,分析一下网站的结构和源码访问: http://kaijiang.zhcw.com/zhcw/inc/3d/3d_wqhg.jsp 跳转到: http://kaijia ...
多自由度系统模态分析与试验
目录 1. 频响函数 2. 频响函数与模态参数之间的关系 2.1 频响函数的任意一行 2.2 频响函数的一列 3. 频响函数的图像 3.1 幅频曲线与相频曲线 3.2 实频曲线与虚频曲线 3.3 频响 ...
分析和解析PHP代码的7大工具
PHP已成为时下最热门的编程语言之一,然而却有许多PHP程序员苦恼找不到合适的工具来帮助自己分析和解析PHP代码.今天小编就为大家介绍几个非常不错的工具,来帮助程序员们提高自己的工作效率,一起来看看吧 ...
多模态理论张德禄_结构动力学中的模态分析(3) —— 模态参数及实验模态分析...
引言前面的文章介绍了模态相关的数学基础及实模态分析. 蒙特遇见卡罗:结构动力学中的模态分析(1) -- 线性系统和频响函数zhuanlan.zhihu.com 蒙特遇见卡罗:结构动力学中的模态分析 ...

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

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

最新文章

热门文章