资源描述：
《结构有限元刚度方程求解基础.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、结构有限元刚度方程求解基础结构有限元法刚度阵的特征：1）一般维数很大2）“0”元素很多3）非“0”元素多集中在矩阵的对角线附近4）一般来说矩阵是对称正定的（特殊的如结构流体耦合问题则刚度阵非对称）有限元的灵魂：求解方程组已知力F求位移δ的情况，则，实际计算中，我们不直接求[K]的逆矩阵，而把[K]转化为上三角或下三角矩阵，最后回代求解方程。线性方程组的解法主要有直接法和迭代法，以下介绍直接法。<1>三角方程组以下矩阵均只考虑方阵。1)向前消去考虑下面方程组如果l11l22≠0,则未知数可确定如下:这就是
2、所谓的向前消去,其一般形式为算法1.1（向前消去：行形式）Ax=b的解覆盖bb(1)=b(1)/L(1,1)fori=2~nb(i)=(b(i)-L(i,1~i-1)b(1~i-1))/L(i,i)end2)向后消去解上三角方程组的类似算法叫向后消去法,算法1.2（向后消去：行形式）Ux=b的解覆盖bb(n)=b(n)/L(n,n)fori=n-1,1,-1b(i)=(b(i)-U(i,i+1~n)b(i+1~n))/U(i,i)end向后消去法的算法实现(Fortran):注:数组du(i)表示向量,
3、whlf(i)表示向量,二维数组whlk(i,j)表示矩阵,变量neqns表示方程数,临时变量tempdu(neqns)=whlf(neqns)/whlk(neqns,neqns)do100i=neqns-1,1,-1temp=0.do200j=i+1,neqnstemp=temp+whlk(i,j)*du(j)200enddodu(i)=(whlf(i)-temp)/whlk(i,i)100enddo基于列的形式：交换循环顺序可得到以上算法的列形式，考虑向前消去，一旦x1解出来，该变量可以从第2~n个
4、方程中去掉，我们可只考虑缩小后的方程组然后我们算出x2，并且从第3~n个方程中去掉x2，依次类推，例如：我们有x1=3，那么我们处理2x2方程组算法1.3（向前消去：列形式）Lx=b的解覆盖bforj=1~n-1b(j)=b(j)/L(j,j)b(j+1~n)=b(j+1~n)-b(j)L(j+1~n,j)endb(n)=b(n)/L(n,n)算法1.4（向后消去：列形式）Ux=b的解覆盖bforj=n,2,-1b(j)=b(j)/U(j,j)b(1~j-1)=b(1~j-1)-b(j)U(1~j-1,
5、j)endb(1)=b(1)/U(1,1)<2>LU分解如前面所见，三角方程比较容易求解，那么我们可以把刚度阵[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/
6、x1,那么更一般的，设且，如果定义：则形如Mk的矩阵的前k个分量都为0时，一般称高斯变换，这样的矩阵是下三角的，元素称为乘子，向量称为高斯向量。上三角化，设，通常可找到高斯变换使得，是上三角矩阵例：如果：则：同样地到n-1步后就可以完成上三角化，注意，分解过程中需要对角线元素不为“0”，有限元刚度阵满足这一条件，如果刚度阵对角线有“0”元素，则说明结构存在刚体位移，无法求解。算法2.1（高斯消去）U存于A的上三角部分fork=1~n-1rows=k+1~nA(rows,k)=A(rows,k)/A(k,
7、k)A(rows,rows~n)=A(rows,rows~n)-A(rows,k)A(k,rows~n)end代码实现：增广矩阵[Kf]的高斯消去do100k=1,neqns-1if(whlk(k,k)==0)exitdo200rows=k+1,neqnsv_gauss=whlk(rows,k)/whlk(k,k)do300j=1,neqnswhlk(rows,j)=whlk(rows,j)-v_gauss*whlk(k,j)300enddowhlf(rows)=whlf(rows)-v_gauss*w
8、hlf(k)200enddo100enddo其他形式：do100k=1,neqns-1if(whlk(k,k)==0)exitdo400rows=k+1，neqnswhlk(rows,k)=whlk(rows,k)/whlk(k,k)!先求高斯向量400enddodo200i=k+1,neqnsdo300j=1,neqnswhlk(i,j)=whlk(i,j)-whlk(i,k)*whlk(k,j)300enddowhlf(i)=whlf(
