韓志強-有限元程序設計

韓志強-有限元程序設計韓志強-有限元程序設計韓志強-有限元程序設計xxx公司韓志強-有限元程序設計文件編號：文件日期：修訂次數：第1.0次更改批準審核制定方案設計，管理制度基于Matlab語言按平面三角形單元劃分的結構有限元程序設計專業(yè)：建筑與土木工程班級：建工研12-2姓名：韓志強學號：0

基于Matlab語言的按平面三角形單元劃分結構有限元程序設計有限單元發(fā)及Matlab語言概述1.有限單元法隨著現代工業(yè)、生產技術的發(fā)展，不斷要求設計高質量、高水平的大型、復雜和精密的機械及工程結構。為此目的，人們必須預先通過有效的計算手段，確切的預測即將誕生的機械和工程結構，在未來工作時所發(fā)生的應力、應變和位移因此，需要尋求一種簡單而又精確的數值分析方法。有限單元法正是適應這種要求而產生和發(fā)展起來的一種十分有效的數值計算方法。有限元法把一個復雜的結構分解成相對簡單的“單元”，各單元之間通過結點相互連接。單元內的物理量由單元結點上的物理量按一定的假設內插得到，這樣就把一個復雜結構從無限多個自由度簡化為有限個單元組成的結構。我們只要分析每個單元的力學特性，然后按照有限元法的規(guī)則把這些單元“拼裝”成整體，就能夠得到整體結構的力學特性。有限單元法基本步驟如下：(1)結構離散：結構離散就是建立結構的有限元模型，又稱為網格劃分或單元劃分，即將結構離散為由有限個單元組成的有限元模型。在該步驟中，需要根據結構的幾何特性、載荷情況等確定單元體內任意一點的位移插值函數。(2)單元分析：根據彈性力學的幾何方程以及物理方程確定單元的剛度矩陣。(3)整體分析:把各個單元按原來的結構重新連接起來，并在單元剛度矩陣的基礎上確定結構的總剛度矩陣，形成如下式所示的整體有限元線性方程：①式中，是載荷矩陣,是整體結構的剛度矩陣,是節(jié)點位移矩陣。(4)載荷移置:根據靜力等效原理，將載荷移置到相應的節(jié)點上，形成節(jié)點載荷矩陣。(5)邊界條件處理:對式①所示的有限元線性方程進行邊界條件處理。(6)求解線性方程:求解式①所示的有限元線性方程，得到節(jié)點的位移。在該步驟中，若有限元模型的節(jié)點越多，則線性方程的數量就越多，隨之有限元分析的計算量也將越大。(7)求解單元應力及應變：根據求出的節(jié)點位移求解單元的應力和應變。(8)結果處理與顯示：在進行程序設計時，還要進入有限元分析的后處理部分，對計算出來的結果進行加工處理，并以各種形式將計算結果顯示出。2.Matlab簡介在用有限元法進行結構分析時，將會遇到大量的數值計算，因而在實用上是一定要借助于計算機和有限元程序，才能完成這些復雜而繁重的數值計算工作。而Matlab是當今國際科學界最具影響力和活力的軟件。它起源于矩陣運算，并已經發(fā)展成一種高度集成的計算機語言。它提供了強大的科學計算，靈活的程序設計流程，高質量的圖形可視化與界面設計，便捷的與其他程序和語言接口的功能。Matlab在各國高校與研究單位起著重大的作用。“工欲善其事，必先利其器”。如果有一種十分有效的工具能解決在教學與研究中遇到的問題，那么Matlab語言正是這樣的一種工具。它可以將使用者從繁瑣、無謂的底層編程中解放出來，把有限的寶貴時間更多地花在解決問題中，這樣無疑會提高工作效率。目前，Matlab已經成為國際上最流行的科學與工程計算的軟件工具，現在的Matlab已經不僅僅是一個“矩陣實驗室”了，它已經成為了一種具有廣泛應用前景的全新的計算機高級編程語言了，有人稱它為“第四代”計算機語言，它在國內外高校和研究部門正扮演著重要的角色。Matlab語言的功能也越來越強大，不斷適應新的要求提出新的解決方法。可以預見，在科學運算、自動控制與科學繪圖領域Matlab語言將長期保持其獨一無二的地位。為此，本例采用Matlab語言編程，以利用其快捷強大的矩陣數值計算功能。問題描述一矩形薄板，一邊固定，承受150kN集中荷載，結構簡圖及按平面三角形單元劃分的有限元模型圖如下所示。材料參數：彈性模量；泊松比：；薄板厚度。在本例中，所取結構模型及數據主要用于程序設計理論分析，與工程實際無關。參數輸入：單元個數NELEM=4節(jié)點個數NNODE=6

