根据MATLAB的平面刚架静力分析

1、基于MATLAB勺平面刚架静力分析本文还利用ANSYS为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应 用MATLAB编制计算程序对以平面刚架结构进行了静力分析。发现两者计算结果大型商用有限元分析软件对矩阵位移法的计算结果进行校核, 相当吻合，验证了计算结果的可靠性。问题描述如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa截面为0.12X 0.2m的实心矩形，现要求解荷载作用下刚架的位移和内力。5ma77图1:4m3m二、矩阵位移法计算程序编制为编制程序方便考虑，本文计算中采用“先处理法”。具体的计算步骤如下。(1) 对结构进行离散化，对结点和单元进行编号，建

2、立结构(整体)坐标系和单元(局部)坐标系，并对结点位移进行编号;(2) 对结点位移分量进行编码，形成单元定位向量(3) 建立按结构整体编码顺序排列的结点位移列向量，计算固端力pe、等效结点荷载Pe及综合结点荷载列向量P ；(4)计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵T形成整体坐标系下的单元刚度矩阵Ke t'K ;0(5)利用单元定位向量形成结构刚度矩阵(6)按式二K 1P求解未知结点位移;(7)计算各单元的杆端力Fe 0根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计EAlEAl0K EAl12EI6EI6EIT4EIlEAl12EI6EIT6EIV2EIl0

3、012EI6EIl3l26EI2EIl2l0012EI6EIl3l26EI4EIl2l转换矩阵则按下式计算。cossinsin 0cos 00000000010 cos0sin0 00 sin00cos 00 1计算程序框图如图2所示，具体的程序代码见附录1。Start1End图2 MATLAB矩阵分析法程序框图三、解题步骤取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图4所示，局部坐标系用单兀中箭头的方向表示，原始数据如下：456刚架结点输入矩阵为,x=0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2;各单元定位向量为,Iocvec1=1 2 3 0 0

4、0;Iocvec2=1 2 3 4 5 6;Iocvec3=4 5 6 7 8 9;Iocvec4=7 8 9 10 11 12; locvec5=10 11 12 0 0 0;输入截面参数,E=2.1e11;%E=210G Paa=0.12; %矩形截面长0.12mb=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

5、/m fe=0.5*q*l(1),0,-q*l(1)A2/12,0.5*q*l(1),0,q*l(1)A2/12;由此计算得到平面刚架整体坐标系下的结点位移(m),d=0.00350.0000-0.00040.0030-0.0005-0.00040.00270.00000.0016-0.00510.0000-0.0006各个单元的杆端力如表1所示,表1各单兀杆端力单元12345 m i端Fx(kN)-17917.04717917.0517917.0517917.05-32083Fy(kN)17507.3731-17507.422492.6322492.6322492.63Mz(kN - m)1

6、897.83076-1897.832092.833-26668.344999.85j端Fx(kN)-32082.953-17917-1791732082.9532082.95Fy(kN)-17507.373-22492.6-22492.6-22492.6-22492.6Mz(kN - m)-37312.596-2092.8326668.34-44999.851249.01四、计算结果校核在ANSY鋪使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力图、剪力图和弯矩图如图 6 7、8所示，命令流见附录2。t i*tL&

7、#39;I'HUI1C图5有限元模型IIILianrAN STSILE -iin“ELey-?MAX -IT 口NTFLPW-7IP-M .Tiii-it .tJE'7-JI. iMT图6轴力图（kN）图7剪力图（kN）图8弯矩图（kN - m）巧剪丄:di“Mgr!-M .昭t-It. LSk."-FHT4ja LhM卯 ptLJ.-1TUt-Jrii屮Hjzr =-5jW =*i.+*从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果比较，得到表2、3。表2各节点位移比较节点号项目矩阵位移法ANSYS误差1Ux(m)000Uy(m)000Rotz

8、(rad)0002Ux(m)0.0034780.00348-2E-06Uy(m)0.00001740.00001740Rotz(rad)-0.00037-0.000365-5E-063Ux(m)0.0030180.00302-2E-06Uy(m)-0.00051-0.0005122E-06Rotz(rad)-0.00038-0.000378-2E-064Ux(m)0.0026870.00269-3E-06Uy(m)0.00003120.00003120Rotz(rad)0.0016240.001624E-065Ux(m)-0.00513-0.0051451.5E-05Uy(m)0.000013

