基于matlab的二维杆单元开发
目录
概述 杆单元为什么要引入局部坐标 局部坐标系下的刚度矩阵 全局坐标系下的刚度矩阵 程序编写&算例 所有的计算代码
概述
杆系结构是工程中应用较为广泛的结构体系，包含了平面或空间形式的梁、钢架和行桁架等等，这些结构的每个单独的杆件都可以视为一个单独的单元，单元与单元之间用节点连接，根据结构体系的不同，节点可以传递轴向力、切向力或者弯矩等等。
帖子以二维平面杆单元为例，采用材料力学更直观的方法建立杆单元的刚度矩阵，详细介绍了具体的推导过程，最后以一个桁架为例讲解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
