概述

杆系结构是工程中应用较为广泛的结构体系，包含了平面或空间形式的梁、钢架和行桁架等等，这些结构的每个单独的杆件都可以视为一个单独的单元，单元与单元之间用节点连接，根据结构体系的不同，节点可以传递轴向力、切向力或者弯矩等等。
帖子以二维平面杆单元为例，采用材料力学更直观的方法建立杆单元的刚度矩阵，详细介绍了具体的推导过程，最后以一个桁架为例讲解matlab程序并采用abaqus做对比计算，文末给出了详细的程序。

杆单元为什么要引入局部坐标系

我们在做一件事情的时候，一定要知道为什么做。
杆单元为什么不直接在全局坐标系中求解，非得人为的引入一个局部坐标系？我想第一次学习杆单元的人都会有这个疑问，然后很遗憾，单单用杆单元部分的知识并不能很好的回答这个问题。一个简单的答案是：为了统一单元、最终是为了简化计算。
我们知道，科学家最喜欢将特殊的东西泛化，然后给出一个简洁美丽的公式来解释广泛的事物。在科学计算中也是如此，就拿这次的杆单元来说，我们脑子里面浮现的是一个结构体系中横七竖八的段杆，就像是体育场上上空的架子，乱糟糟的，用数学的语言描述就是他们的空间位置都不相同，这在直觉上意味着极大的工作量。
这时候，天马行空的想法随着工作量的增加会自然而然的从脑子里面冒出来（很多发明一开始都是天马行空的偷懒想法）：如果我们能只求一次某个量，然后乘以一个系数把这个量映射到所有的杆单元上面就好了，这个想法就是引入局部坐标系的火苗，有了这个思想的指导，后面的数学家就开始干活了，最终的杆单元计算过程就是上面那句加粗的文字那样：首先在局部坐标系下计算刚度矩阵，然后乘以一个坐标转换矩阵，将据局部坐标系下的刚度矩阵映射到全局坐标系中参与组装。下面列出一些引入局部坐标系具体方便了哪些

简化计算
局部坐标系可以简化杆单元的几何和材料属性的描述，使得在进行刚度矩阵和质量矩阵的计算时更为方便。通过将杆单元的方向与局部坐标系对齐，可以使得计算过程中的数学表达更为直观。
方便载荷施加
在局部坐标系中，可以更方便地定义和施加载荷。比如，集中载荷或分布载荷可以直接在杆单元的局部坐标系中施加，避免了转换成全球坐标系的复杂性。
更好地处理变形和应力分析
杆单元在变形时会有不同方向的应力和应变，局部坐标系可以帮助明确杆单元的拉伸、压缩及剪切行为，进而提高计算精度。
支持多方向的连接
在复杂的结构中，杆单元可能以不同的方向连接。在局部坐标系中，每个单元的方向都是独立定义的，这使得处理复杂结构的连接关系更为灵活。
易于与其他单元类型结合
在实际工程应用中，结构通常由多种单元组成。引入局部坐标系可以更方便地将杆单元与其他类型的单元（如面单元或体单元）结合，确保它们之间的力学性质和边界条件协调一致。
我在这里只是从定性的角度给出了一些解释，要想彻底的从底层理解这一套工作思路，还需要学习等参元，等参元的引入把简化计算的思想发挥的淋漓尽致，在等参元中，科学家通过坐标变换解决了积分难题，

局部坐标系下的刚度矩阵

材料力学给出了二维杆单元的受力特点

每个杆单元的节点只有两个平动自由度
二维杆单元的节点不传递弯矩
外荷载均为作用在节点上的集中力
基于此，我们从一个结构体系中取出一个杆件，并引入与杆轴向重合的局部坐标系，设杆单元的截面积为，长度为，两端的节点分别用和表示，并且从节点到节点的方向为局部有坐标轴，垂直于杆轴向为，具体的示意图为根据材料力学的知识，对两端节点列方程

整理成矩阵形式

以上就是二维杆单元在局部坐标系中的刚度矩阵，下面推导局部坐标系与全局坐标系的转化矩阵，即全局坐标系的刚度矩阵。

全局坐标系下的刚度矩阵

推导全局坐标系下的刚度矩阵的重点就是推导转换矩阵，即轴向量与全局坐标系之间的转换关系，如下图展示了轴向数据与全局坐标系数据通过上面图片可知，整体坐标系下，杆单元两端节点的位移矩阵为

局部坐标系下，杆单元两端节点的位移矩阵为

根据三角函数关系，写出全局坐标系位移与局部坐标系位移之间的关系为

整理成矩阵形式

其中，转换矩阵为

其中

