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

PAGE 三角形常應(yīng)變單元程序的編制與使用有限元法是求解微分方程邊值問題的一種通用數(shù)值方法，該方法是一種基于變分法（或變分里茲法）而發(fā)展起來的求解微分方程的數(shù)值計(jì)算方法，以計(jì)算機(jī)為手段，采用分片近似，進(jìn)而逼近整體的研究思想求解物理問題。有限元分析的基本步驟可歸納為三大步：結(jié)構(gòu)離散、單元分析和整體分析。開始輸入初始數(shù)據(jù)生成單剛集成總剛開始輸入初始數(shù)據(jù)生成單剛集成總剛施加約束信息生成荷載向量邊界條件處理計(jì)算結(jié)點(diǎn)位移計(jì)算單元應(yīng)力計(jì)算結(jié)果整理結(jié)束Matlab語(yǔ)言是進(jìn)行矩陣運(yùn)算的強(qiáng)大工具，因此，用Matlab語(yǔ)言編寫有限元中平面問題的程序有優(yōu)越性。本章將詳細(xì)介紹如何利用Matlab語(yǔ)言編制三角形常應(yīng)變單元的計(jì)算程序，程序流程圖見圖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程序說明%*******************************************************************%三角形常應(yīng)變單元求解結(jié)構(gòu)主程序%*******************************************************************功能：運(yùn)用有限元法中三角形常應(yīng)變單元解平面問題的計(jì)算主程序?；舅枷耄?jiǎn)卧Y(jié)點(diǎn)按右手法則順序編號(hào)。荷載類型：可計(jì)算結(jié)點(diǎn)荷載。說明：主程序的作用是通過賦值語(yǔ)句、讀取和寫入文件、函數(shù)調(diào)用等完成算法的全過程，即實(shí)現(xiàn)程序流程圖的程序表達(dá)。%1程序準(zhǔn)備 formatshorte %設(shè)定輸出類型 clearall %清除所有已定義變量clc %清屏說明：formatshorte －設(shè)定計(jì)算過程中顯示在屏幕上的數(shù)字類型為短格式、科學(xué)計(jì)數(shù)法；clearall －清除所有已定義變量，目的是在本程序的運(yùn)行過程中，不會(huì)發(fā)生變量名相同等可能使計(jì)算出錯(cuò)的情況；clc －清屏，使屏幕在本程序運(yùn)行開始時(shí)%2全局變量定義globalNNODENPIONNELEMNVFIXNFORCECOORDLNODSYOUNGPOISSTHICKglobalFORCEFIXEDglobalBMATXDMATXSMATXAREAglobalASTIFASLODASDISPglobalFP1說明：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說明：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說明：后處理法中置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;說明：平面應(yīng)力問題的彈性矩陣,其中，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;說明：det表示求矩陣行列式的值，，其中分別表示一個(gè)三角形單元的三個(gè)節(jié)點(diǎn)坐標(biāo)。MATLAB中若一行中無法寫下一個(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);說明：應(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向量說明：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說明：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說明：置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)位移向量說明：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說明：通過循環(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%說明：第一行：讀入程序控制信息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列為