有限元編程大作業(yè)報告

1、本科生實驗報告書四節(jié)點等參單元有限元分析的FORTRAN程序目錄1 問題概述 (1)2 四節(jié)點四邊形等參單元介紹 (1)3 單元應(yīng)力磨平方法介紹 (4)4 程序流程設(shè)計 (6)4.1 程序設(shè)計概述 4.2 流程圖5 程序結(jié)構(gòu)及程序說明 (8)6 程序應(yīng)用及算例分析 (9)6.1 算例概述6.2 算例ANSYS求解 6.3 算例程序數(shù)值解6.4 算例分析7. 總結(jié) (15)7. 附錄 ()1. 問題概述等參單元是有限元方法中使用最廣泛的單元類型。等參單元的位移模式和坐標變換均采用相同的形函數(shù)，這種坐標變換叫做等參變換。通過等參變換可以將自然（局部）坐標中幾何形狀規(guī)則的單元轉(zhuǎn)換成總體（笛卡爾）坐標

2、中形狀扭曲的單元，因而使得單元有較好的適應(yīng)性。本問題首先對平面四節(jié)點四邊形等參單元的形函數(shù)、應(yīng)力矩陣和等效節(jié)點力矩陣、應(yīng)力磨平公式等的推導(dǎo)和計算求解。并通過設(shè)計FORTRAN求解程序進行編程求解，最后給出算例（受集中荷載的懸臂梁）并進行求解，將解與ANSYS的解進行比較。在這個過程中，采用了高斯三點積分和高斯兩點積分，這種積分方法的求解效率較高而且精度也較好。在問題的最后，嘗試去分析引起數(shù)值解誤差的原因，并分析四節(jié)點等參單元的若干特性。2. 四節(jié)點四邊形等參單元介紹邊長為2的正方形單元（如下圖所示），在其形心處安置一個局部坐標。單元角結(jié)點i的坐標（i,i）分別為1,因此單元四條邊界的方程可用簡

3、單公式1=0和1=0逐一給出。12340 圖2-1 母單元 圖2-2 四邊形單元形函數(shù)Ni的表達式：Ni=1+i(1+i)/4位移函數(shù)： 坐標變換式： ，單元應(yīng)變矩陣 式中單元節(jié)點的位移列陣； 這里記號，。根據(jù)復(fù)合函數(shù)求導(dǎo)規(guī)則，有 從而有 因此，單元內(nèi)的應(yīng)力可以表示為 應(yīng)力矩陣可以寫成分塊形式對于平面應(yīng)力情形， 單元剛度矩陣是其中積分采用三點高斯積分，（高斯積分點的總數(shù)），和是加權(quán)系數(shù)，和是單元內(nèi)的坐標.。對于三點高斯積分，高斯積分點的位置： ,，。結(jié)構(gòu)剛度矩陣結(jié)構(gòu)結(jié)點荷載列陣 注意，對于上兩式中的理解不是簡單的疊加而是按照對應(yīng)的自由度集成?？倓偲胶夥匠虖氖缴鲜角蟪觯簩⒒卮牒椭?，得到和。等效

4、節(jié)點力(1) 集中力 將集中荷載作用點取為結(jié)點，集中荷載就可以直接轉(zhuǎn)化為等效結(jié)點力。(2) 體積力 等效結(jié)點力按下式計算Fvie=FvxiFvyie=-11-11NiPvxPvyhJdd (i=1,2,3,4)其中積分采用高斯2點積分，（高斯積分點的總數(shù)），和或是加權(quán)系數(shù)，和是單元內(nèi)的坐標。對于2點高斯積分，高斯積分點的位置： 1=1=-1/3 ; 1=1 ; 2=2=1/3 ; 2=1(3) 表面力 單元的ij邊上兩個結(jié)點的等效結(jié)點力按下式計算 Fsie=-11Nixi,+yi,yi,-xi,hd (i=3,4)3. 應(yīng)力磨平方法介紹 為了減少改進應(yīng)力結(jié)果的工作量，可以采用單元應(yīng)力的局部磨平

5、。當單元足夠小時，磨平可以在各個單元內(nèi)進行。對于等參元來說，有了單元內(nèi)應(yīng)力局部磨平的方法，可以方便地利用精度較高的高斯積分點(最佳應(yīng)力點)的應(yīng)力值來改進等參元結(jié)點應(yīng)力的近似性質(zhì)。以二維單元為例，如果是二次等參元，插值函數(shù)中的完全多項式是二次，即p2。對于C0型單元m1。這時p - m+1=2,則在22個高斯積分點上有限元的應(yīng)力計算值 有較高的精度。如果進行應(yīng)力磨平時只要求得到4各角點的改進應(yīng)力值，即使在計算位移時是8結(jié)點二次等參元，但進行應(yīng)力磨平時，單元結(jié)點數(shù)可只取4個，即用二維雙線性單元進行應(yīng)力磨平。磨平式各應(yīng)力分量分別進行，這時應(yīng)力磨平插值函數(shù)Ni應(yīng)采用雙線性函數(shù)，即 （3.1）在此情況下

