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

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

为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应用MATLAB 编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS 大型商用有限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果相当吻合，验证了计算结果的可靠性。

一、 问题描述

如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa ，截面为0.12×0.2m 的实心矩形，现要求解荷载作用下刚架的位移和内力。

5m

4m

3m

图1

二、矩阵位移法计算程序编制

为编制程序方便考虑，本文计算中采用“先处理法”。具体的计算步骤如下。

（1） 对结构进行离散化，对结点和单元进行编号，建立结构（整体）坐标系

和单元（局部）坐标系，并对结点位移进行编号； （2） 对结点位移分量进行编码，形成单元定位向量e λ；

（3） 建立按结构整体编码顺序排列的结点位移列向量δ，计算固端力e F P 、等

效结点荷载E P 及综合结点荷载列向量P ；

（4） 计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵T 形成整体坐标

系下的单元刚度矩阵e T e K T K T = ； （5） 利用单元定位向量形成结构刚度矩阵K ； （6） 按式1=K P δ- 求解未知结点位移； （7） 计算各单元的杆端力e F 。

根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计算。

32322

23

2

32

22

0000

1261260

064620

00001261260062640

EA

EA l

l EI EI

EI EI l l l l EI EI EI EI l l l l K EA EA l l EI EI EI EI l l l l EI EI EI EI l l

l l ⎡⎤-

⎢⎥⎢

⎥⎢⎥-

⎢⎥

⎢

⎥⎢⎥-

⎢⎥⎢

⎥=⎢⎥-⎢⎥

⎢

⎥⎢⎥---⎢⎥⎢

⎥⎢⎥-⎢⎥⎣

⎦

^.

转换矩阵则按下式计算。

cos sin 0000sin cos 0000001000000cos sin 0000sin cos 00

1T ααααααα

α⎡⎤⎢⎥-⎢⎥⎢⎥=⎢

⎥⎢⎥⎢⎥-⎢⎥⎣⎦

计算程序框图如图2所示，具体的程序代码见附录1。

图2 MATLAB矩阵分析法程序框图

三、解题步骤

取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图4所示，局部坐标系用单元中箭头的方向表示，原始数据如下：

1

3

4

5

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所示，

四、计算结果校核

在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。

图5 有限元模型图6 轴力图(kN)

图7 剪力图（kN）

图8 弯矩图（kN·m）

从ANSYS计算结果中提取各结点位移、内力，并与矩阵位移法分析的结果比较，得到表2、3。

节点号项目矩阵位移法ANSYS 误差

1

Ux(m) 0 0 0 Uy(m) 0 0 0 Rotz(rad) 0 0 0

2

Ux(m) 0.003478 0.00348 -2E-06 Uy(m) 0.0000174 0.0000174 0 Rotz(rad) -0.00037 -0.000365 -5E-06