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

1、圖1程序流程圖三角形常應(yīng)變單元程序的編制與使用有限元法是求解微分方程邊值問題的一種通用數(shù)值方法， 該方法是一種基于 變分法或變分里茲法而開展起來的求解微分方程的數(shù)值計算方法， 以電腦為 手段，采用分片近似，進(jìn)而逼近整體的研究思想求解物理問題。有限元分析的根本步驟可歸納為三大步：結(jié)構(gòu)離散、單元分析和整體分析。 對于平面問題，結(jié)構(gòu)離散常用的網(wǎng)格形狀有三角形、矩形、任意四邊形，以三個 頂點(diǎn)為節(jié)點(diǎn)的三角形單元是最簡單的平面單元，它較矩形或四邊形對曲邊邊界有 更好的適應(yīng)性，而矩形或四邊形單元較三節(jié)點(diǎn)三角 形有更高的計算精度。Matlab語言是進(jìn)行矩陣運(yùn)算的強(qiáng)大工具，因 此，用Matlab語言編寫有限兀中

2、平面冋題的程序 有優(yōu)越性。本章將詳細(xì)介紹如何利用 Matlab語言 編制三角形常應(yīng)變單元的計算程序，程序流程圖見 圖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計算結(jié)果整理。計算結(jié)果整理包括位移和應(yīng) 力兩個方面；位移計算結(jié)果一般不需要特別 的處理，利用計算出的節(jié)點(diǎn)位移分量，就可 畫出結(jié)構(gòu)任意方向的位移云圖；而應(yīng)力解的 誤

3、差表現(xiàn)在單元內(nèi)部不滿足平衡方程，單元 與單元邊界處應(yīng)力一般不連續(xù)，在邊界上應(yīng) 力解一般與力的邊界條件不相符合。1.1 程序說明%*三角形常應(yīng)變單元求解結(jié)構(gòu)主程序%* 功能：運(yùn)用有限元法中三角形常應(yīng)變單元解平面問題的計算主程序。 根本思想：單元結(jié)點(diǎn)按右手法那么順序編號。 荷載類型：可計算結(jié)點(diǎn)荷載。說明：主程序的作用是通過賦值語句、讀取和寫入文件、函數(shù)調(diào)用等完成算 法的全過程，即實現(xiàn)程序流程圖的程序表達(dá)。%1 程序準(zhǔn)備format short e%設(shè)定輸出類型clear all%去除所有已定義變量clc%清屏說明：format short e 設(shè)定計算過程中顯示在屏幕上的數(shù)字類型為短格式、科 學(xué)計

4、數(shù)法；clear all 去除所有已定義變量，目的是在本程序的運(yùn)行過程中，不會 發(fā)生變量名相同等可能使計算出錯的情況；clc 清屏，使屏幕在本程序運(yùn)行開始時%2 全局變量定義globalNNODENPIONNELEMNVFIXNFORCE COORDLNODSYOUNGPOISSTHICKglobalFORCEFIXEDglobalBMATXDMATXSMATXAREAglobalASTIFASLODASDISPglobalFP1說明：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

5、單元定義數(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)與仇3) 分別為約束點(diǎn) x 方向和 y 方向的約束情況 ,受約束為 1 否那么為 0BMATX 單元應(yīng)變矩陣 (3*6)， DMATX 單元彈性矩陣 (3*3)，SMATX 單元應(yīng)力矩陣 (3*6)，AREA 單元面積ASTIF 總體剛度矩陣，ASLOD 總體荷載向量，ASDISP結(jié)點(diǎn)位移向量FP1數(shù)據(jù)文件指針3 翻開文件FP1=

