《计算力学》课程大作业
小组成员：成 曦（080681） 宋世伟（080683） 张璘琳（080684） 冯东明（080745）
计算力学——课程设计
目录
一、任务概要 ..................................... 3 二、程序说明 ..................................... 4 三、程序流程图 ................................... 5 四、程序源代码 .................................. 12 五、计算结果比较与分析 .......................... 13 六、程序优化 .................................... 27 七、存在问题 .................................... 30 八、编程体会 .................................... 30 九、各人工作量 .................................. 31
第 2 页 共 31 页
计算力学——课程设计
一、任务概要
图1 结构概图
本次《计算力学》课程大作业，要求我们利用所学有限元知识，编制程序分析图1所示的带斜梁的平面刚架结构。
为了便于分析，我们首先规定了结构的部分尺寸与构件截面的相关特征，详见图2。并对结构作了适当简化，将剪力墙上部变平，这样可能与实际情况更于接近。在程序优化部分中，我们添加了平面四节点的等参元程序，可以解决如图1所示的剪力墙上部为斜线的情况。
图2 平面框架尺寸布置以及节点与单元编号
梁、柱我们规定了统一尺寸，均为0.3m?0.6mm,A?0.18m2,IZ?0.054m4。
第 3 页 共 31 页
计算力学——课程设计
二、程序说明
程序功能：实现带斜梁的平面刚架结构的静力分析。 使用语言：C++。
数据输入：通过C++语言中的输入输出流实现原始数据的输入，即通过框架、剪力墙与耦合模块中的数据输入文件input.txt来实现结构的节点、单元、材料、截面与荷载等基本数据的输入。
数据输出：在C++运行过程中，我们同样通过输入输出流来实现计算数据的输出，即在框架、剪力墙与耦合模块中的数据输出文件output.txt中写入节点的力与位移等相关参量。
对于该有限元程序，基于完全原创的理念，结合计算力学相关知识，拟定的具体编制思路如下：
1. 首先，从总体上把握程序的模块。访程序可分为框架模块，剪力墙模块，耦合模块和其他的模块（如：四节点单元等参元、荷载等效）等等；
2. 然后，确定了各个模块之间的关系，剪力墙与框架平行，然后共同派生出耦合模块； 3. 具体到各个模块的编制：考虑到对此程序以前没有接触，所以笔者从一个解题者的角度出发，思考具体的方案乃至细节；
4. 细节上具体有许多，如何将它们统筹规划，并协调合理的呈现出来是关键，而这其中的要点就是对于数据的处理，数据是否合理成为了我们修改程序的核心；
5. 思考合理的数据输入，要使其既能够尽量简单，又能够便于处理和程序的编制，还要能提供足够的计算需要的数据；
6. 着手解决每一个模块，从剪力墙模块开始，关于其自有变量的定义，这很难一开始就定义完整，而是在以后编制的过程中，将它们前后对应，补充修改，最后才能最终确定，有输入数据函数inputdata()，获取单元刚度矩阵函数getKelement()，获取总刚矩阵函数getKwall()，获取结构荷载向量函数getload()，计算处理函数caculate()和打印函数print()等等；
7. 采用大体上同样的思路解决了框架模块的编制，但其中也有其特有的部分，如数据的输入并不是完全相同，框架单元有角度特征，而剪力墙单元没有。剪力墙单元刚度矩阵是8×8的矩阵，而框架是6×6的矩阵，而且荷载的具体含义也不是完全相同，但是理念上是相同的；
8. 最后在完成了两个基本模块的编制和调试后，开始着手编制耦合模块，具体有如下思路：首先共同划分并编号，然后对于相交点将其自由度统一取为3，将原来剪力墙中的单刚加0扩展，最后得到总刚和总的荷载列向量，然和联合求解，最后通过单刚求出每一个节点上的力和位移，同时对于公共点，都采用了应力均化的思想，使得结果精度更高。
9. 对于改程序的具体细节的优化，在初步工作完成后，我们想到了对于改程序的某些细节的优化，比如：存储方式，计算求解线性方程组的方法，截面的优化，和结果的输出优化等等。
10. 完成编制，得出结果，并与当下流行的有限元软件的计算结果进行比较。最后思考并对于结果的偏差进行合理的猜测！
第 4 页 共 31 页
计算力学——课程设计
三、程序流程图
（1） 框架模块的具体流程图如图3所示。 Frame类由以下几部分组成： 1、 成员变量：
const int Nelement=10; //定义单元数； const int Nnode=9; //定义节点数； double gp[Nelement][2]; //定义单元包含的节点； double element[Nelement][4]; //定义单元属性(长度，角度，面积，惯性矩)； double E; //定义弹性模量；
int dof[Nelement][6][2]; //获取自由度的整体编码和局部编码； double gk[Nelement][6][6]; //得到单元刚度的数组； double Kframe[Nnode*3][Nnode*3]; //得到总体刚度的数组； double gl[Nnode][3]; //获得节点荷载列阵；
double loadframe[Nnode*3]; //获得框架的荷载列阵(单位为kN,m)； double displacement[Nnode*3]; //位移列阵(单位为m)；
double Fpoint[Nnode][3]; //算出节点上的力(单位为kN,m) ； 2、 成员函数
①void inputdata() //用于输入相关数据，包括节点编号、坐标与单元属性；
②void getKelement() //获取单元刚度矩阵函数；
这里我们采用的是基于经典梁单元理论的刚度矩阵，程序中输入的是一个6?6的二维数组。
③void KT(int a,double b) //单刚从单元坐标向总体坐标转换的函数，第一个是转置矩阵，第二个是角度，其中a-单元编号，b-杆件的弧度；
之前我们定义的单元刚度矩阵是建立在局部坐标系基础之上的，我们需要将其转换成总体坐标下。由计算力学知识我们可知：[K]?[T][K][T]，故在该函数中我们首先由初始数据中杆件的弧度来求得坐标转换矩阵[T]，再由矩阵的转置与乘法来实现单元刚度矩阵由局部坐标系到总体坐标系下的转换。
④ void getKframe() //单元刚度合成整体刚度矩阵函数；
在输入数据中，我们定义了节点编码与自由度。其中自由度是按一定顺序输入的，故其本质就是每个节点三个自由度的定位向量。在定位向量的基础上，我们可以进行总体刚度矩阵的合成。具体思路是：每个单元刚度矩阵中的元素的行标与列标均有具体的定位向量与其对应。在合成总刚的时候，总体刚度矩阵中的行标与列阵就是定位向量，然后我们利用for循环在各个单元刚度矩阵中寻找与某定位向量对应的元素，然后相加得到最终的总刚。
关键函数如下： for(i=0;i 第 5 页 共 31 页 Te