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

1、三角形常應(yīng)變單元程序的編制與使用有限元法是求解微分方程邊值問題的一種通用數(shù)值方法，該方法是一種基于變分法（或變分里茲法）而發(fā)展起來的求解微分方程的數(shù)值計(jì)算方法，以計(jì)算機(jī)為手段，采用分片近似，進(jìn)而逼近整體的研究思想求解物理問題。有限元分析的基本步驟可歸納為三大步：結(jié)構(gòu)離散、單元分析和整體分析。開始輸入初始數(shù)據(jù)生成單剛集成總剛施加約束信息生成荷載向量邊界條件處理計(jì)算結(jié)點(diǎn)位移計(jì)算單元應(yīng)力計(jì)算結(jié)果整理結(jié)束對于平面問題，結(jié)構(gòu)離散常用的網(wǎng)格形狀有三角形、矩形、任意四邊形，以三個頂點(diǎn)為節(jié)點(diǎn)的三角形單元是最簡單的平面單元，它較矩形或四邊形對曲邊邊界有更好的適應(yīng)性，而矩形或四邊形單元較三節(jié)點(diǎn)三角形有更高的計(jì)算精

2、度。Matlab語言是進(jìn)行矩陣運(yùn)算的強(qiáng)大工具，因此，用Matlab語言編寫有限元中平面問題的程序有優(yōu)越性。本章將詳細(xì)介紹如何利用Matlab語言編制三角形常應(yīng)變單元的計(jì)算程序，程序流程圖見圖1。有限元法中三節(jié)點(diǎn)三角形分析結(jié)構(gòu)的步驟如下：1） 整理原始數(shù)據(jù)，如材料性質(zhì)、荷載條件、約束條件等，離散結(jié)構(gòu)并進(jìn)行單元編碼、結(jié)點(diǎn)編碼、結(jié)點(diǎn)位移編碼、選取坐標(biāo)系。2） 單元分析，建立單元剛度矩陣。3） 整體分析，建立總剛矩陣。4） 建立整體結(jié)構(gòu)的等效節(jié)點(diǎn)荷載和總荷載矩陣5） 邊界條件處理。6） 解方程，求出節(jié)點(diǎn)位移。7） 求出各單元的單元應(yīng)力。8） 計(jì)算結(jié)果整理。計(jì)算結(jié)果整理包括位移和應(yīng)力兩個方面；位移計(jì)算結(jié)

3、果一般不需要特別的處理，利用計(jì)算出的節(jié)點(diǎn)位移分量，就可畫出結(jié)構(gòu)任意方向的位移云圖；而應(yīng)力解的圖1 程序流程圖誤差表現(xiàn)在單元內(nèi)部不滿足平衡方程，單元與單元邊界處應(yīng)力一般不連續(xù)，在邊界上應(yīng)力解一般與力的邊界條件不相符合。1.1 程序說明%*% 三角形常應(yīng)變單元求解結(jié)構(gòu)主程序%* l 功能：運(yùn)用有限元法中三角形常應(yīng)變單元解平面問題的計(jì)算主程序。l 基本思想：單元結(jié)點(diǎn)按右手法則順序編號。l 荷載類型：可計(jì)算結(jié)點(diǎn)荷載。l 說明：主程序的作用是通過賦值語句、讀取和寫入文件、函數(shù)調(diào)用等完成算法的全過程，即實(shí)現(xiàn)程序流程圖的程序表達(dá)。%-1 程序準(zhǔn)備format short e %設(shè)定輸出類型clear all

4、%清除所有已定義變量clc%清屏l 說明：format short e 設(shè)定計(jì)算過程中顯示在屏幕上的數(shù)字類型為短格式、科學(xué)計(jì)數(shù)法； clear all 清除所有已定義變量，目的是在本程序的運(yùn)行過程中，不會發(fā)生變量名相同等可能使計(jì)算出錯的情況； clc 清屏，使屏幕在本程序運(yùn)行開始時%-2 全局變量定義global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXED global BMATX DMATX SMATX AREA global ASTIF ASLOD ASDISPglobal F

