基于Matlab語言的按平面三角形單元劃分的結(jié)構(gòu)有限元程序設(shè)計說明書

可編輯版13/13基于Matlab語言的按平面三角形單元劃分的結(jié)構(gòu)有限元程序設(shè)計專業(yè)：建筑與土木工程班級：建工研12-2姓名：韓志強學(xué)號：471220580基于Matlab語言的按平面三角形單元劃分結(jié)構(gòu)有限元程序設(shè)計有限單元發(fā)及Matlab語言概述1.有限單元法隨著現(xiàn)代工業(yè)、生產(chǎn)技術(shù)的發(fā)展,不斷要求設(shè)計高質(zhì)量、高水平的大型、復(fù)雜和精密的機械及工程結(jié)構(gòu)。為此目的,人們必須預(yù)先通過有效的計算手段,確切的預(yù)測即將誕生的機械和工程結(jié)構(gòu),在未來工作時所發(fā)生的應(yīng)力、應(yīng)變和位移因此,需要尋求一種簡單而又精確的數(shù)值分析方法。有限單元法正是適應(yīng)這種要求而產(chǎn)生和發(fā)展起來的一種十分有效的數(shù)值計算方法。有限元法把一個復(fù)雜的結(jié)構(gòu)分解成相對簡單的"單元",各單元之間通過結(jié)點相互連接。單元內(nèi)的物理量由單元結(jié)點上的物理量按一定的假設(shè)內(nèi)插得到,這樣就把一個復(fù)雜結(jié)構(gòu)從無限多個自由度簡化為有限個單元組成的結(jié)構(gòu)。我們只要分析每個單元的力學(xué)特性,然后按照有限元法的規(guī)則把這些單元"拼裝"成整體,就能夠得到整體結(jié)構(gòu)的力學(xué)特性。有限單元法基本步驟如下：<1>結(jié)構(gòu)離散：結(jié)構(gòu)離散就是建立結(jié)構(gòu)的有限元模型,又稱為網(wǎng)格劃分或單元劃分,即將結(jié)構(gòu)離散為由有限個單元組成的有限元模型。在該步驟中,需要根據(jù)結(jié)構(gòu)的幾何特性、載荷情況等確定單元體內(nèi)任意一點的位移插值函數(shù)。<2>單元分析：根據(jù)彈性力學(xué)的幾何方程以及物理方程確定單元的剛度矩陣。<3>整體分析:把各個單元按原來的結(jié)構(gòu)重新連接起來,并在單元剛度矩陣的基礎(chǔ)上確定結(jié)構(gòu)的總剛度矩陣,形成如下式所示的整體有限元線性方程：①式中,是載荷矩陣,是整體結(jié)構(gòu)的剛度矩陣,是節(jié)點位移矩陣。<4>載荷移置:根據(jù)靜力等效原理,將載荷移置到相應(yīng)的節(jié)點上,形成節(jié)點載荷矩陣。<5>邊界條件處理:對式①所示的有限元線性方程進行邊界條件處理。<6>求解線性方程:求解式①所示的有限元線性方程,得到節(jié)點的位移。在該步驟中,若有限元模型的節(jié)點越多,則線性方程的數(shù)量就越多,隨之有限元分析的計算量也將越大。<7>求解單元應(yīng)力及應(yīng)變根據(jù)求出的節(jié)點位移求解單元的應(yīng)力和應(yīng)變。<8>結(jié)果處理與顯示進入有限元分析的后處理部分,對計算出來的結(jié)果進行加工處理,并以各種形式將計算結(jié)果顯示出。2.Matlab簡介在用有限元法進行結(jié)構(gòu)分析時,將會遇到大量的數(shù)值計算,因而在實用上是一定要借助于計算機和有限元程序,才能完成這些復(fù)雜而繁重的數(shù)值計算工作。而Matlab是當(dāng)今國際科學(xué)界最具影響力和活力的軟件。它起源于矩陣運算,并已經(jīng)發(fā)展成一種高度集成的計算機語言。它提供了強大的科學(xué)計算,靈活的程序設(shè)計流程,高質(zhì)量的圖形可視化與界面設(shè)計,便捷的與其他程序和語言接口的功能。Matlab在各國高校與研究單位起著重大的作用。"工欲善其事,必先利其器"。如果有一種十分有效的工具能解決在教學(xué)與研究中遇到的問題,那么Matlab語言正是這樣的一種工具。它可以將使用者從繁瑣、無謂的底層編程中解放出來,把有限的寶貴時間更多地花在解決問題中,這樣無疑會提高工作效率。目前,Matlab已經(jīng)成為國際上最流行的科學(xué)與工程計算的軟件工具,現(xiàn)在的Matlab已經(jīng)不僅僅是一個"矩陣實驗室"了,它已經(jīng)成為了一種具有廣泛應(yīng)用前景的全新的計算機高級編程語言了,有人稱它為"第四代"計算機語言,它在國內(nèi)外高校和研究部門正扮演著重要的角色。Matlab語言的功能也越來越強大,不斷適應(yīng)新的要求提出新的解決方法?？梢灶A(yù)見,在科學(xué)運算、自動控制與科學(xué)繪圖領(lǐng)域Matlab語言將長期保持其獨一無二的地位。為此,本例采用Matlab語言編程,以利用其快捷強大的矩陣數(shù)值計算功能。問題描述一矩形薄板,一邊固定,承受150kN集中荷載,結(jié)構(gòu)簡圖及按平面三角形單元劃分的有限元模型圖如下所示。材料參數(shù)：彈性模量；泊松比：；薄板厚度。在本例中,所取結(jié)構(gòu)模型及數(shù)據(jù)主要用于程序設(shè)計理論分析,與工程實際無關(guān)。參數(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é)點坐標(biāo)數(shù)組COORD=節(jié)點力數(shù)組FORCE=[40-150]約束信息數(shù)組FIXED=以上數(shù)值數(shù)據(jù)為程序運行前輸入的初始數(shù)據(jù),存為"471220580.txt"文本格式,同時必須放在Matlab工作目錄下,路徑不對程序不能自動讀取指定初始文件,運行出錯。初始數(shù)據(jù)文本文件輸入格式如下圖：Matlab語言程序源代碼：程序中變量說明NNODE單元節(jié)點數(shù)NPION總結(jié)點數(shù)NELEM單元數(shù)NVFIX受約束邊界點數(shù)FIXED約束信息數(shù)組NFORCE節(jié)點力數(shù)FORCE節(jié)點力數(shù)組COORD結(jié)構(gòu)節(jié)點坐標(biāo)數(shù)組LNODS單元定義數(shù)組YOUNG彈性模量POISS泊松比THICK厚度B單元應(yīng)變矩陣<3*6>D單元彈性矩陣<3*3>S單元應(yīng)力矩陣<3*6>A單元面積ESTIF單元剛度矩陣ASTIF總體剛度矩陣ASLOD總體荷載向量ASDISP節(jié)點位移向量ELEDISP單元節(jié)點位移向量STRESS單元應(yīng)力FP1數(shù)據(jù)文件指針Matlab語言程序代碼%******************************************************************************%初始化及數(shù)據(jù)調(diào)用formatshorte %設(shè)定輸出類型clear %清除內(nèi)存變量FP1=fopen<'471220580.txt','rt'>;%打開輸入數(shù)據(jù)文件,讀入控制數(shù)據(jù)NELEM=fscanf<FP1,'%d',1>,%單元個數(shù)〔單元編碼總數(shù)NPION=fscanf<FP1,'%d',1>,%結(jié)點個數(shù)〔結(jié)點編碼總數(shù)NVFIX=fscanf<FP1,'%d',1>%受約束邊界點數(shù)NFORCE=fscanf<FP1,'%d',1>,%結(jié)點荷載個數(shù)YOUNG=fscanf<FP1,'%e',1>,%彈性模量POISS=fscanf<FP1,'%f',1>,%泊松比THICK=fscanf<FP1,'%f',1>%厚度LNODS=fscanf<FP1,'%d',[3,NELEM]>'%單元定義數(shù)組〔單元結(jié)點號%相應(yīng)為單元結(jié)點號〔編碼、按逆時針順序輸入COORD=fscanf<FP1,'%f',[2,NPION]>'%結(jié)點坐標(biāo)數(shù)組%坐標(biāo)：x,y坐標(biāo)〔共NPOIN組FORCE=fscanf<FP1,'%f',[3,NFORCE]>'%結(jié)點力數(shù)組〔受力結(jié)點編號,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>%*****************************************************************************%計算當(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)變矩陣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;%*****************************************************************************%求應(yīng)力矩陣S=D*BS=D*B;ESTIF=B'*S*THICK*A;%求解單元剛度矩陣a=LNODS<i,:>;%臨時向量,用來記錄當(dāng)前單元的節(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>;%根據(jù)節(jié)點編號對應(yīng)關(guān)系將單元剛度分塊疊加到總剛%度矩陣中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%*****************************************************************************%求解內(nèi)力ASDISP=ASTIF\ASLOD'%計算節(jié)點位移向量ELEDISP<1:6>=0;%當(dāng)前單元節(jié)點位移向量fori=1:NELEMforj=1:3ELEDISP<j*2-1:j*2>=ASDISP<LNODS<i,j>*2-1:LNODS<i,j>*2>;%取出當(dāng)前單元的節(jié)點位移向量endiSTRESS=D*B1<:,:,i>*ELEDISP'%求內(nèi)力endfclose<FP1>;%關(guān)閉數(shù)據(jù)文件程序運行結(jié)果NELEM=4NPION=6NVFIX=2NFORCE=1YOUNG=200000000POISS=2.0000e-001THICK=2.0000e-003LNODS=126234245256COORD=001020211101FORCE=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.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〔說明：以上12個ASDISP輸出結(jié)果數(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輸出結(jié)果數(shù)據(jù)表示單元應(yīng)力SX,SY,SXY,i為單元號。ANSYS建模比較分析利用ANSYS有限元分析軟件,完全按照前面Matlab程序設(shè)計的各項條件參數(shù)以及單元劃分方式建立ANSYS模型,其求解結(jié)果與以上程序計算結(jié)果比較,驗證程序是否正確。ANSYS模型節(jié)點單元如下圖所示：ANSYS模型求解變形圖如下所示：ANSYS求解節(jié)點位移結(jié)果列表顯示如下：〔單位：mPRINTUNODALSOLUTIONPERNODE*****POST1NODALDEGREEOFFREEDOMLISTING*****THEFOLLOWINGDEGREEOFFREEDOMRESULTSAREINTHEGLOBALCOORDINATESYSTEMNODEUXUYUZUSUM10.00000.00000.00000.00002-0.80607E-03-0.15848E-020.00000.17780E-023-0.10281E-02-0.44727E-020.00000.45893E-0240.11937E-02-0.46947E-020.00000.48441E-0250.86670E-03-0.18880E-020.00000.20774E-0260.00000.00000.00000.0000MAXIMUMABSOLUTEVALUESNODE4404VALUE0.11937E-02-0.46947E-020.00000.48441E-02ANSYS求解單元應(yīng)力結(jié)果列表顯示如下：〔單位：KN/m2PRINTSELEMENTSOLUTIONPERELEMENT*****POST1ELEMENTNODALSTRESSLISTING*****THEFOLLOWINGX,Y,ZVALUESAREINGLOBALCOORDINATESELEMENT=1PLANE182NODESXSYSZSXYSYZ1-0.16793E+06-33586.0.0000-0.13207E+060.00002-0.16793E+06-33586.0.0000-0.13207E+060.00006-0.16793E+06-33586.0.0000-0.13207E+060.0000ELEMENT=2PLANE182NODESXSYSZSXYSYZ2-55503.-55503.0.0000-55503.0.00003-55503.-55503.0.0000-55503.0.00004-55503.-55503.0.0000-55503.0.0000ELEMENT=3PLANE182NODESXSYSZSXYSYZ255503.-49526.0.0000-94497.0.0000455503.-49526.0.0000-94497.0.0000555503.-49526.0.0000-94497.0.0000ELEMENT=4PLANE182NODESXSYSZSXYSYZ20.16793E+06-27040.0.0000-17932.0.000050.16793E+06-27040.0.0000-17932.0.000060.16793E+06-27040.0.0000-17932.0.0000結(jié)論通過比較Matlab語言設(shè)計程序運行結(jié)果和ANSYS建模分析結(jié)果可知,兩種方式算出的結(jié)果完全一致,說程序設(shè)計正確。所以本程序適用于按三角形單元劃分的平面結(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=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>;%生成特定大小總體剛度矩陣并置0fori=1:NELEM%生成彈性矩陣DD=YOUNG/<1-POISS^2>*[1POISS0;POISS10;00<1-POISS>/2]%計算當(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)變矩陣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;%求應(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<:,<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=ASTIF\ASLOD'ELEDISP<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>*ELEDISP'endfclose<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>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>;%生成特定大小總體