上面就是局部坐标系向全局坐标系转换数据的矩阵，局部和全局量尊遵循上述转换关系。对于全局坐标系下的刚度矩阵，则有

结合上文中的局部坐标系下的刚度矩阵，代入转换矩阵有

好了，现在有了全局坐标系下的杆单元刚度矩阵，下面就可以写程序了。

程序编写&算例

这里直接给出一个算例，以算例的形式详细讲解程序。二维杆单元的弹性模量，截面面积，具体的网格和边界条件示意图为上面的桁架一共15个杆单元，单元编号用斜体字表示；9个节点，节点用有背景色的方块表示，在节点1和6上施加竖直向的集中力，力的大小，节点9和4为的两个平动自由度被限制。
下面开始介绍关键的代码，全部的代码在文末给出。
首先是杆单元的材料属性，杆单元只需要弹性模量和截面积。

E=210e3; % 弹性模量 (GPa=1000 N/mm^2)
A=500; % 杆单元截面积 (mm^2)

然后给出模型的信息。

node=9; % 模型节点总数
element=15; % 杆单元的数目
nodecord=[  1.5,  0.866025388;
  2.5,  0.866025388;
   2.,           0.;
   3.,           0.;
   1.,           0.;
  0.5,  0.866025388;
   0.,           0.;
 -0.5,  0.866025388;
  -1.,           0.];
xcord=nodecord(:,1); % 模型的x坐标
ycord=nodecord(:,2); % 模型的y坐标
elementcon=[1, 2   ; 
2, 3    ;
3, 4    ;
4, 2    ;
1, 5    ;
6, 1    ;
6, 7    ;
7, 5    ;
 3, 1   ;
 5, 3   ;
 5, 6   ;
 8, 6   ;
 7, 8   ;
 9, 7   ;
 8, 9   ];
Wd =2*node; % 模型的整体自由度

下面计算杆单元刚度矩阵和施加节点力。计算杆单元的全部代码见文末。

% 计算杆单元刚度矩阵
[S]=planetruss(E,A, Wd,element,elementcon,xcord,ycord);
% 施加荷载 (units: N)
F=zeros(Wd,1);
% 给第二个自由度施加45e3N的力，即给节点1的y方向施加45e3N的力
F(2)=-45e3; 
% 给第十二个自由度施加45e3N的力，即给节点6的y方向施加45e3N的力
F(12)=-45e3;

下面给节点4和9施加固定边界条件。

%% 施加固定位移条件
gdof=1: Wd;
cdof=[7 8 17 18];

求解方程，这里是静力计算，实质上的计算量就是矩阵求逆。

