勾都 2010021212
摘要：随着现代科技的快速发展，人们在工程领域的研究也越来越深入，对科学的研究要求快速化，简单化，精确化，实用化。所以市场上出现了一大批针对工程分析与运用的软件，matlab就以实用，简单，精确而为广大用户推重。ATLAB 的名称源自 Matrix Laboratory ，它是一种科学计算软件，专门以矩阵的形式处理数据。 MATLAB 将高性能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。目前 MATLAB 产品族可以用来进行：数值分析 数值和符号计算 工程与科学绘图 数字图像处理 数字信号处理 通讯系统设计与仿真 财务与金融工程。本文是基于MATLAB的对工程分析，主要介绍了用MATLAB对梁元的分析与计算，和相关图的绘制。
关键词：工程分析 梁元 图的绘制
引言：1有限元法的步骤
（1）离散化域
（2）写出单元刚度矩阵
（3）集成整体刚度矩阵
（4）引入边界条件
（5）解方程
（6）后处理
由上步骤可看出，结决问题的过程结合使用了matlab和某些有限的手动操作（步骤1、4、5）。可以看出，所有冗长、反复的计算都可由matlab完成。
2用于有限元分析的函数
1 BeamElementStiffness(E,I,L)——该函数用于计算弹性模量E、转动惯量I、长度L的梁的单元刚度矩阵。返回4x4的单元刚度矩阵k。
2 BeamAssemble（K,k,i,j）——该函数连接节点i和节点j的梁元的单元刚度矩阵k集成到整体刚度矩阵K。每集成一个单元，该函数都返回2nX2n的整体刚度矩阵K.
3 BeamElementForces(k,u)——该函数用单元刚度矩阵k和单元节点位移矢量u计算单元节点矢量。返回4x1的单元节点力矢量f。
4 BeamElementShearDiagram(f,L)——该函数绘制节点力矢量为f和长度为L的单元剪力图。
5 BeamElementMomentDiagram(f,L)——该函数绘制节点力矢量f和长度L的单元弯矩曲线图。
基础知识:梁元是总体坐标和局部坐标一致的二维有限元，用线性函数描述。梁元的系数有弹性模量 E、惯性矩 I、长度 L。如下图1-1。每个梁元有2个节点，并且假定他是水平的。忽略轴向的形变，元刚度矩阵如下
梁有4个自由度——每个节点有2个自由度（横位移和转角）。约定位移向上为正，转角逆时针为正。所以，有n个节点的结构其整体刚度矩阵K是2nX2n。
更据其整体刚度矩阵k，就可求出以下方程组：
U为结构点位移矢量，F是结构点载荷矢量。边界条件被手动赋值给矢量U和F。然后用分解和高斯消去法解上方程组。一旦求出未知的位移和支反力，就可用下式求出每个单元的节点力矢量：
f是4x1的单元节点力矢量，u是4x1的单元节点位移矢量。每个u矢量的第一个和第2个分量分别是第一个节点的横位移和转角，第3个和第4个分量则分别是第2个节点的横位移和转角。
实际运用：如下图的梁结构。假设求：
（1）该结构的整体刚度矩阵
（2）节点2的位移
（3）节点2和节点3的转角
（4）节点1和节点3的支反力
（5）每个单元的力（剪力和弯矩）
（6）每个单元的剪力图
（7）每个单元的弯矩图
离散化域
将一个节点放置在集中载荷作用点位置以便求出改点的待求量（位移、转角、剪力、弯矩）。所以我们将定义域分为2个单元3个节点。
| 单元编号 | 节点i | 节点j |
| 1 2 | 1 2 | 2 3 |
>> E=210e6
E =
>> I=60e-6
I =
6.0000e-005
>> L=2
L =
>> k1=BeamElementStiffness(E,I,L)
>> k1=BeamElementStiffness(E,I,L)
k1 =
100 100 -100 100
100 25200 -100 12600
-100 -100 100 -100
100 12600 -100 25200
>> k2=BeamElementStiffness(E,I,L)
k2 =
100 100 -100 100
100 25200 -100 12600
-100 -100 100 -100
100 12600 -100 25200
集成整体刚度矩阵：
该结构有3个节点，所以整体刚度矩阵是6x6的。因此，为了得到整体刚度矩阵我们要生成一个6x6的零矩阵。由于该结构只有2个梁元，所以我们只需两次调用matlab的BeamAssemble函数就可以得到整体刚度矩阵K。每次对函数的调用都会集成一个单元。
>> K=zeros(6,6)
K =
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
>> K=BeamAssemble(K,k1,1,2)
K =
100 100 -100 100 0 0
100 25200 -100 12600 0 0
-100 -100 100 -100 0 0
100 12600 -100 25200 0 0
0 0 0 0 0 0
0 0 0 0 0 0
>> K=BeamAssemble(K,k2,2,3)
K =
100 100 -100 100 0 0
100 25200 -100 12600 0 0
-100 -100 37800 0 -100 100
100 12600 0 50400 -100 12600
0 0 -100 -100 100 -100
0 0 100 12600 -100 25200
引入边界条件：
由上得到的整体刚度矩阵，可得到该结构的方程组。
边界条件：
将边界条件代入方程组，得：
解方程：
用分解法和高斯消去法求解方程组。对方程组进行分解，提取整体刚度矩阵K的第3行、第4行的第3列、第4列、第6列，第6行的第3列、第4列，以及第6行的第6列作为子矩阵。得到：
用matlab得
>> k=[K(3:4,3:4) K(3:4,6);K(6,3:4) K(6,6)]
k =
37800 0 100
0 50400 12600
100 12600 25200
>> f=[-20;0;0]
f =
>> u=k\\f
u =
1.0e-003 *
由结果可看出：节点2的垂直位移是0.9259m（竖直向下），节点2和节点3的转角分别为0.1984rad（顺时针方向）0.7937rad（逆时针方向）。
后处理：
我们用matlab的命令求出节点1的支反力和内力（剪力和弯矩）。建立结构节点位移矢量U，计算节点力矢量F。
>> U=[0;0;u(1);u(2);0;u(3)]
U =
1.0e-003 *
0
0
-0.9259
-0.1984
0
0.7937
>> F=K*U
F =
13.7500
15.0000
-20.0000
0.0000
6.2500
0.0000
由上结果可知：节点1的支反力13.75KN(方向向上)，弯矩15KN.m(逆时针方向)。节点3的支反力6.25KN(方向向上)。
建立单元节点位移矢量u1和u2，调用matlab函数计算单元力矢量f1和f2.
>> f1=BeamElementForces(k1, u1)
f1 =
13.7500
15.0000
-13.7500
12.5000
>> f2=BeamElementForces(k2, u2)
f2 =
-6.2500
-12.5000
6.2500
0.0000
由上可得每个单元的每一端的剪力和弯矩。单元1左端的剪力是13.75KN，弯矩是-15KN；右端剪力是-13.75KN，弯矩是12.5KN。单元2左端的剪力是-6.25KN，弯矩是-12.5KN；右端剪力是6.25KN，弯矩是0KN。
绘制图形：
调用matlab函数绘制剪力曲线图和弯矩曲线图。
>> BeamElementShearDiagram(f1,L)
>> BeamElementShearDiagram(f2,L)
>> BeamElementMomentDiagram(f1,L)
>> BeamElementMomentDiagram(f2,L)
结论：
由上可看出，matlab编的计算函数可以精确的计算工程中复杂的，有大量计算的问题。
并且修改matlab函数的个别参数，就可运用到其他的相似工程问题中，如对网格元，二次四边形元…….的计算。
参考文件：
Matlab论坛上相关资料，
Matab有限元分析与运用下载本文
