平面/空间杆系结构有限元编程计算（MATLAB）

GitHub链接：MATLAB源码
博主QQ：915339719，有问题可以随时交流

程序简介

本程序可以对绝大多数的平面/空间杆系结构进行静力学计算、自振频率计算，以及平面杆系结构的稳定性计算，具有通用性

程序m文件

程序有6个.m文件：1、Link.m（主文件，包含静力计算/稳定性计算/自振频率计算） 2、 eLength.m（计算单元长度） 3、Transformation.m（计算静力问题的坐标转换阵） 4、stabTransformation.m（计算稳定性问题的坐标转换阵） 5、Cholesky.m（乔里斯基分解计算全局节点位移） 6、Delete.m（约束处理-划行划列方法）

程序变量

程序里面的代码都有详细的注释，如何实现不作赘述，这里对程序需要输入参数的变量进行说明，了解了如何输入参数，程序就完全可以当做黑盒子使用。

程序定义的大部分变量如下图所示，这里只需要注意6个变量：

1、AE（单元特征数）：存储材料信息。第一行存储弹性模量，第二行存储截面面积，第三行为质量密度，一列代表一种材料号，比如图中就只有1种材料，弹性模量为200∗109200*10^9200∗109，截面面积为10−410^{-4}10−4，质量密度为10310^3103。

2、NMN（单元特征数类信息）：存储材料编号。列号就是单元号(第一列就代表1号单元的信息)，里面填该单元的材料编号，计算时按照材料编号取对应编号下的弹性模量、截面面积和质量密度。图中代表所有单元都是1号材料，即它们的弹性模量均取200∗109200*10^9200∗109，截截面面积均取10−410^{-4}10−4，质量密度均取10310^3103。

3、ME（单元信息）：存储单元的节点编号。列号就是单元号(第一列就代表1号单元的信息)，里面填该单元的两个节点编号(杆单元只有两个节点)，第一行为第一个节点，第二行为第二个节点，顺序任意。

4、Coordinates（节点坐标）：存储节点坐标。这里行号是节点号(第一行代表1号节点的坐标信息)，平面问题填两个坐标即可，空间问题填三个坐标，坐标值由自己建立坐标系确定。

5、NRR（节点约束信息）：存储节点约束信息。列号是节点号，第一、二、三行分别依次存储x、y、z方向的约束。约定有该方向的约束就填1。平面问题只需要填x、y方向的约束信息。

6、P（节点载荷信息）：存储节点载荷值。列号是节点号，第一、二、三行分别依次存储x、y、z方向的载荷值，按照实际值填写，按照自己建立的坐标系规定正负。
这里只能填集中力载荷，有些杆系问题会有均布载荷和弯矩的情况，可以先自行把均布载荷转化为集中力、把弯矩转化为等大反向的平行力作用在节点上再写入P中即可。

细节说明

稳定性计算只能针对结构只有一个外力载荷的情况，它计算的就是该力在该方向下的临界值。它的理论方法就是在某个节点上施加竖直或者水平方向的单位力，得到单位力下的杆内力，代入初应力矩阵（几何刚度阵）进行计算，具体的不深入介绍。程序计算稳定性问题要借助静力计算时施加的外载，通过线性关系可以计算出单位力下的杆内力。也就是说外载的方向就是单位力的方向，计算出来的临界力就是外载方向的临界力。因此每次计算稳定性问题时下图方框内要改成外力载荷坐标号。
建议不进行稳定性计算或者自振频率计算时把整个模块的代码都注释掉，以免出错。

案例分析

给了一个桁架的案例进行说明。

题目：平面桁架的支座和载荷如下图所示，F=1KNF=1KNF=1KN，a=1ma=1ma=1m，E=200∗109E=200*10^9E=200∗109，A=1cm2A=1cm^2A=1cm2，试求支座反力，所有杆的内力，求FFF的临界力大小，求桁架结构的固有频率（自振频率）。

先给杆和节点任意编号，这里作者取的杆编号是：
CD−1CD-1CD−1 BC−2BC-2BC−2 EF−3EF-3EF−3 AD−4AD-4AD−4 AB−5AB-5AB−5 AE−6AE-6AE−6 BE−7BE-7BE−7 DF−8DF-8DF−8 CF−9CF-9CF−9
给A−FA-FA−F六个点按顺序编号1−61-61−6
建立坐标系在AAA点（图中未画出），ABABAB方向为x轴，ADADAD方向为y轴。写出坐标
A(0,0)A(0,0)A(0,0)，B(1,0)B(1,0)B(1,0)，C(1,1)C(1,1)C(1,1)，D(0,1)D(0,1)D(0,1)

现在填写参数：

改为添加的外载

直接求解

节点1x方向约束反力计算精度问题，视为0即可

经过简单的理论计算，可以知道程序解答与理论解答是一致的（正负代表杆的拉压状态）

临界力一般相差比较大，命令行窗口可能会显示0，打开变量直接看值即可，这里不过多解释结果，只说下取值方法：正数取最小值，负数取绝对值最小的那个值。这里即取
8.3452e6，-2.6106e7，说明在F点施加竖直方向的力，方向向下的临界值为-2.6106e7，方向向上的临界值为8.3452e6。超过此值就会有杆出现稳定破坏。

自振频率个数与结构自由度个数相等，平面杆单元的节点有两个自由度，此题有6个节点就有12个自由度，但有3个自由度被约束，因此就有9个自由度即9个自振频率。在这里面的最小值就是基频。

补充

作者基于Winform写了一个能够计算杆系结构和梁系结构的应用程序，除了具有上述求解杆系结构的功能外，还能进行平面/空间梁系结构的静力计算和平面稳定性计算（压杆稳定问题）。

框架参照了ANSYS APDL，模型基于开源的ANYCAD三维控件，缺点是不能对模型进行操作，只能输入点的数据，然后给线段编号，最后绘制线段图形。不过已经对于杆梁计算已经足够。

现在主体功能已经开发完毕，若有兴趣，可以联系作者，相互交流。