6、fopen(input.txt,rt);%翻開輸入數(shù)據(jù)文件存放初始數(shù)據(jù)說明：FP1=fopen(input.txt,rt);-翻開已存在的輸入數(shù)據(jù)文件 input.txt，且設(shè)置其 為只讀格式，使程序在執(zhí)行過程中不能改變輸入文件中的數(shù)值，并用文件句柄FP1 來執(zhí)行FP2=fopen(output.txt,wt); -翻開輸出數(shù)據(jù)文件， 該文件不存在時， 通過此 命令創(chuàng)立新文件，該文件存在時那么將原有內(nèi)容全部刪除。 該文件設(shè)置為可寫格式， 可在程序執(zhí)行過程中向輸出文件寫入數(shù)據(jù)。% 4 讀入程序控制信息%結(jié)點(diǎn)個數(shù)結(jié)點(diǎn)編碼總數(shù)%單元個數(shù)單元編碼總數(shù)結(jié)點(diǎn)荷載個數(shù)彈性模量泊松比%厚度%單元定義數(shù)組單元結(jié)

7、點(diǎn)號NPION=fscanf(FP1,%d,1)NELEM=fscanf(FP1,%d,1)NFORCE=fscanf(FP1,%d,1)NVFIX=fscanf(FP1,%d,1)YOUNG=fscanf(FP1,%e,1)POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%d,1)LNODS=fscanf(FP1,%d,3,NELEM)說明： 建立 LNODS 矩陣，該矩陣指出了每一單元的連接信息。矩陣的每一行針對每一單元，共計NELEM ;每一列相應(yīng)為單元結(jié)點(diǎn)號編 碼、按逆時針順序輸入。命令中，3,NELEM表示讀取NELEM行3列數(shù)據(jù)賦值給LNODS矩陣 顯

8、然，LNODS(i,1:3)依次表示i單元的i， j，k結(jié)點(diǎn)號。COORD=fscanf(FP1,%f,2,NPION)%結(jié)點(diǎn)坐標(biāo)數(shù)組說明：建立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=fsca nf(FP1,%f,3,NFORCE)% 結(jié)點(diǎn)力數(shù)組說明：(n,3) n:受力結(jié)點(diǎn)數(shù)目,(n,1):作用點(diǎn),(n,2):x方向,(n,3):y方向FIXED=fsca nf(FP1,%d,3,i nf)% 約束信息數(shù)組說明：(n,3) n:受約束

9、節(jié)點(diǎn)數(shù)目，(n,1):約束點(diǎn)(n,2)與(n,3)分別為約束點(diǎn)x方向和y 方向的約束情況，受約束為1否那么為0總體說明：從輸入文件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)。采用了命令fscanf,其中：d表示讀入整數(shù)格式，f 表示讀入浮點(diǎn)數(shù)；1 表示讀取1個數(shù)，A,B形式表示讀A行B列數(shù)組，A,B表示將A,B轉(zhuǎn)置，inf 表示正無窮。%-nteJ5調(diào)用子程生成單剛，

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

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

12、2);%跟據(jù)節(jié)點(diǎn)編號對應(yīng)關(guān)系將單元 剛度分塊疊加到總剛矩陣中endendend說明：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)到整體的位移編碼。%

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

