弹性平面问题的有限元电算求解.docVIP

有限元方法 上机报告 班级：040712 学号姓名：尹鹏一道弹性平面问题的有限元电算求解 班级：040712 学号 姓名：尹鹏 问题重述： 图示为一个悬臂深梁，水平方向的均匀分布力q作用在自由端截面上。若采用图示简单网格，去两个单元，求悬臂梁的位移。其中弹性模量，泊松比，梁的厚度，均布载荷，几何尺寸如图。 有限元法基本原理和程序功能介绍： 此程序为求解弹性平面问题而设计，算法的基本思路为有限元方法，只适合用三角形单元对弹性体进行离散化，以求解平面应力问题。 (一)，离散化： 对于一个具体的问题，需要人工进行电算前的准备工作，即离散化；具体需要进行诸如弹性问题的分类，物理模型的简化，基本单元类型的选取，单元划分，约束简化及单元参数的准备等工作。其次还需要做部分单元分析的工作，如整体结点编号，单元结点编号，结点坐标的确定等。综上，下面列出需要人工准备并输入程序的数据： 材料相关性质参数：弹性模量，泊松比，板厚度； 单元结点编号矩阵（布尔矩阵）； 结点坐标矩阵； 已知的等效结点载荷矩阵； 已知的结点约束矩阵。 (二)，单元分析和整体分析： 将这些工作结束后就可以将准备好的数据输入计算机，让程序进行单元刚度矩阵的计算，整体分析，位移和应力的求解等后续工作了。 本程序中的整体分析过程较为特殊，并没用采用课本上介绍的扩大阶数法和对号集成法，而采用了基于单刚子块的对号集成法。程序中没用形成单元刚度矩阵和各个单元的贡献矩阵，而只计算了每个单元的九个单刚子块，并依据布尔矩阵，将各个单刚子块按照对应关系“搬家”到整刚中。这种方法与上述两种方法相比更为简单，更容易编程实现。另外，利用整刚的对称性，程序采用上三角矩阵的形式储存整刚，这样可以提高运算速度，并节省内存。（但为了满足输出的直观性，上三角形式的整刚将以满阵的形式输出，实际上并没用利用上三角矩阵的优点。） (三)，结点位移的求解： 在整体分析之后，需要结合整体所受载荷和约束情况对问题做最后的求解。 第一，形成整体载荷列向量。与手工计算不同的是不需要分析各个结点的载荷情况，程序中只要求计算者利用载荷移置法将外载荷移置到对应结点上，并形成已知的结点等效载荷矩阵。矩阵的第一列为已知结点载荷的自由度编号，第二列为等效结点载荷值。 其二，形成约束条件。此方法与手工计算的方法相同，计算者需要辨别各结点约束的方向和类型，最后形成已知的结点约束矩阵。矩阵的第一列为已知结点约束的自由度编号，第二列为约束值。 其三，将约束引入方程组，并形成整体载荷列阵。这个部分是整个程序中最精彩的部分之一。采用乘“大数”法将已知的结点约束条件引入方程组中，即联系已知结点约束矩阵，将整刚中对角线上对应元素乘以罚因子（通常取~），同时对应修改该结点的载荷项，即联系已知结点约束矩阵，将该自由度已知的约束值乘以罚因子后赋给该自由度的载荷项。最后将已知结点载荷矩阵和修改后的载荷项按照自由度对应关系集成到整体载荷列阵中。 对于原方程组含有未知的结点载荷项与未知的结点位移项分立在等式的两端，为高斯消元法解方程组带来不便。乘“大数”法的意义在于将结点约束的相关已知条件等价的引入方程组右端的整体载荷列阵中，使得原方程组变为 其中，表示置“大数”后的整刚，表示结点位移列向量，是方程组的未知项，表示将已知结点约束条件引入载荷列阵后的整体载荷列阵。 而目前求解此方程组已不再困难。 第四，方程组的求解。联系乘“大数”后的总体刚度矩阵和整体载荷列阵，使用高斯消元法，最终将各个自由度的位移解出，以结点位移列向量的形式表现出来。 主要符号说明： (一)，常量符号： NE：单元总数； NJ：结点总数； NJF：自由度总数()； NG：已知结点载荷的自由度数目； NZ：支杆总数(或约束自由度数目)； E：材料弹性模量； VM：材料泊松比； T：板厚度； (二)，矩阵和向量符号： ie[NE][3]：布尔矩阵； z[NJ][2]：结点坐标矩阵，即； b[NE][3],c[NE][3]：形函数参数矩阵，由，(轮换)计算得到； aa[NE]：单元面积向量； a[NE]：单元刚度矩阵公因子向量，计算得到； g[2][2]：单元刚度矩阵子块； zk[NJF][NJF]：总体刚度矩阵； PJ[NG][2]：已知结点载荷矩阵； ZJ[NZ][2]：已知结点约束矩阵； P[NJF]：整体载荷列向量； (三)，其它变量： je, i1, i2, i3, i4, i5, i6, i7, i8, i, j, j4, j7, j8, k, k1, ii, m, n, p, q, I, cc为程序各部分的控制变量。 源程序及运行结果： (一)，源程序代码： #include stdlib.h #include stdio.h #define