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

2、理原始數(shù)據(jù)，如材料性質、荷載條件、約束條件等，離散結構并進行單元編碼、結點編碼、結點位移編碼、選取坐標系。2）單元分析，建立單元剛度矩陣。3）整體分析，建立總剛矩陣。4）建立整體結構的等效節(jié)點荷載和總荷載矩陣5）邊界條件處理。6）解方程，求出節(jié)點位移。7）求出各單元的單元應力。8）計算結果整理。一、程序設計網(wǎng)格劃分如圖，將薄板如圖劃分為6行，并建立坐標系，則PPXXYY單元編號節(jié)點編號剛度矩陣的集成建立與總剛度矩陣等維數(shù)的空矩陣，已變單元剛度矩陣的集成。由單元分析已知節(jié)點、單元的排布規(guī)律，繼而通過循環(huán)計算求得每個單元對應的節(jié)點序號。通過循環(huán)逐個計算：（1）每個單元對應2種單元剛度矩陣中的哪一種

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

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

5、數(shù) NPION 總結點數(shù) NELEM 單元數(shù) NVFIX 受約束邊界點數(shù) FIXED 約束信息數(shù)組 NFORCE 節(jié)點力數(shù) FORCE 節(jié)點力數(shù)組 COORD 結構節(jié)點坐標數(shù)組 LNODS 單元定義數(shù)組 YOUNG 彈性模量 POISS 泊松比 THICK 厚度 B 單元應變矩陣(3*6) D 單元彈性矩陣(3*3) S 單元應力矩陣(3*6) A 單元面積 ESTIF 單元剛度矩陣 ASTIF 總體剛度矩陣 ASLOD 總體荷載向量 ASDISP 節(jié)點位移向量 ELEDISP 單元節(jié)點位移向量 STRESS 單元應力 %*%初始化clearformat short e %設定輸出類型 cle

6、ar %清除內存變量 NELEM=36 %單元個數(shù)（單元編碼總數(shù)） NPION=28 %結點個數(shù)（結點編碼總數(shù)） NVFIX=2 %受約束邊界點數(shù) NFORCE=1 %結點荷載個數(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ù)組（單元結點號） %相應為單元結點號（編碼）、按逆時針順序輸入 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 %結點坐標數(shù)組 %坐標：x，y 坐標（共 NPOIN 組） FORCE=1 0 -15 %結點力數(shù)組（受力結點編號, x 方向,y 方向） FIXED=22 1 1;28 1 1 %約束信息（約束點，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) %*%計算當前單元的面積 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 %*%生成應變矩陣 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;%* %求應力矩陣 S=D*B S=D*B; ESTIF=B'*S*THICK*A; %求解單元剛度矩陣 a=LNODS(i,:); %臨時向量,用來記錄當前單元的節(jié)點編號 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é)點編號對應關系將單元剛度分塊疊加到總剛 %度矩陣中 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 %* %求解內力 ASDISP=ASTIFASLOD' %計算節(jié)點位移向量 ELEDISP(1:6)=0; %當前單元節(jié)點位移向量 for i=1:NELEM for j=1:3 ELEDISP(