9、40.00001340Rotz(rad)-0.00056-0.000557-3E-066Ux(m)000Uy(m)000Rotz(rad)000表3各结点内力比较节点号项目矩阵位移法ANS YS误差1Fx(kN)-32.083-32.080-0.003Fy(kN)-17.507-17.503-0.004Mz(kN m)-37.313-37.307-0.0062Fx(kN)-17.917-17.9200.003Fy(kN)17.50717.5030.004Mz(kN m)1.8981.909-0.0113Fx(kN)-17.917-17.9200.003Fy(kN)-22.493-22.4970

10、.004Mz(kN m)-2.093-2.1100.0184Fx(kN)-17.917-17.9200.003Fy(kN)-22.493-22.4970.004Mz(kN m)26.66826.682-0.0145Fx(kN)-32.083-32.080-0.003Fy(kN)-22.493-22.4970.004Mz(kN m)-45.000-44.999-0.0016Fx(kN)32.08332.0800.003Fy(kN)-22.493-22.4970.004Mz(kN m)51.24951.2390.010由表2、表3的结果对比可知，两种方法的计算结果十分接近，误差均可以忽略不计，从而

11、验证了计算结果的可靠性与准确性。四、结论通过对一个平面刚架静力分析问题的求解， 本文编制的平面刚架静力分析程 序计算结果与有限元软件 ANSYS+算结果校核后，发现两者计算结果十分接近，误差可忽略不计，说明该程序计算结果的可靠性与精确性。附录1矩阵位移法计算程序pmgj.m计算主程序clearcIc%结点坐标x=0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2;% x1=0;y1=0;% x2=0;y2=-5;% x3=1.63;y3=-6.37;% x4=4;y4=-5;% x5=4;y5=-1;% x6=4;y6=2;%单元定位向量Iocvec1=1 2 3 0 0 0;

12、Iocvec2=1 2 3 4 5 6; locvec3=4 5 6 7 8 9;Iocvec4=7 8 9 10 11 12; Iocvec5=10 11 12 0 0 0;%刚架的材料特性截面特性E=2.1e11;%E=210G Pa a=0.12; %矩形截面长0.12m b=0.2; %矩形截面宽0.2m A=a*b;1=*"3)/12;clear a b%单元长度计算for i=1:5dx=x(i,1)-x(i+1,1);dy=x(i,2)-x(i+1,2);I(i)=(dx2+dy2)0.5;end%生成转换矩阵t1=coortra ns(x(2,1),x(2,2),x(

