這次程序自己編的一個(gè)簡單的門式剛架的程序，可以和上一個(gè)梁單元的例子對(duì)比學(xué)習(xí)。程序中所用的單元相當(dāng)于ansys中的beam3。根據(jù)旁邊的說明語句，可以很容易的讀懂這個(gè)程序。
%simple finite element program,beam element,分析門式剛架，采用梁單元，考慮軸向變形,相當(dāng)于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
%總長度為1ength，劃分的份數(shù)為elnum份，ellen是每個(gè)單元的長度
length=100
elnum=10
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; z
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] .
totalelnum=3*elnum
totalnonum=3*elnum+1
totalfreedom=3*totalnonum
%將其剛度矩陣由單元局部坐標(biāo)變化為整體坐標(biāo),變化的是左右兩邊的柱子
rr=pi/2
te=[cos(rr),sin(rr),0,0,0,0;
-sin(rr),cos(rr),0,0,0,0;
0,0,1,0,0,0;
0,0,0,cos(rr),sin(rr),0;
0,0,0,-sin(rr),cos(rr),0;
0,0,0,0,0,1]
kewhole=te'*ke*te
%進(jìn)行整體剛度矩陣的組裝
k=zeros(totalfreedom,totalfreedom)
%生成單元定位向量
dir(1,=1:6
for i=2:totalelnum
dir(i,=dir(i-1,+3
end
%按照單元定位向量進(jìn)行定位,實(shí)行分塊定位的方法,定位的次序?yàn)樽笾瑱M梁，右柱
%左柱
k(dir(1,,dir(1,)=kewhole
for i=2:elnum
k(dir(i,,dir(i,)=k(dir(i,,dir(i,)+kewhole
end
%橫梁
for i=elnum+1:elnum*2
k(dir(i,,dir(i,)=k(dir(i,,dir(i,)+ke
end
%右柱
for i=elnum*2+1:totalelnum
k(dir(i,,dir(i,)=k(dir(i,,dir(i,)+kewhole
end
%荷載為垂直方向的均布荷載，只加在橫梁上，在局部坐標(biāo)系中屬于第二個(gè)自由度的方向,均布荷載值為20，均布于橫梁上
p=zeros(totalfreedom,1)
p(3*elnum+2:3:6*elnum+2)=-ones(elnum+1,1)*20*ellen
%邊界條件為兩柱末端固定，對(duì)應(yīng)的節(jié)點(diǎn)為第一個(gè)節(jié)點(diǎn)和最后一個(gè)節(jié)點(diǎn)，約束的自由度為全部的自由度，采用劃0置1法換算剛度矩陣
for i=1:3
k(i,=0
k(:,i)=0
end
for i=1:3
k(i,i)=1
end
for i=totalfreedom:-1:totalfreedom-2
k(i,=0
k(:,i)=0
end
for i=totalfreedom-2:1:totalfreedom
k(i,i)=1
end
%對(duì)于節(jié)點(diǎn)力進(jìn)行置0處理
p(1:3)=0
p(totalfreedom-2:totalfreedom)=0
%進(jìn)行求解 求出節(jié)點(diǎn)位移
q=k\p
%從節(jié)點(diǎn)位移中提取橫梁撓度，即縱向位移
naodu=q(3*elnum+2:3:6*elnum+2)
plot(0:ellen:length,naodu)
grid
1 I! _( G% R3 w0 T9 v3 M6 g

計(jì)算結(jié)果的撓度曲線

平面梁單元的說明

局部坐標(biāo)和整體坐標(biāo)變換的圖示