6、，4個結(jié)點上的應(yīng)力可由高斯點上的應(yīng)力外推得到，即令在22高斯積分點上有。4個高斯積分點的座標(見圖3-1)如下: 圖3-1 二維等參單元將高斯點坐標代入(3.1)式得到下面的等式： (3.2） 式中等號左邊的應(yīng)力列陣是有限元中已求出的4個高斯點相應(yīng)的應(yīng)力分量；是磨平后應(yīng)力的結(jié)點值；轉(zhuǎn)換矩陣由(3.22)式的第二式代入高斯點坐標后的插值函數(shù)值構(gòu)成。由（3.2）式求逆可得 （3.3）其中 （3.4）各應(yīng)力分量均可用（3.3）求解。這種改進結(jié)點應(yīng)力的方法亦稱之為應(yīng)力插值外推。求得改進的應(yīng)力結(jié)點值后，如需要求單元內(nèi)部的應(yīng)力值仍可按(3.1)第一式進行計算。 采用單元局部應(yīng)力磨平的方法，對于同一結(jié)點，由

7、不同相鄰單元求得的應(yīng)力改進值通常是不相同的?？砂严嚓P(guān)單元求得的改進結(jié)點值再取平均作為最后的結(jié)點應(yīng)力值。4. 程序流程設(shè)計4.1 程序設(shè)計概述本程序的設(shè)計以P3單元的FORTRAN程序為基礎(chǔ)，在其架構(gòu)之上修改而來。整體的架構(gòu)與P3單元相似，但是在應(yīng)力應(yīng)變矩陣、主單元剛度矩陣、整體剛度矩陣、等效結(jié)點力以及應(yīng)力磨平的算法方面有著較大的差異。具體的算法方面，在剛度矩陣的計算用了3點高斯積分，在等效結(jié)點力的計算用了2點高斯積分。應(yīng)力磨平采用了單元內(nèi)應(yīng)力局部磨平的方法，對于同一結(jié)點，由不同相鄰單元求得的應(yīng)力改進值通常是不相同的。把相關(guān)單元求得的改進結(jié)點值再取平均作為最后的結(jié)點應(yīng)力值。4.2 流程圖5. 程

8、序結(jié)構(gòu)5.1 程序結(jié)構(gòu)5.2 子程序功能說明P4INPUT 讀入單元數(shù)據(jù)、節(jié)點數(shù)據(jù)、節(jié)點約束條件及各類單元屬性DEA 形成主對角元素地址P4SSM 形成總體剛度矩陣P4STIFB 計算單元剛度矩陣P4BMATB 計算當前單元的應(yīng)變矩陣P4MODB 計算當前單元的彈性矩陣P4DBE 計算應(yīng)力矩陣P4XJACM 計算當前單元的形函數(shù)和雅可比行列式BOUNDARY 邊界條件處理LDLT LDLT分解結(jié)構(gòu)剛度矩陣P4LOAK 計算單元荷載矩陣形成節(jié)點總荷載向量用置0置1法處理邊界條件時產(chǎn)生的荷載向量修正項，修正節(jié)點總荷載向量SOLV 處理荷載向量,回代求結(jié)點位移P4STREB 計算單元應(yīng)力和主應(yīng)力P4

9、GAUSTR 單元內(nèi)應(yīng)力局部磨平GSTREB 計算磨平后的單元應(yīng)力P4SUMSTRS 計算磨平后各節(jié)點的應(yīng)力6. 程序應(yīng)用及算例分析6.1 算例概述考慮一個右端固定的懸臂梁如下圖5-1，其材料參數(shù)為：=0.24kg/mm3 ，E=2e5MPa，=0.2；幾何參數(shù)為：L=8mm, H=2mm, b=1mm。在左端上方有集中荷載F=-100N，不考慮懸臂梁的自重，求解懸臂梁的撓度以及各截面的應(yīng)力。 圖6-16.2 算例ANSYS求解撓度圖圖6-2X方向應(yīng)力圖6-3Y方向應(yīng)力圖6-46.3 算例程序數(shù)值解將懸臂梁進行單元網(wǎng)格劃分如下圖，劃分為8個單元15個結(jié)點。 圖6-5將劃分結(jié)果及相關(guān)參數(shù)輸入程序

10、的input文檔，啟動程序進行求解。求解結(jié)果如下： 表6-1 節(jié)點應(yīng)力表6-2 節(jié)點位移6.4 算例分析為了驗證P4求解程序的準確性，我們將懸臂梁最大受拉區(qū)（頂邊）的位移解和應(yīng)力解在ANSYS下的和在本程序下得到的進行比較如下圖：圖6-6 位移解的比較由圖可以看出，本程序的位移數(shù)值解與ANSYS的解十分接近，最大的節(jié)點誤差在11.417%左右。圖6-7 受拉區(qū)x向應(yīng)力解比較由圖可以看出，本程序的應(yīng)力數(shù)值解與ANSYS的解十分接近，最大的節(jié)點誤差在15.39%左右。應(yīng)力磨平結(jié)果分析可見對于四邊形等參單元應(yīng)力磨平前后，局部坐標為（0，0）的點的應(yīng)力值不變。7. 總結(jié)本程序采用的四結(jié)點平面等參單元具有較高的求解精度，對于平面的問題也有很好的適應(yīng)性。P4等參單元的求解還是和P3三角形常應(yīng)變單元相比還是比較復(fù)雜的。雖然求解的流程基本一致，但是等參單元的求解需要解決剛度陣和等效節(jié)點力求解時的數(shù)值積分問題，本程序中采用應(yīng)用比較廣泛的高斯求積法。關(guān)于高斯積分的階數(shù)，在李人憲所著的有限元法基礎(chǔ)中指出，“用數(shù)值積分代替精確積分, 無疑計算時會產(chǎn)生誤差. 為盡量減少這一誤差,希望采用盡可能高的數(shù)值積分階次” 。文獻亦同時指出:鑒于N 很大時計算費時, 建議在二維情況下, 四節(jié)點單元的N 取為2 , 八節(jié)點單元的N 取為3 較好。本程序則混合采用二階