计算力学——课程设计
表12 SAP2000结果比较（节点力） 节点力比较（剪力墙部分） X Y X Y 节点编号 13（1） 15（2） C++运算结果 389.744 1378.31 326.848 875.989 ANSYS分析结果 443.978 2155.506 421.24 2028.65 绝对误差 54.23 777.20 94.39 1152.66 相对误差 12.22% 36.06% 22.41% 56.82% 注：节点编号中括号内的数字为SAP2000分析中的节点编号，括号外的数字为C++计算中的节点编号。
可能是由于ANSYS与SAP2000中的算法存在差异，两种电算结果本身就存在差距。但SAP2000结果与C++结果比较之后，结果变化趋势大致相同；此外，我们可以发现节点位移差距稍大，节点力差距稍小。但均在同一数量级之上。
从C++结果与ANSYS 10、SAP2000两种电算结果比较中，我们可以发现C++结果从变化趋势上是比较合理的，具体数值差距也不是很大，绝大多数均在一个数量级上。这充分说明了我们编制的C++程序从算法上来说是比较合理的，具有一定可行性。
其次，误差分析：
反观所编C++程序，经初步分析，我们认为两种计算结果产生差距的主要原因有以下几点：
1、 输入数据中存在近似数
在框架模块的输入数据中，存在梁单元的角度，即element[Nelement][1]，在初始数据中我们采用了弧度的近似值而非精确值，从而导致坐标转换矩阵中相关元素不准确，故直接会影响到总刚的集成与方程的求解。 2、 单元刚度矩阵存在区别
在C++源程序编制中，梁单元的单元刚度矩阵采用的基于经典梁单元的单元刚度矩阵，平面四节点单元的单元刚度矩阵则是我们通过自己推导而得。故C++程序中的单元刚度矩阵可能与ANSYS中Beam3和Plane 42单元的单元刚度矩阵有所出入，从而导致计算结果存在差距。 3、 线性方程组的解法
在C++源程序编制中，线性方程组的解法采用列主元高斯消去法。这种方法在低自由度时效率较高。而对于自由度较多的复杂问题来说，总体刚度矩阵具有较大的稀疏性，采用列主元高斯消去法会产生较大的截断误差。ANSYS中采用的应该是精确度较高的迭代法。故方程组解法的不同也会造成计算结果存在一定差距。 4、 荷载方向的定义有区别
在框架模块中，计算所得的节点力与ANSYS的分析结果比较，可以发现弯矩存在一个正负号的区别。经分析可知，在C++编制程序中，我们将弯矩按逆时针旋转为正方向；而在ANSYS中，应该与材料力学中规定相同，弯矩按顺时针旋转为正。 5、 编程思路存在区别
我们是以一个初学者的角度去看待问题、分析问题，只能根据所学计算力学知识逐步编写程序；而ANSYS与SAP2000程序则是两种经过多年完善并得到大众认可的程序。故编程思路会存在一定区别。ANSYS与SAP2000中存在许多优化的地方，而我们所编写的程序中未能得以体现；而且ANSYS与SAP2000两种程序中考虑问题应该比我们全面。故结果之间存在一定差距，我们认为是在情理之中的。
第 26 页 共 31 页
计算力学——课程设计
六、程序优化
由于此次是我们第一次尝试较大规模的有限元程序编制，故在程序编写过程中遇到了许多问题，考虑问题也不全面。纵观源程序，我们发现还有许多值得改进的方面。由于时间问题，我们只将以下几部分作为独立程序进行编写，未能其溶入到主程序之中，望尹老师予以原谅。
1、 等参元
在剪力墙模块中，我们使用的是四节点平面单元。起初单元刚度矩阵，我们是通过手算，之后将规则四节点单元的刚度矩阵在程序中加以定义，通用性比较差，且在剪力墙上部是斜向的，网格划分后的四节点单元是不规则的，故须用等参元加以计算。回顾有限元相关知识后，我们认为利用2*2高斯积分，可以精确得出相应单元的刚度矩阵。
高斯积分点可取为(?13,?13)，(13,?13)，(13,13)，(?13,13)。
图12 四节点等参元
我们在优化过程中，定义了一个等参元类，可以用于计算任意四节点等参元的刚度矩阵。 Isoparametric类，由以下几部分组成： （7） 成员变量 double GaussPoint[4][2]; //存储高斯积分点 double Coordinate[4][2]; //存储自然坐标 double J[4][2][2]; //存储Jacobi矩阵 double InversJ[4][2][2]; //存储Jacobi的逆矩阵 double Diff[4][2][4]; //存储四个高斯积分点形函数的微分值 double dN[4][2][4]; //存储自然坐标下的微分值 double B[4][3][8]; //存储各高斯点的B矩阵 double BT[4][8][3]; //存储B的逆矩阵 double D[4][3][3]; //存储刚度矩阵 double BTDB[4][8][8]; //存储各高斯点的刚度矩阵 double K[8][8]; //存储单元的刚度矩阵 （8） 成员函数
① Isoparametric(double integral[4][2]) //integral[4][2]为等参元的整体坐标
在该函数中，定义高斯积分点的坐标。由于我们利用的是2*2高斯积分，故存在四个高斯积分点。在这个程序中，我们可以输入高斯积分点的坐标。
② void getDiff() //获得在正则坐标下的形函数的偏微分矩阵；
第 27 页 共 31 页
计算力学——课程设计
在该函数中，我们对各个节点的形函数对正则坐标求偏微分，求得
?Ni?Ni与（其????中i=1,2,3,4）。
③ void getJ() //获得Jacobi矩阵;
形函数的偏微分函数（Diff矩阵）与高斯积分点的坐标矩阵相乘，便可求得Jacobi
矩阵。
④ void getiverseJ() //获得Jacobi矩阵的逆矩阵；
⑤ void getdN() //获得形函数在整体坐标下的偏微分方程 形函数对整体坐标的偏微分矩阵（
?Ni?Ni与），可由InversJ矩阵与Diff矩阵相?x?y乘而得。
⑥ void getB() //获得B矩阵；
⑦ void getD(double E,double mu) //获取单元的D矩阵；
⑧ void getBTDB() //通过计算获得BTDB矩阵；
在该函数中，由矩阵转置与矩阵乘法可求得BDB矩阵。再通过四节点的BDB矩阵相加而得单元刚度矩阵。
2、 荷载等效
在程序编写过程中，荷载数据我们是通过荷载等效原则，通过手算得到节点荷载，再加 添加到荷载输入文件中。在实际问题中，由于外部与内部荷载相当复杂，故手算过程相当复杂，且忽略了计算机的实效性。
在源程序中，我们列举了以下几种荷载类型，包括集中力与分布荷载。
TT图 几种荷载类型示意图
结合相关结构力学知识，我们按荷载等效原则，将各种荷载类型下梁单元两端的荷载以C++语言的形式写入到源程序中。
其主要程序如下： void equivalent() { double a,b; a=Position/element; b=1-a;
第 28 页 共 31 页
计算力学——课程设计
switch (LoadType) { case 0: gl[int(a)][Direct]+=Value;break; case 1:
gl[0][1]+=Value*(element-Position)*(element-Position)*(1+2*a)/(element*element); gl[1][1]+=Value*Position*Position*(1+2*b)/(element*element); gl[0][2]+=-Value*Position*(element-Position)*(element-Position)/(element*element);
gl[1][2]+=-Value*Position*Position*(element-Position)/(element*element);break; case 2: gl[0][1]+=Value*element/2; gl[1][1]+=Value*element/2; gl[0][2]+=-Value*element*element/12; gl[1][2]+=-Value*element*element/12;break; case 3: gl[0][1]+=-6*Value*a*b/(element*element*element); gl[1][1]+=-gl[0][1]; gl[0][2]+=Value*(element-Position)*(2-3*b)/element; gl[1][2]+=-Value*Position*(2-3*a)/element;break; case 4: gl[0][1]+=3*Value*element/20; gl[1][1]+=7*Value*element/20; gl[0][2]+=-Value*element*element/30; gl[1][2]+=-Value*element*element/20;break; default: break; } }
具体运算程序可见equiload.cpp文件。
3、 一维变带宽存储
框架与剪力墙的总体刚度矩阵，我们在源程序中是按二维数组进行存储。其中二维数组的行数与列数是按所有节点的自由度总和来确定的。由计算力学知识，我们可知总体刚度矩阵具有对称性、稀疏性和带形分布等特点。故总体刚度矩阵中绝大部分元素为零，故按上述方法存储总刚，将浪费大量的存储空间。为达到高效率，我们可以采用一维变带宽存储总刚中非零元素。
我们通过寻找各个自由度的相关自由度来确定总刚中各行元素的半带宽，从而确定各主元在一维变带宽存储的总刚中的位置。
在函数中，我们先寻找每个单元自由度（梁单元6个、平面四节点单元8个）中的
最小自由度，再根据最小自由度来确定单元中各个自由度的半带宽。 具体程序可见：bandai.cpp文件
4、 线性方程组求解方法的优化
在源代码中，我们利用列主元高斯消去法求解线性方程组。由于总体刚度矩阵具有较大的稀疏性，故在求解过程中会存在截断误差与舍入误差，从而造成结果的精确度不高。
相比之下，迭代法适用高自由度、精度要求高的计算模型。特别适合实体单元的分析。
故方程组的解法也是本程序有待进一步完善的方面。
第 29 页 共 31 页
计算力学——课程设计
七、存在问题
在此次编程过程中，我们也发现了很多方面。其中有些问题经过小组讨论与资料查询得到了解决，但还有较多问题仍存疑惑。在此提出，望尹老师可以帮助我们共同解决。
1、 两种单元的耦合方法：
在编制的程序中，我们是将平面四节点单元中每个节点由原来的2个自由度扩充为3个自由度，增辑一个转角自由度?，且将?对应的刚度系数设定为零。然后将框架的总体刚度矩阵与剪力墙的总体刚度矩阵进行合成，对耦合点自由度对应的刚度系数进行叠加，从而得到总体刚度矩阵。这种方法借鉴于计算力学中介绍的薄壳单元耦合方法。不知这种方法是否合适？
且这种方法在耦合点处只能传递Fx与Fy，而不能传递M。 2、 单元刚度矩阵： 在C++程序中，我们对梁单元采用的是经典梁单元的单元刚度矩阵，平面四节点是根据高斯积分得到的单元刚度矩阵。不知这两种刚度矩阵是否与ANSYS中BEAM3与PLANE42两种单元的单元刚度矩阵相对应？
3、 框架位移中角度的误差： 对于框架，C++程序运算结果与ANSYS分析结果比较中，我们发现X方向与Y方向的位移以及节点力差距均不大，唯独角度的差距较大，多数达到90%左右。不知是不是由于单位不同的原因？又为何角度差距如此之大，反而根据位移计算所得的节点力差距均不是很大？这个问题困扰了我们很久。
八、编程体会
在整个编程的过程中，我们有很多的感触，最深刻的就是从简单做起。这个程序虽然不是很难，但是如果想一口气就把它写完，且很完善，没有漏缺还是很不容易的。就比如说其中的一些函数，如耦合模块中总体刚度矩阵的形成函数getKfw（），在没有编制时，是很难把当中的所有细节都思考完整的，这需要你在编制的过程中不停的修改和思考，才能最后得到合适的满意的程序代码。当然，一开始的总体把握是非常重要的，它能够为以后努力确定明确的方向，并且能够更好的思考每一个程序块的编制！相信很多组都会遇到无从下手的问题，我们组开始也遇到了同样的问题，当时没有谁敢保证可以把程序编得很好，但是我们相信只要我们努力学习思考，就一定可以有所突破，我们把问题尽可能的简单化，实在不行，我们把所有的数据都编进去，就是一个一对一的计算器，没有任何通用性可言！对于这样的问题，大家还是很有信心的，于是在大家的共同的努力下，终于完成了剪力墙模块的编制，后来又完成了框架模块的编制，最后耦合完成，优化程序一个个的完成！就这样，大家在不知不觉中，完成了整个作业！所以，就像俗话说的，万事开头难，而开始最难的就是让自己有信心完成，并且找到思路，而我们一开始很好的做到了这些，所以以后就顺利多了，当然还是有一定的难度的。我想，最难的就是程序的调试了。
好不容易完成了程序的编制，这个时候，大家已经很疲倦了，还带着点自满，但是发现程序中又没能运行出比较理想的结果，浮躁开始蔓延！让大家再去修改自己负责的程序时，效率很低，两天过了几乎没有任何进展！于是，大家决定休息一天，好好调整后，第二天让大家交换程序进行修改，这样每个人看到的都是全新的代码，不会有视觉疲劳，也不会有心理上的疲劳，反而更有尽头了，最后终于发现了一些细小的错误，虽然细小，不过是致命的，他会让你的结果谬以千里！我们调试修改程序的策略是：先从简单问题入手，比如框架我们便分析了单层单跨的平面框架。分别从单刚矩阵、坐标转换矩阵以及总体刚度矩阵集成三部
第 30 页 共 31 页
计算力学——课程设计
分入手，这些步骤我们可以通过手算加以推导。经手算结果与C++结果比较，我们发现了问题的所在。改完程序，得到了比较理想的结果输出，大家终于可以放松了，而且都很有感慨，于是写下来，得到了上面这段话，作为对自己多日来辛勤劳动的一种总结吧！
总之，计算力学是一门多学科综合的课程，所要求的掌握内容多、理解层次高，是结构工程专业不可或缺的专业课。通过此次大作业，我们既熟悉了C++编程的相关的知识，更加深了对计算力学这门课程的理解，可谓受益匪浅。
九、各人工作量
成曦：部分调试，部分优化（耦合），模块的框架设计，荷载的处理和单刚到总刚的合成处理和计算矩阵的前处理等等），SAP2000的分析，及部分报告整理； 联系方式：13655165465
宋世伟：部分ANSYS分析；部分调试，部分优化（耦合中部分程序，等参部分，荷载等效，友好界面的尝试编制，算法的改进，一维变带宽存储等），及部分报告整理； 联系方式：13675131759
张璘琳：部分ANSYS分析；模块中部分函数的编制（列主元高斯消去，单刚的形成，打印函数及其数据的分析处理等等）部分报告整理，及查找资料； 联系方式：13851874856
冯东明：部分ANSYS分析；模块中部分函数的编制（节点编码的方式及其原则，数据的输入输出等）部分报告整理，及查找资料。 联系方式：13675110317
第 31 页 共 31 页