受約束邊界點數NVFIX=2節(jié)點荷載個數NFORCE=1彈性模量YOUNG=2e8泊松比POISS=厚度THICK=單元節(jié)點編碼數組LNODS=節(jié)點坐標數組COORD=節(jié)點力數組FORCE=[40-150]約束信息數組FIXED=以上數值數據為程序運行前輸入的初始數據，存為“”文本格式，同時必須放在Matlab工作目錄下，路徑不對程序不能自動讀取指定初始文件，運行出錯。初始數據文本文件輸入格式如下圖：Matlab語言程序源代碼：程序中變量說明NNODE單元節(jié)點數NPION總結點數NELEM單元數NVFIX受約束邊界點數FIXED約束信息數組NFORCE節(jié)點力數FORCE節(jié)點力數組COORD結構節(jié)點坐標數組LNODS單元定義數組YOUNG彈性模量POISS泊松比THICK厚度B單元應變矩陣(3*6)D單元彈性矩陣(3*3)S單元應力矩陣(3*6)A單元面積ESTIF單元剛度矩陣ASTIF總體剛度矩陣ASLOD總體荷載向量ASDISP節(jié)點位移向量ELEDISP單元節(jié)點位移向量STRESS單元應力FP1數據文件指針Matlab語言程序代碼%******************************************************************************%初始化及數據調用formatshorte %設定輸出類型clear %清除內存變量FP1=fopen('','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])'%結點坐標數組%坐標：x，y坐標（共NPOIN組）FORCE=fscanf(FP1,'%f',[3,NFORCE])'%結點力數組（受力結點編號,x方向,y方向）FIXED=fscanf(FP1,'%d',[3,NVFIX])'%約束信息（約束點，x約束，y約束）%有約束為1，無約束為0%*****************************************************************************%生成單元剛度矩陣并組成總體剛度矩陣ASTIF=zeros(2*NPION,2*NPION);%生成特定大小總體剛度矩陣并置0%*****************************************************************************fori=1:NELEM%生成彈性矩陣DD=[1POISS0;POISS10;00(1-POISS)/2]*YOUNG/(1-POISS^2)%*****************************************************************************%計算當前單元的面積A=det([1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)])/2%*****************************************************************************%生成應變矩陣Bforj=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)0b(2)0b(3)0;...0c(1)0c(2)0c(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é)點編號forj=1:3fork=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);%根據節(jié)點編號對應關系將單元剛度分塊疊加到總剛%度矩陣中endendend%*****************************************************************************%將約束信息加入總體剛度矩陣（對角元素改一法）fori=1:NVFIXifFIXED(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;%對角元素為1endifFIXED(i,3)==1ASTIF(:,FIXED(i,1)*2)=0;%一列為零ASTIF(FIXED(i,1)*2,:)=0;%一行為零ASTIF(FIXED(i,1)*2,FIXED(i,1)*2)=1;%對角元素為1endend%*****************************************************************************%生成荷載向量ASLOD(1:2*NPION)=0;%總體荷載向量置零fori=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%*****************************************************************************%求解內力ASDISP=ASTIF\ASLOD'%計算節(jié)點位移向量ELEDISP(1:6)=0;%當前單元節(jié)點位移向量fori=1:NELEMforj=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出當前單元的節(jié)點位移向量endiSTRESS=D*B1(:,:,i)*ELEDISP'%求內力endfclose(FP1);%關閉數據文件程序運行結果NELEM=4NPION=6NVFIX=2NFORCE=1YOUNG=0POISS=THICK=LNODS=126234245256COORD=001020211101FORCE=40-150FIXED=111611D=+008+0070+007+008000+007A=D=+008+0070+007+008000+007A=D=+008+0070+007+008000+007A=D=+008+0070+007+008000+007A=ASDISP=0000（說明：以上12個ASDISP輸出結果數據表示節(jié)點位移向量，即各節(jié)點分別在X，Y方向上產生的位移。單位：m）i=1STRESS=+005+004+005i=2STRESS=+004+004+004i=3STRESS=+004+004+004i=4STRESS=+005+004+004（說明：以上四組STRESS輸出結果數據表示單元應力SX,SY,SXY，i為單元號。單位：KN/m2）ANSYS建模比較分析利用ANSYS有限元分析軟件，完全按照前面Matlab程序設計的各項條件參數以及單元劃分方式建立ANSYS模型，其求解結果與以上程序計算結果比較，驗證程序是否正確。ANSYS模型節(jié)點單元如下圖所示：ANSYS模型求解變形圖如下所示：ANSYS求解節(jié)點位移結果列表顯示如下：（單位：m）PRINTUNODALSOLUTIONPERNODE*****POST1NODALDEGREEOFFREEDOMLISTING*****THEFOLLOWINGDEGREEOFFREEDOMRESULTSAREINTHEGLOBALCOORDINATESYSTEMNODEUXUYUZUSUM123456MAXIMUMABSOLUTEVALUESNODE4404VALUEANSYS求解單元應力結果列表顯示如下：（單位：KN/m2）PRINTSELEMENTSOLUTIONPERELEMENT*****POST1ELEMENTNODALSTRESSLISTING*****THEFOLLOWINGX,Y,ZVALUESAREINGLOBALCOORDINATESELEMENT=1PLANE182NODESXSYSZSXYSYZ1+06-33586.+062+06-33586.+0