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

结构有限元刚度方程求解基础 结构有限元法刚度阵的特征： 1）一般维数很大 2）“0”元素很多 3）非“0”元素多集中在矩阵的对角线附近 4）一般来说矩阵是对称正定的（特殊的如结构流体耦合问题则刚度阵非对称） 有限元的灵魂： 求解方程组 已知力 F 求位移 δ 的情况，则，实际计算中，我们不直接求[K]的逆矩阵，而把[K]转化为上三角或下三角矩阵，最后回代求解方程 线性方程组的解法主要有直接法和迭代法，以下介绍直接法 三角方程组三角方程组以下矩阵均只考虑方阵 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) LULU 分解分解如前面所见，三角方程比较容易求解，那么我们可以把刚度阵[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 时，一般称高斯变换，这样的矩阵是下三角的，元素称为乘子，向量称为高斯向量。

上三角化，设，通常可找到高斯变换使得，是上三角矩阵 例： 如果： 则： 同样地 到 n-1 步后就可以完成上三角化， 注意，分解过程中需要对角线元素不为“0”，有限元刚度阵满足这一条件，如 果刚度阵对角线有“0”元素，则说明结构存在刚体位移，无法求解 算法 2.1（高斯消去）U 存于 A 的上三角部分 for k=1~n-1 rows=k+1~n A(rows,k)=A(rows,k)/A(k,k) A(rows,rows~n)=A(rows,rows~n)-A(rows,k)A(k,rows~n) end 代码实现：增广矩阵[K f]的高斯消去 do 100 k=1, neqns-1 if(whlk(k,k)==0) exit do 200 rows=k+1, neqns v_gauss=whlk(rows,k)/whlk(k,k) do 300 j=1,neqns whlk(rows,j)=whlk(rows,j)-v_gauss*whlk(k,j) 300 enddo whlf(rows)=whlf(rows)- v_gauss*whlf(k) 200 enddo 100 enddo 其他形式： do 100 k=1, neqns-1 if(whlk(k,k)==0) exit do 400 rows=k+1，neqns whlk(rows,k)=whlk(rows,k)/whlk(k,k) !先求高斯向量 400 enddo do 200 i=k+1, neqns do 300 j=1,neqns whlk(i,j)=whlk(i,j)-whlk(i,k)*whlk(k,j) 300 enddo whlf(i)=whlf(i)- whlk(i,k)*whlf(k) 200 enddo 100 enddo 使用后原矩阵被覆盖 例： [U]存于[K]的上三角. 特殊的特殊的 LULU 分解分解当问题中出现如对称性,稀疏性等特性时,需要改进一般矩阵算法以提高效率. 1) 分解 可将一般矩阵 A 分解成 ,其中[L],[M]为下三角矩阵,[D]为对角矩阵. 假设已知,L 的前 j-1 列,D 的前 j-1 个对角元素,M 的前 j-1 行,为求 L 的第 j 列 L(j+1~n,j),M 的第 j 行 M(j,1~j-1)和 D 的对角元 d(j),取出方程第 j 列等式,即: A(1~n,j)=Lv 其中（表示 v 是[U]中的第 j 列），现将 v 的上半部 v(1~j)定义为已知的下三角方程组: 的解,一旦求得 v,则可以计算: d(j)=v(j) M(j,i)=v(i)/d(i), i=1~j-1 v 的下半部 v(j+1~n)有关系式 重新组织后可用于求 L 的第 j 列： 算法 3.1: for j=1~n 从 L(1~j,1~j)v(1~j)=A(1~j,j)解出 v(1~j) for i=1~j-1 M(j,i)=v(i)/d(i) end d(j)=v(j) L(j+1~n,j)=(A(j+1~n,j)-L(j+1~n,1~j-1)v(1~j-1))/v(j) end 2) 分解(LDU) 如果[A]是对称的，则分解中的[M]=[L]. 向量 v(1~j)是由的前 j 个分量定义的，由于[M]=[L]，得 由 L（1~j，1~j）v（1~j）=A（1~j，j）的第 j 个方程: L（j，1~j）v（1~j）=A（j，j） 则有关系式 L（j，1~j-1）v（1~j-1）+ L（j，j）v（j） =A（j，j） ð v（j）=A（j，j）- L（j，1~j-1）*v(1~j-1) 故有: 算法 3.2： for j=1~n for i=1~j-1 v(i)=A(j,j)-L(j,1~j-1)v(1~j-1) end v(j)=A(j,j)-L(j,1~j-1)v(1~j-1) d(j)=v(j) L(j+1~n,j)=(A(j+1~n,j)-L(j+1~n,1~j-1)v(1~j-1))/v(j) end 代码实现，[K]被下三角 L，对角 D 覆盖 program main implicit none parameter neqns=3 real::whlk(neqns,neqns)=(/10,20,30,20,45,80,30,80,171/) real::v(neqns)=0，aa integer::i,j,k,l !j=1 时单独计算，否则出错 aa=1.0/whlk(1,1) do 1000 i=2,neqns whlk(i,1)=whlk(i,1)*aa 1000 enddo do 100 j=2,neqns !计算 v(1~j-1) do 200 i=1,j-1 v(i)=whlk(j,i)*whlk(i,i) 200 enddo !计算 v(j) do 300 k=1,j-1 v(j)=whlk(j,j)-whlk(j,k)*v(k) ！储存 d(j) whlk(j,j)=v(j) 300 enddo !计算 L(j+1~n,j) do 400 k=j+1,neqns do 500 l=1,j-1 whlk(k,j)=(whlk(k,j)-whlk(k,l)*v(l))/v(j) 500 enddo 400 enddo 100 enddo 例： 执行完后 分解后可用 Ly=b（向前消去），Dz=y，LTx=Z（向后消去）得到 Ax=b 的解. 改进，不用 v 数组,200 和 300 循环可以合并，500 循环中的 V 也可以用上面计 算的刚度阵中的元素相乘得到。

