beamdeflection.m

可以计算两端固定的平面梁的挠度

%写了一个平面梁单元的例子，相当于ansys里面的beam3单元。当划分的份数逐渐增大时，计算得到的梁的挠度值逐渐趋于稳定，不在变化。 
%simple finite element program,beam element,采用梁单元，考虑轴向变形,相当于ansys中的beam3 
%xiao yong, erc lab, hanyang university 
clear 
clc 
elmu=1*10^5 ;
widh=5 ;
height=10 ;
area=widh*height; 
iner=(widh*height^3)/12 ;
%总长度为length，划分的份数为elnum份，ellen是每个单元的长度 
length=100 ;
elnum=15 ;
ellen=length/elnum ;
nonum=elnum+1 ;
nofreedom=3 ;
freedom=nonum*nofreedom 
ke=[elmu*area/ellen,0,0,-elmu*area/ellen,0,0; 
    0,12*elmu*iner/(ellen^3),6*elmu*iner/(ellen^2),0,-12*elmu*iner/(ellen^3),6*elmu*iner/(ellen^2); 
    0,6*elmu*iner/(ellen^2),4*elmu*iner/ellen,0,-6*elmu*iner/(ellen^2),2*elmu*iner/ellen; 
    -elmu*area/ellen,0,0,elmu*area/ellen,0,0; 
    0,-12*elmu*iner/(ellen^3),-6*elmu*iner/(ellen^2),0,12*elmu*iner/(ellen^3),-6*elmu*iner/(ellen^2); 
    0,6*elmu*iner/(ellen^2),2*elmu*iner/ellen,0,-6*elmu*iner/(ellen^2),4*elmu*iner/ellen] ;
%进行整体刚度矩阵的组装 
k=zeros(freedom,freedom) ;
%生成单元定位向量 
dir(1,:)=1:6 ;
for i=2:elnum 
    dir(i,:)=dir(i-1,:)+3 ;
end 
%按照单元定位向量进行定位,实行分块定位的方法 
k(dir(1,:),dir(1,:))=ke ;
for i=2:elnum 
    k(dir(i,:),dir(i,:))=k(dir(i,:),dir(i,:))+ke ;
end 
%荷载为垂直方向的均布荷载，属于第二个自由度的方向 
p=zeros(freedom,1) ;
p(2:3:3*nonum)=-ones(nonum,1)*20*ellen 
%边界条件为两端固定，对应的节点为第一个节点和最后一个节点，约束的自由度为全部的自由度，采用划0置1法换算刚度矩阵 
for i=1:3 
    k(i,:)=0 
    k(:,i)=0 
end 
for i=1:3 
   k(i,i)=1 
end 
for i=freedom:-1:freedom-2 
    k(i,:)=0 
    k(:,i)=0 
end 
for i=freedom-2:1:freedom 
    k(i,i)=1 
end 
%对于节点力进行置0处理 
p(1:3)=0 
p(freedom-2:freedom)=0 
%进行求解 求出节点位移 
q=k\p 
%从节点位移中提取挠度，即纵向位移 
naodu=q(2:3:freedom) 
plot(0:ellen:length,naodu) 
grid