基于MATLAB的平面刚架静力分析.doc

基于MATLAB的平面刚架静力分析 为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应用MATLAB编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS大型商用有限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果相当吻合，验证了计算结果的可靠性。一、 问题描述如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa，截面为0.120.2m的实心矩形，现要求解荷载作用下刚架的位移和内力。图1二、矩阵位移法计算程序编制为编制程序方便考虑，本文计算中采用“先处理法”。具体的计算步骤如下。（1） 对结构进行离散化，对结点和单元进行编号，建立结构（整体）坐标系和单元（局部）坐标系，并对结点位移进行编号；（2） 对结点位移分量进行编码，形成单元定位向量；（3） 建立按结构整体编码顺序排列的结点位移列向量，计算固端力、等效结点荷载 及综合结点荷载列向量；（4） 计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵形成整体坐标系下的单元刚度矩阵 ；（5） 利用单元定位向量形成结构刚度矩阵；（6） 按式 求解未知结点位移；（7） 计算各单元的杆端力。根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计算。转换矩阵则按下式计算。计算程序框图如图2所示，具体的程序代码见附录1。图2 MATLAB矩阵分析法程序框图三、解题步骤取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图4所示，局部坐标系用单元中箭头的方向表示，原始数据如下： 图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=210GPaa=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/mfe=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 各单元杆端力单元12345i端Fx(kN)-17917.04717917.0517917.0517917.05-32083Fy(kN)17507.3731-17507.422492.6322492.6322492.63Mz(kNm)1897.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(kNm)-37312.596-2092.8326668.34-44999.851249.01四、计算结果校核在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。图5 有限元模型图6 轴力图(kN)图7 剪力图（kN）图8 弯矩图（kNm）从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果比较，得到表2、3。表2 各节点位移比较节点号项目矩阵位移法ANSYS误差1Ux(m)000Uy(m)000Rotz(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.00001340.00001340Rotz(rad)-0.00056-0.000557-3E-066Ux(m)000Uy(m)000Rotz(rad)000表3 各结点内力比较节点号项目矩阵位移法ANSYS误差1Fx(kN)-32.083-32.080-0.003Fy(kN)-17.507-17.503-0.004Mz(kNm)-37.313-37.307-0.0062Fx(kN)-17.917-17.9200.003Fy(kN)17.50717.5030.004Mz(kNm)1.8981.909-0.0113Fx(kN)-17.917-17.9200.003Fy(kN)-22.493-22.4970.004Mz(kNm)-2.093-2.1100.0184Fx(kN)-17.917-17.9200.003Fy(kN)-22.493-22.4970.004Mz(kNm)26.66826.682-0.0145Fx(kN)-32.083-32.080-0.003Fy(kN)-22.493-22.4970.004Mz(kNm)-45.000-44.999-0.0016Fx(kN)32.08332.0800.003Fy(kN)-22.493-22.4970.004Mz(kNm)51.24951.2390.010由表2、表3的结果对比可知，两种方法的计算结果十分接近，误差均可以忽略不计，从而验证了计算结果的可靠性与准确性。四、结论通过对一个平面刚架静力分析问题的求解，本文编制的平面刚架静力分析程序计算结果与有限元软件ANSYS计算结果校核后，发现两者计算结果十分接近，误差可忽略不计，说明该程序计算结果的可靠性与精确性。附录1 矩阵位移法计算程序pmgj.m 计算主程序clear clc%结点坐标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;%单元定位向量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=210GPaa=0.12; %矩形截面长0.12mb=0.2; %矩形截面宽0.2mA=a*b;I=(a*b3)/12;clear a b% 单元长度计算for i=1:5 dx=x(i,1)-x(i+1,1); dy=x(i,2)-x(i+1,2); l(i)=(dx2+dy2)0.5;end%生成转换矩阵t1=coortrans(x(2,1),x(2,2),x(1,1),x(1,2);t2=coortrans(x(2,1),x(2,2),x(3,1),x(3,2);t3=coortrans(x(3,1),x(3,2),x(4,1),x(4,2);t4=coortrans(x(4,1),x(4,2),x(5,1),x(5,2);t5=coortrans(x(5,1),x(5,2),x(6,1),x(6,2);%结点荷载（结构坐标系下）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);%单元等效节点荷载（结构坐标系下）q=10e3; %10kN/mfe=0.5*q*l(1),0,-q*l(1)2/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;%生成单元刚度矩阵(单元坐标系）k1=elemstiffm(E,A,l(1),I);k2=elemstiffm(E,A,l(2),I);k3=elemstiffm(E,A,l(3),I);k4=elemstiffm(E,A,l(4),I);k5=elemstiffm(E,A,l(5),I);%将单元坐标系下的单元刚度矩阵转换到结构坐标系下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);NonLog=(locvec1=0);ij=find(NonLog);IJ=locvec1(NonLog);K(IJ,IJ)=K(IJ,IJ)+K1(ij,ij);F(IJ)=F(IJ)+f(1,ij)+fe(ij);clear NonLog ij IJNonLog=(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 NonLog ij IJNonLog=(locvec3=0);ij=find(NonLog);IJ=locvec3(NonLog);K(IJ,IJ)=K(IJ,IJ)+K3(ij,ij);F(IJ)=F(IJ)+f(3,ij);clear NonLog ij IJNonLog=(locvec4=0);ij=find(NonLog);IJ=locvec4(NonLog);K(IJ,IJ)=K(IJ,IJ)+K4(ij,ij);F(IJ)=F(IJ)+f(4,ij);clear NonLog ij IJNonLog=(locvec5=0);ij=find(NonLog);IJ=locvec5(NonLog);K(IJ,IJ)=K(IJ,IJ)+K5(ij,ij);F(IJ)=F(IJ)+f(5,ij);clear NonLog ij IJ%节点位移d=KF;%单元1杆端力（结构坐标系下）d1=zeros(6,1);NonLog=(locvec1=0);ij=find(NonLog);ij1=find(NonLog);IJ=locvec1(NonLog);d1=d(IJ);d1(ij1)=0;F1=K1*d1-fe;%单元2杆端力（结构坐标系下）d2=zeros(6,1);NonLog=(locvec2=0);ij=find(NonLog);ij1=find(NonLog);IJ=locvec2(NonLog);d2=d(IJ);d2(ij1)=0;F2=K2*d2-f(2,:);%单元3杆端力（结构坐标系下）d3=zeros(6,1);NonLog=(locvec3=0);ij=find(NonLog);ij1=find(NonLog);IJ=locvec3(NonLog);d3=d(IJ);d3(ij1)=0;F3=K3*d3-f(3,:);%单元4杆端力（结构坐标系下）d4=zeros(6,1);NonLog=(locvec4=0);ij=find(NonLog);ij1=find(NonLog);IJ=locvec4(NonLog);d4=d(IJ);d4(ij1)=0;F4=K4*d4-f(4,:);%单元5杆端力（结构坐标系下）d5=zeros(6,1);NonLog=(locvec5=0);ij=find(NonLog);ij1=find(NonLog);IJ=locvec5(NonLog);d5=d(IJ);d5(ij1)=0;F5=K5*d5-f(5,:);coortrans.m 转换矩阵生成函数% 转换矩阵生成函数（单元坐标系-结构坐标系）%y=coortrans(x1,y1,x2,y2),x1,y1,x2,y2为单元两端结点在结构坐标系下的坐标，单位都为mfunction y=coortrans(x1,y1,x2,y2)a=atan(y2-y1)/(x2-x1);c=cos(a);s=sin(a);t=zeros(6);t(1,1)=c;t(1,2)=s;t(2,1)=-s;t(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;endelemstiffm.m 单元刚度矩阵生成函数% 单元刚度矩阵生成函数% y=elemstiffm(E,A,l,I),E,A,l,I均采用国际单位制 Pa m2 m m4function y=elemstiffm(E,A,l,I)i1=E*A/l;i2=12*E*I/(l3);i3=6*E*I/(l2);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);k(4,1)=k(1,4);k(5,2)=k(2,5);k(5,3)=k(3,5);k(6,2)