结构动力学大作业
对于如下结构，是研究质量块的质量变化和在简支梁上位置的变化对整个系统模态的影响。
1
以上为一个简支梁结构。集中质量块放于梁上，质量块距简支梁的左端点距离为L.
将该简支梁简化为欧拉伯努利梁，并离散为N 个单元。每个单元有两个节点，四个自由度。
单元的节点位移可表示为：
]1122,,,e v v δθθ⎡=⎣
则单元内一点的挠度可计作：
带入边界条件：
1
3
32210)(x a x a x a a x v +++=0
1)0(a v x v ===3
322102)(L a L a L a a v L x v +++===1
10
d d a x v
x ===θ2
321232d d L a L a a x
v L
x ++===θ1
0v a =
[]12
3
4N N N N N =
建立了单元位移模式后，其动能势能均可用节点位移表示。单元的动能为：
00111()222
l l T T T
ke e e e e y E dx q
N Ndxq q mq t ρρ∂===∂⎰⎰ 其中m 为单元质量阵，并有：
l
T m N Ndx ρ=⎰
带入公式后积分可得：
222215622541322413354
1315622420133224l l l l l l l m l l l l l l ρ-⎡⎤
⎢⎥-⎢⎥=
⎢⎥
-⎢⎥
---⎣
⎦
单元势能可表示为
22
200
11()()22
2
T
l l
T
T e pe e e e q y E EI dx EI N N dxq q Kq x ∂''''==
=∂⎰⎰ 其中K 为单元刚度矩阵，并有
()l
T K EI N N dx ''''=⎰
2
23
2212
612664621261266264l l l l l l EI k l l l l l l l -⎡⎤⎢⎥-⎢⎥=⎢⎥---⎢⎥-⎣⎦
以上为单元类型矩阵，通过定义全局位移矩阵，可以得到系统刚度矩阵和系统质量矩
1
1θ=a )2(1)(3211222θθ+--=L
v v L a )(1)(22122133θθ++-=
L
v v L a 1232133222231)(θ⎪⎪⎭⎫
⎝⎛+-+⎪⎪⎭⎫ ⎝⎛+-=L x L x x v L x L x x v 2
2232332223θ⎪⎪⎭
⎫ ⎝⎛-+⎪⎪⎭⎫ ⎝⎛-+L x L x v L x L x 2
4231211)()()()()(θθx N v x N x N v x N x v +++=
阵。
当集中质量块加到简支梁上时，可以认为系统的刚度矩阵不变，而质量矩阵在 有基础上有所增加。
当集中质量块位于第i 个单元内，距单元左节点的位移为0x ，则其动能可表示
为：
000()()2
T T
e e m V N x N x δδ=
则质量块附加到总体刚度矩阵的子矩阵尾为：
000()()T m N x N x
例如当质量为0m 的质量块位于单元的右节点上时，其质量阵可表示为
000000
00000100
00
0m m ⎡⎤⎢⎥⎢
⎥'=⎢⎥⎢⎥⎣⎦
根据以上理论，编写matlab 程序。（相关参数：密度
7860KG/M3， 长度
L=1m ，截面尺寸0.02m*0.02m ）
以下为matlab 源代码： 主程序
clear; clc;
Beam_InputData541; %输入相关参数，划分了40个单 %元
k=zeros(No_nel*No_dof,No_nel*No_dof); %初始化单元刚度阵和质量阵 m=zeros(No_nel*No_dof,No_nel*No_dof); kk=zeros(Sys_dof,Sys_dof);
mm=zeros(Sys_dof,Sys_dof); %初始化总体质量阵和刚度阵
index=zeros(No_nel*No_dof,1); for
iel=1:No_el nd(1)=iel; nd(2)=iel+1; leng=0.025;
[k,m]=BeamElement11(prop,leng); %调用子函数形成单元刚度矩阵 %质量矩阵
index=femEldof(nd,No_nel,No_dof);
kk=femAssemble1(kk,k,index);
mm=femAssemble1(mm,m,index); %调用子函数形成总体刚度矩阵 end %和质量矩阵
Mzhiliangkuai=1/5*prop(3)*prop(6)*1; %将集中质量块加进总体质量矩
Ms=zeros(82,82); %之中
Ms(41,41)=Mzhiliangkuai; %此为质量块为于梁中央时的处mm=Ms+mm; 理方法
kk(1,:)=[]; %引入边界条件，划去相应的kk(80,:)=[]; %行和列
kk(:,1)=[];
kk(:,80)=[];
mm(1,:)=[];
mm(80,:)=[];
mm(:,1)=[];
mm(:,80)=[];
[V,D]=eig(kk,mm); %求解特征值和特征向量
%或使用雅克比迭代求解特征值与特征向量，以下为代码
%function[aa,v]=jac(a)
%此程序用jacobi方法求实对称矩阵的全部特征值和特征向量
%输入x：n x n矩阵
%----------------------------------------------------%
% n=length(a);
% aa=zeros(1,n);
% u=zeros(n,n);
% l=0;
% v是n阶单位矩阵
% v=eye(n,n);
% while (1)
% l=l+1;
% fm=0.0;
% for i=1:n
% for j=1:n
% if i==j
% continue;
% else
% d=abs(a(i,j));
% if d>fm
% fm=d;
% p=i;
% q=j;
% end
% end
% end
