结构有限元刚度方程求解基础.pdf

结构有限元刚度方程求解基础 结构有限元法刚度阵的特征： 1）一般维数很大 2 ）“0”元素很多 3 ）非“0”元素多集中在矩阵的对角线附近 4 ）一般来说矩阵是对称正定的（特殊的如结构流体耦合问题则刚度阵非对称） 有限元的灵魂： 求解 方程组 已知力 F 求位移 δ 的情况，则 ，实际计算中，我们不直接求 [K] 的逆矩阵，而把 [K] 转化为上三角或下三角矩阵，最后回代求解方程。 线性方程组的解法主要有直接法和迭代法，以下介绍直接法。 1三角方程组 以下矩阵均只考虑方阵。 1) 向前消去 考虑下面方程组 如果 l11l22 ≠0, 则未知数可确定如下 : 这就是所谓的向前消去 , 其一般形式为 算法 1.1 （向前消去：行形式） Ax=b 的解覆盖 b b(1)=b(1)/L(1,1) for i=2~n b(i)=(b(i)-L(i,1~i-1)b(1~i-1))/L(i,i) end 2) 向后消去 解上三角方程组 的类似算法叫向后消去法 , 算法 1.2 （向后消去：行形式） Ux=b 的解覆盖 b b(n)=b(n)/L(n,n) for i=n-1,1,-1 b(i)=(b(i)-U(i,i+1~n)b(i+1~n))/U(i,i) end 向后消去法的算法实现 (Fortran): 注 : 数组 du(i) 表示向量 ,whlf(i) 表示向量 , 二维数组 whlk(i,j) 表示矩 阵 , 变量 neqns 表示方程数 , 临时变量 temp du(neqns)=whlf(neqns)/whlk(neqns,neqns) do 100 i=neqns-1, 1, -1 temp=0. do 200 j=i+1, neqns temp = temp + whlk(i,j)*du(j) 200 enddo du(i) = (whlf(i)-temp)/whlk(i,i) 100 enddo 基于列的形式：交换循环顺序可得到以上算法的列形式，考虑向前消去，一旦 x1 解出来，该变量可以从第 2~n 个方程中去掉，我们可只考虑缩小后的方程组 然后我们算出 x2 ，并且从第 3~n 个方程中去掉 x2 ，依次类推，例如： 我们有 x1=3，那么我们处理 2x2 方程组 算法 1.3 （向前消去：列形式） Lx=b 的解覆盖 b for j=1~n-1 b(j)=b(j)/L(j,j) b(j+1~n)=b(j+1~n)-b(j)L(j+1~n,j) end b(n)=b(n)/L(n,n) 算法 1.4 （向后消去：列形式） Ux=b 的解覆盖 b for j=n,2,-1 b(j)=b(j)/U(j,j) b(1~j-1)=b(1~j-1)-b(j)U(1~j-1,j) end b(1)=b(1)/U(1,1) 2LU 分解 如前面所见，三角方程比较容易求解，那么我们可以把刚度阵 [K] 经过适当的线 性变换等价成一个三角方程组— Gauss消去法 例：在方程组 中，将第一个方程乘以 2，并且在第二个方程中减去它则： 矩阵形式： 这就是 n=2 的 Gauss 消去法 此算法线性变换过程也可写成矩阵形式， 则最终求出一个下三角矩阵 L 和一个上 三角矩阵 U，使得 A=LU， 然后原始问题 Ax=b 可通过两个三角求解过程求得： Ly=b,Ux=y = Ax=LUx=Ly=b 在 n=2 的水平上，如果 x1 ≠0且 τ=x2/x1, 那么 更一般的，设 且 ，如果 定义： 则 形如 Mk 的矩阵的前 k 个分量都为 0 时，一般称高斯变换，这样的矩阵是下三角 的，元素 称为乘子，向量 称为高斯向量。 上三角化，设