《有限元方法与MATLAB程序设计》 课件 11自由振动与弹性稳定

第11章自由振动与弹性稳定§11.1刚架自由振动§11.2刚架弹性稳定§11.1自由振动§11.1.1质量矩阵的推导特征值问题非零解条件（频率方程）外力虚功内力虚功虚功原理惯性力结点力质量矩阵弹性矩阵假设自由振动211.1.1质量矩阵的推导质量矩阵3例11-1

求图示刚架结构自振频率。单元①和②局部坐标下刚度矩阵为单位长度质量4例11-14

3②1

2

①124323232323③单元③刚度矩阵5结构刚度矩阵例11-1

求图示刚架结构自振频率。①和②单元质量矩阵为单位长度质量6例11-14

3③1

2

①124323232323②③单元质量矩阵7结构质量矩阵例11-1带入参数求解出自振频率求图示刚架结构自振频率。为单位长度质量8§11.1.2质量矩阵集中质量矩阵9例11-2function

FramDyna

%刚架模态分析gxy=[0,0;

0,2;

0,4;

3,4;

6,4;

6,2;

6,0];ndel=[1,2;2,3;3,4;4,5;5,6;6,7];nd=size(gxy,1);ne=size(ndel,1);A1=0.01;I1=833.33e-8;

A2=0.015;I2=2812.5e-8;AI=[A1,I1;A1,I1;A2,I2;A2,I2;A1,I1;A1,I1];Em=2.5e10;r=2.5e3;dofix=[1:3,19:21];dofree=setdiff(1:3*nd,dofix);4m6m等截面匀质钢架，弹性模量20GPa，密度2.5kg/m3。其他尺寸如图和表所示。横截面面积/m2横截面惯性矩/m4竖杆0.01833.33e-8横杆0.0152812.5e-8③①%弹性模量和质量密度%位移约束自由度1%无约束自由度24④

5⑤6⑥7x3m2m%结点坐标%单元信息y%

结点总数

3%

单元总数

2m

②3m10形成质量阵和刚度阵K=zeros(3*nd,3*nd);M=zeros(3*nd,3*nd);for

el=1:nexy=gxy(ndel(el,:),:);%刚度矩阵%质量矩阵%单元结点号和结点坐标N=kron(3*ndel(el,:),[1,1,1])-[2,1,0,2,1,0];%单元自由度%单元刚度矩阵%组装结构刚度矩阵%单元质量矩阵%组装结构质量矩阵[ke,T]=FramStif(xy,AI(el,:));K(N,N)=K(N,N)+Em*T’*ke*T;me=FramMass(xy,AI(el,1));M(N,N)=M(N,N)+r*T"*me*T;end%求自振频率w=sqrt(1./diag(D))";disp(sprintf("w=%7.3f,%7.3f,%7.3f",w));w=

14.858,

46.481,117.884求解和输出[V,D]=eigs(M(dofree,dofree),K(dofree,dofree),

5);11刚架质量矩阵function

m=FramMass(xy,A) %

