基于某MATLAB的平面刚架静力分析报告.docx

标准 基于 MATLAB的平面刚架静力分析 为了进一步理解有限元方法计算的过程， 本文根据矩阵位移法的基本原理应 用 MATLAB编制计算程序对以平面刚架结构进行了静力分析。 本文还利用 ANSYS 大型商用有限元分析软件对矩阵位移法的计算结果进行校核， 发现两者计算结果 相当吻合，验证了计算结果的可靠性。 一、 问题描述 如图 1 所示的平面刚架，各杆件的材料及截面均相同， E=210GPa，截面为 0.12 × 0.2m 的实心矩形，现要求解荷载作用下刚架的位移和内力。 40° 40kN 30° 10kN/m 4m 5m 4m 50kN 3m 图 1 二、矩阵位移法计算程序编制 为编制程序方便考虑， 本文计算中采用 “先处理法”。具体的计算步骤如下。 文案 标准 （1） 对结构进行离散化，对结点和单元进行编号，建立结构（整体）坐标系和单元（局部）坐标系，并对结点位移进行编号； （2） 对结点位移分量进行编码，形成单元定位向量 e ； （3） 建立按结构整体编码顺序排列的结点位移列向量 ，计算固端力 PFe 、等 效结点荷载 PE 及综合结点荷载列向量 P ； （4） 计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵 T 形成整体坐标 系下的单元刚度矩阵 K e T T K eT ； （5） 利用单元定位向量形成结构刚度矩阵 K ； （6） 按式 =K 1P 求解未知结点位移； （7） 计算各单元的杆端力 F e 。 根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计 算。 EA 0 0 EA 0 l 0 l 0 12EI 6EI 12EI 6EI l 3 l 2 0 l 3 l 2 6EI 4EI 6EI 2EI 0 l 2 l 0 l 2 l K EA EA 0 0 0 l 0 l 0 12EI 6EI 12EI 6EI l 3 l 2 0 3 l 2 l 0 6EI 2EI 6EI 4EI l 2 l 0 l 2 l 文案 标准 转换矩阵则按下式计算。 cos sin 0 0 0 0 sin cos 0 0 0 0 T 0 0 1 0 0 0 0 0 0 cos sin 0 0 0 0 sin cos 0 0 0 0 0 0 1 计算程序框图如图 2 所示，具体的程序代码见附录 1。 文案 标准 Start 输入刚架结点坐标矩阵 x 输入各单元定位向量 locvec 输入单元截面参数 E,A,I 计算单元长度向量 l(i) 生成转换矩阵 T 输入结点荷载列阵 f （整体坐标 系） 输入结点等效荷载列阵 fe （整体坐 标系） 生成单元刚度矩阵 ki （局部坐标 系） 将ki( 局部坐标系 ) 转换 Ki( 整体坐 标系 ) 组装总体刚度矩阵 K 生成综合结点荷载列阵 F 计算结点位移 d=F\K 计算各单元杆端力 Fi End 图 2 MATLAB矩阵分析法程序框图 文案 标准 三、解题步骤 取整体坐标系如图 3 所示，对结构进行离散化， 对结点和单元进行编号如图 所示，局部坐标系用单元中箭头的方向表示，原始数据如下： 3 ③ ② 2 4 ④ ① x 5 1 θ ⑤ y 6 图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 Fx(kN) -17917.047 17917.05 17917.