基于Matlab語(yǔ)言的按平面三角形單元?jiǎng)澐值慕Y(jié)構(gòu)有限元程序設(shè)計(jì)模板

1、基于Matlab語(yǔ)言的按平面三角形單元?jiǎng)澐值慕Y(jié)構(gòu)有限元程序設(shè)計(jì)專業(yè)：建筑與土木工程班級(jí)：姓名：學(xué)號(hào):基于Matlab語(yǔ)言的按平面三角形單元?jiǎng)澐纸Y(jié)構(gòu)有限元程序設(shè)計(jì)有限單元發(fā)及Matlab語(yǔ)言概述1.有限單元法隨著現(xiàn)代工業(yè)、生產(chǎn)技術(shù)的發(fā)展，不斷要求設(shè)計(jì)高質(zhì)量、高水平的大型、復(fù)雜和精密的機(jī)械及工程結(jié)構(gòu)。為此目的，人們必須預(yù)先通過(guò)有效的計(jì)算手段，確切的預(yù)測(cè)即將誕生的機(jī)械和工程結(jié)構(gòu)，在未來(lái)工作時(shí)所發(fā)生的應(yīng)力、應(yīng)變和位移因此，需要尋求一種簡(jiǎn)單而又精確的數(shù)值分析方法。有限單元法正是適應(yīng)這種要求而產(chǎn)生和發(fā)展起來(lái)的一種十分有效的數(shù)值計(jì)算方法。有限元法把一個(gè)復(fù)雜的結(jié)構(gòu)分解成相對(duì)簡(jiǎn)單的“單元”，各單元之間通過(guò)結(jié)點(diǎn)相

2、互連接。單元內(nèi)的物理量由單元結(jié)點(diǎn)上的物理量按一定的假設(shè)內(nèi)插得到，這樣就把一個(gè)復(fù)雜結(jié)構(gòu)從無(wú)限多個(gè)自由度簡(jiǎn)化為有限個(gè)單元組成的結(jié)構(gòu)。我們只要分析每個(gè)單元的力學(xué)特性，然后按照有限元法的規(guī)則把這些單元“拼裝”成整體，就能夠得到整體結(jié)構(gòu)的力學(xué)特性。有限單元法基本步驟如下：(1)結(jié)構(gòu)離散：結(jié)構(gòu)離散就是建立結(jié)構(gòu)的有限元模型，又稱為網(wǎng)格劃分或單元?jiǎng)澐?，即將結(jié)構(gòu)離散為由有限個(gè)單元組成的有限元模型。在該步驟中，需要根據(jù)結(jié)構(gòu)的幾何特性、載荷情況等確定單元體內(nèi)任意一點(diǎn)的位移插值函數(shù)。(2)單元分析：根據(jù)彈性力學(xué)的幾何方程以及物理方程確定單元的剛度矩陣。(3)整體分析：把各個(gè)單元按原來(lái)的結(jié)構(gòu)重新連接起來(lái)，并在單元?jiǎng)偠染?/p>

3、陣的基礎(chǔ)上確定結(jié)構(gòu)的總剛度矩陣，形成如下式所示的整體有限元線性方程：1=K肪式中，作是載荷矩陣，k】是整體結(jié)構(gòu)的剛度矩陣，心是節(jié)點(diǎn)位移矩陣。(4)載荷移置：根據(jù)靜力等效原理，將載荷移置到相應(yīng)的節(jié)點(diǎn)上，形成節(jié)點(diǎn)載荷矩陣。(5)邊界條件處理：對(duì)式所示的有限元線性方程進(jìn)行邊界條件處理。(6)求解線性方程：求解式所示的有限元線性方程，得到節(jié)點(diǎn)的位移。在該步驟中，若有限元模型的節(jié)點(diǎn)越多，則線性方程的數(shù)量就越多，隨之有限元分析的計(jì)算量也將越大。(7)求解單元應(yīng)力及應(yīng)變根據(jù)求出的節(jié)點(diǎn)位移求解單元的應(yīng)力和應(yīng)變。(8)結(jié)果處理與顯示進(jìn)入有限元分析的后處理部分，對(duì)計(jì)算出來(lái)的結(jié)果進(jìn)行加工處理，并以各種形式將計(jì)算結(jié)果