5、P1l 說明：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厚度 FORCE 節(jié)點(diǎn)力數(shù)組(n,3) n:受力節(jié)點(diǎn)數(shù)目,(n,1):作用點(diǎn),(n,2):x方向,(n,3):y方向； FIXED 約束信息數(shù)組(n,3) n:受約束節(jié)點(diǎn)數(shù)目, (n,1):約束點(diǎn) (n,2)與(n,3)分別為約束點(diǎn)x方向和y方向的約束情況,受約束為1否則為0 BMATX單元應(yīng)變矩陣(3*6)， DMATX單元彈性矩陣(3*3)，SMATX單元應(yīng)力矩陣(3*

6、6)，AREA單元面積 ASTIF總體剛度矩陣，ASLOD總體荷載向量，ASDISP結(jié)點(diǎn)位移向量FP1數(shù)據(jù)文件指針3 打開文件FP1=fopen('input.txt','rt'); %打開輸入數(shù)據(jù)文件 存放初始數(shù)據(jù)l 說明：FP1=fopen('input.txt','rt'); 打開已存在的輸入數(shù)據(jù)文件input.txt，且設(shè)置其為只讀格式，使程序在執(zhí)行過程中不能改變輸入文件中的數(shù)值，并用文件句柄FP1來執(zhí)行FP2=fopen('output.txt','wt');打開輸出數(shù)據(jù)文件，該文件不存在

7、時，通過此命令創(chuàng)建新文件，該文件存在時則將原有內(nèi)容全部刪除。該文件設(shè)置為可寫格式，可在程序執(zhí)行過程中向輸出文件寫入數(shù)據(jù)。%-4 讀入程序控制信息NPION=fscanf(FP1,'%d',1) %結(jié)點(diǎn)個數(shù)（結(jié)點(diǎn)編碼總數(shù)）NELEM=fscanf(FP1,'%d',1) %單元個數(shù)（單元編碼總數(shù)）NFORCE=fscanf(FP1,'%d',1) 結(jié)點(diǎn)荷載個數(shù)NVFIX=fscanf(FP1,'%d',1) YOUNG=fscanf(FP1,'%e',1) 彈性模量POISS=fscanf(FP1,'%f&#

8、39;,1) 泊松比 THICK=fscanf(FP1,'%d',1) %厚度LNODS=fscanf(FP1,'%d',3,NELEM)' %單元定義數(shù)組（單元結(jié)點(diǎn)號）l 說明：建立LNODS矩陣，該矩陣指出了每一單元的連接信息。矩陣的每一行針對每一單元，共計(jì) NELEM；每一列相應(yīng)為單元結(jié)點(diǎn)號（編碼）、按逆時針順序輸入。命令中，3,NELEM表示讀取NELEM行3列數(shù)據(jù)賦值給LNODS矩陣。顯然，LNODS(i,1:3)依次表示i單元的i，j，k結(jié)點(diǎn)號。COORD=fscanf(FP1,'%f',2,NPION)' %結(jié)點(diǎn)坐標(biāo)

9、數(shù)組l 說明：建立COORD矩陣，該矩陣用來存儲各結(jié)點(diǎn)x，y方向的坐標(biāo)值。從FP1文件中讀取全部結(jié)點(diǎn)個數(shù)NPOIN的坐標(biāo)值，從1開始按順序讀取。COORD(i，1：2)表示第i個結(jié)點(diǎn)的x，y坐標(biāo)。FORCE=fscanf(FP1,'%f',3,NFORCE)' %結(jié)點(diǎn)力數(shù)組l 說明： (n,3) n:受力結(jié)點(diǎn)數(shù)目,(n,1):作用點(diǎn),(n,2):x方向,(n,3):y方向FIXED=fscanf(FP1,'%d',3,inf)' %約束信息數(shù)組l 說明： (n,3) n:受約束節(jié)點(diǎn)數(shù)目, (n,1):約束點(diǎn) (n,2)與(n,3)分別為約束點(diǎn)x方

10、向和y方向的約束情況,受約束為1否則為0l 總體說明：從輸入文件FP1中讀入結(jié)點(diǎn)個數(shù)，單元個數(shù)，結(jié)點(diǎn)荷載個數(shù)，受約束邊界點(diǎn)數(shù)，彈性模量，泊松比，厚度，單元定義數(shù)組，結(jié)點(diǎn)坐標(biāo)數(shù)組，結(jié)點(diǎn)力數(shù)組，約束信息數(shù)組；程序中彈性模量僅輸入了一個值，表明本程序僅能求解一種材料構(gòu)成的結(jié)構(gòu)，如：鋼筋混凝土結(jié)構(gòu)、鋼結(jié)構(gòu)，不能求解鋼筋混凝土鋼組合結(jié)構(gòu)。infinite采用了命令fscanf，其中：%d表示讀入整數(shù)格式，%f'表示讀入浮點(diǎn)數(shù)；1表示讀取1個數(shù)，A,B形式表示讀A行B列數(shù)組，A,B表示將A,B轉(zhuǎn)置，inf表示正無窮。%-5 調(diào)用子程生成單剛，組成總剛并加入約束信息 ASSEMBLE()%-6 調(diào)用