[U]=solve(Wd,gdof,cdof,S,F);
displacements=[gdof' U]

下面是一些后处理的程序，以注释的形式给出解释说明。

%% 统计变形后的节点坐标
ux=1:2:2*node-1;
uy=2:2:2*node;
UX=U(ux);
UY=U(uy);
newcordx=xcord+UX;
newcordy=ycord+UY;
newcord=[newcordx,newcordy]
%% 计算杆单元的力
[p]=trussforce(E,A,element,elementcon,xcord,ycord,U);
A=1:size(p,1);
Elemental_Force=[A' p]
%% 绘制杆单元
SC=100; %缩放因子，便以观察变形
Snewcordx=xcord+SC*UX;
Snewcordy=ycord+SC*UY;
Snewcord=[Snewcordx,Snewcordy];
drawplanetruss(element,elementcon,nodecord,newcord,Snewcord);

运行上面的matlab程序，并将变形前后的桁架绘制在同一个窗口中观察变形，图形为可以看到，变形的形状符合直觉，为了进一步验证程序的计算效果，采用abaqus创建同样的模型并对比计算结果，abaqus的模型为 abaqus计算完毕后，统计所有节点的竖直向位移与matlab程序对比，绘制表格为

matlab	abaqus
-0.003214286	-0.00321429
-0.001071429	-0.00107143
-0.002285714	-0.00228571
0	-4.50E-32
-0.003428571	-0.00342857
-0.003214286	-0.00321429
-0.002285714	-0.00228571
-0.001071429	-0.00107143
0	-4.50E-32

可以看到，上面的计算结果一模一样，跟造数据了一样，事实上就算是我造数据了，你们也不知道，但是按我的习惯，我必然会把所的计算文件放在文末。

所有的计算代码

下面是matlab的计算主程序，主程序调用的函数我就不放在这里了，有需要的可以私聊我，私聊必发！

clc ; clear ; % 清空变量
%%
% 杆单元材料参数
E=210e3; % 弹性模量 (GPa=1000 N/mm^2)
A=500; %杆单元截面积 (mm^2)
node=9; % 模型节点总数
element=15; % 模型单元总数
nodecord=[  1.5,  0.866025388;
  2.5,  0.866025388;
   2.,           0.;
   3.,           0.;
   1.,           0.;
  0.5,  0.866025388;
   0.,           0.;
 -0.5,  0.866025388;
  -1.,           0.];
xcord=nodecord(:,1); ycord=nodecord(:,2);
elementcon=[1, 2   ; 
2, 3    ;
3, 4    ;
4, 2    ;
1, 5    ;
6, 1    ;
6, 7    ;
7, 5    ;
 3, 1   ;
 5, 3   ;
 5, 6   ;
 8, 6   ;
 7, 8   ;
 9, 7   ;
 8, 9   ];
Wd =2*node; % Whole DOFs
% 计算杆单元刚度矩阵
[S]=planetruss(E,A, Wd,element,elementcon,xcord,ycord);
% 施加节点力 (units: N)
F=zeros(Wd,1);
F(2)=-45e3;
F(12)=-45e3;
%% 施加固定位移边界条件
gdof=1: Wd;
cdof=[7 8 17 18];
%% 求解位移
[U]=solve(Wd,gdof,cdof,S,F);
displacements=[gdof' U]
%% 统计变形后的坐标
ux=1:2:2*node-1;
uy=2:2:2*node;
UX=U(ux);
UY=U(uy);
newcordx=xcord+UX;
newcordy=ycord+UY;
newcord=[newcordx,newcordy]
%% 计算单元力
[p]=trussforce(E,A,element,elementcon,xcord,ycord,U);
A=1:size(p,1);
Elemental_Force=[A' p]
%% 绘制变形前后的对比图
SC=100; % 缩放因子，便于观察
Snewcordx=xcord+SC*UX;
Snewcordy=ycord+SC*UY;
Snewcord=[Snewcordx,Snewcordy];
drawplanetruss(element,elementcon,nodecord,newcord,Snewcord);

下面是abaqus的计算inp文件，这个inp文件拿过来就能直接算，算完的结果跟上面的matlab一样，曲线吻合的跟假的一样，非常能提升自己搞科研的自信心，建议保存起来，每当其他程序出bug的时候运行一下这个程序，以此欺骗和麻痹自己。

*Heading
** Job name: Job-1 Model name: Model-1
** Generated by: Abaqus/CAE 6.14-2
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*Node
      1,          1.5,  0.866025388
      2,          2.5,  0.866025388
      3,           2.,           0.
      4,           3.,           0.
      5,           1.,           0.
      6,          0.5,  0.866025388
      7,           0.,           0.
      8,         -0.5,  0.866025388
      9,          -1.,           0.
*Element, type=T2D2
1, 1, 2
2, 2, 3
3, 3, 4
4, 4, 2
5, 1, 5
6, 6, 1
7, 6, 7
8, 7, 5
 9, 3, 1
10, 5, 3
11, 5, 6
12, 8, 6
13, 7, 8
14, 9, 7
15, 8, 9
*Nset, nset=Set-1
 1, 5, 6, 7, 8
*Elset, elset=Set-1
  6,  7,  8, 11, 13
*Nset, nset=Set-6, generate
 1,  9,  1
*Elset, elset=Set-6, generate
  1,  15,   1
*Nset, nset=fixed
 4, 9
*Nset, nset=middle
 1, 6
** Section: Section-1
*Solid Section, elset=Set-6, material=Material-1
500.,
*End Part
**  
**
** ASSEMBLY
**
*Assembly, name=Assembly
**  
*Instance, name=Part-1-1, part=Part-1
*End Instance
**  
*End Assembly
** 
** MATERIALS
** 
*Material, name=Material-1
*Elastic
210000.,0.
** 
** BOUNDARY CONDITIONS
** 
** Name: fixed Type: Displacement/Rotation
*Boundary
Part-1-1.fixed, 1, 1
Part-1-1.fixed, 2, 2
** ----------------------------------------------------------------
** 
** STEP: Step-1
** 
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
** 
** LOADS
** 
** Name: middle   Type: Concentrated force
*Cload
Part-1-1.middle, 2, -45000.
** 
** OUTPUT REQUESTS
** 
*Restart, write, frequency=0
** 
** FIELD OUTPUT: F-Output-1
** 
*Output, field, variable=PRESELECT
** 
** HISTORY OUTPUT: H-Output-1
** 
*Output, history, variable=PRESELECT
*End Step

来源：有限元先生

目录

概述

杆单元为什么要引入局部坐标系

局部坐标系下的刚度矩阵

全局坐标系下的刚度矩阵

程序编写&算例

所有的计算代码