基于MATLAB的平面刚架有限元分析,基于MATLAB的平面刚架静力分析

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

基于MATLAB的平面刚架静力分析 为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应用MATLAB编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS大型商用有限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果相当吻合，验证了计算结果的可靠性。 一、 问题描述 如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa，截面为0.12×0.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=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 各单元杆端力 单元 1 2 3 4 5 i端 Fx(kN) -17917.047 17917.05 17917.05 17917.05 -32083 Fy(kN) 17507.3731 -17507.4 22492.63 22492.63 22492.63 Mz(kN·m) 1897.83076 -1897.83 2092.833 -26668.3 44999.85 j端 Fx(kN) -32082.953 -17917 -17917 32082.95 32082.95 Fy(kN) -17507.373 -22492.6 -22492.6 -22492.6 -22492.6 Mz(kN·m) -37312.596 -2092.83 26668.34 -44999.8 51249.01 四、计算结果校核 在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。 图5 有限元模型 图6 轴力图(kN) 图7 剪力图(kN) 图8 弯矩图(kN·m) 从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果比较，得到表2、3。 表2 各节点位移比较 节点号 项目 矩阵位移法 ANSYS 误差 1 Ux(m) 0 0 0 Uy(m) 0 0 0 Rotz(rad) 0 0 0 2 Ux(m) 0. 0.00348 -2E-06 Uy(m) 0. 0. 0 Rotz(rad) -0.00037 -0. -5E-06 3 Ux(m) 0. 0.00302 -2E-06 Uy(m) -0.00051 -0. 2E-06 Rotz(rad) -0.00038 -0. -2E-06 4 Ux(m) 0. 0.00269 -3E-06 Uy(m) 0. 0. 0 Rotz(rad) 0. 0.00162 4E-06 5 Ux(m) -0.00513 -0. 1.5E-05 Uy(m) 0. 0. 0 Rotz(rad) -0.00056 -0. -3E-06 6 Ux(m) 0 0 0 Uy(m) 0 0 0 Rotz(rad) 0 0 0 表3 各结点内力比较 节点号 项目 矩阵位移法 ANSYS 误差 1 Fx(kN) -32.083 -32.080 -0.003 Fy(kN) -17.507 -17.503 -0.004 Mz(kN·m) -37.313 -37.307 -0.006 2 Fx(kN) -17.917 -17.920 0.003 Fy(kN) 17.507 17.503 0.004 Mz(kN·m) 1.898 1.909 -0.011 3 Fx(kN) -17.917 -17.920 0.003 Fy(kN) -22.493 -22.497 0.004 Mz(kN·m) -2.093 -2.110 0.018 4 Fx(kN) -17.917 -17.920 0.003