1、一、題目如圖1所示，一個(gè)厚度均勻的三角形薄板，在頂點(diǎn)作用沿板厚方向均勻分布的豎向載荷。已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。試計(jì)算薄板的位移及應(yīng)力分布。要求：1. 編寫有限元計(jì)算機(jī)程序，計(jì)算節(jié)點(diǎn)位移及單元應(yīng)力。（劃分三角形單元，單元數(shù)不得少于30個(gè)）；2. 采用有限元軟件分析該問題（有限元軟件網(wǎng)格與程序設(shè)計(jì)網(wǎng)格必須一致），詳細(xì)給出有限元軟件每一步的操作過程，并將結(jié)果與程序計(jì)算結(jié)果進(jìn)行對比（任選取三個(gè)點(diǎn)，對比位移值）；3. 提交程序編寫過程的詳細(xì)報(bào)告及計(jì)算機(jī)程序；4. 所有同學(xué)參加答辯，并演示有限元計(jì)算程序。有限元法中三節(jié)點(diǎn)三角形分析結(jié)構(gòu)的步驟如下：1）整

2、理原始數(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é)果整理。一、程序設(shè)計(jì)網(wǎng)格劃分如圖，將薄板如圖劃分為6行，并建立坐標(biāo)系，則PPXXYY單元編號節(jié)點(diǎn)編號剛度矩陣的集成建立與總剛度矩陣等維數(shù)的空矩陣，已變單元剛度矩陣的集成。由單元分析已知節(jié)點(diǎn)、單元的排布規(guī)律，繼而通過循環(huán)計(jì)算求得每個(gè)單元對應(yīng)的節(jié)點(diǎn)序號。通過循環(huán)逐個(gè)計(jì)算：（1）每個(gè)單元對應(yīng)2種單元剛度矩陣中的哪一種

3、； （2）該單元對應(yīng)總剛度矩陣的那幾行哪幾列 （3）將該單元的單元剛度矩陣加入總剛度矩陣的對應(yīng)行列循環(huán)又分為3層循環(huán)：（1）最外層：逐行計(jì)算 （2）中間層：該行逐個(gè)計(jì)算 （3）最里層：區(qū)分為第 奇/偶 數(shù)個(gè)計(jì)算單元剛度的集成：邊界約束的處理：劃0置1法適用：這種方法適用于邊界節(jié)點(diǎn)位移分量為已知(含為0)的各種約束。  做法：  （1） 將總剛矩陣K中相應(yīng)于已知位移行主對角線元素置1，其他元素改為零；同時(shí)將載荷列陣R中相應(yīng)元素用已知位移置換。  這樣，由該方程求得的此位移值一定等于已知量。 （2） 將K中已

4、知位移相應(yīng)的列的非主對角成元素也置0，以保持K的對稱性。  當(dāng)然，在已知位移分量不為零的情況下，這樣做就改變了方程左端的數(shù)值，為保證方程成立，須在方程右端減去已知位移對該方程的貢獻(xiàn)已知位移和相應(yīng)總剛元素的乘積。  若約束為零位移約束時(shí)，此步則可省去。 特點(diǎn)： （1） 經(jīng)以上處理同樣可以消除剛性位移(約束足夠的前提下)，去掉未知約束反力。 （2） 但這種方法不改變方程階數(shù)，利于存貯。 （3） 不過，若是要求出約束反力，仍要重新計(jì)算各個(gè)劃去的總剛元素。程序如下：變量說明NNODE 單元節(jié)點(diǎn)

5、數(shù) NPION 總結(jié)點(diǎn)數(shù) NELEM 單元數(shù) NVFIX 受約束邊界點(diǎn)數(shù) FIXED 約束信息數(shù)組 NFORCE 節(jié)點(diǎn)力數(shù) FORCE 節(jié)點(diǎn)力數(shù)組 COORD 結(jié)構(gòu)節(jié)點(diǎn)坐標(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é)點(diǎn)位移向量 ELEDISP 單元節(jié)點(diǎn)位移向量 STRESS 單元應(yīng)力 %*%初始化clearformat short e %設(shè)定輸出類型 cle

