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

工程计算实践作业1基于MATLAB的平面刚架静力分析为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应用MATLAB编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS大型商用有限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果相当吻合，验证了计算结果的可靠性。一、问题描述如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa，截面为0.12×0.2m的实心矩形，现要求解荷载作用下刚架的位移和内力。5m4m3m4m10kN/m40kN50kN30°40°图1二、矩阵位移法计算程序编制为编制程序方便考虑，本文计算中采用“先处理法”。具体的计算步骤如下。工程计算实践作业2（1）对结构进行离散化，对结点和单元进行编号，建立结构（整体）坐标系和单元（局部）坐标系，并对结点位移进行编号；（2）对结点位移分量进行编码，形成单元定位向量e；（3）建立按结构整体编码顺序排列的结点位移列向量，计算固端力eFP、等效结点荷载EP及综合结点荷载列向量P；（4）计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵T形成整体坐标系下的单元刚度矩阵eTeKTKT；（5）利用单元定位向量形成结构刚度矩阵K；（6）按式1=KP求解未知结点位移；（7）计算各单元的杆端力eF。根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计算。323222323222000012612600646200000012612600626400EAEAllEIEIEIEIllllEIEIEIEIllllKEAEAllEIEIEIEIllllEIEIEIEIllll工程计算实践作业3转换矩阵则按下式计算。cossin0000sincos0000001000000cossin0000sincos0000001T计算程序框图如图2所示，具体的程序代码见附录1。工程计算实践作业4Start输入刚架结点坐标矩阵x输入各单元定位向量locvec输入单元截面参数E,A,I生成转换矩阵T计算单元长度向量l(i)输入结点荷载列阵f（整体坐标系）输入结点等效荷载列阵fe（整体坐标系）生成单元刚度矩阵ki（局部坐标系）将ki(局部坐标系)转换Ki(整体坐标系)组装总体刚度矩阵K生成综合结点荷载列阵F计算结点位移d=F\K计算各单元杆端力FiEnd图2MATLAB矩阵分析法程序框图工程计算实践作业5三、解题步骤取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图4所示，局部坐标系用单元中箭头的方向表示，原始数据如下：xθy123456①②③④⑤图3图4刚架结点输入矩阵为，x=[00;0-5;1.63-6.37;4-5;4-1;42];各单元定位向量为，locvec1=[123000];locvec2=[123456];locvec3=[456789];locvec4=[789101112];locvec5=[101112000];输入截面参数，E=2.1e11;%E=210GPaa=0.12;%矩形截面长0.12mb=0.2;%矩形截面宽0.2m输入整体坐标系下各单元结点荷载列阵，f(1,:)=zeros(1,6);f(2,:)=[000040e30];f(3,:)=zeros(1,6);f(4,:)=[000-50e300];f(5,:)=zeros(1,6);输入整体坐标系下单元1等效节点荷载工程计算实践作业6q=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.00350.0000-0.00040.0030-0.0005-0.00040.00270.00000.0016-0.00510.0000-0.0006各个单元的杆端力如表1所示，表1各单元杆端力单元12345i端Fx(kN)-17917.04717917.0517917.0517917.05-32083Fy(kN)17507.3731-17507.422492.6322492.6322492.63Mz(kN·m)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(kN·m)-37312.596-2092.8326668.34-44999.851249.01四、计算结果校核在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。工程计算实践作业7图5有限元模型图6轴力图(kN)工程计算实践作业8图7剪力图（kN）图8弯矩图（kN·m）从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果比较，得到表2、3。表2各节点位移比较节点号项目矩阵位移法ANSYS误差1Ux(m)000Uy(m)000Rotz(rad)0002Ux(m)0.0034780.00348-2E-06Uy(m)0.00001740.00001740工程计算实践作业9Rotz(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(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.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的结果对比可知，两种方法的计算结果十分接近，误差均可以忽略不计，从而验证了计算结果的可靠性与准确性。工程计算实践作业10四、结论通过对一个平面刚架静力分析问题的求解，本文编制的平面刚架静力分析程序计算结果与有限元软件ANSYS计算结果校核后，发现两者计算结果十分接近，误差可忽略不计，说明该程序计算结果的可靠性与精确性。工程计算实践作业11附录1矩阵位移法计算程序pmgj.m计算主程序clearclc%结点坐标x=[00;0-5;1.63-6.37;4-5;4-1;42];%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=[123000];locvec2=[123456];locvec3=[456789];locvec4=[789101112];locvec5=[101112000];%刚架的材料特性截面特性E=2.1e11;%E=210GPaa=0.12;%矩形截面长0.12mb=0.2;%矩形截面宽0.2mA=a*b;I=(a*b^3)/12;clearab%单元长度计算fori=1:5dx=x(i,1)-x(i+1,1);dy=x(i,2)-x(i+1,2);l(i)=(dx^2+dy^2)^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));工程计算实践作业12t3=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,:)=[000040e30];f(3,:)=zeros(1,6);f(4,:)=[000-50e300];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)';clearNonLogijIJ工程计算实践作业13NonLog=(locvec2~=0);ij=find(NonLog);IJ=locvec2(NonLog);K(IJ,IJ)=K(IJ,IJ)+K2(ij,ij);F(IJ)=F(IJ)+f(2,ij)';clearNonLogijIJNonLog=(locvec3~=0);ij=find(NonLog);IJ=locvec3(NonLog);K(IJ,IJ)=K(IJ,IJ)+K3(ij,ij);F(IJ)=F(IJ)+f(3,ij)';clearNonLogijIJNonLog=(locvec4~=0);ij=find(NonLog);IJ=locvec4(NonLog);K(IJ,IJ)=K(IJ,IJ)+K4(ij,ij);F(IJ)=F(IJ)+f(4,ij)';clearNonLogijIJNonLog=(locvec5~=0);ij=find(NonLog);IJ=locvec