MATLAB在结构力学分析中的应用

吴小希 周小波 李彬

摘 要:传统的手算方法解超静定结构工作量繁重，有时甚至是不可能，运用结构有限元编程的一般方法，通过两个实例的对照，展示MATLAB在结构力学分析中的应用，MATLAB具有高性能，方法具有普遍的适用性,实现弯矩图自动绘制。

关键词: MATLAB结构有限元弯矩图

Abstract：While using the traditional manual method to resolve complex statically indeterminate structures, it is heavy workloads, sometimes even impossible，using finite element programming of the general method, Based on two examples, This paper introduces a method of application of MATLAB in structure mechanics, MATLAB has the advantages of high performance, it can be applied to many kinds of structures, realization of automatic drawing bending moment diagram.

Key words: MATLAB; Finite element; Bend moment diagram

引言

結构力学[3]中，常利用传统的力法与位移法求解超静定结构，力法是几何问题，位移法把复杂的几何图乘转化为代数运算，但它们基本未知量很多时，系数构成的矩阵计算巨大，两者都不能满足科研工作者的需要。应用MATLAB软件丰富可靠的矩阵运算、数据处理、图形绘制等便利工具，可使得计算和图象一体化。对于结构力学计算是十分有利的工具。

1基本方法

MATLAB结构有限元编程的基本思路是先分后合，即将结构分成各个单元和节点，桁架与刚架已经离散化，对于连续系统这一步极其重要，然后进行单元分析，集成整体刚度矩阵，引入边界条件，最后解方程。在求解平面桁架结构，虽然结构简单，用手算可得各杆件的轴力，但重复的过程太多，现在使用MATLAB语言来编制有限元位移法的程序时，则编程的难度明显降低，对有限元位移法的概念的理解更加深入，编程所需时间也大大减少。

图1为一平面桁架，各杆E=70GPaA=0.004，试用矩阵位移法求解各杆轴力

图1

解：平面桁架元是既有局部坐标又有总体坐标的二维有限元；对各结点和单元进行编号，建立结构坐标系( 图1 )

第一步,利用MATLAB函数

y=Plane Truss Element Length(x1, y1, x2, y2)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); ％ 局部坐标中杆件长度

第二步, MATLAB函y=Plane Truss Element Stiffness(E ,A ,L ,theta)

x=theta*pi/180; C= cos (x); S=sin(x);

y=E*A/L*[C*C C*S -C*C -C*S; C*S S*S -C*S -S*S;-C*C -C*S C*C C*S;-C*S -S*S C*S S*S];% 总体坐标中建立各单元的刚度 矩阵

第三步,建立整体刚度阵。该结构有4个节点,每个节点有两个自由度(可考虑支座沉降),为了得到整体刚度阵K，首先利用生成一个8×8的0矩阵，因为该结构有4个单元，所以4次调用M a t lab的Plane Truss Assemble函数;其中K为整体刚度阵, k为单元刚度阵, i j为单元两端在整体节点上的编号。

y=Plane Truss Assemble (K, k, i , j)

K (2*i-1, 2*i-1) =K (2*i-1, 2*i-1) +k (1, 1);

K (2*i-1, 2*i) =K (2*i-1, 2*i) + k (1, 2);

K (2*i-1, 2*j-1) = K (2*i-1, 2*j-1)+ k (1,3);

K (2*i-1, 2*j) =K (2*i-1, 2*j) +k (1, 4);

K (2* i , 2*i-1) =K (2* i, 2*i-1) +k (2, 1);

K (2*i, 2*i) =K (2*i, 2*i) +k (2, 2);

K (2*i , 2*j-1)=K(2*i,2*j-1)+k(2,3);

K (2*i, 2*j) =K (2*i, 2*j) + k (2, 4);

K (2*j-1, 2*i-1) =K (2*j, 2*i-1) +k (3, 1);

K (2*j-1, 2*i) =K (2*j-1, 2*i) + k (3, 2);

K (2*j-1, 2*j-1) =K (2*j-1, 2*j-1) +k (3, 3);

