计算力学程序分析报告(成曦、宋世伟、张璘琳、冯东明)

《计算力学》课程大作业

小组成员：成 曦（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

计算力学——课程设计

withini=1; within1=n; } if(dof[m][n][0]==j) { withinj=1; within2=n; } if(withini*withinj==1) { Kframe[i][j]+=gk[m][within1][within2]; withini=0; withinj=0; } } } }

⑤void getload() //编出荷载列阵获取函数； 在成员变量中，我们定义了double gl[Nnode][3]为节点荷载列阵，在初始数据中我

们输入了该变量，故可直接获得节点的荷载列阵。

实际情况中，我们往往遇到的荷载并非直接施加于节点上，出现概率较大的往往

是梁上集中荷载、集中弯矩、均布荷载或线性分布荷载。但我们可以利用荷载等效原则，将各种荷载状况转换为节点荷载。 具体函数可见：equiload.cpp文件。

⑥ void caculate() //利用先处理法对总体刚度矩阵进

行处理，再求得各节点位移；

在求解平衡方程之前，我们先对刚度矩阵前处理。由于框架与剪力墙部分边界结

点是固定支承，它们的所有位移分量已知为零。故位移分量已知为零的结点，其结点编码均可设为“0”,且编码为“0”的结点一律不设置地址，只有位移分量未知的结点才可进入提供的内存空间。

按此原则，我们可以形成计算矩阵Ab[i][j]，然后根据高斯消去法求解方程，求

得各节点的位移。最后将各节点的位移反带到各单元的单元刚度矩阵中，便可求得节点力。

⑦ void print() //打印出计算结果力和位移

在此函数中，我们将各节点的位移与荷载输入到output.txt文件中，以便查询。

第 6 页 共 31 页

计算力学——课程设计

开始 输入原始数据 提取数据生成梁单元 读取节点自由度的的整体编码和局部编码 计算坐标转换矩阵 组装总刚矩阵 计算刚架单元单刚 组装总荷载向量zonghezaixiangli计算单元节点荷载列阵 解方程求解节点位移（高斯消去法） 计算节点力 输出相关计算结果 图3 框架模块流程图

（2） 剪力墙模块的具体流程图如图4所示。 Wall类由以下几部分组成： 1、成员变量：

const int Nmeshx=2; //定义剪力墙X方向上的划分数 const int Nmeshy=4; //定义剪力墙X方向上的划分数 double point[4][2]; //定义剪力墙结构的四个角点 double E,v; //定义弹性模量和泊松比 double gp[Nmeshy+1][Nmeshx+1][2]; //得到内部节点的坐标数组 int dof[Nmeshx*Nmeshy][8][2]; //获取自由度的整体编码和局部编码 double gk[Nmeshy][Nmeshx][8][8]; //得到单元刚度的数组 double Kwall[(Nmeshx+1)*(Nmeshy+1)*2][(Nmeshx+1)*(Nmeshy+1)*2]; //得到总体刚度的数组 double gl[(Nmeshx+1)*(Nmeshy+1)][2]; //获得节点荷载列阵 double loadwall[(Nmeshx+1)*(Nmeshy+1)*2]; //获得剪力墙的荷载列阵(单位为kN,m) double displacement[(Nmeshx+1)*(Nmeshy+1)*2]; //位移列阵(单位为m) double Fpoint[(Nmeshx+1)*(Nmeshy+1)*2][2]; //算出节点上的力(单位为kN,m)

第 7 页 共 31 页

计算力学——课程设计

定义剪力墙的角点坐标 划分网格，获得内部节点的坐标数组 计算节点自由度的的整体编码和局部编码 等参元 集成剪力墙总刚 计算四节点平面单元刚度矩阵 组装总荷载向量zonghezaixiangli解方程求解节点位移 计算单元节点荷载列阵 计算节点力 输出相关计算结果 图4 剪力墙模块流程图

3、 成员函数

①void inputdata() //用于输入相关数据，包括剪力墙四个角点的坐标、单元属性与节点荷载；

②void getpoint() //获取节点各自由度的局部编码与整体编码；

在程序中，我们定义了常变量Nmeshx与Nmeshy，即两个方向的划分数。根据网格划分数与单元节点数，我们可以确定各节点的坐标。其中，gp[i][j][0]存储横坐标，而gp[i][j][1]存储纵坐标。

③void getKelement() //获取单元刚度矩阵函数；

根据网格划分情况与单元节点数，我们可以得到各节点自由度的整体编码与局部编码。其中，dof[i][j][0]存储整体编码，而dof[i][j][1]存储局部编码。根据计算力学中介绍 平面四节点单元特性，我们可以计算得出其单元刚度矩阵，是一个8?8的矩阵。 ④void getKwall() //单元刚度合成整体刚度矩阵函数；

在输入数据中，我们定义了节点编码与自由度。其中自由度是按一定顺序输入的，故其本质就是每个节点两个自由度的定位向量。在定位向量的基础上，我们可以进行总体刚度矩阵的合成。具体思路是：每个单元刚度矩阵中的元素的行标与列标均有具体的定位向量与其对应。在合成总刚的时候，总体刚度矩阵中的行标与列阵就是定位向量，然后我们利用for循环在各个单元刚度矩阵中寻找与某定位向量对应的元素，然后相加得到最终的总刚。

关键函数如下： int withini,withinj,within1,within2; //判断单元是否包含两个作用的自由度 for(i=0;i<(Nmeshx+1)*(Nmeshy+1)*2;i++) {

第 8 页 共 31 页

计算力学——课程设计