4、顯示出。2.Matlab簡(jiǎn)介在用有限元法進(jìn)行結(jié)構(gòu)分析時(shí)，將會(huì)遇到大量的數(shù)值計(jì)算，因而在實(shí)用上是一定要借助于計(jì)算機(jī)和有限元程序，才能完成這些復(fù)雜而繁重的數(shù)值計(jì)算工作。而Matlab是當(dāng)今國(guó)際科學(xué)界最具影響力和活力的軟件。它起源于矩陣運(yùn)算，并已經(jīng)發(fā)展成一種高度集成的計(jì)算機(jī)語(yǔ)言。它提供了強(qiáng)大的科學(xué)計(jì)算，靈活的程序設(shè)計(jì)流程，高質(zhì)量的圖形可視化與界面設(shè)計(jì)，便捷的與其他程序和語(yǔ)言接口的功能。Matlab在各國(guó)高校與研究單位起著重大的作用?！肮び破涫拢叵壤淦鳌?。如果有一種十分有效的工具能解決在教學(xué)與研究中遇到的問(wèn)題，那么Matlab語(yǔ)言正是這樣的一種工具。它可以將使用者從繁瑣、無(wú)謂的底層編程中解放出來(lái)

5、，把有限的寶貴時(shí)間更多地花在解決問(wèn)題中，這樣無(wú)疑會(huì)提高工作效率。目前，Matlab已經(jīng)成為國(guó)際上最流行的科學(xué)與工程計(jì)算的軟件工具，現(xiàn)在的Matlab已經(jīng)不僅僅是一個(gè)“矩陣實(shí)驗(yàn)室”了，它已經(jīng)成為了一種具有廣泛應(yīng)用前景的全新的計(jì)算機(jī)高級(jí)編程語(yǔ)言了，有人稱它為“第四代”計(jì)算機(jī)語(yǔ)言，它在國(guó)內(nèi)外高校和研究部門(mén)正扮演著重要的角色。Matlab語(yǔ)言的功能也越來(lái)越強(qiáng)大，不斷適應(yīng)新的要求提出新的解決方法?？梢灶A(yù)見(jiàn)，在科學(xué)運(yùn)算、自動(dòng)控制與科學(xué)繪圖領(lǐng)域Matlab語(yǔ)言將長(zhǎng)期保持其獨(dú)一無(wú)二的地位。為此，本例采用Matlab語(yǔ)言編程，以利用其快捷強(qiáng)大的矩陣數(shù)值計(jì)算功能。二、問(wèn)題描述一矩形薄板，一邊固定，承受150kN集

6、中荷載，結(jié)構(gòu)簡(jiǎn)圖及按平面三角形單元?jiǎng)澐值挠邢拊Ｐ蛨D如下所示。材料參數(shù)：彈性模量E=2父108KN/m2;泊松比：卜=0.2;薄板厚度2mm。在本例中，所取結(jié)構(gòu)模型及數(shù)據(jù)主要用于程序設(shè)計(jì)理論分析，與工程實(shí)際無(wú)關(guān)。F=150kN三、參數(shù)輸入：?jiǎn)卧獋€(gè)數(shù)NELEM=4節(jié)點(diǎn)個(gè)數(shù)NNODE=6受約束邊界點(diǎn)數(shù)NVFIX=2節(jié)點(diǎn)荷載個(gè)數(shù)NFORCE=1彈性模量YOUNG=2e8泊松比POISS=0.2厚度THICK=0.002彳261234單元節(jié)點(diǎn)編碼數(shù)組LNODS=2452561節(jié)點(diǎn)坐標(biāo)數(shù)組COORD =-01 I2210 100111節(jié)點(diǎn)力數(shù)組FORCE =4 0-150 11約束信息數(shù)組FIXED =

7、 6 111 1以上數(shù)值數(shù)據(jù)為程序運(yùn)行前輸入的初始數(shù)據(jù)，存為471220580.txt ” 文本格式，同時(shí)必須放在Matlab工作目錄下，路徑不對(duì)程序不能自動(dòng)讀取指定初始文件，運(yùn)行出錯(cuò)。初始數(shù)據(jù)文本文件輸入格式如下圖：四、Matlab語(yǔ)言程序源代碼:1.程序中變量說(shuō)明NNODE單元節(jié)點(diǎn)數(shù)B單兀應(yīng)及矩陣(3*6)NPION總結(jié)點(diǎn)數(shù)D單元彈性矩陣(3*3)NELEM單元數(shù)S單元應(yīng)力矩陣(3*6)NVFIX受約束邊界點(diǎn)數(shù)A單兀回積FIXED約束信息數(shù)組ESTIF單元?jiǎng)偠染仃嘚FORCE下點(diǎn)力數(shù)ASTIF總體剛度矩陣FORCE節(jié)點(diǎn)力數(shù)組ASLOD總體荷載向量COORD結(jié)構(gòu)節(jié)點(diǎn)坐標(biāo)數(shù)組ASDISP節(jié)點(diǎn)位