11、子程生成荷載向量 FORMLOAD()%-7 計(jì)算結(jié)點(diǎn)位移向量 ASDISP=ASTIFASLOD'%-8調(diào)用子程計(jì)算單元應(yīng)力 WRITESTRESS()%*9 關(guān)閉輸出數(shù)據(jù)文件 fclose(FP2);%*讀取ASSEMBLE 子程%*function ASSEMBLE()% 所引用的全局變量：global NPION NELEM NVFIX LNODS ASTIF THICK global BMATX SMATX AREA FIXED%- -% 計(jì)算單剛并生成總剛 ASTIF(1:2*NPION,1:2*NPION)=0; %張成特定大小總剛矩陣并置0l 說明：Array stif

12、fness建立單元剛度矩陣ASTIF，該矩陣的行列數(shù)均為2*NPION ，NPION表示結(jié)點(diǎn)數(shù)，每個結(jié)點(diǎn)有兩個方向的力和位移。%- for i=1:NELEM FORMSMATX(i) %調(diào)用應(yīng)力子程序 ESTIF=BMATX'*SMATX*THICK*AREA; %求解單元剛度矩陣 a=LNODS(i,:); %臨時向量,用來記錄當(dāng)前單元的節(jié)點(diǎn)編號 for j=1:3 for k=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,

13、k*2-1:k*2); %跟據(jù)節(jié)點(diǎn)編號對應(yīng)關(guān)系將單元剛度分塊疊加到總剛矩陣中 end endendl 說明：FORMSMATX(i)調(diào)用應(yīng)力子程序，提取i單元的應(yīng)力矩陣SMATX；a=LNODS(i,:)記錄i單元的三個結(jié)點(diǎn)編號；forend循環(huán)語句表示行從1到3循環(huán)，列從1到3循環(huán)，將單剛中的元素疊加至總剛中：ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)表示總剛中第a(j)*2-1到:a(j)*2行，第a(k)*2-1到a(k)*2列的元素由單剛中第j*2-1到j(luò)*2行，第k*2-1到k*2列的元素疊加而得，a(j)*2即將單元中的位移編碼對應(yīng)到整體的位

14、移編碼。%-%將約束信息加入總剛（置0置1法）NUM=1; %計(jì)數(shù)器(當(dāng)前已分析的節(jié)點(diǎn)數(shù))NVFIX：受約束邊界點(diǎn)數(shù) i=0; %計(jì)數(shù)器(當(dāng)前已處理的約束數(shù)) tmp(NVFIX)=0; %臨時存被約束的序號while i<NVFIX for j=-1:0 if FIXED(NUM,j+3)=1 %若發(fā)現(xiàn)約束temporary i=i+1; %計(jì)數(shù)器自增 tmp(i)=FIXED(NUM)*2+j; %求約束序號 end end NUM=NUM+1;endl 說明： tmp(NVFIX)=0，形成一個元素值均為0的一行，NVFIX列的行向量，執(zhí)行while語句，首先判斷i是否小于控制數(shù)據(jù)

15、NVFIX，若小于則往下進(jìn)行，若不小于則退出。執(zhí)行for語句，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ì)算整體約束序號，將序號放入tmp行向量中的i列。%-for i=1:NVFIX ASTIF(1:2*NPION,tmp(i)=0; %將一約束序號處的總剛列向量清0 ASTIF(tmp(i),1:2*NPION

16、)=0; %.將一約束序號處的總剛行向量清0 ASTIF(tmp(i),tmp(i)=1; %將行列交叉處的元素置為1 endl 說明：后處理法中置0置1法，設(shè)（包括），則將總剛中的主元素K jj換為1，j行和j列的其他元素均改為零。%*% 讀取FORMSMATX子程%* function FORMSMATX(ELEMENT) %計(jì)算應(yīng)力矩陣 %引用所需的全局變量global NPION NELEM COORD LNODS YOUNG POISSglobal BMATX DMATX SMATX AREA%-%生成彈性矩陣D a=YOUNG/(1-POISS2); DMATX(1,1)=1*a;

17、DMATX(1,2)=POISS*a;DMATX(2,1)=POISS*a;DMATX(2,2)=1*a;DMATX(3,3)=(1-POISS)*a/2;l 說明： 平面應(yīng)力問題的彈性矩陣,其中，E為彈性模量，為泊松比。%-i=ELEMENT; %i為當(dāng)前所計(jì)算的單元號 %計(jì)算當(dāng)前單元的面積 AREA=det(1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);. 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);. 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);)/2;l 說明：det表

