2D四杆桁架结构的有限元分析实例

实例：2D四杆桁架结构的有限元分析学习有限元方法的一个最佳途径，就是在充分掌握基本概念的基础上亲自编写有限元分析程序，这就需要一个良好的编程环境或平台。

MATLAB软件就是这样一个平台，它以功能强大、编程逻辑直观、使用方便见长。

将提供有限元分析中主要单元完整的MATLAB程序，并给出详细的说明。

1D杆单元的有限元分析MATLAB程序(Bar1D2Node)最简单的线性杆单元的程序应该包括单元刚度矩阵、单元组装、单元应力等几个基本计算程序。

下面给出编写的线性杆单元的四个MATLAB函数。

Bar1D2Node _Stiffness(E,A,L)该函数计算单元的刚度矩阵，输入弹性模量E，横截面积A和长度L，输出单元刚度矩阵k(2×2)。

Bar1D2Node _Assembly(KK,k,i,j)该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j，输出整体刚度矩阵KK。

Bar1D2Node _Stress(k,u,A)该函数计算单元的应力，输入单元刚度矩阵k、单元的位移列阵u(2×1)以及横截面积A计算单元应力矢量，输出单元应力stress。

Bar1D2Node_Force(k,u)该函数计算单元节点力矢量，输入单元刚度矩阵k和单元的位移列阵u(2×1)，输出2×1的单元节点力矢量forces。

基于1D杆单元的有限元分析的基本公式，写出具体实现以上每个函数的MATLAB程序如下。

%%%%%%%%%%% Bar1D2Node %% begin %%%%%%%%%function k=Bar1D2Node_Stiffness(E, A, L)%该函数计算单元的刚度矩阵%输入弹性模量E，横截面积A和长度L%输出单元刚度矩阵k(2×2)%---------------------------------------k=[E*A/L -E*A/L; -E*A/L E*A/L];%%%%%%%%%%%%%%%%%%%%%%%%%%function z=Bar1D2Node_Assembly(KK,k,i,j)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k，单元的节点编号i、j%输出整体刚度矩阵KK%-----------------------------------DOF(1)=i;DOF(2)=j;for n1=1:2for n2=1:2KK(DOF(n1), DOF(n2))= KK(DOF(n1), DOF(n2))+k(n1, n2);endendz=KK;%------------------------------------------------------------function stress=Bar1D2Node_Stress(k, u, A)%该函数计算单元的应力%输入单元刚度矩阵k, 单元的位移列阵u(2×1)%输入横截面积A计算单元应力矢量%输出单元应力stress%-----------------------------------stress=k*u/A;%-----------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%function forces=Bar1D2Node_Force(k, u)%该函数计算单元节点力矢量%输入单元刚度矩阵k和单元的位移列阵u(2×1)%输出2×1的单元节点力分量forces%-----------------------------------------forces=k*u;%%%%%%%%%%% Bar1D2Node %% end %%%%%%%%%【四杆桁架结构的有限元分析—数学推导】如图所示的结构，各杆的弹性模量和横截面积都为E=29.54×10N/mm2, A=100mm 2，试求解该结构的节点位移、单元应力以及支反力。

图1 四杆桁架结构解答：对该问题进行有限元分析的过程如下。

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

表1 节点及坐标表2 单元编号及对应节点表3 各单元的长度及轴线方向余弦节点x y 单元节点1 节点2 单元l x n y n1 0 0 ① 12 ①400 1 02 400 0 ②3 2 ②300 0 -13 400 300 ③ 1 3 ③500 0.8 0.64 0 300 ④ 4 3 ④400 1 0（2）各个单元的矩阵描述由于所分析的结构包括有斜杆，所以必须在总体坐标下对节点位移进行表达，所推导的单元刚度矩阵也要进行变换，各单元经坐标变换后的刚度矩阵如下。

（3）建立整体刚度方程将所得到的各个单元刚度矩阵按节点编号进行组装，可以形成整体刚度矩阵，同时将所有节点载荷也进行组装。

刚度矩阵：K= K(1) +K(2)+K(3)+K(4)节点位移：q= [u1v1u2v2u3v3u4v4]T节点力：P=R+F=[ R x1R y12×104R y20 2.5×104 R x4R y4]T其中(R x1, R y1)为节点1处沿x和y方向的支反力，R y2为节点2处y方向的支反力，(R x4, R y4) 为节点4处沿x和y方向的支反力。

整体刚度方程为(4) 边界条件的处理及刚度方程求解边界条件BC(u)为：u1=v1=v2=u4=v4=0，代入整体刚度方程中，经化简后有对该方程进行求解，有则所有的节点位移为(5) 各单元应力的计算其中T为坐标转换矩阵；同理，可求出其它单元的应力。

(6) 支反力的计算将节点位移的结果代入整体刚度方程中，可求出2D杆单元的有限元分析程序(Bar2D2Node)编写平面桁架单元的单元刚度矩阵、单元组装、单元应力的计算程序。

编写的平面桁架单元的四个MATLAB 函数如下。

Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)该函数计算单元的刚度矩阵，输入弹性模量E，横截面积A，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）和角度alpha（单位是度），输出单元刚度矩阵k(4×4)。

Bar2D2Node_Assembly(KK,k,i,j)该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j，输出整体刚度矩阵KK。

Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)该函数计算单元的应力，输入弹性模量E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha （单位是度）和单位节点位移矢量u，返回单元应力标量。

Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)该函数计算单元的应力，输入弹性模量E，横截面积A，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）和单元节点位移矢量u，返回单元节点力。

基于2D杆单元的基本公式，可以编写出具体实现以上每个函数的MATLAB程序如下。

%%%%%%%%%%% Bar2D2Node %% begin %%%%%%%%%%%%%%function k=Bar2D2Node_Stiffness(E, A, x1,y1, x2, y2, alpha)%该函数计算单元的刚度矩阵%输入弹性模量E，横截面积A%输入第一个节点坐标（x1, y1），第二个节点坐标（x2, y2），角度alpha（单位是度）%输出单元刚度矩阵k(4×4)%-------------------------------------------------L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=alpha*pi/180;C=cos(x);S=sin(x);k=E*A/L*[C*C C*S -C*C -C*S;C*S S*S -C*S -S*S;-C*C -C*S C*C C*S;-C*S -S*S C*S S*S]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function z = Bar2D2Node_Assembly(KK, k, i, j)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k，单元的节点编号i、j%输出整体刚度矩阵KK%--------------------------------------------------------DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;for n1=1:4for n2=1:4KK(DOF(n1), DOF(n2))= KK(DOF(n1), DOF(n2))+k(n1, n2);endendz=KK;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%--function stress= Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)%该函数计算单元的应力%输入弹性模量E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）%输入角度alpha（单位是度）和单位节点位移矢量u%返回单元应力标量stress%------------------------------------------------L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=alpha*pi/180;C=cos(x);S=sin(x);stress=E/L*[-C -S C S]*u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function forces= Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)%该函数计算单元的应力%输入弹性模量E，横截面积A%输入第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）%输入单元节点位移矢量u%返回单元节点力forces%-------------------------------------------------------------L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=alpha*pi/180;C=cos(x);S=sin(x);forces= E*A/L*[-C -S C S]*u;%%%%%%%%%%% Bar2D2Node %% end %%%%%%%%%%%%%%【四杆桁架结构的有限元分析—MATLAB—(Bar2D2Node)】仍就图1所示结构，基于MATLAB平台求解该结构的节点位移、单元应力以及支反力。