1. 一维杆单元有限元分析程序
一维刚单元的局部坐标系（单元坐标系）与全局坐标系相同。
1.1 线性杆单元
如图所示是一个杆单元，由两个节点i和j，局部坐标系的X轴沿着杆的方向，由i节点指向j节点，每个节点有一个自由度，在使用时需要定义杆长L、截面积A和弹性模型E。
单元刚度矩阵为：
线性杆单元有限元分析MATLAB程序主要包括单元刚度矩阵（Bar1D2Node_Stiffness）、单元组装(Bar1D2Node_Assembly)、单元应力（Bar1D2Node_Stress）和节点力（Bar1D2Node_Forces）四个函数，分别如下：
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 y = Bar1D2Node_Assemble(K,k,i,j)
% 该函数进行单元刚度矩阵的组装
% 输入：单元刚度矩阵k和单元节点编号i，j
% 输出：整体刚度矩阵K
DOF(1) = i;
DOF(2) = j;
for n1 = 1:2for n2 = 1:2K(DOF(n1),DOF(n2)) = K(DOF(n1),DOF(n2))+k(n1,n2);end
end
y = K;function forces = Bar1D2Node_Forces(k,u)
% 该函数计算线性杆单元的节点力
% 输入：单元刚度矩阵k，单元的位移列阵u（2*1）
% 输出：单元节点力forces
forces = k*u;function stress = Bar1D2Node_Stress(k,u,A)
% 该函数计算线性杆单元的应力
% 输入：单元刚度矩阵k，单元的位移列阵u（2*1），单元横截面积A
% 输出：单元应力stress
stress = k * u/A;1.2 二阶杆单元
二阶杆单元是带中间节点的杆单元，如下图所示，单元节点为i,j,m，杆单元界面为A，弹性模量为E，坐标系为i指向j的x方向。
单元刚度矩阵为：
二阶杆单元有限元分析MATLAB程序主要包括单元刚度矩阵（Bar1D3Node_Stiffness）、单元组装(Bar1D3Node_Assembly)、单元应力（Bar1D3Node_Stress）和节点力（Bar1D3Node_Forces）四个函数，分别如下：
function k = Bar1D3Node_Stiffness(E,A,L)
% 该函数计算二阶杆单元的刚度矩阵
% 输入：弹性模量E，杆的横截面积A和长度L
% 输出：单元刚度矩阵k（3*3）
k = E*A/(3*L)*[7 1 -8;1 7 -8;-8 -8 16];function y = Bar1D3Node_Assemble(K,k,i,j,m)
% 该函数进行二阶杆单元刚度矩阵的组装
% 输入：单元刚度矩阵k和单元节点编号i，j, m
% 输出：整体刚度矩阵K
DOF(1) = i;
DOF(2) = j;
DOF(3) = m;
for n1 = 1:3for n2 = 1:3K(DOF(n1),DOF(n2)) = K(DOF(n1),DOF(n2))+k(n1,n2);end
end
y = K;function forces = Bar1D3Node_Forces(k,u)
% 该函数计算二阶杆单元的节点力
% 输入：单元刚度矩阵k，单元的位移列阵u（3*1）
% 输出：单元节点力forces
forces = k*u;function stress = Bar1D3Node_Stress(k,u,A)
% 该函数计算二阶杆单元的应力
% 输入：单元刚度矩阵k，单元的位移列阵u（3*1），单元横截面积A
% 输出：单元应力stress
stress = k * u/A;2. 二维杆单元有限元分析程序
未完待续....
3. 三维杆单元有限元分析程序
参考文献
曾攀《有限元基础教程》；