8、移向量LNODS單元定義數(shù)組ELEDISP單元節(jié)點(diǎn)位移向量YOUNG彈性模量STRESS單元應(yīng)力POISS泊松比FP1數(shù)據(jù)文件指針THICK厚度2.Matlab語(yǔ)言程序代碼%*%初始化及數(shù)據(jù)調(diào)用format short eclearFP1=fopen( 471220580.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=

9、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) %設(shè)定輸出類(lèi)型%清除內(nèi)存變量打開(kāi)輸入數(shù)據(jù)文件，讀入控制數(shù)據(jù)單元個(gè)數(shù)（單元編碼總數(shù)）結(jié)點(diǎn)個(gè)數(shù)（結(jié)點(diǎn)編碼總數(shù)）受約束邊界點(diǎn)數(shù)結(jié)點(diǎn)荷載個(gè)數(shù)彈性模量泊松比厚度單元定義數(shù)組（單元結(jié)點(diǎn)號(hào)）%目應(yīng)為單元結(jié)點(diǎn)號(hào)（編碼）、按逆時(shí)針順序輸入結(jié)點(diǎn)坐標(biāo)數(shù)組%座標(biāo)：x,y坐標(biāo)（共NPOIN組）結(jié)點(diǎn)力數(shù)組（受力結(jié)點(diǎn)編號(hào)，x方向,y方向）約束信息（約束點(diǎn)

10、，x約束，y約束）有約束為1,無(wú)約束為0%*成單元?jiǎng)偠染仃嚥⒔M成總體剛度矩陣ASTIF=zeros(2*NPION,2*NPION);%生成特定大小總體剛度矩陣并置0%*fori=1:NELEM%生成彈性矩陣DD=1POISS0;POISS10;00(1-POISS)/2*YOUNG/(1-POISSA2)%*%計(jì)算當(dāng)前單元的面積AA=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%*%生成應(yīng)

11、變矩陣Bforj=0:2b(j+1)=COORDL(NODS(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;%*%求應(yīng)力矩陣S=D*BS=D*B;ESTIF=B*S*THICK*A;%求解單元?jiǎng)偠染仃嘺=LNODS(i,:);

12、%臨時(shí)向量,用來(lái)記錄當(dāng)前單元的節(jié)點(diǎn)編號(hào)forj=1:3fork=1:3ASTIF(aO)*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é)點(diǎn)編號(hào)對(duì)應(yīng)關(guān)系將單元?jiǎng)偠确謮K疊加到總剛%度矩陣中endendend%*%將約束信息加入總體剛度矩陣(對(duì)角元素改一法)fori=1:NVFIXifFIXED(i,2)=1ASTIF(:,(FIXED(i,1)*2-1)=0;%一列為零ASTIF(FIXED(i,1)*2-1),:)=0;%一行為零AST

13、IF(FIXED(i,1)*2-1),(FIXED(i,1)*2-1)=1;%對(duì)角元素為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;%對(duì)角元素為1endend%*%生成荷載向量ASLOD(1:2*NPION)=0;%總體荷載向量置零fori=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%*%求解內(nèi)力ASDISP=ASTIFASLOD%計(jì)算節(jié)點(diǎn)位

14、移向量ELEDISP(1:6)=0;%當(dāng)前單元節(jié)點(diǎn)位移向量fori=1:NELEMforj=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出當(dāng)前單元的節(jié)點(diǎn)位移向量endiSTRESS=D*B1(:,:,i)*ELEDISP%求內(nèi)力endfclose(FP1);%關(guān)閉數(shù)據(jù)文件五、程序運(yùn)行結(jié)果NELEM=4NPION=6NVFIX=2NFORCE=1YOUNG=200000000POISS=2.0000e-001THICK=2.0000e-003LNODS=1 262 34245256COORD=001 02 03 11 1

15、01FORCE=40-150FIXED=111611D=2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.3333e+007A=5.0000e-001D=2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.3333e+007A=5.0000e-001D=2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.3333e+007A=5.0000e-001D=2.0833e+0084.1667e+00704.1667e+0072.0833e+0080008.

16、3333e+007A=5.0000e-001ASDISP=00-8.0607e-004-1.5848e-003-1.0281e-003-4.4727e-0031.1937e-003-4.6947e-0038.6670e-004-1.8880e-00300（說(shuō)明：以上12個(gè)ASDISP輸出結(jié)果數(shù)據(jù)表示節(jié)點(diǎn)位移向量，即各節(jié)點(diǎn)分別在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+0

