【MATLAB算例】3.2.5(2) 四杆桁架结构的有限元分析(Bar2D2Node)

如图3-8所示的结构，各个杆的弹性模量和横截面积都为E?29.5?10N/mm,

42A?100mm2。试基于MATLAB平台求解该结构的节点位移、单元应力以及支反力。

图3-8 四杆桁架结构

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

对该结构进行自然离散，节点编号和单元编号如图3-8所示，有关节点和单元的信息见表3-1~表3-3。

（2）计算各单元的刚度矩阵(基于国际标准单位)

建立一个工作目录，将所编制的用于平面桁架单元分析的4个MATLAB函数放置于该工作目录中，分别以各自函数的名称给出文件名，即：Bar2D2Node_Stiffness，Bar2D2Node_Assembly，Bar2D2Node_Stress，Bar2D2Node_Forces。然后启动MATLAB，将工作目录设置到已建立的目录中，在MATLAB环境中，输入弹性模量E、横截面积A，各点坐标x1,y1,x2,y2,x3,y3,x4,y4,角度alpha 1, alpha 2和alpha 3,然后分别针对单元1，2，3和4，调用4次Bar2D2Node_Stiffness，就可以得到单元的刚度矩阵。相关的计算流程如下。

>>E=2.95e11; >>A=0.0001; >>x1=0; >>y1=0; >>x2=0.4; >>y2=0; >>x3=0.4; >>y3=0.3; >>x4=0; >>y4=0.3; >> alpha1=0; >> alpha2=90;

>> alpha3=atan(0.75)*180/pi;

>> k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2,alpha1)

k1 = 73750000 0 -73750000 0 0 0 0 0 -73750000 0 73750000 0 0 0 0 0 >> k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3,alpha2) k2 = 1.0e+007 *

0.0000 0.0000 -0.0000 -0.0000 0.0000 9.8333 -0.0000 -9.8333 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -9.8333 0.0000 9.8333

>> k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3) k3 = 1.0e+007 *

3.7760 2.8320 -3.7760 -2.8320 2.8320 2.1240 -2.8320 -2.1240 -3.7760 -2.8320 3.7760 2.8320 -2.8320 -2.1240 2.8320 2.1240

>> k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3,alpha1) k4 = 73750000 0 -73750000 0 0 0 0 0 -73750000 0 73750000 0 0 0 0 0

（3） 建立整体刚度方程

由于该结构共有4个节点，因此，设置结构总的刚度矩阵为KK(8×8)，先对KK清零，然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。相关的计算流程如下。

>>KK=zeros(8,8);

>>KK=Bar2D2Node_Assembly (KK,k1,1,2); >>KK=Bar2D2Node_Assembly (KK,k2,2,3); >>KK=Bar2D2Node_Assembly (KK,k3,1,3); >>KK=Bar2D2Node_Assembly (KK,k4,4,3) KK= 1.0e+008 *

1.1151 0.2832 -0.7375 0 -0.3776 -0.2832 0 0 0.2832 0.2124 0 0 -0.2832 -0.2124 0 0 -0.7375 0 0.7375 0.0000 -0.0000 -0.0000 0 0 0 0 0.0000 0.9833 -0.0000 -0.9833 0 0 -0.3776 -0.2832 -0.0000 -0.0000 1.1151 0.2832 -0.7375 0 -0.2832 -0.2124 -0.0000 -0.9833 0.2832 1.1957 0 0 0 0 0 0 -0.7375 0 0.7375 0 0 0 0 0 0 0 0 0

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

由图3-8可以看出，节点1的位移将为零，即u1?0, v1?0，节点2的位移v2?0，节点4的u4?0，v4?0。节点载荷F3＝10N。采用高斯消去法进行求解，注意：MATLAB中的反斜线符号“\\”就是采用高斯消去法。

该结构的节点位移为：q?u1而节点力为：P?R?F?Rx1?v1u2Ry1v2u3v3u42?104Ry2v4?

T?0?2.5?104Rx4Ry4

?T其中，(Rx1,Ry1)为节点1处沿x和y方向的支反力，Ry2为节点2处y方向的支反力，

(Rx4,Ry4)为节点4处沿x和y方向的支反力。相关的计算流程如下。

>>k=KK([3,5,6],[3,5,6]) k =1.0e+008 *

0.7375 -0.0000 -0.0000 -0.0000 1.1151 0.2832 -0.0000 0.2832 1.1957

>> p=[20000;0;-25000]; >>u=k\\p

u = 1.0e-003 *

0.2712 0.0565 -0.2225 [这里将列排成了一行，以节省篇幅]

由此可以看出，所求得的结果u2?0.271 2mm,u3?0.056 5mm,v3??0.222 5mm，则所有节点位移为

q??000.271 200.056 5?0.222 500?mm (3-75)

与前面通过数学推导得到的结果相同，见式(3-72)。

（5）支反力的计算

在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。将整体的位移列阵q（采用国际单位）代回原整体刚度方程，计算出所有的节点力P，按上面的对应关系就可以找到对应的支反力。相关的计算流程如下。

>> q=[0 0 0.0002712 0 0.0000565 -0.0002225 0 0]' q = 1.0e-003 *

0 0 0.2712 0 0.0565 -0.2225 0 0 [这里将列排成了一行，以节省篇幅] >>P=KK*q P = 1.0e+004 *

-1.5833 0.3126 2.0001 2.1879 -0.0001 -2.5005 -0.4167 0 [这里将列排成了行]

T

按对应关系，可以得到对应的支反力为

(3-76)

（6）各单元的应力计算

先从整体位移列阵q中提取出单元的位移列阵，然后，调用计算单元应力的函数Bar2D2Node_ElementStress，就可以得到各个单元的应力分量。当然也可以调用上面的Bar2D2Node_ElementForces(E,A,x1,y1,x2,y2,alpha,u)函数来计算单元的集中力，然后除以面积求得单元应力。相关的计算流程如下。

>>u1=[q(1);q(2);q(3);q(4)] u1 = 1.0e-003 *

0 0 0.2712 0 [这里将列排成了一行，以节省篇幅] >> stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha1,u1)