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:2
for n2=1:2
收集于网络，如有侵权请联系管理员删除
KK(DOF(n1), DOF(n2))= KK(DOF(n1), DOF(n2))+k(n1, n2);
end
end
z=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 n
1 0 0 ① 1
2 ①400 1 0
2 400 0 ②
3 2 ②300 0 -1
3 400 300 ③ 1 3 ③500 0.8 0.6
4 0 300 ④ 4 3 ④400 1 0
（2）各个单元的矩阵描述
由于所分析的结构包括有斜杆，所以必须在总体坐标下对节点位移进行表达，所推导的单元刚度矩阵也要进行变换，各单元经坐标变换后的刚度矩阵如下。
收集于网络，如有侵权请联系管理员删除
收集于网络，如有侵权请联系管理员删除
（3）建立整体刚度方程
将所得到的各个单元刚度矩阵按节点编号进行组装，可以形成整体刚度矩阵，同时将所有节点载荷也进行组装。
刚度矩阵： K = K (1) +K (2)+K (3)+K (4) 节点位移：q = [u 1 v 1 u 2 v 2 u 3 v 3 u 4 v 4]T
节点力： P=R+F=[ R x 1 R y 1 2×104 R y 2 0 2.5×104 R x 4 R y 4]T
其中(R x 1 , R y 1)为节点1处沿x 和y 方向的支反力，R y 2为节点2处y 方向的支反力，(R x 4 , R y 4) 为节点4处沿x 和y 方向的支反力。 整体刚度方程为
