平面三角形單元常應(yīng)變單元matlab程序的編制

PAGE 三角形常應(yīng)變單元程序的編制與使用有限元法是求解微分方程邊值問(wèn)題的一種通用數(shù)值方法，該方法是一種基于變分法（或變分里茲法）而發(fā)展起來(lái)的求解微分方程的數(shù)值計(jì)算方法，以計(jì)算機(jī)為手段，采用分片近似，進(jìn)而逼近整體的研究思想求解物理問(wèn)題。有限元分析的基本步驟可歸納為三大步：結(jié)構(gòu)離散、單元分析和整體分析。開(kāi)始輸入初始數(shù)據(jù)生成單剛集成總剛開(kāi)始輸入初始數(shù)據(jù)生成單剛集成總剛施加約束信息生成荷載向量邊界條件處理計(jì)算結(jié)點(diǎn)位移計(jì)算單元應(yīng)力計(jì)算結(jié)果整理結(jié)束Matlab語(yǔ)言是進(jìn)行矩陣運(yùn)算的強(qiáng)大工具，因此，用Matlab語(yǔ)言編寫有限元中平面問(wèn)題的程序有優(yōu)越性。本章將詳細(xì)介紹如何利用Matlab語(yǔ)言編制三角形常應(yīng)變單元的計(jì)算程序，程序流程圖見(jiàn)圖1。有限元法中三節(jié)點(diǎn)三角形分析結(jié)構(gòu)的步驟如下：整理原始數(shù)據(jù)，如材料性質(zhì)、荷載條件、約束條件等，離散結(jié)構(gòu)并進(jìn)行單元編碼、結(jié)點(diǎn)編碼、結(jié)點(diǎn)位移編碼、選取坐標(biāo)系。單元分析，建立單元?jiǎng)偠染仃?。整體分析，建立總剛矩陣。建立整體結(jié)構(gòu)的等效節(jié)點(diǎn)荷載和總荷載矩陣邊界條件處理。解方程，求出節(jié)點(diǎn)位移。求出各單元的單元應(yīng)力。計(jì)算結(jié)果整理。計(jì)算結(jié)果整理包括位移和應(yīng)力兩個(gè)方面；位移計(jì)算結(jié)果一般不需要特別的處理，利用計(jì)算出的節(jié)點(diǎn)位移分量，就可畫出結(jié)構(gòu)任意方向的位移云圖；而應(yīng)力解的圖1程序流程圖誤差表現(xiàn)在單元內(nèi)部不滿足平衡方程，單元圖1程序流程圖與單元邊界處應(yīng)力一般不連續(xù)，在邊界上應(yīng)力解一般與力的邊界條件不相符合。1.1程序說(shuō)明%*******************************************************************%三角形常應(yīng)變單元求解結(jié)構(gòu)主程序%*******************************************************************功能：運(yùn)用有限元法中三角形常應(yīng)變單元解平面問(wèn)題的計(jì)算主程序?；舅枷耄?jiǎn)卧Y(jié)點(diǎn)按右手法則順序編號(hào)。荷載類型：可計(jì)算結(jié)點(diǎn)荷載。說(shuō)明：主程序的作用是通過(guò)賦值語(yǔ)句、讀取和寫入文件、函數(shù)調(diào)用等完成算法的全過(guò)程，即實(shí)現(xiàn)程序流程圖的程序表達(dá)。%1程序準(zhǔn)備 formatshorte %設(shè)定輸出類型 clearall %清除所有已定義變量clc %清屏說(shuō)明：formatshorte －設(shè)定計(jì)算過(guò)程中顯示在屏幕上的數(shù)字類型為短格式、科學(xué)計(jì)數(shù)法；clearall －清除所有已定義變量，目的是在本程序的運(yùn)行過(guò)程中，不會(huì)發(fā)生變量名相同等可能使計(jì)算出錯(cuò)的情況；clc －清屏，使屏幕在本程序運(yùn)行開(kāi)始時(shí)%2全局變量定義globalNNODENPIONNELEMNVFIXNFORCECOORDLNODSYOUNGPOISSTHICKglobalFORCEFIXEDglobalBMATXDMATXSMATXAREAglobalASTIFASLODASDISPglobalFP1說(shuō)明：NNODE—單元結(jié)點(diǎn)數(shù)，NPION—總結(jié)點(diǎn)數(shù)，NELEM—單元數(shù)，NVFIX—受約束邊界點(diǎn)數(shù)，NFORCE—結(jié)點(diǎn)力數(shù)，COORD—結(jié)構(gòu)結(jié)點(diǎn)坐標(biāo)數(shù)組，LNODS—單元定義數(shù)組，YOUNG—彈性模量，POISS—泊松比，THICK—厚度第a(k)*2-1到a(k)*2列的元素由單剛中第j*2-1到j(luò)*2行，第k*2-1到k*2列的元素疊加而得，a(j)*2即將單元中的位移編碼對(duì)應(yīng)到整體的位移編碼。%%將約束信息加入總剛（置0置1法） NUM=1;%計(jì)數(shù)器(當(dāng)前已分析的節(jié)點(diǎn)數(shù))NVFIX：受約束邊界點(diǎn)NVFIX：受約束邊界點(diǎn)數(shù)tmp(NVFIX)=0;%臨時(shí)存被約束的序號(hào) whilei<NVFIXforj=-1:0ifFIXED(NUM,j+3)==1%若發(fā)現(xiàn)約束temporaryi=i+1;%計(jì)數(shù)器自增temporarytmp(i)=FIXED(NUM)*2+j;%求約束序號(hào)endendNUM=NUM+1; end說(shuō)明：tmp(NVFIX)=0，形成一個(gè)元素值均為0的一行，NVFIX列的行向量，執(zhí)行while語(yǔ)句，首先判斷i是否小于控制數(shù)據(jù)NVFIX，若小于則往下進(jìn)行，若不小于則退出。執(zhí)行for語(yǔ)句，F(xiàn)IXED(NUM,j+3)表示約束信息數(shù)組中第NUM行，第j+3列的元素，j從-1到0，即j+3表示2到3列，即約束信息數(shù)組中描述結(jié)點(diǎn)x和y方向受約束的情況，判斷FIXED(NUM,j+3)若等于1，則約束數(shù)自增，若不等于1則跳出。FIXED(NUM)表示FIXED(NUM，1)，tmp(i)=FIXED(NUM)*2+j計(jì)算整體約束序號(hào)，將序號(hào)放入tmp行向量中的i列。% fori=1:NVFIXASTIF(1:2*NPION,tmp(i))=0;%將一約束序號(hào)處的總剛列向量清0ASTIF(tmp(i),1:2*NPION)=0;%.將一約束序號(hào)處的總剛行向量清0ASTIF(tmp(i),tmp(i))=1;%將行列交叉處的元素置為1end說(shuō)明：后處理法中置0置1法，設(shè)（包括），則將總剛中的主元素Kjj換為1，j行和j列的其他元素均改為零。%*******************************************************************%讀取FORMSMATX子程%*******************************************************************functionFORMSMATX(ELEMENT)%計(jì)算應(yīng)力矩陣 %引用所需的全局變量 globalNPIONNELEMCOORDLNODSYOUNGPOISS globalBMATXDMATXSMATXAREA% %生成彈性矩陣D a=YOUNG/(1-POISS^2); DMATX(1,1)=1*a;DMATX(1,2)=POISS*a;DMATX(2,1)=POISS*a;DMATX(2,2)=1*a;DMATX(3,3)=(1-POISS)*a/2;說(shuō)明：平面應(yīng)力問(wèn)題的彈性矩陣,其中，E為彈性模量，為泊松比。% i=ELEMENT;%i為當(dāng)前所計(jì)算的單元號(hào)%計(jì)算當(dāng)前單元的面積AREA=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;說(shuō)明：det表示求矩陣行列式的值，，其中分別表示一個(gè)三角形單元的三個(gè)節(jié)點(diǎn)坐標(biāo)。MATLAB中若一行中無(wú)法寫下一個(gè)完整的命令，則可以在行尾加入3個(gè)連續(xù)的英文句號(hào)，表示命令余下的部分在下一行出現(xiàn)。%%生成應(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); endBMATX=[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*AREA);說(shuō)明：應(yīng)變矩陣rem表示求余函數(shù)，rem（x，y）命令返回的是x-n.*y，當(dāng)y0時(shí)，n=fix(x./y)，其中fix為最近的整數(shù)取整。% SMATX=DMATX*BMATX;%求應(yīng)力矩陣S=D*B%*******************************************************************%讀取FORMLOAD子程%*******************************************************************functionFORMLOAD()%本子程生成荷載向量%%所需引用的全局變量globalASLODNPIONNFORCEFORCE%ASLOD(1:2*NPION)=0;%張成特定大小的0向量說(shuō)明：ASLOD為總體荷載向量，每個(gè)結(jié)點(diǎn)有x，y兩個(gè)方向的結(jié)點(diǎn)力。%fori=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end說(shuō)明：FORCE(i,1)為作用點(diǎn),FORCE(i,2:3)為x,y方向的結(jié)點(diǎn)力。%%將有約束處的荷載置為0NUM=1;i=0;tmp(NVFIX)=0; whilei<NVFIXforj=-1:0ifFIXED(NUM,j+3)==1i=i+1;tmp(i)=FIXED(NUM)*2+j;endendNUM=NUM+1; end fori=1:NVFIXASLOD(tmp(i))=0;end說(shuō)明：置0置1法，同上。%*******************************************************************ASDISP=ASTIF\ASLOD'%計(jì)算結(jié)點(diǎn)位移向量%*******************************************************************%讀取WRITESTRESS子程%*******************************************************************functionWRITESTRESS()%求解內(nèi)力，即兩個(gè)方向的正應(yīng)力與一個(gè)剪應(yīng)力（)%%所引用的全局變量globalNELEMNPIONSMATXASDISPLNODS%ELEDISP(1:6)=0;%當(dāng)前單元的結(jié)點(diǎn)位移向量說(shuō)明：ELEDIS為單元的結(jié)點(diǎn)位移%fori=1:NELEMforj=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endFORMSMATX(i)%%%調(diào)用子程求應(yīng)力矩陣iSTRESS=SMATX*ELEDISP'%求內(nèi)力end說(shuō)明：通過(guò)循環(huán)，依次從結(jié)構(gòu)位移列陣中對(duì)號(hào)，賦值給第i個(gè)單元的單元位移向量ELEDISP。%1.2程序應(yīng)用舉例【例1】利用三角形三節(jié)點(diǎn)常應(yīng)變單元程序計(jì)算圖2所示的結(jié)構(gòu)。設(shè)彈性模量，泊松比為0，厚度為1。%輸入數(shù)據(jù)文件input.txt為：64161.0e00.01312524325635圖20.02.00.01.01.01.00.00.01.00.02.00.010-1110210411501601%說(shuō)明：第一行：讀入程序控制信息NPION=fscanf(FP1,'%d',1)%結(jié)點(diǎn)個(gè)數(shù)（結(jié)點(diǎn)編碼總數(shù)）NELEM=fscanf(FP1,'%d',1)%單元個(gè)數(shù)（單元編碼總數(shù)）NFORCE=fscanf(FP1,'%d',1)％結(jié)點(diǎn)荷載個(gè)數(shù)NVFIX=fscanf(FP1,'%d',1)％受約束邊界點(diǎn)數(shù)YOUNG=fscanf(FP1,'%e',1)％彈性模量POISS=fscanf(FP1,'%f',1)％泊松比THICK=fscanf(FP1,'%d',1)%厚度第二、三、四、五行：讀入單元連接信息：LNODS=fscanf(FP1,'%d',[3,NELEM])';%單元定義數(shù)組，單元結(jié)點(diǎn)號(hào)，逆時(shí)針輸入第六、七、八、九、十、十一行：讀入結(jié)點(diǎn)坐標(biāo)COORD=fscanf(FP1,'%f',[2，NPION])';%結(jié)點(diǎn)坐標(biāo)值，第1列為