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

结构有限元刚度方程求解基础 结构有限元法刚度阵的特征： 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时，一般称高斯变换，这样的矩阵是下三角的，元素称为乘子，向量称为高斯向量。 上三角化，设，通常可找到高斯变换使得，是上三角矩阵 例： 如果： 则： 同样地 到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 end