matlab 对梁元的分析与运用

                          勾都  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个节点。

>> 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有限元分析与运用