基于Matlab語言的按平面三角形單元劃分的結構有限元程序設計

基于Matlab語言的按平面三角形單元劃分的結構有限元程序設計專 業(yè)： 建筑與土木工程班 級： 建工研12-2姓 名： 韓志強學 號： 基于Matlab語言的按平面三角形單元劃分結構有限元程序設計一、 有限單元發(fā)及Matlab語言概述1. 有限單元法隨著現(xiàn)代工業(yè)、生產(chǎn)技術的發(fā)展，不斷要求設計高質(zhì)量、高水平的大型、復雜和精密的機械及工程結構。為此目的，人們必須預先通過有效的計算手段，確切的預測即將誕生的機械和工程結構，在未來工作時所發(fā)生的應力、應變和位移因此，需要尋求一種簡單而又精確的數(shù)值分析方法。有限單元法正是適應這種要求而產(chǎn)生和發(fā)展起來的一種十分有效的數(shù)值計算方法。有限元法把一個復雜的結構分解成相對簡單的“單元”，各單元之間通過結點相互連接。單元內(nèi)的物理量由單元結點上的物理量按一定的假設內(nèi)插得到，這樣就把一個復雜結構從無限多個自由度簡化為有限個單元組成的結構。我們只要分析每個單元的力學特性，然后按照有限元法的規(guī)則把這些單元“拼裝”成整體，就能夠得到整體結構的力學特性。有限單元法基本步驟如下：(1)結構離散：結構離散就是建立結構的有限元模型，又稱為網(wǎng)格劃分或單元劃分，即將結構離散為由有限個單元組成的有限元模型。在該步驟中，需要根據(jù)結構的幾何特性、載荷情況等確定單元體內(nèi)任意一點的位移插值函數(shù)。 (2)單元分析：根據(jù)彈性力學的幾何方程以及物理方程確定單元的剛度矩陣。(3)整體分析:把各個單元按原來的結構重新連接起來，并在單元剛度矩陣的基礎上確定結構的總剛度矩陣，形成如下式所示的整體有限元線性方程： 式中，是載荷矩陣, 是整體結構的剛度矩陣, 是節(jié)點位移矩陣。 (4)載荷移置:根據(jù)靜力等效原理，將載荷移置到相應的節(jié)點上，形成節(jié)點載荷矩陣。(5)邊界條件處理:對式所示的有限元線性方程進行邊界條件處理。 (6)求解線性方程:求解式所示的有限元線性方程，得到節(jié)點的位移。在該步驟中，若有限元模型的節(jié)點越多，則線性方程的數(shù)量就越多，隨之有限元分析的計算量也將越大。(7)求解單元應力及應變 根據(jù)求出的節(jié)點位移求解單元的應力和應變。 (8)結果處理與顯示 進入有限元分析的后處理部分，對計算出來的結果進行加工處理，并以各種形式將計算結果顯示出。2. Matlab簡介在用有限元法進行結構分析時，將會遇到大量的數(shù)值計算，因而在實用上是一定要借助于計算機和有限元程序，才能完成這些復雜而繁重的數(shù)值計算工作。而Matlab是當今國際科學界最具影響力和活力的軟件。它起源于矩陣運算，并已經(jīng)發(fā)展成一種高度集成的計算機語言。它提供了強大的科學計算，靈活的程序設計流程，高質(zhì)量的圖形可視化與界面設計，便捷的與其他程序和語言接口的功能。Matlab在各國高校與研究單位起著重大的作用?！肮び破涫拢叵壤淦鳌?。如果有一種十分有效的工具能解決在教學與研究中遇到的問題，那么Matlab語言正是這樣的一種工具。它可以將使用者從繁瑣、無謂的底層編程中解放出來，把有限的寶貴時間更多地花在解決問題中，這樣無疑會提高工作效率。 目前，Matlab已經(jīng)成為國際上最流行的科學與工程計算的軟件工具，現(xiàn)在的Matlab已經(jīng)不僅僅是一個“矩陣實驗室”了，它已經(jīng)成為了一種具有廣泛應用前景的全新的計算機高級編程語言了，有人稱它為“第四代”計算機語言，它在國內(nèi)外高校和研究部門正扮演著重要的角色。Matlab語言的功能也越來越強大，不斷適應新的要求提出新的解決方法?？梢灶A見，在科學運算、自動控制與科學繪圖領域Matlab語言將長期保持其獨一無二的地位。為此，本例采用Matlab語言編程，以利用其快捷強大的矩陣數(shù)值計算功能。二、 問題描述一矩形薄板，一邊固定，承受150kN集中荷載，結構簡圖及按平面三角形單元劃分的有限元模型圖如下所示。材料參數(shù)：彈性模量；泊松比：；薄板厚度。在本例中，所取結構模型及數(shù)據(jù)主要用于程序設計理論分析，與工程實際無關。三、 參數(shù)輸入：單元個數(shù)NELEM=4 節(jié)點個數(shù)NNODE=6受約束邊界點數(shù)NVFIX=2 節(jié)點荷載個數(shù)NFORCE=1彈性模量YOUNG=2e8 泊松比POISS=0.2 厚度THICK=0.002 單元節(jié)點編碼數(shù)組LNODS =節(jié)點坐標數(shù)組COORD =節(jié)點力數(shù)組FORCE =4 0 -150約束信息數(shù)組FIXED =以上數(shù)值數(shù)據(jù)為程序運行前輸入的初始數(shù)據(jù)，存為“.txt”文本格式，同時必須放在Matlab工作目錄下，路徑不對程序不能自動讀取指定初始文件，運行出錯。初始數(shù)據(jù)文本文件輸入格式如下圖：四、 Matlab語言程序源代碼：1. 程序中變量說明NNODE 單元節(jié)點數(shù) NPION 總結點數(shù)NELEM 單元數(shù) NVFIX 受約束邊界點數(shù)FIXED 約束信息數(shù)組 NFORCE 節(jié)點力數(shù)FORCE 節(jié)點力數(shù)組 COORD 結構節(jié)點坐標數(shù)組 LNODS 單元定義數(shù)組YOUNG 彈性模量 POISS 泊松比THICK 厚度 B 單元應變矩陣(3*6) D 單元彈性矩陣(3*3) S 單元應力矩陣(3*6)A 單元面積 ESTIF 單元剛度矩陣ASTIF 總體剛度矩陣 ASLOD 總體荷載向量 ASDISP 節(jié)點位移向量 ELEDISP 單元節(jié)點位移向量STRESS 單元應力 FP1 數(shù)據(jù)文件指針2. Matlab語言程序代碼%*%初始化及數(shù)據(jù)調(diào)用format short e %設定輸出類型clear %清除內(nèi)存變量FP1=fopen(.txt,rt); %打開輸入數(shù)據(jù)文件，讀入控制數(shù)據(jù)NELEM=fscanf(FP1,%d,1), %單元個數(shù)（單元編碼總數(shù)）NPION=fscanf(FP1,%d,1), %結點個數(shù)（結點編碼總數(shù)）NVFIX=fscanf(FP1,%d,1) %受約束邊界點數(shù)NFORCE=fscanf(FP1,%d,1), %結點荷載個數(shù)YOUNG=fscanf(FP1,%e,1), %彈性模量POISS=fscanf(FP1,%f,1), %泊松比 THICK=fscanf(FP1,%f,1) %厚度LNODS=fscanf(FP1,%d,3,NELEM) %單元定義數(shù)組（單元結點號）%相應為單元結點號（編碼）、按逆時針順序輸入COORD=fscanf(FP1,%f,2,NPION) %結點坐標數(shù)組 %坐標：x，y坐標（共NPOIN組）FORCE=fscanf(FP1,%f,3,NFORCE) %結點力數(shù)組（受力結點編號, x方向,y方向）FIXED=fscanf(FP1,%d,3,NVFIX) %約束信息（約束點，x約束，y約束） %有約束為1，無約束為0%*%生成單元剛度矩陣并組成總體剛度矩陣ASTIF=zeros(2*NPION,2*NPION); %生成特定大小總體剛度矩陣并置0%*for i=1:NELEM%生成彈性矩陣D D= 1 POISS 0; POISS 1 0; 0 0 (1-POISS)/2*YOUNG/(1-POISS2)%*%計算當前單元的面積A A=det(1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2); 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2); 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)/2%*%生成應變矩陣B for j=0:2 b(j+1)= COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNODS(i,(rem(j+2),3)+1),2); c(j+1)=-COORD(LNODS(i,(rem(j+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1); end B=b(1) 0 b(2) 0 b(3) 0;. 0 c(1) 0 c(2) 0 c(3);. c(1) b(1) c(2) b(2) c(3) b(3)/(2*A);B1( :,:,i)=B;%*%求應力矩陣S=D*BS=D*B; ESTIF=B*S*THICK*A; %求解單元剛度矩陣 a=LNODS(i,:); %臨時向量,用來記錄當前單元的節(jié)點編號 for j=1:3 for k=1:3 ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2) =ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2); %根據(jù)節(jié)點編號對應關系將單元剛度分塊疊加到總剛%度矩陣中 end end end%*%將約束信息加入總體剛度矩陣（對角元素改一法） for i=1:NVFIX if FIXED(i,2)=1 ASTIF(:,(FIXED(i,1)*2-1)=0; %一列為零 ASTIF(FIXED(i,1)*2-1),:)=0; %一行為零 ASTIF(FIXED(i,1)*2-1),(FIXED(i,1)*2-1)=1; %對角元素為1 end if FIXED(i,3)=1 ASTIF( :,FIXED(i,1)*2)=0; %一列為零 ASTIF(FIXED(i,1)*2,:)=0; %一行為零 ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %對角元素為1 end end%*%生成荷載向量 ASLOD(1:2*NPION)=0; %總體荷載向量置零 for i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end%*%求解內(nèi)力 ASDISP=ASTIFASLOD %計算節(jié)點位移向量 ELEDISP(1:6)=0; %當前單元節(jié)點位移向量 for i=1:NELEM for j=1:3 ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); %取出當前單元的節(jié)點位移向量 end i STRESS=D*B1(:, :, i)*ELEDISP %求內(nèi)力 end fclose(FP1); %關閉數(shù)據(jù)文件五、 程序運行結果NELEM = 4NPION = 6NVFIX = 2NFORCE = 1YOUNG = POISS = 2.0000e-001THICK = 2.0000e-003LNODS = 1 2 6 2 3 4 2 4 5 2 5 6COORD = 0 0 1 0 2 0 2 1 1 1 0 1FORCE = 4 0 -150FIXED = 1 1 1 6 1 1D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A = 5.0000e-001D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A = 5.0000e-001D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A = 5.0000e-001D = 2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A = 5.0000e-001ASDISP = 0 0 -8.0607e-004 -1.5848e-003 -1.0281e-003 -4.4727e-003 1.1937e-003 -4.6947e-003 8.6670e-004 -1.8880e-003 0 0（說明：以上12個ASDISP輸出結果數(shù)據(jù)表示節(jié)點位移向量，即各節(jié)點分別在X，Y方向上產(chǎn)生的位移。）i = 1STRESS = -1.6793e+005 -3.3586e+004 -1.3207e+005i = 2STRESS = -5.5503e+004 -5.5503e+004 -5.5503e+004i = 3STRESS = 5.5503e+004 -4.9526e+004 -9.4497e+004i = 4STRESS = 1.6793e+005 -2.7040e+004 -1.7932e+004（說明：以上四組STRESS輸出結果數(shù)據(jù)表示單元應力SX,SY,SXY，i為單元號。）六、 ANSYS建模比較分析利用ANSYS有限元分析軟件，完全按照前面Matlab程序設計的各項條件參數(shù)以及單元劃分方式建立ANSYS模型，其求解結果與以上程序計算結果比較，驗證程序是否正確。ANSYS模型節(jié)點單元如下圖所示：ANSYS模型求解變形圖如下所示：ANSYS求解節(jié)點位移結果列表顯示如下：（單位：m）PRINT U NODAL SOLUTION PER NODE* POST1 NODAL DEGREE OF FREEDOM LISTING * THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEM NODE UX UY UZ USUM 1 0.0000 0.0000 0.0000 0.0000 2 -0.80607E-03-0.15848E-02 0.0000 0.17780E-02 3 -0.10281E-02-0.44727E-02 0.0000 0.45893E-02 4 0.11937E-02-0.46947E-02 0.0000 0.48441E-02 5 0.86670E-03-0.18880E-02 0.0000 0.20774E-02 6 0.0000 0.0000 0.0000 0.0000 MAXIMUM ABSOLUTE VALUES NODE 4 4 0 4 VALUE 0.11937E-02-0.46947E-02 0.0000 0.48441E-02ANSYS求解單元應力結果列表顯示如下：（單位：KN/m2）PRINT S ELEMENT SOLUTION PER ELEMENT* POST1 ELEMENT NODAL STRESS LISTING * THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATES ELEMENT= 1 PLANE182 NODE SX SY SZ SXY SYZ 1 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 2 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182 NODE SX SY SZ SXY SYZ 2 -55503. -55503. 0.0000 -55503. 0.0000 3 -55503. -55503. 0.0000 -55503. 0.0000 4 -55503. -55503. 0.0000 -55503. 0.0000 ELEMENT= 3 PLANE182 NODE SX SY SZ SXY SYZ 2 55503. -49526. 0.0000 -94497. 0.0000 4 55503. -49526. 0.0000 -94497. 0.0000 5 55503. -49526. 0.0000 -94497. 0.0000 ELEMENT= 4 PLANE182 NODE SX SY SZ SXY SYZ 2 0.16793E+06 -27040. 0.0000 -17932. 0.0000 5 0.16793E+06 -27040. 0.0000 -17932. 0.0000 6 0.16793E+06 -27040. 0.0000 -17932. 0.0000 結論通過比較Matlab語言設計程序運行結果和ANSYS建模分析結果可知，兩種方式算出的結果完全一致，說程序設計正確。所以本程序適用于按三角形單元劃分的平面結構有限元分析。format short eclearFP1=fopen(LinearTriangleElement of Thin plate.txt,rt);NELEM=fscanf(FP1,%d,1)NPION=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)YOUNG=fscanf(FP1,%e,1)POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%f,1)LNODS=fscanf(FP1,%d,3,NELEM)COORD=fscanf(FP1,%f,2,NPION)FORCE=fscanf(FP1,%f,3,NFORCE)FIXED=fscanf(FP1,%d,3,NVFIX)ASTIF=zeros(2*NPION,2*NPION); %生成特定大小總體剛度矩陣并置0for i=1:NELEM%生成彈性矩陣DD= YOUNG/(1-POISS2)*1 POISS 0;POISS 1 0;0 0 (1-POISS)/2%計算當前單元的面積AA=det(1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)/2%生成應變矩陣Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNODS(i,(rem(j+2),3)+1),2);c(j+1)=-COORD(LNODS(i,(rem(j+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1);endB=b(1) 0 b(2) 0 b(3) 0;.0 c(1) 0 c(2) 0 c(3);.c(1) b(1) c(2) b(2) c(3) b(3)/(2*A);B1( :,:,i)=B;%求應力矩陣S=D*BS=D*B;ESTIF=B*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);endendend%將約束信息加入總體剛度矩陣for i=1:NVFIXif FIXED(i,2)=1ASTIF(:,(FIXED(i,1)*2-1)=0;ASTIF(FIXED(i,1)*2-1),:)=0;ASTIF(FIXED(i,1)*2-1),(FIXED(i,1)*2-1)=1;endif FIXED(i,3)=1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷載向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解內(nèi)力ASDISP=ASTIFASLODELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:, :, i)*ELEDISPendfclose(FP1);format short eclearFP1=fopen(LinearTriangleElement of Thin plate.txt,rt);NELEM=fscanf(FP1,%d,1)NPION=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)YOUNG=fscanf(FP1,%e,1)POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%f,1)LNODS=fscanf(FP1,%d, NELEM,3)COORD=fscanf(FP1,%f, NPION,2)FORCE=fscanf(FP1,%f, NFORCE,3)FIXED=fscanf(FP1,%d, NVFIX,3)ASTIF=zeros(2*NPION,2*NPION); %生成特定大小總體剛度矩陣并置0for i=1:NELEM%生成彈性矩陣DD= YOUNG/(1-POISS2)*1 POISS 0;POISS 1 0;0 0 (1-POISS)/2%計算當前單元的面積AA=det(1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(