《有限元基础教程》_【MATLAB算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)

如图4-20所示为一矩形薄平板，在右端部受集中力F100 000N作用，材料常数为：

7弹性模量E110Pa，泊松比1，板的厚度t0.1m。基于MATLAB平台求解该

3结构的节点位移、支反力以及单元应力。

图4-20

解答：对该问题进行有限元分析的过程如下。 （1）结构的离散化与编号

将结构离散为二个3节点三角形单元，单元编号及节点编号如图4-20(b)所示。 （2）计算各单元的刚度矩阵(以国际标准单位)

首先在MATLAB环境下，输入弹性模量E、泊松比NU、薄板厚度t和平面应力问题性质指示参数ID，然后针对单元1和单元2，分别两次调用函数Triangle2D3Node_Stiffness，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)。

>> E=1e7; >> NU=1/3; >> t=0.1; >> ID=1;

>> k1=Triangle2D3Node_Stiffness(E,NU,t,2,0,0,1,0,0,ID) k1 = 1.0e+006 *

0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750

-0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188

>>k2=Triangle2D3Node_Stiffness(E,NU,t,0,1,2,0,2,1,ID) k2 = 1.0e+006 *

0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188

（3） 建立整体刚度方程

由于该结构共有4个节点，则总共的自由度数为8，因此，结构总的刚度矩阵为KK(8×8)，先对KK清零，然后两次调用函数Triangle2D3Node_Assembly进行刚度矩阵的组装。

>>KK = zeros(8,8);

>>KK=Triangle2D3Node_Assembly(KK,k1,2,3,4); >>KK=Triangle2D3Node_Assembly(KK,k2,3,2,1) KK = 1.0e+006 *

Columns 1 through 6

0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 -0.3750 -0.1875 0.6563 0 0 0.3750 -0.1875 -1.1250 0 1.2188 0.3750 0 -0.2813 -0.1875 0 0.3750 0.6563 0 -0.1875 -0.0938 0.3750 0 0 1.2188

0 0 -0.2813 -0.1875 -0.3750 -0.1875 0 0 -0.1875 -0.0938 -0.1875 -1.1250

Columns 7 through 8 0 0 0 0 -0.2813 -0.1875 -0.1875 -0.0938 -0.3750 -0.1875 -0.1875 -1.1250 0.6563 0.3750 0.3750 1.2188

（4） 边界条件的处理及刚度方程求解

由图4-20(b)可以看出，节点3和节点4的两个方向的位移将为零，即

u30,v30,u40,v40，因此，将针对节点1和节点2的位移进行求解，节点1和节点

2的位移将对应KK矩阵中的前4行和前4列，则需从KK(8×8)中提出，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解。注意：MATLAB中的反斜线符号“\\”就是采用高斯消去法。

>>k=KK(1:4,1:4) k = 1.0e+006 *

0.6563 0.3750 -0.3750 -0.1875 0.3750 1.2188 -0.1875 -1.1250 -0.3750 -0.1875 0.6563 0 -0.1875 -1.1250 0 1.2188

>>p=[0;-5000;0;-5000]

p = 0 -5000 0 -5000 [这里将列排成了一行，以节省篇幅]

>>u=k\\p

u = 0.0188 -0.09 -0.0150 -0.0842 [这里将列排成了一行，以节省篇幅]

由此可以看出，所求得的结果为：u10.018 8,v10.0 9,u20.015 1,v20.084 2。 （5）支反力的计算

由方程(4-183)可知，在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。先将上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵U(8×1)，再代回原整体刚度方程，计算出所有的节点力P(8×1)，按式(4-179)的对应关系就可以找到对应的支反力。

>>U=[u;0;0;0;0]; >>P=KK*U P = 1.0e+004 *

-0.0000 -0.5000 0 -0.5000 -2.0000 -0.0702 2.0000 1.0702 [这里将列排成了一行]

由式(4-179)的对应关系，可以得到对应的支反力为Rx320 000,Ry3702,Rx420 000,Ry410 702。

（6）各单元的应力计算

先从整体位移列阵U(8×1)中提取出单元的位移列阵，然后，调用计算单元应力的函数Triangle2D3Node_Stress，就可以得到各个单元的应力分量。

>>u1=[U(3);U(4);U(5);U(6);U(7);U(8)]

u1 = -0.0150 -0.0842 0 0 0 0

>>stress1=Triangle2D3Node_Stress(E,NU,2,0,0,1,0,0,u1,ID) stress1 = 1.0e+005 *

-0.8419 -0.2806 -1.5791 [这里将列排成了一行，以节省篇幅]

>>u2=[U(5);U(6);U(3);U(4);U(1);U(2)]

u2 = 0 0 -0.0150 -0.0842 0.0188 -0.09 [这里将列排成了一行，以节省篇幅]

>>stress2=Triangle2D3Node_Stress(E,NU,0,1,2,0,2,1,u2,ID) stress2 =1.0e+004 *

8.4187 -2.53 -4.2094 [这里将列排成了一行，以节省篇幅]

可以看出：计算得到的单元1的应力分量为x84 190Pa,y28 060Pa,

xy157 910Pa；单元2的应力分量为

x84 187Pa,

y28 953Pa,xy42 094Pa。