刚架质量矩阵dl=xy(2,:)-xy(1,:);l=sqrt(dl*dl");

%杆长m=[140,0,0,70,0,0;0,156,22*l,0,

54,-13*l;0,

22*l,

4*l*l,0,

13*l,-3*l*l;..70,0,0,140,0,0;0,

54,13*l,0,156,-22*l;0,-13*l,-3*l*l,0,-22*l,

4*l*l];m=m*A*l/420;

%单元质量阵12§11.2弹性稳定§11.2.1弹性稳定方程特征值问题分支点失稳增量形式的刚度方程失稳切线刚度矩阵小变形有非零解的条件临界力方程无干扰时有干扰时13§11.2弹性稳定§11.2.1弹性稳定方程临界力方程刚架刚度矩阵刚架刚度矩阵（忽略轴向变形）刚架几何刚度矩阵忽略轴向变形14例11-3求临界力结构标识结点位移向量(3)单元刚度矩阵单元②③单元①(4)结构刚度矩阵15例11-3求临界力结构标识结点位移向量(5)单元几何刚度矩阵(6)结构几何刚度矩阵16(7)结构稳定方程临界力17例11-4求临界力③①124yxF3②0.5l0.5l0.5l0.5l④

518function

FramBuckgxy=[0,0;0,0.5;0,1;0.5,1;1,1];%结点坐标%单元信息%结点总数%单元总数ndel=[1,2;2,3;3,4;4,5];nd=size(gxy,1);ne=size(ndel,1);AI=repmat([0.04;1.2e-4],1,ne);dofix=[1:3,13:15];dofree=setdiff(1:3*nd,dofix);Em=1;F=

zeros(3*nd,1);F(8)=-1;%杆件横截面面积和惯性矩%位移约束自由度%无约束自由度%弹性模量%结点力列向量%置结点力杆件截面相同A=400cm2，I=1200cm2。19刚度矩阵，静力分析K=

zeros(3*nd,3*nd);for

el=1:ne[ke,T]=FramStif(gxy(ndel(el,:),:),AI(:,el));N=kron(3*ndel(el,:),[1,1,1])-[2,1,0,2,1,0];%单元自由度%组装总刚度矩阵K(N,N)=K(N,N)+Em*T’*ke*T;endU=zeros(3*nd,1);U(dofree)=K(dofree,dofree)\F(dofree);%结点位移几何刚度矩阵和求解临界力Kg=

zeros(3*nd,3*nd);Fe=zeros(6,1);for

el=1:neN=kron(3*ndel(el,:),[1,1,1])-[2,1,0,2,1,0];%单元自由度[ke,T]=FramStif(gxy(ndel(el,:),:),AI(:,el));Fe=Em*ke*T*U(N);

%杆端力fsprintf("%4i%4i%14.4g%14.4g%14.4g\n%8i%14.4g%14.4g%14.4g\n",i,ndel(el,1),Fe(1:3),ndel(el,2),Fe(4:6)));%单元几何刚度矩阵%组装结构刚度矩阵kg=FramGeSt(gxy(ndel(el,:),:));Kg(N,N)=Kg(N,N)+Fe(4)*

T"*

kg*T;end[V,D]=eigs(K(dofree,dofree),-Kg(dofree,dofree),1,"SM");fprintf("Fcr=%7.3f\n",D));20刚架几何矩阵function

ks=FramGeSt(xy)dl=xy(2,:)-xy(1,:);l=sqrt(dl*dl");%杆长a=atan2(dl(2),dl(1));%杆轴线方向角kg=zeros(6,6);kg(2:3,2:3)=[

12/l,

1;

1,4*l/3];kg(2:3,5:6)=[-12/l,

1;-1,

-l/3];kg(5:6,2:3)=[-12/l,-1;

1,

-l/3];kg(5:6,5:6)=[

12/l,-1;-1,4*l/3];21特征值问题MATLAB代码[V,D]

=

eigs(A,B)MATLAB代码[V,D]=eigs(K,-Kg)22桁架结构稳定方程单元刚度矩阵在第3章式(3.28)给出几何刚度矩阵在第10章式(10.56)给出23算例-桁架function

TrusBuckgxy=[0,0;1,0;0,1];ndel=[1,3;3,2];nd=size(gxy,1);ne=size(ndel,1);%桁架的整体失稳临界力%结点坐标*%单元信息，即每个单元的结点号*%结点总数%单元总数EA=[1,1];dofix=1:4;%抗拉刚度EA*%约束自由度*dofree=setdiff(1:2*nd,dofix);

%

无约束自由度F=zeros(2*nd,1); %

结点力列向量F(6)=-1;

%置结点力K=sparse(2*nd,2*nd);kg=zeros(4,4);kg([2,4],[2,4])=[1,-1;-1,1];k0=zeros(4,4);k0([1,3],[1,3])=[1,-1;-1,1];%几何刚度矩阵%单元刚度矩阵②①13yxFl2l24算例-桁架for

el=1:ne[T,l]=TrusRota(gxy(ndel(el,:),:));N(2:2:4)=2*ndel(el,:);N(1:2:4)=N(2:2:4)-1;%单元自由度%组装总刚度矩阵K(N,N)=K(N,N)+EA(el)/l*T"*k0*T;endU=zeros(2*nd,1);U(dofree)=K(dofree,dofree)\F(dofree);Kg=sparse(2*nd,2*nd);for

el=1:ne%结点位移N(2:2:4)=2*ndel(el,:);N(1:2:4)=N(2:2:4)-1;%单元自由度[T,l]=TrusRota(