基于MATLAB 的平面刚架静力分析
为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应用MATLAB 编制计算程序对以平面刚架结构进行了静力分析。
本文还利用ANSYS 大型商用有限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果相当吻合，验证了计算结果的可靠性。
一、 问题描述
如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa ，截面为0.12×0.2m 的实心矩形，现要求解荷载作用下刚架的位移和内力。
5m
4m
3m
图1
二、矩阵位移法计算程序编制
为编制程序方便考虑，本文计算中采用“先处理法”。
具体的计算步骤如下。
（1） 对结构进行离散化，对结点和单元进行编号，建立结构（整体）坐标系
和单元（局部）坐标系，并对结点位移进行编号； （2） 对结点位移分量进行编码，形成单元定位向量e λ；
（3） 建立按结构整体编码顺序排列的结点位移列向量δ，计算固端力e F P 、等
效结点荷载E P 及综合结点荷载列向量P ；
（4） 计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵T 形成整体坐标
系下的单元刚度矩阵e T e K T K T = ； （5） 利用单元定位向量形成结构刚度矩阵K ； （6） 按式1=K P δ- 求解未知结点位移； （7） 计算各单元的杆端力e F 。
根据上述步骤编制了平面刚架的分析程序。
程序中单元刚度矩阵按下式计算。
32322
23
2
32
22
0000
1261260
064620
00001261260062640
EA
EA l
l EI EI
EI EI l l l l EI EI EI EI l l l l K EA EA l l EI EI EI EI l l l l EI EI EI EI l l
l l ⎡⎤-
⎢⎥⎢
⎥⎢⎥-
⎢⎥
⎢
⎥⎢⎥-
⎢⎥⎢
⎥=⎢⎥-⎢⎥
⎢
⎥⎢⎥---⎢⎥⎢
⎥⎢⎥-⎢⎥⎣
⎦
转换矩阵则按下式计算。
cos sin 0000sin cos 0000001000000cos sin 0000sin cos 00
1T ααααααα
α⎡⎤⎢⎥-⎢⎥⎢⎥=⎢
⎥⎢⎥⎢⎥-⎢⎥⎣⎦
计算程序框图如图2所示，具体的程序代码见附录1。
图2 MATLAB矩阵分析法程序框图
三、解题步骤
取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图4所示，局部坐标系用单元中箭头的方向表示，原始数据如下：
1
3
4
5
6
⑤
图3 图4 刚架结点输入矩阵为，
x=[0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2];
各单元定位向量为，
locvec1=[1 2 3 0 0 0];
locvec2=[1 2 3 4 5 6];
locvec3=[4 5 6 7 8 9];
locvec4=[7 8 9 10 11 12];
locvec5=[10 11 12 0 0 0];
输入截面参数，
E=2.1e11;%E=210GPa
a=0.12; %矩形截面长0.12m
b=0.2; %矩形截面宽0.2m
输入整体坐标系下各单元结点荷载列阵，
f(1,:)=zeros(1,6);
f(2,:)=[0 0 0 0 40e3 0];
f(3,:)=zeros(1,6);
f(4,:)=[0 0 0 -50e3 0 0];
f(5,:)=zeros(1,6);
输入整体坐标系下单元1等效节点荷载
q=10e3; %10kN/m
fe=[0.5*q*l(1),0,-q*l(1)^2/12,0.5*q*l(1),0,q*l(1)^2/12];
由此计算得到平面刚架整体坐标系下的结点位移(m),
d=
0.0035
0.0000
-0.0004
0.0030
-0.0005
-0.0004
0.0027
0.0000
0.0016
-0.0051
0.0000
-0.0006
各个单元的杆端力如表1所示，
表1 各单元杆端力
四、计算结果校核
在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。
计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。
图5 有限元模型图6 轴力图(kN)
图7 剪力图（kN）
图8 弯矩图（kN·m）
从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果比较，得到表2、3。
节点号项目矩阵位移法ANSYS 误差
1
Ux(m) 0 0 0 Uy(m) 0 0 0 Rotz(rad) 0 0 0
2 Ux(m) 0.003478 0.00348 -2E-06
Uy(m) 0.0000174 0.0000174 0
表3 各结点内力比较
由表2、表3的结果对比可知，两种方法的计算结果十分接近，误差均可以忽略不计，从而验证了计算结果的可靠性与准确性。
四、结论
通过对一个平面刚架静力分析问题的求解，本文编制的平面刚架静力分析程序计算结果与有限元软件ANSYS计算结果校核后，发现两者计算结果十分接近，误差可忽略不计，说明该程序计算结果的可靠性与精确性。
附录1 矩阵位移法计算程序
附录2 ANSYS建模计算命令流finish
/clear
/prep7
B=0.120$H=0.200$E=210000000
et,1,beam3
mp,ex,1,E
mp,prxy,1,0.3
r,1,B*H,B*H*H*H/12,H
k,1
k,2,0,5
k,3,1.6304,6.3681
k,4,4,5
k,5,4,1
k,6,4,-2
*set,i
*do,i,1,5
l,i,i+1
*enddo
lesize,all,0.5
latt,1,1,1
lmesh,all
dk,1,all
dk,6,all
fk,3,fy,-40
fk,5,fx,-50
lsel,s,,,1
esll,s
sfbeam,all,1,pres,10
allsel,all
dtran
ftran
sftran
/pbc,all,,2
/psf,pres,norm,2,0,1
eplot
/solu
solve
/post1
/pbc,u,,1 !显示支座约束符号，并图形显示变形pldisp,1 !将当前主要结果列表显示
presol,elem
!/pnum,sval,1
etable,mi,smisc,6
etable,mj,smisc,12
plls,mi,mj,-1 !弯矩图kN.m
etable,qi,smisc,2
etable,qj,smisc,8
plls,qi,qj,-1 !剪力图kN
etable,ni,smisc,1
etable,nj,smisc,7
plls,ni,nj,1 !轴力图kN。