顺便调整下循环编号 program main ！code LDU implicit none parameter neqns=3 real::whlk(neqns,neqns)=(/10,20,30,20,45,80,30,80,171/) real::aa,whlf(neqns)=(/1,2,3/),u(neqns) integer::i,j,k,l aa=1.0/whlk(1,1) do 1000 i=2,neqns whlk(i,1)=whlk(i,1)*aa 1000 enddo do 100 i=2,neqns do 200 k=1,i-1 whlk(i,i)=whlk(i,i)-whlk(i,k)*whlk(i,k)*whlk(k,k) 200 enddo do 400 j=i+1,neqns do 500 k=1,i-1 whlk(j,i)=(whlk(j,i)-whlk(j,k)*whlk(i,k)*whlk(k,k)) 500 enddo whlk(j,i)= whlk(j,i)/ whlk(i,i) 400 enddo 100 enddo write(*,110)((whlk(i,j),j=1,3),i=1,3) 110 format(3e12.4) u=whlf 把 whlf 赋值给 u do 2000 i=1,neqns ！下三角向前消去解 Ly=whlf do 2100 k=1,i-1 u(i)=u(i)-whlk(i,k)*u(k) 2100 enddo 2000 enddo do 2200 i=1,neqns !解 Dz=y u(i)=u(i)/whlk(i,i) 2200 enddo do 2300 i=neqns-1,1,-1 !解上三角 U=LTx=z do 2400 j=i+1,n u(i)=u(i)-whlk(j,i)*u(j) 2400 enddo 2300 enddo end P．S． 以上代码只用到了结构刚度阵的对称特性，关于利用刚度阵对称正定的特性， 则存在 Cholesky 分解。

如果 A 是对称正定的方阵，则存在唯一一个对角元全部大于 0 的下三角矩阵G，满足 比较等式的第 j 列可得： 也就是说： 如果 G 的前 j-1 列已知，则可以计算出 v，由，L 是单位上三角，D是对角阵，所以 得出： （Cholesky）算法 3.3 for j=1~n v（j~n）=A(j~n,j) for k=1~j-1 v(j~n)=v(j~n)-G(j,k)G(j~n,k) end G(j~n,j)=v(j~n)/sqrt(v(j)) program main implicit none parameter neqns=2 real::whlk。