基于MATLAB 的平面刚架静力分析
为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应 用
MATLAB编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS 大型商用有
限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果 相当吻合，验证了计算结果的可靠性。
一、问题描述
如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa,截面为0.22 X0.2m的实心矩形，现要求解荷载作用下刚架的位移和内力。
5m
4m
图1
二、矩阵位移法计算程序编制
为编制程序方便考虑，本文计算中采用“先处理法”。具体的计算步骤如下。
（1） 对结构进行离散化，对结点和单元进行编号，建立结构（整体）坐标系 和单元（局部）坐标系，并对结点位移进行编号；
（2） 对结点位移分量进行编码，形成单元定位向量e;
（3）建立按结构整体编码顺序排列的结点位移列向量，计算固端力耳、等 F
效结点荷载PE及综合结点荷载列向量P；
(4)
计算个单元局部坐标系的刚度矩阵, 通过坐标变换矩阵T形成整体坐标
系下的单元刚度矩阵Ke TTKeT ;
（5） 利用单元定位向量形成结构刚度矩阵K； （6） 按式二K iP求解未知结点位移；
（7）计算各单元的杆端力Fe
根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计
•EA
I
0
0 ■ — ■
1
0
0
oti
0 \"1
IZtl
orzi h I2
0 I3 I2
oEi
州
tzl oEi
zEi
0 0
k
1
I2
1
-EA
rA
I
0 0
0 0
1
■
5ET «■■■» •
12 Ei
0
0
IZ tl
k I2
I3 I2
Otzl■ ■0
Ztzl■ •
Otzl■ ■
4tZl I2
1
0
I2
1
算。
转换矩阵则按下式计算。
cos sin 0
T
sin cos 1 cos 0 0
0
0 0
0 0 00
0 00 0
0 0 0 0
0 0 0
八
sin 0 0
sin 0
cos 0
0 1
计算程序框图如图2所示， 具体的程序代码见附录1。
End
2 MATLAB矩阵分析法程序框图图
三、解题步骤
取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图
4所示，局部坐标系用单元中箭头的方向表示,
原始数据如下:
图3
刚架结点输入矩阵为,
x=[0 0;0 -5;1.63 -6・37;4 -5;4 -1;4 2];
各单元定位向量为，
locvecl=[1 locvec2=[1 4 5 6]; locvec3=[4
locvec4=[7 10 11 12]; locvec5=[10 11
12 0 0 0];
输入截面参数,
E=2 ・1ell;%E=210GPa a=0.12; %矩形截面长0.12m b=0.2;令矩形截面宽0・2m
输入整体坐标系下各单元结点荷载列阵,
f (1, :)=zeros(lz 6); f (2, :) = [0 0 0 0 40e3 0]; f (3, :)=zeros(lz 6); f (4, :) = [0 0 0 -50e3 0 0]; f (5,
:)=zeros(lz 6); 输入整体坐标系下单元1等效节点荷载 q=10e3;
%10kN/m
fe=[O・5*q*l(1),Oz-q*l(1)2/12,0.5*q*l(1),OAq*l(1)2/12]; 由此计算得到平面刚架整体坐
AA
标系下的结点位移(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各单元杆端力
单元 Fx (kN) ・ A11I 1 -17917.047 17507.3731 1897.83076 -32082.953 -17507.373 -37312.596 2 17917.05 -17507.4 -1897.83 -17917 -22492.6 -2092.83 3 17917.05 22492.63 2092.833 -17917 -22492.6 26668.34 4 17917.05 22492.63 -26668.3 32082.95 -22492.6 一5 -32083 22492.63 44999.85 32082.95 -22492.6 51249.01 丄瑞 Fv(kN) Mz(kN • m) Fx (kN) ・亠山 肿而 Fy(kN) Mz(kN • m)
44999.8 四、计算结果校核
在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模 型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力 图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。
图5有限元模型
LIIJE 空 17F巳33
AN5YS
=a TIWEaa HD 时
USM =-2&・7^83
EL£M€=21
WAZ =-1 丁. -SOST Ei£^=a II
-16 ..9m
-J4..?8Si
^..650
-20J5942
-12-«$342
一药 J 38 幺
一23 “67卩
TL・・
一1 ^..5617
-17..5027
图6轴力图(kN)
I/IMZ 3TPE.33
AN$Y$
KH5
•SUB
TH厢勺 Cl 3
HDM =-32・*7亨坊 BLffMPi
MAX =32・Q7S6 EI£W=32 if
-15 “07*
-2「餡g
-17..的2
-10”务笛2
-3..564 34-
3M5W?>
i&..«»U
17.^2
乳.•第处
MnO7J«
图7剪力图(kN)
LBPS 3VPE3S
SUB »J •rnwe-j IMB MJ
ANSYS
KIM
wnu =-51 ・：2：3E2 EU?^=3S
宓工=Z焙5
ELFA*=3O
TO/轻
-4.4CS4S
3.r32«?3
U..3SS1
2^.4131
图8弯矩图(kN・m)
从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果 比较，得到表2、3o
表2各节点位移比较
节点号 1 项目 Ux (m) Uy(m) Rotz(rad) Ux(m) Uy (m) 矩阵位移法 ANSYS 误差 0 0 0 0. 003478 0.0000174 0 0 0 0. 00348 0.0000174 0 0 0 -2E-06 0 2
Rotz (rad) -0. 00037 -0. 000365 -5E-06
Ux (m) 3 Uy(m) Rotz (rad) Ux(m) 4 Uy (m) Rotz (rad) Ux (m) 3 Uy (m) Rotz(rad) 0. 003018 -0. 00051 -0. 00038 0. 002687 0.0000312 0.001624 -0. 00513 0.0000134 -0. 00056 0 0 0 0. 00302 -0. 000512 -0. 000378 0. 00269 0.0000312 0. 00162 -0. 005145 0.0000134 -0. 000557 0 0 0 -2E-06 2E-06 -2E-06 -3E-06 0 4E-06 1. 5E-05 0 -3E-06 0 0 0 Ux(m) 6
Uy(m) Rotz (rad) 表3各结点内力比较
节点号 1 项冃 Fx(kN) Fy(kN) Mz (kN • m) Fx (kN) 矩阵位移法 -32. 083 -17. 507 -37. 313 -17.917 17. 507 1.898 -17.917 -22. 493 -2. 093 -17.917 -22. 493 26. 668 -32. 083 -22. 493 -45. 000 32. 083 -22. 493 51.249 ANSYS -32. 080 -17. 503 -37. 307 -17. 920 17. 503 1.909 -17. 920 -22. 497 -2. 110 -17. 920 -22. 497 26. 682 -32. 080 -22. 497 -44. 999 32. 080 -22. 497 51.239 误差 -0. 003 -0. 004 -0. 006 0. 003 0. 004 -0.011 0. 003 0. 004 0.018 0. 003 0. 004 -0. 014 -0. 003 0. 004 -0. 001 0. 003 0. 004 0.010 2 Fy(kN) Mz(kN • m) Fx(kN) 3 Fy(kN) Mz (kN • m) Fx (kN) Fy(kN) Mz(k\\ • m) Fx(kN) Fy(kN) Mz (kN • m) Fx (kN) 4 5 6 Fy(kN) Mz(kN • m) 由表2、表3的结果对比可知，两种方法的计算结果十分接近，误差均可以 忽略不计，从而验证了计算结果的可靠性与准确性。
四、结论
通过对一个平面刚架静力分析问题的求解，本文编制的平面刚架静力分析程 序计算结果与有限元软件ANSYS计算结果校核后，发现两者计算结果十分接近， 误差可忽
略不计，说明该程序计算结果的可靠性与精确性。
附录1矩阵位移法计算程序
pmgj.m 计算主程序 clear clc
%结点坐标
x=[0 0;0 -5;1.63 ・6・37;4 -5;4 -1;4 2]; % xl=O;yl=O; % x2=0;y2=-5; %x3=1.63;y3=-6.37; % x4=4;y4=-5; % x5=4;y5=-l; % x6=4;y6=2;
%单元定位向量 locvecl=[1 2 3 000]; Iocvec2=[l 2 345 6];
Iocvec3=[4 5 6 7 8 9]; Iocvec4=[7 8 9 10 11 12]; Iocvec5=[10 11 12 0 0 0];
%刚架的材料特性截面特性
E=2.1ell;%E=210GPa a=0.12; %矩形截面长0.12m b=0.2; %矩形截面宽0.2m A=a*b; l=(a*bA3)/12; clear a b
%单元长度计算
for i=l:5
dx=x(i,l)-x(i+l,l);
dy=x(i,2)-x(i+l,2)； l(i)=(dxA2+dyA2)A0.5; end
%生成转换矩阵
tl=coortrans(x(乙 l),x(2,2),x(l,l),x(l,2)); t2=coortrans(x(2,l),x(2/2),x(3,l),x(3,2)); t3=coortra
ns(x(3J),x(3,2),x(4JLx(4,2));
t4=coortrans(x(4,l),x(4z2),x(5,l)/X(5,2)); t5=coortrans(x(5,l)/x(5/2),x(6,l),x(6,2));
%结点荷载(结构坐标系下)
f(l,:)=zeros(l?6);
f(2,:)=[0 0 0 040e3 0]; f(3/)=zeros(l,6);
f(4/)=[0 0 0 -50e3 0 0]; f(5/)=zeros(l,6);
%单元等效节点荷载(结构坐标系下)
q=10e3;
%10kN/m
fe=[0.5*q*l(l),0,-q*l(l)A2/12,0.5*q*l(lL0,q*l(l)A2/12]; %
单元坐标系下的值
fel=[0/0.5*q*l(l),q*l(l)A2/12/0,0.5*q*l(l),-q*l(l)A2/12];
%生成单元刚度矩阵(单元坐标系)
kl=elemstiffm(E,AJ(l)/l)； k2=elemstiffm(E,AJ(2)J); k3=elemstiffm(E/A,l(3)J); k4=elemstiffm(E,A/l(4)J); k5=elemstiffm(E,AJ(5)J);
%将单元坐标系下的单元刚度矩阵转换到结构坐标系下
Kl=tl'*kl*tl; K2=t2**k2*t2; K3=t3'*k3*t3; K4=t4'*k4*t4; K5=t5,*k5*t5;
%组装总体刚度矩阵
K=zeros(12); F=zeros(12,l); NonLog=(locvecl~=0); ij=find(NonLog); IJ=locvecl(NonLog); K(IJ/IJ)=K(IJ,IJ)+Kl(ij/ij); F(IJ)=F(IJ)+f(l/ij)'+fe(ij)1; clear Non Log ij IJ
NonLog=(locvec2~=0); ij=find(NonLog); IJ=locvec2(NonLog); K(IJ/IJ)=K(IJ,IJ)+K2(ij/ij); F(IJ)=F(IJ)+f(2/ij)-; clear Non Log ij IJ
NonLog=(locvec3~=0); ij=fi nd(N on Log);
IJ=locvec3(NonLog); K(IJ/IJ)=K(IJ,IJ)+K3(ij/ij); F(IJ)=F(IJ)+f(3/ij)1; clear Non Log ij IJ
NonLo 沪(Iocvec4~=0); ij=fi nd(N on Log); IJ=locvec4(NonLog); K(IJ,IJ)=K(IJ,IJ)+K4(ijzij); F(IJ)=F(IJ)+f(4Jj)'; clear Non Log ij IJ NonLog=(locvec5~=0); ij=fi nd(NonLog); IJ=locvec5(NonLog); K(IJ/IJ)=K(IJ,IJ)+K5(ij/ij); F(IJ)=F(IJ)+f(5/ij)1; clear Non Log ij IJ
%节点位移
d=K\\F;
%单元1杆端力(结构坐标系下) dl=zeros(6,l);
NonLog=(locvecl~=0); ij=fi nd(N on Log); ijl=find(~NonLog); IJ=locvecl(NonLog); dl=d(IJ); dl(ijl)=O; Fl=Kl*dl-fe,;
%单元2杆端力(结构坐标系下) d2=zeros(6,l);
NonLog=(locvec2~=0); ij=find(NonLog); ijl=find(~NonLog); IJ=locvec2(NonLog); d2=d(IJ); d2(ijl)=0; F2=K2*d2-f(2/)';
%单元3杆端力(结构坐标系下) d3=zeros(6,l)；
NonLog=(locvec3~=0); ij=fi nd(N on Log); ijl=find(~NonLog); IJ=locvec3(NonLog); d3=d(IJ); d3(ijl)=O; F3=K3*d3-f(3/:)1;
%单元4杆端力(结构坐标系下) d4=zeros(6,l);
NonLog=(locvec4~=0); ij=fi nd(N on Log); ijl=find(~NonLog); IJ=locvec4(NonLog); d4=d(IJ); d4(ijl)=0; F4=K4*d4-f(4,:)1;
%单元5杆端力(结构坐标系下) d5=zeros(6,l);
NonLog=(locvec5~=0); ij=find(NonLog); ijl=find(~NonLog); IJ=locvec5(NonLog); d5=d(IJ); d5(ijl)=O; F5=K5*d5-f(5,:)1;
coortrans.m 转换矩阵生成函数 %转换矩阵生成函数(单元坐标系・＞结构坐标
系)%y=coortrans(xl/yl,x2,y2),xl,yl/x2/y2为单元两端结点在结构坐标系下的坐 标，单位都为
m
function y二coortrans(xbyl.,x2,y2) a=atan((y2-yl)/(x2-xl)); c=cos(a); s二sin(a); t=zeros(6); t(l,l)=c;t(l,2)=s; t(2,l)=-s;t(2,2)=c; t(3,3)=l;
t(4,4)=c;t(4z5)=s; t(5/4)=-s;t(5,5)=c; t(6,6)=l; y=t； end
单元冈ij度矩阵生成函数 %单元刚度矩阵生成函数~~I
% y=elemstiffm(E/A/IJ),E,AJ,l 均采用国际单位制 Pa m2 m m4 function y二elemstiffm(E,A」」)
elemstiffm.m il=E*A/l; i2=12*E*l/(lA3); i3=6*E*l/(lA2); i4=4*E*l/l; i5=2*E*l/l; k=zeros(6); k(lj)=il;k ⑴ 4)“2;
k(2,2)=i2;k(2,3)=i3;k(2,5)=-i2;k(2,6)=i3; k(3,3)=i4;k(3,5)=-i3;k(3,6)=i5; k(4,4)=il;
k(5,5)=i2;k(5,6)=-i3;
k(6,6)=i4; k(3,2)=k(2,3); k(4zl)=k(l,4);
k(5,2)=k(2,5);k(5,3)=k(3,5); k(6/2)=k(2,6);k(6,3)=k(3,6);k(6,5)=k(5,6); y 二 k; end
附录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
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 IJJ+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,1O allsel,all
文案
dtran ftran sftran
/pbc,all,,2 /psf,pres, norm20J eplot /solu solve /postl /pbc,u„1
!显示支座约束符号，并图形显示变形
pldisp,1 !将当前主要结果列表显示presol,elem !/pnum,sval,1 etable,mi,smisc,6 etable 耐 smisc,12
plls5mi,mj,-1 !弯矩图 kN.m etable,qi,smisc,2 etable,qj,smisc,8 plls5qi5qj,-1 !剪力图 kN etable,ni,smisc,1 etable,nj,smisc,7 plls,ni,nj,1 !轴力图 kN
文案
因篇幅问题不能全部显示，请点此查看更多更全内容