13、1,1),x(1,2); t2=coortra ns(x(2,1),x(2,2),x(3,1),x(3,2);t3=coortra ns(x(3,1),x(3,2),x(4,1),x(4,2);t4=coortra ns(x(4,1),x(4,2),x(5,1),x(5,2); t5=coortra ns(x(5,1),x(5,2),x(6,1),x(6,2);zeros(1,6);0 0 0 0 40e3 0;zeros(1,6);0 0 0 -50e3 0 0;zeros(1,6);%结点荷载(结构坐标系下) f(1,:) f(2,:) f(3,:) f(4,:) f(5,:)%单元等效节

14、点荷载(结构坐标系下)q=10e3;%10kN/mfe=0.5*q*l(1),0,-q*l(1)A2 /12,0.5*q*l(1),0,q*l(1)2 /12;%单元坐标系下的值fe1=0,0.5*q*l(1),q*l(1)2 /12,0,0.5*q*l(1),-q*l(1)2 /12;%生成单元刚度矩阵(单元坐标系) k仁elemstiffm(E,A,l(1),l); k2=elemstiffm(E,A,l(2),l); k3=elemstiffm(E,A,l(3),l); k4=elemstiffm(E,A,l(4),I); k5=elemstiffm(E,A,l(5),I);%将单元坐标

15、系下的单元刚度矩阵转换到结构坐标系下K1=t1'*k1*t1;K2=t2'*k2*t2;K3=t3'*k3*t3;K4=t4'*k4*t4;K5=t5'*k5*t5;%组装总体刚度矩阵K=zeros(12);F=zeros(12,1);Non Log=(locvec1=0); ij=fi nd(N on Log); IJ=locvec1(No nLog); K(IJ,IJ)=K(IJ,IJ)+K1(ij,ij); F(IJ)=F(IJ)+f(1,ij)'+fe(ij)'clear Non Log ij IJ Non Log=(locvec

16、2=0); ij=fi nd(N on Log); IJ=locvec2(No nLog); K(IJ,IJ)=K(IJ,IJ)+K2(ij,ij); F(IJ)=F(IJ)+f(2,ij)'clear Non Log ij IJNon Log=(locvec3=0); ij=fi nd(N on Log); IJ=locvec3(No nLog); K(IJ,IJ)=K(IJ,IJ)+K3(ij,ij); F(IJ)=F(IJ)+f(3,ij)'clear Non Log ij IJNon Log=(locvec4=0); ij=fi nd(N on Log); IJ=loc

17、vec4( Non Log); K(IJ,IJ)=K(IJ,IJ)+K4(ij,ij); F(IJ)=F(IJ)+f(4,ij)' clear Non Log ij IJNon Log=(locvec5=0); ij=fi nd(N on Log); IJ=locvec5(N on Log); K(IJ,IJ)=K(IJ,IJ)+K5(ij,ij); F(IJ)=F(IJ)+f(5,ij)' clear Non Log ij IJ %节点位移 d=KF;%单元1杆端力(结构坐标系下) d1=zeros(6,1);Non Log=(locvec1=0); ij=fi nd(N o

18、n Log);ij1=fi nd(No nLog);IJ=locvec1(No nLog); d1=d(IJ);d1(ij1)=0;F1=K1*d1-fe'%单元2杆端力(结构坐标系下)d2=zeros(6,1):No nLog=(locvec2=0); ij=fi nd(N on Log); ij1=fi nd(No nLog); IJ=locvec2(No nLog); d2=d(IJ);d2(ij1)=0; F2=K2*d2-f(2,:)'%单元3杆端力(结构坐标系下)d3=zeros(6,1);Non Log=(locvec3=0);ij=fi nd(N on Log)

19、;ij1=fi nd(No nLog);IJ=locvec3(No nLog); d3=d(IJ);d3(ij1)=0;F3=K3*d3-f(3,:)'%单元4杆端力(结构坐标系下)d4=zeros(6,1);Non Log=(locvec4=0); ij=fi nd(N on Log);ij1=fi nd(No nLog);IJ=locvec4( Non Log); d4=d(IJ);d4(ij1)=0;F4=K4*d4-f(4,:)'%单元5杆端力(结构坐标系下)d5=zeros(6,1);Non Log=(locvec5=0);ij=fi nd(N on Log);ij1

20、=fi nd(No nLog);IJ=locvec5(N on Log); d5=d(IJ);d5(ij1)=0;F5=K5*d5-f(5,:)'coortra ns.m转换矩阵生成函数%转换矩阵生成函数(单元坐标系-结构坐标系)%y=coortra ns(x1,y1,x2,y2),x1,y1,x2,y2为单元两端结点在结构坐标系下的坐标，单位都为mfun cti on y=coortra ns(x1,y1,x2,y2)a=ata n(y2-y1)/(x2-x1);c=cos(a);s=s in( a);t=zeros(6);t(1,1)=c;t(1,2)=s;t(2,1)=-s;t(

21、2,2)=c;t(3,3)=1;t(4,4)=c;t(4,5)=s;t(5,4)=-s;t(5,5)=c;t(6,6)=1;y=t; end单元刚度矩阵生成函数elemstiffm.m%单元刚度矩阵生成函数% y=elemstiffm(EAI,l),E,A,l,均采用国际单位制Pa m2 m m4fun cti on y=elemstiffm(E,A,l,l)i1=E*A/l;i2=12*E*l/(|A3); i3=6*E*l/(|A2);i4=4*E*I/l;i5=2*E*I/l;k=zeros(6);k(1,1)=i1;k(1,4)=-i1; 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)=i1;k(5,5)=i2;k(5,6)=-i3;k(6,6)=i4;k(3,2)=k(2,3