14、，假設(shè)不小于那么退出。執(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)假設(shè)等于1,那么約束數(shù)自增，假設(shè)不 等于1那么跳出。FIXED(NUM)表示 FIXED(NUM，1)，tmp(i)=FIXED(NUM)*2+j 計算整體約 束序號，將序號放入tmp行向量中的i列。%for i=1:NVFIXASTIF(1:2*NPION,tmp(i)=0;%將一約束序號處的總剛列向量清0ASTIF(tmp(i),1:2*NPION)=0;

15、%.將一約束序號處的總剛行向量清0ASTIF(tmp(i),tmp(i)=1;%將行列交叉處的元素置為1end說明：0丨，那么將總剛中的主元素后處理法中置0置1法，設(shè)j Cj包括CjKjj換為1, j行和j列的其他元素均改為零。%(*%讀取FORMSMATX子程%(*fun ctio n FORMSMATX(ELEMENT)% 計算應(yīng)力矩陣%引用所需的全局變量global NPION NELEM COORD LNODS YOUNG POISS global BMATX DMATX SMATX AREA%生成彈性矩陣Da=YOUNG/(1-POISSA2);DMATX(1,1)=1*a;DMAT

16、X(1,2)=POISS*a;DMATX(2,1)=POISS*a;DMATX(2,2)=1*a;DMATX(3,3)=(1-POISS)*a/2;說明：1, ,0 0,0平面應(yīng)力問題的彈性矩陣D,1,0，其中，E為彈性模量,為泊松比%i=ELEMENT;%i為當(dāng)前所計算的單元號%計算當(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;說明：1, xi,

17、yidet表示求矩陣行列式的值,1,Xj ,yj其中xi(i, j,m)分別表示一個三1, Xm, ym角形單元的三個節(jié)點(diǎn)坐標(biāo)。MATLAB中假設(shè)一行中無法寫下一個完整的命令，那么可以在行尾參加3個連續(xù)的英文句號，表示命令余下的局部在下一行出現(xiàn)。%生成應(yīng)變矩陣Bfor j=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

18、) 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);說明：bl,0、 1應(yīng)變矩陣 Bl0, c (I i, j, m)2ACl ,blrem表示求余函數(shù)，remx, y命令返回的是x-n.*y，當(dāng)y 0時，n=fix(x./y)， 其中fix為最近的整數(shù)取整。%SMATX=DMATX*BMATX;%求應(yīng)力矩陣 S=D*BQ%*%讀取FORMLOAD子程Q%*fun ctio n FORMLOAD()%本子程生成荷載向量%所需引用的全局變量global ASLOD NPION NFORCE FORCEASLOD

19、(1:2*NPION)=0;%張成特定大小的 0 向量說明：ASLOD 為總體荷載向量，每個結(jié)點(diǎn)有 x， y 兩個方向的結(jié)點(diǎn)力。%for i=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end 說明：F0RCE(i,1)為作用點(diǎn)，F(xiàn)ORCE(i,2:3)為x,y方向的結(jié)點(diǎn)力。%將有約束處的荷載置為 0NUM=1;i=0;tmp(NVFIX)=0;while iNVFIXfor j=-1:0if FIXED(NUM,j+3)=1i=i+1; tmp(i)=FIXED(NUM)*2+j;endend NUM=NUM+1;endf

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

21、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)位移列陣中對號， 賦值給第 i 個單元的單元位移向量ELEDISP。%1.2程序應(yīng)用舉例-;1【例1】利用三角形三節(jié)點(diǎn)常應(yīng)變單元程序計算 圖2所示的結(jié)構(gòu)。設(shè)彈性模量E 1，泊松比為0,厚度為1。%輸入數(shù)據(jù)文件input.txt為:6 4 1 6 1.0e0 0.0 13 1 25 2 43 2 56 3 51 0 -11 1 02 1 04 1 15 0 16 0 1%說明：第一行：讀入程序控制信息%結(jié)點(diǎn)個數(shù)結(jié)點(diǎn)編碼總數(shù)%單元個數(shù)單元編碼總數(shù)%結(jié)點(diǎn)荷載個數(shù)%受約束邊界點(diǎn)數(shù)%彈性模量NPION=fsca nf(FP1,%d,1)NELEM=fsca nf(FP1,%d,1)NFORCE=fsca nf(FP1,%d,1)NVFIX=fsca nf(FP1,%d,1)YOUNG=fsca nf(FP1,%e,1)泊松比%厚度POISS=fscanf(FP1,%f,1)THICK=fscanf(FP1,%d,1) 第二、三、四、五行：讀入單元連接信息：LNODS=fscanf(FP1,%d,3,NELEM);%單元定義數(shù)組，單元結(jié)點(diǎn)號，逆時