结构总刚方程的特性及其求解方法.doc

第九章 结构总刚方程的特性及其求解方法一、 结构总刚矩阵的特性l 对称性由： l 稀疏性当结构的离散单元较多时，对一个节点p的变形能提供刚度仅与该节点相连的那些单元有关（称这些单元上的节点为节点p的相关联节点）。相关联节点总比总节点数要少的多，故总体刚度矩阵中每一行的非零元素只占该行元素总数的一少部分。这样的矩阵称为稀疏矩阵。可利用总刚阵的稀疏性，设法只存储总刚阵中与非零元素相关的那些矩阵元素，可大大节省计算机的内存。l 带状性 带状性是指稀疏矩阵中的元素比较集中地位于对角线元附近。当我们进行结构离散化网格，只要对节点编号稍加注意，就总能做到使每一个节点与相关联节点的编号比较接近。这样，就可使总刚阵中的非零元素较集中于对角线的两侧。l 奇异性 有限元的分析计算步骤中，最初形成单元刚阵和组集总刚阵时，未考虑分片差值曲线（面）对边界条件的满足。从能量原理讲，此时无法获得最小能量解。由此，必须置入适当的边界条件。边界条件的置入是在离散化单元组集后的能量变分后进行的（因为此时获得了明确的结构刚度矩阵方程）。从物理上分析，在置入边界条件前，结构总刚阵是奇异的。二、结构总刚的处理方案l 奇异消除的置大数法及其原理 以前学过强加边界条件的划行划列及置“1”法。虽然这些方法是准确的，但都要对矩阵结构做很多处理，显然不十分方便。置大数法是强加边界条件的一种近似方法，该方法对刚阵的处理却十分简单。/ 置大数法步骤 设给定的边界条件： （为一个或多个节点位移）在结构总刚对应行的对角线上，取： （D为一大数，102040 ） 在载荷列阵相应的行上，取：/ 验证置大数法的近似性 对应第i行刚阵，写出其相应刚度方程： D / 置大数法可用来近似求解支反力在给定节点位移的点上，节点力是未知的（注意在位移与载荷边界同点上的处理原则，在这些点上载荷不再考虑）。在结构总刚度方程中，即对应给定节点位移的那一行刚度方程的右端是未知数。即： 由置大数我们获得了全部位移解，此时的方程应为： 用D遍除： 上式近似于： 但应注意：一定是给定零位移条件下有较好的精度，否则精度十分差。如： 显然差之千里。 / 罚函数约束变分原理 置大数法可归入罚函数约束变分原理当中，以下简介之。 考虑对一个弹性体的无约束变分问题，增加一组约束方程组： (研究域) 该约束条件的置入，可采用乘子法。 但拉氏乘子法变分原理是以扩大计算变量为代价的，当然，可通过对约束方程的了解和推导，建立起拉氏乘子的物理形式，从而获得广义变分原理。但总的来说，计算量扩大或变得难以处理。罚函数约束变分原理没有上述缺点，也不过多改变刚度方程性 质，但解的性质是随罚函数的增大而趋向好的近似性。 思路及方法： 及 罚函数约束变分原理是把直接约束条件转换成积分约束，即新泛函为： （ 求解时再不考虑的条件） 其中，是一个罚数，如果本身是要取极小值的，则应当是一个正数（因为后面的积分大于或等于零），使泛函取驻值而得到的解只是近似地满足约束条件，值越大，约束条件满足得越好。 直观上看，要使有极小值，若当很大时，如：为有限值，则不可能取极小值，故只有当很大时，（这不是证明，只能说明这种解的趋势）。 举一个例子： （引入单个约束时不需积分） 关于参数的极小值，即： 可以看到随的增大， 2 6 10 100a1=-12.00 -12.00-12.00-12.00-12.00a2=-13.5-13.00 -12.43-12.27-12.03 在有限元中，若引入一强制边界条件： 能量泛函修正为： ； 由 这正是上面所讲的置大数法。v/ 结构总刚的存贮方案之一：变带宽一维存贮法 因为结构总刚是稀疏对称的，形如右图。为节省内存，可采用存贮每一行从第一个非“0” 元到对角线元的办法。各行顺序接起来，并用一管理数组记录每一行从第一个非“0”元到对角线元的个数即可。这就是常用的一种有限元结构总刚的存贮方案。从第一个非“0”元到对角线元结构的元素个数叫做每行的带宽，记为。1 n1n1+1 n2n2+1 n3. 叫作矩阵的包络。 三、结构总刚方程的计算方法 （1）分解法由对称满秩矩阵的三角分解定理，可作如下唯一分解：其中 下三角阵 总刚方程的计算： 记：记： . 得： . 由方程是下三角阵，可采用向下回代过程计算 即： 由方程中是上三角单位阵，可采用向上回代过程求得 剩余的问题是如何将进行分解，采用对等法： （可只计算下三角阵部分元素即可） 由此， ；仔细研究上述数值计算，可以发现：在每行的第一个非“0”元素前，分解过程不产生非“0”元。即元素仅处于原总刚的带宽之间，这保证了一维存贮方式的有效性。在分解K阵过程中，每一个元素仅用到过一次，这保证了分解好的元素可以放在总刚的原位置上。这样不必开设新的内存空间，仅用原总刚存贮分解好的总刚。仔细研究以下的分解计算步骤： ikm1由公式：（从上图中，只能看到第一个非零元）看求和号内元素乘积规律，可以发现：a. 从1到列，元素都为“0”，又为“0”。故不会在带宽之外产生非“0”元。b. (第一个非“0”元)c. d. 再细加程序判断，即可以形成算法程序。作业： 研究中间变量及的生成算法细图。（2）子结构方法我们以前多次谈到内节点和内自由度问题，其实是指不与其它节点或自由度相关联的那些，或仅在局部单元内部，而不与外部发生联系的那些节点或自由度。如图a中的点i或图b梁元两端的转角自由度。这个概念可以延伸到一个复杂的结构上。把任何数量毗连的元素拼在一起成为一个新的结构单元，称之为“超单元”或“子结构”，用它外边界上的节点和其它的子结构以及一般的元素共同组集成整个结构。这就是“子结构法”的思想。梁元子结构的方法的关键点是把“超单元”的那些内部节点或自由度通过等价变换（称为缩象），用外边界上的自由度来表征（包括节点或载荷列阵）。这对于求解大型复杂结构（例如整架飞机）的优点是可大大降低对机器的存贮量。具体子结构法步骤： 设：内部节点的编号为，外边界上点的编号为由此可形成该超单元的刚度矩阵，即： 应当注意，这里的Si是已知的, 即作用在内节点上的等效载荷。用矩阵块乘法得： 由式 代入式得： 由方程与其它子结构刚度方程组集，通过计算线性方程组可获得解，再通过式获得解。具体计算思路：证： ， 这样可以利用方法对进行分解，从而利用矩阵分解法获得及，而不必计算的逆。这是因为的稀疏性较好，而则可能是满秩的，而利用分解法求解方程，可保持矩阵的稀疏性，而不必用求逆方法计算。四、关于应力计算在位移模式的有限元方法中，第一步获得的解是结点位移，由结点位移通过单元的几何矩阵求得应变、应力，即对于应力计算（在位移有限元方法中）应当有几点概念：应力的连续性程度较位移的连续性程度要低一个阶次，特别在单元的边界上，可能发生应力的不连续的状态。如最明显的常应力单元导致问题的最主要根源在于形状函数的选择，特别是分片插值函数的选择。应力计算的精度较位移计算的精度要低一个阶次（一个数量级），主要误差来源有二，一是计算过程的数值误差的传递；二是插值函数的导数精确逼近。在等参元当中，可以证明Gauss点上的应力与其他点上的应力计算精度要高一个阶次。因此，在等参元中一般仅计算Gauss点上的应力，代表单元的应力水平。对于应力不连续的情形可采用多种修正方案来处理，使得应力计算趋于合理化。 第十章 现代有限元数值计算技术的发展 各种能量原理的分区组合应用，以求解关心域的高精