for(int j=0;j<(Nmeshx+1)*(Nmeshy+1)*2;j++) { for(int m1=0;m1

⑤void getload() //编出荷载列阵获取函数； 在成员变量中，我们定义了double gl[Nnode][3]为节点荷载列阵，在初始数据中我

们输入了该变量，故可直接获得节点的荷载列阵。

实际情况中，我们往往遇到的荷载并非直接施加于节点上，出现概率较大的往往

是梁上集中荷载、集中弯矩、均布荷载或线性分布荷载。但我们可以利用荷载等效原则，将各种荷载状况转换为节点荷载。 具体函数可见：equiload.cpp文件。

⑥ void caculate() //利用先处理法对总体刚度矩阵进行处

理，再求得各节点位移；

在求解平衡方程之前，我们先对刚度矩阵前处理。由于框架与剪力墙部分边界结

点是固定支承，它们的所有位移分量已知为零。故位移分量已知为零的结点，其结点编码均可设为“0”,且编码为“0”的结点一律不设置地址，只有位移分量未知的结点才可进入提供的内存空间。

按此原则，我们可以形成计算矩阵Ab[i][j]，然后根据高斯消去法求解方程，求

得各节点的位移。最后将各节点的位移反带到各单元的单元刚度矩阵中，便可求得节点力。

⑦ void print() //打印出计算结果力和位移

在此函数中，我们将各节点的位移与荷载输入到output.txt文件中，以便查询。

第 9 页 共 31 页

计算力学——课程设计

（3） 耦合模块具体流程 fw类由以下几部分组成： 1、 成员变量：

const int Njoint=4; //定义节点数； int joint[Njoint][2]; //第一个元素表示剪力墙的整体坐标，第二个元素表示框架的整体编码 double Kframewall[(Nmeshx+1)*(Nmeshy+1)*3+Nnode*3][(Nmeshx+1)*(Nmeshy+1)*3+ Nnode*3]; //所有点的刚度矩阵 double Kwallch[(Nmeshx+1)*(Nmeshy+1)*3][(Nmeshx+1)*(Nmeshy+1)*3]; //剪力墙的变换后的刚度矩阵 double loadfw[(Nmeshx+1)*(Nmeshy+1)+Nnode][3]; //所有点的荷载 double displacementfw[(Nmeshx+1)*(Nmeshy+1)*3+Nnode*3]; //所有点的位移

2、成员函数

①void inputdata() //读取框架与剪力墙耦合点的剪力墙节点编号与框架的节点编号； ② Kwallchange() //剪力墙总刚的扩充；

由于剪力墙模块中，我们使用的是平面四节点单元，每个节点只有两个自由度。而框架梁单元中每个节点有三个自由度。故如果对两个单元进行耦合，我们认为需对剪力墙的单元刚度矩阵进行扩充，即：每个节点由原来的两个自由度扩充为三个自由度，单元刚度矩阵也随之由原来的8?8矩阵扩充为12?12矩阵。扩充的那两个自由度对应刚度矩阵中的元素均置0。 ③getKfw() //所有点的总刚处理，对所有点进行新的编码，剪力墙原编码统一加上一个数16形成新的编码；

在这里我们重新定义了耦合后的刚度矩阵double Kframewall[(Nmeshx+1)*(Nmeshy +1)*3+Nnode*3][(Nmeshx+1)*(Nmeshy+1)*3+Nnode*3]。其数组大小是根据框架的总结点数与剪力墙的划分情况而确定的。

在框架总体刚度矩阵与剪力墙刚度矩阵耦合时，首先根据节点自由度判断是否为耦合点自由度，若是则剪力墙总体刚度矩阵与框架总体刚度矩阵中的相关项进行叠加。 for(i=0;i<(Nmeshx+1)*(Nmeshy+1)*3;i++) //判断耦合节点，从而将

剪力墙总刚与框架总刚进行集成 { for(int j=0;j<(Nmeshx+1)*(Nmeshy+1)*3;j++) { int mi1,mi2,mi3,mj1,mj2,mj3; mi1=i/3; mi2=i%3; mi3=0; mj1=j/3; mj2=j%3; mj3=0; for(int n=0;n

第 10 页 共 31 页

计算力学——课程设计

表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 页

计算力学——课程设计

分入手，这些步骤我们可以通过手算加以推导。经手算结果与C++结果比较，我们发现了问题的所在。改完程序，得到了比较理想的结果输出，大家终于可以放松了，而且都很有感慨，于是写下来，得到了上面这段话，作为对自己多日来辛勤劳动的一种总结吧！

总之，计算力学是一门多学科综合的课程，所要求的掌握内容多、理解层次高，是结构工程专业不可或缺的专业课。通过此次大作业，我们既熟悉了C++编程的相关的知识，更加深了对计算力学这门课程的理解，可谓受益匪浅。

九、各人工作量

成曦：部分调试，部分优化（耦合），模块的框架设计，荷载的处理和单刚到总刚的合成处理和计算矩阵的前处理等等），SAP2000的分析，及部分报告整理； 联系方式：13655165465

宋世伟：部分ANSYS分析；部分调试，部分优化（耦合中部分程序，等参部分，荷载等效，友好界面的尝试编制，算法的改进，一维变带宽存储等），及部分报告整理； 联系方式：13675131759

张璘琳：部分ANSYS分析；模块中部分函数的编制（列主元高斯消去，单刚的形成，打印函数及其数据的分析处理等等）部分报告整理，及查找资料； 联系方式：13851874856

冯东明：部分ANSYS分析；模块中部分函数的编制（节点编码的方式及其原则，数据的输入输出等）部分报告整理，及查找资料。 联系方式：13675110317

第 31 页 共 31 页

本文来源：https://www.bwwdw.com/article/562.html