結(jié)構(gòu)有限元剛度方程求解基礎(chǔ)

《結(jié)構(gòu)有限元剛度方程求解基礎(chǔ)》由會員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫。

1、結(jié)構(gòu)有限元剛度方程求解基礎(chǔ)結(jié)構(gòu)有限元法剛度陣的特征：1）一般維數(shù)很大2）“0”元素很多3）非“0”元素多集中在矩陣的對角線附近4）一般來說矩陣是對稱正定的（特殊的如結(jié)構(gòu)流體耦合問題則剛度陣非對稱）有限元的靈魂：求解方程組已知力F求位移δ的情況，則，實(shí)際計(jì)算中，我們不直接求[K]的逆矩陣，而把[K]轉(zhuǎn)化為上三角或下三角矩陣，最后回代求解方程。線性方程組的解法主要有直接法和迭代法，以下介紹直接法。<1>三角方程組以下矩陣均只考慮方陣。1)向前消去考慮下面方程組如果l11l22≠0,則未知數(shù)可確定如下:這就是所謂的

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向后消去法的算法實(shí)現(xiàn)(Fortran):注:數(shù)組du(i)表示向量,whlf(i

3、)表示向量,二維數(shù)組whlk(i,j)表示矩陣,變量neqns表示方程數(shù),臨時(shí)變量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基于列的形式：交換循環(huán)順序可得到以上算法的列形式，考慮向前消去，一旦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,j)endb(1)=b(

5、1)/U(1,1)<2>LU分解如前面所見，三角方程比較容易求解，那么我們可以把剛度陣[K]經(jīng)過適當(dāng)?shù)木€性變換等價(jià)成一個三角方程組—>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,那么更一般的，設(shè)且，如果

6、定義：則形如Mk的矩陣的前k個分量都為0時(shí)，一般稱高斯變換，這樣的矩陣是下三角的，元素稱為乘子，向量稱為高斯向量。上三角化，設(shè)，通?？烧业礁咚棺儞Q使得，是上三角矩陣?yán)喝绻簞t：同樣地到n-1步后就可以完成上三角化，注意，分解過程中需要對角線元素不為“0”，有限元剛度陣滿足這一條件，如果剛度陣對角線有“0”元素，則說明結(jié)構(gòu)存在剛體位移，無法求解。算法2.1（高斯消去）U存于A的上三角部分fork=1~n-1rows=k+1~nA(rows,k)=A(rows,k)/A(k,k)A(rows,rows~n)=A

7、(rows,rows~n)-A(rows,k)A(k,rows~n)end代碼實(shí)現(xiàn)：增廣矩陣[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*whlf(k)200enddo100endd

8、o其他形式：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(