17、04-9.4497e+004i=4STRESS=1.6793e+005-2.7040e+004-1.7932e+004（說(shuō)明：以上四組STRESS出結(jié)果數(shù)據(jù)表示單元應(yīng)力SX,SY,SXYi為單元號(hào)。）六、ANSYS建模比較分析利用ANSYSt限元分析軟件，完全按照前面Matlab程序設(shè)計(jì)的各項(xiàng)條件參數(shù)以及單元?jiǎng)澐址绞浇NSYS1型，其求解結(jié)果與以上程序計(jì)算結(jié)果比較，驗(yàn)證程序是否正確。ANSYS模型節(jié)點(diǎn)單元如下圖所示:ELEMEMT5ELEMNUMANAPR201301:33:58PLOTNO.1ANAPR 26 2013 13:24:06PLOT NO. 1hamhiqiangANSYS模

18、型求解變形圖如下所示:DISELAOEMEldTSTEP-1SUB=1TIME=1DMX=.004844hanzhiqiangANSYS求解節(jié)點(diǎn)位移結(jié)果列表顯示如下：（單位：m）PRINTUNODALSOLUTIONPERNODE*POST1NODALDEGREEOFFREEDOMLISTING*THEFOLLOWINGDEGREEOFFREEDOMRESULTSAREINTHEGLOBALCOORDINATESYSTEMNODEUXUYUZUSUM10.00000.00000.00000.00000.17780E-020.45893E-020.48441E-020.20774E-022-0.

19、80607E-03-0.15848E-020.00003-0.10281E-02-0.44727E-020.000040.11937E-02-0.46947E-020.000050.86670E-03-0.18880E-020.000060.00000.00000.00000.0000MAXIMUMABSOLUTEVALUESNODE4404VALUE0.11937E-02-0.46947E-020.00000.48441E-02ANSYS求解單元應(yīng)力結(jié)果列表顯示如下:（單位：KN/m2）PRINTSELEMENTSOLUTIONPERELEMENT*POST1ELEMENTNODALSTRE

20、SSLISTING*THEFOLLOWINGX,Y,ZVALUESAREINGLOBALCOORDINATESELEMENT= 1NODE SX SYPLANE1821 -0.16793E+06 -33586.2 -0.16793E+06 -33586.6 -0.16793E+06 -33586.SZ : 0.0000 0.0000 0.0000SXYSYZ-0.13207E+06-0.13207E+06-0.13207E+060.00000.00000.0000NODESXSYSZSXYSYZ2-55503.-55503.0.0000-55503.0.00003-55503.-55503.0

21、.0000-55503.0.00004-55503.ELEMENT=-55503.0.00003PLANE182-55503.0.0000NODESXSYSZSXYSYZ255503.-49526.0.0000-94497.0.0000455503.-49526.0.0000-94497.0.0000555503.ELEMENT=-49526.0.00004PLANE182-94497.0.0000NODESXSYSZSXYSYZELEMENT=2PLANE18220.16793E+06-27040.0.0000-17932.0.000050.16793E+06-27040.0.0000-17

22、932.0.000060.16793E+06-27040.0.0000-17932.0.0000結(jié)論通過(guò)比較Matlab語(yǔ)言設(shè)計(jì)程序運(yùn)行結(jié)果和ANSYSt模分析結(jié)果可知，兩種方式算出的結(jié)果完全一致，說(shuō)程序設(shè)計(jì)正確。所以本程序適用于按三角形單元?jiǎng)澐值钠矫娼Y(jié)構(gòu)有限元分析。formatshorteclearFP1=fopen(LinearTriangleElementofThinplate.txt,rt);NELEM=fscanf(FP1,%d,1)NPION=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)YOUNG=fsc

23、anf(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);%生成特定大小總體剛度矩陣并置0fori=1:NELEM%生成彈性矩陣DD=YOUNG/(1-POISSA2)*1POISS0;POISS10;00(1-POISS)/2%計(jì)算當(dāng)前單元的面積AA=det

24、(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%生成應(yīng)變矩陣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)0

25、b(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;%求應(yīng)力矩陣S=D*BS=D*B;ESTIF=B*S*THICK*A;a=LNODS(i,:);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);endendend%將約束信息加入總體剛度矩陣fori=1:NVFIXifFIXED(i,2)=1ASTIF(:,(

26、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;endifFIXED(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;fori=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解內(nèi)力ASDISP=ASTIFA

27、SLODELEDISP(1:6)=0;fori=1:NELEMforj=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:,:,i)*ELEDISPendfclose(FP1);formatshorteclearFP1=fopen(LinearTriangleElementofThinplate.txt,rt);NELEM=fscanf(FP1,%d,1)NPION=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)YOUN

28、G=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);%生成特定大小總體剛度矩陣并置0fori=1:NELEM%生成彈性矩陣DD=YOUNG/(1-POISSA2)*1POISS0;POISS10;00(1-POISS)/2%計(jì)算當(dāng)前單元的面積AA=det(1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);1COORD(LNODS(i