18、示求矩陣行列式的值，其中分別表示一個三角形單元的三個節(jié)點(diǎn)坐標(biāo)。MATLAB中若一行中無法寫下一個完整的命令，則可以在行尾加入3個連續(xù)的英文句號，表示命令余下的部分在下一行出現(xiàn)。%- %生成應(yīng)變矩陣B for j=0:2 b(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);end BMATX=b(1) 0 b(2) 0 b(3) 0;. 0 c(1) 0

19、 c(2) 0 c(3);. c(1) b(1) c(2) b(2) c(3) b(3)/(2*AREA);l 說明：應(yīng)變矩陣rem表示求余函數(shù)，rem（x，y）命令返回的是x-n.*y，當(dāng)y0時，n=fix(x./y)，其中fix為最近的整數(shù)取整。%- SMATX=DMATX*BMATX; %求應(yīng)力矩陣S=D*B%*% 讀取FORMLOAD子程%* function FORMLOAD() %本子程生成荷載向量%-%所需引用的全局變量global ASLOD NPION NFORCE FORCE%-ASLOD(1:2*NPION)=0; %張成特定大小的0向量l 說明：ASLOD為總體荷載向量

20、，每個結(jié)點(diǎn)有x，y兩個方向的結(jié)點(diǎn)力。%- for i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); endl 說明：FORCE(i,1)為作用點(diǎn),FORCE(i,2:3)為x,y方向的結(jié)點(diǎn)力。%-%將有約束處的荷載置為0 NUM=1; i=0; tmp(NVFIX)=0; while i<NVFIX for j=-1:0 if FIXED(NUM,j+3)=1 i=i+1; tmp(i)=FIXED(NUM)*2+j; end end NUM=NUM+1; endfor i=1:NVFIX ASLOD(tmp(i)=

21、0; endl 說明：置0置1法，同上。%*ASDISP=ASTIFASLOD' %計(jì)算結(jié)點(diǎn)位移向量%*% 讀取WRITESTRESS子程%* function WRITESTRESS()%求解內(nèi)力，即兩個方向的正應(yīng)力與一個剪應(yīng)力（)%- %所引用的全局變量 global NELEM NPION SMATX ASDISP LNODS%-ELEDISP(1:6)=0; %當(dāng)前單元的結(jié)點(diǎn)位移向量l 說明：ELEDIS為單元的結(jié)點(diǎn)位移%- for i=1:NELEM for j=1:3 ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2

22、); end FORMSMATX(i) % %調(diào)用子程求應(yīng)力矩陣 i STRESS=SMATX*ELEDISP' %求內(nèi)力 endl 說明：通過循環(huán)，依次從結(jié)構(gòu)位移列陣中對號，賦值給第i個單元的單元位移向量ELEDISP。%- 1.2 程序應(yīng)用舉例【例1】利用三角形三節(jié)點(diǎn)常應(yīng)變單元程序計(jì)算圖2所示的結(jié)構(gòu)。設(shè)彈性模量，泊松比為0，厚度為1。 %-輸入數(shù)據(jù)文件input.txt為：6 4 1 6 1.0e0 0.0 13 1 25 2 43 2 5 6 3 5 圖2 0.0 2.00.0 1.01.0 1.00.0 0.01.0 0.02.0 0.01 0 -11 1 02 1 04 1

23、15 0 16 0 1%-說明：第一行：讀入程序控制信息NPION=fscanf(FP1,'%d',1) %結(jié)點(diǎn)個數(shù)（結(jié)點(diǎn)編碼總數(shù)）NELEM=fscanf(FP1,'%d',1) %單元個數(shù)（單元編碼總數(shù)）NFORCE=fscanf(FP1,'%d',1) 結(jié)點(diǎn)荷載個數(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