6、ar %清除內(nèi)存變量 NELEM=36 %單元個(gè)數(shù)（單元編碼總數(shù)） NPION=28 %結(jié)點(diǎn)個(gè)數(shù)（結(jié)點(diǎn)編碼總數(shù)） NVFIX=2 %受約束邊界點(diǎn)數(shù) NFORCE=1 %結(jié)點(diǎn)荷載個(gè)數(shù) YOUNG=2e11 %彈性模量 POISS=0.25 %泊松比 THICK=0.1 %厚度 LNODS=1 2 3;2 4 5;2 5 3;3 5 6; 4 7 8;4 8 5;5 8 9;5 9 6; 6 9 10;7 11 12;7 12 8;8 12 13; 8 13 9;9 13 14;9 14 10;10 14 15; 11 16 17;11 17 12; 12 17 18; 12 18 13; 13

7、18 19; 13 19 14;14 19 20;14 20 15; 15 20 21;16 22 23;16 23 17;17 23 24; 17 24 18;18 24 25;18 25 19;25 19 26; 19 26 20;20 26 27;20 27 21;21 27 28 %單元定義數(shù)組（單元結(jié)點(diǎn)號） %相應(yīng)為單元結(jié)點(diǎn)號（編碼）、按逆時(shí)針順序輸入 COORD=0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3; 1.5 3;-2.25 4.5;-0.75 4.5;0.75 4.5; 2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6; -3.7

8、5 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5; 2.25 7.5;3.75 7.5;-4.5 9;-3 9; -1.5 9;0 9;1.5 9;3 9;4.5 9 %結(jié)點(diǎn)坐標(biāo)數(shù)組 %坐標(biāo)：x，y 坐標(biāo)（共 NPOIN 組） FORCE=1 0 -15 %結(jié)點(diǎn)力數(shù)組（受力結(jié)點(diǎn)編號, x 方向,y 方向） FIXED=22 1 1;28 1 1 %約束信息（約束點(diǎn)，x 約束，y 約束） %有約束為 1，無約束為 0 %*%生成單元剛度矩陣并組成總體剛度矩陣 ASTIF=zeros(2*NPION,2*NPION); %生成特定大小總體剛度矩陣并置 0 %*for i=1:

9、NELEM %生成彈性矩陣 D D= 1 POISS 0; POISS 1 0; 0 0 (1-POISS)/2*YOUNG/(1-POISS2) %*%計(jì)算當(dāng)前單元的面積 A=-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 %*%生成應(yīng)變矩陣 B for j=0:2 b(j+1)=COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNODS(i,

10、(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 B=b(1) 0 b(2) 0 b(3) 0; 0 c(1) 0 c(2) 0 c(3); c(1) b(1) c(2) b(2) c(3) b(3)/(2*A); B1( :,:,i)=B;%* %求應(yīng)力矩陣 S=D*B S=D*B; ESTIF=B'*S*THICK*A; %求解單元剛度矩陣 a=LNODS(i,:); %臨時(shí)向量,用來記錄當(dāng)前單元的節(jié)點(diǎn)編號 for j=1:3 for

11、k=1:3 ASTIF(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é)點(diǎn)編號對應(yīng)關(guān)系將單元剛度分塊疊加到總剛 %度矩陣中 end end end %*%將約束信息加入總體剛度矩陣（對角元素改一法） for i=1:NVFIX if FIXED(i,2)=1 ASTIF(:,(FIXED(i,1)*2-1)=0; %一列為零 ASTIF(FIXED(i,1)*2-1),:)=0; %一行為零 ASTIF(FIXED(i,1)

12、*2-1),(FIXED(i,1)*2-1)=1; %對角元素為 1 end %*%生成單元剛度矩陣并組成總體剛度矩陣 %* if FIXED(i,3)=1 ASTIF( :,FIXED(i,1)*2)=0; %一列為零 ASTIF(FIXED(i,1)*2,:)=0; %一行為零 ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %對角元素為 1 end end %*%生成荷載向量 ASLOD(1:2*NPION)=0; %總體荷載向量置零 for i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end %* %求解內(nèi)力 ASDISP=ASTIFASLOD' %計(jì)算節(jié)點(diǎn)位移向量 ELEDISP(1:6)=0; %當(dāng)前單元節(jié)點(diǎn)位移向量 for i=1:NELEM for j=1:3 ELEDISP(