《有限元基础教程》_【MATLAB算例】3.2.5四杆桁架结构的有限元分析(Bar2D2Node)

【MATLAB算例】3.2.5(2)四杆桁架结构的有限元分析(Bar2D2Node)如图3-8所示的结构，各个杆的弹性模量和横截面积都为4229.510/ENmm,2100Amm。试基于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,角度alpha1,alpha2和alpha3,然后分别针对单元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=737500000-7375000000000-7375000007375000000000k2=Bar2D2Node_Stiffness(E,A,x2,y2,x3,y3,alpha2)k2=1.0e+007*0.00000.0000-0.0000-0.00000.00009.8333-0.0000-9.8333-0.0000-0.00000.00000.0000-0.0000-9.83330.00009.8333k3=Bar2D2Node_Stiffness(E,A,x1,y1,x3,y3,alpha3)k3=1.0e+007*3.77602.8320-3.7760-2.83202.83202.1240-2.8320-2.1240-3.7760-2.83203.77602.8320-2.8320-2.12402.83202.1240k4=Bar2D2Node_Stiffness(E,A,x4,y4,x3,y3,alpha1)k4=737500000-7375000000000-7375000007375000000000（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.11510.2832-0.73750-0.3776-0.2832000.28320.212400-0.2832-0.212400-0.737500.73750.0000-0.0000-0.000000000.00000.9833-0.0000-0.983300-0.3776-0.2832-0.0000-0.00001.11510.2832-0.73750-0.2832-0.2124-0.0000-0.98330.28321.1957000000-0.737500.7375000000000（4）边界条件的处理及刚度方程求解由图3-8可以看出，节点1的位移将为零，即10u,10v，节点2的位移20v，节点4的40u，40v。节点载荷3F＝10N。采用高斯消去法进行求解，注意：MATLAB中的反斜线符号“\”就是采用高斯消去法。该结构的节点位移为：Tvuvuvuvu44332211q而节点力为：TyxyyxRRRRR4442411105.20102FRP其中，11(,)xyRR为节点1处沿x和y方向的支反力，2yR为节点2处y方向的支反力，44(,)xyRR为节点4处沿x和y方向的支反力。相关的计算流程如下。k=KK([3,5,6],[3,5,6])k=1.0e+008*0.7375-0.0000-0.0000-0.00001.11510.2832-0.00000.28321.1957p=[20000;0;-25000];u=k\pu=1.0e-003*0.27120.0565-0.2225[这里将列排成了一行，以节省篇幅]由此可以看出，所求得的结果2330.2712,0.0565,0.2225ummummvmm，则所有节点位移为000.271200.05650.222500Tmmq(3-75)与前面通过数学推导得到的结果相同，见式(3-72)。（5）支反力的计算在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。将整体的位移列阵q（采用国际单位）代回原整体刚度方程，计算出所有的节点力P，按上面的对应关系就可以找到对应的支反力。相关的计算流程如下。q=[000.000271200.0000565-0.000222500]'q=1.0e-003*000.271200.0565-0.222500[这里将列排成了一行，以节省篇幅]P=KK*qP=1.0e+004*-1.58330.31262.00012.1879-0.0001-2.5005-0.41670[这里将列排成了行]按对应关系，可以得到对应的支反力为(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*000.27120[这里将列排成了一行，以节省篇幅]stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha1,u1)stress1=2.0001e+008u2=[q(3);q(4);q(5);q(6)]u2=1.0e-003*0.271200.0565-0.2225[这里将列排成了一行，以节省篇幅]stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,alpha2,u2)stress2=-2.1879e+008u3=[q(1);q(2);q(5);q(6)]u3=1.0e-003*000.0565-0.2225[这里将列排成了一行，以节省篇幅]stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,alpha3,u3)stress3=-52097000u4=[q(7);q(8);q(5);q(6)]u4=1.0e-003*000.0565-0.2225[这里将列排成了一行，以节省篇幅]stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,alpha1,u4)stress4=41668750可以看出：计算得到的单元1的应力为82.000110Pae；单元2的应力为82.187910Pa，单元3的应力为75.209710Pa，单元4的应力为74.16710Pa。与前面通过数学推导得到的结果相同，见式(3-73)。