K (2*j-1, 2*j) =K (2*j-1, 2*j) + k (3, 4);

K (2*j, 2*i-1) =K (2*j, 2*i-1) + k (4, 1);

K (2*j, 2*i) =K (2*j, 2*i) +k (4, 2);

K (2*j, 2 *j-1) =K (2*j,2*j-1)+k (4,3);

K (2*j, 2*j) = K (2*j, 2*j) +k (4, 4);

y=K;

第四步k=K(3:6,3:6);％边界条件下刚度矩阵

f=[0;30;30;0];％形成荷载向量

u=kf;％分解法和高斯消去法,得到结点位移

u = [0.00100.0006 ；0.0011-0.0003]

％结点2、3的结点位移

U=[0;0;u;0;0]; ％结构各节点位移矢量

第五步，M a t lab函数

Plane Truss Element Force (E, A ,L ,theta ,u)

x=theta*pi/180;C=cos(x);S=sin(x); y=E*A/L*[-C -S C S]*u;

可得：F1 =39.8018F2=0F3 = 28.5646 F4 = -13.8618 F5 = 9.801F6 = -20.1982

%各杆件的轴力

图2E=210GPa, I=5*10-6 q=7KN/M，绘制弯矩图。

图2

解：对连续结构单元进行编号十分重要，梁单元是既有局部坐标又有总体坐标的二维有限元，用线性函数表示，主程序根据交互输入的原始数据形成单元刚度矩阵，再根据整体刚度矩阵集成规则，将单元刚度矩阵形成整体刚度矩阵。通过引人支承条件，然后分解和高斯消去法解方程，得到结点位移，进而求出各单元杆端弯矩。

第一步，MATLAB函数k=Beam Element Stiffness(E,I,L)

y=E*I/(L*L*L)*[12 6*L -12 6*L; 6*L 4*L*L -6*L 2*L*L;-12 -6*L 12 -6*L; 6*L 2*L*L -6*L 4*L*L]; %单元刚度

第二步，整体刚度的建立，两者都是二维有限元，程序相同，根据划分单元数，多次调用函数。

第三步，计算等效节点载荷。按照结构力学的方法可以求得

M = [-9.333 9.333]

第四步，引入边界条件，节点2，3 的转角为a1 a2；其余为0；得出边界条件下结构刚度矩阵k

第五步，MATLAB函数

f=Beam Element Forces（k ,u）求得杆端弯矩

M1 (3.111,-6.222)M2= (-6.222, -6.222) M3= (-6.222,3.111);%杆端弯矩

绘制弯矩图：在命令窗口输入如下命令：

q=7

l1=linspace (0, 4);

M1=linspace (3.111,-6.2 22);

l= [4, 6, 8];

mid=q*4^2/8-(6.222+6.222)/2

M= [ -6.222 , mid,-6.222]

P= polyfit (l, M, 2)

l2=4:0.1:8

M2=P (1)*l2.^2+P(2)*l2+P(3);

l3=linspace (8, 12);

M3=linspace (-6.222, 3.111);

plot (l1,M1,l2,M2,l3,M3) 自動生成弯矩图3:

2．结束语

通过两个例子表明, 在结构力学中引入MATLAB，简单的分步编程，即可完成有关问题的过程分析、大量计算和绘图，平面桁架与连续梁单元调用同一函数，即可求出整体刚度矩阵，大大提高了效率。更多地了解和掌握MATLAB,对于我们的教学和科研工作将是十分有益的。

参考文献：

[1]P. I. Kattan著, 韩来彬译. MATLAB有限元分析与应用[M].北京:清华大学出版社,2004.

[2]马晓光,于国清．MATLAB在结构力学中的应用．白城师范学院学报[J]2006，

20（4）99-102．

[3]龙驭球，包世华．结构力学I (第二版)［M］．北京：高等教育出版社，2006．

作者简介：

吴小希(1987-)，男（汉族），湖南新化,湖南科技大学研究生，主要研究领域为桥梁的振动控制。单位：湖南科技大学土木工程学院桥梁研究所

注：文章内所有公式及图表请用PDF形式查看。