有限元編程講義

1、有限元編程示例題目描述題目描述: :如下圖所示的平面桁架，桿件長(zhǎng)度、彈性模量、截面積以如下圖所示的平面桁架，桿件長(zhǎng)度、彈性模量、截面積以及所受節(jié)點(diǎn)力及所受節(jié)點(diǎn)力P的大小可以自行定義。求節(jié)點(diǎn)位移及桿件軸的大小可以自行定義。求節(jié)點(diǎn)位移及桿件軸力。力。例一：桁架例一：桁架解題思路：解題思路： 建立模型建立模型 集成總剛集成總剛 求解位移求解位移 求解桿件軸力求解桿件軸力 輸出結(jié)果輸出結(jié)果建立模型：建立模型：定義節(jié)點(diǎn)坐標(biāo)定義節(jié)點(diǎn)坐標(biāo)Node = zeros(10,2) ;x=-1*L; %L為橫桿長(zhǎng)度為橫桿長(zhǎng)度f(wàn)or i=1:2:10 x=x+L; Node(i,:)=x 0;endx=-1*L;fo

2、r i=2:2:10 x=x+L; Node(i,:)=x H;%H為豎桿長(zhǎng)度為豎桿長(zhǎng)度end模型相關(guān)參數(shù)輸入模型相關(guān)參數(shù)輸入H=input(豎桿長(zhǎng)度(m):); L=input(水平桿長(zhǎng)度(m):);E=input(桿件彈性模量(Gpa):); A=input(桿件截面積(m2):);a=input(節(jié)點(diǎn)力P(kN)：);節(jié)點(diǎn)編號(hào)方式節(jié)點(diǎn)編號(hào)方式定義單元，即儲(chǔ)存單元兩端的節(jié)點(diǎn)號(hào)定義單元，即儲(chǔ)存單元兩端的節(jié)點(diǎn)號(hào)Element=zeros(21,2);for i=1:2:7 Element(5/2*i-3/2,:)=i,i+1; Element(5/2*i-1/2,:)=i,i+2; Eleme

3、nt(5/2*i+1/2,:)=i,i+3; endfor i=2:2:8 Element(5*i/2-1,:)=i,i+1; Element(5*i/2,:)=i,i+2;endElement(21,:)=9,10;加下劃線的為單元編號(hào)加下劃線的為單元編號(hào)集成總剛：集成總剛：xi = Node( Element( ie, 1 ), 1 ) ;%ie為單元號(hào)，以下相同為單元號(hào)，以下相同yi = Node( Element( ie, 1 ), 2 ) ;xj = Node( Element( ie, 2 ), 1 ) ;yj = Node( Element( ie, 2 ), 2 ) ;獲取單元

4、兩端節(jié)點(diǎn)坐標(biāo)獲取單元兩端節(jié)點(diǎn)坐標(biāo)L = ( (xj-xi)2 + (yj-yi)2 )(1/2) ;計(jì)算桿件長(zhǎng)度計(jì)算桿件長(zhǎng)度形成等效荷載列陣形成等效荷載列陣f=0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;%每個(gè)節(jié)點(diǎn)兩個(gè)自由度，每個(gè)節(jié)點(diǎn)兩個(gè)自由度，a為之前輸入的節(jié)點(diǎn)力為之前輸入的節(jié)點(diǎn)力計(jì)算從局部坐標(biāo)到整體坐標(biāo)的坐標(biāo)轉(zhuǎn)換矩陣計(jì)算從局部坐標(biāo)到整體坐標(biāo)的坐標(biāo)轉(zhuǎn)換矩陣Tfunction T = TransformMatrix( ie )%ie為單元號(hào)為單元號(hào)c = (xj-xi)/L ;s = (yj-yi)/L ;T= c -s 0 0 s c 0 0 0 0 c

5、-s 0 0 s c ;計(jì)算單元?jiǎng)偠染仃囉?jì)算單元?jiǎng)偠染仃噆k = E*A/L 0 -E*A/L 0 0 0 0 0 -E*A/L 0 E*A/L 0 0 0 0 0 ; T = TransformMatrix( ie ) ; k = T*k*transpose(T) ;% transpose(T) 為為T的轉(zhuǎn)置矩陣的轉(zhuǎn)置矩陣2集成整體剛度矩陣集成整體剛度矩陣Kfor ie=1:1:21 %按單元順序進(jìn)行循環(huán)按單元順序進(jìn)行循環(huán) k=PlaneTrussElementStiffness(ie); %計(jì)算第計(jì)算第ie個(gè)單元的單剛個(gè)單元的單剛 m=Element(ie,1); %ie單元的首節(jié)點(diǎn)號(hào)單元

6、的首節(jié)點(diǎn)號(hào) n=Element(ie,2); %ie單元的末節(jié)點(diǎn)號(hào)單元的末節(jié)點(diǎn)號(hào) K(2*m-1,2*n-1)=k(1,3); K(2*m-1,2*n)=k(1,4); K(2*m,2*n-1)=k(2,3); K(2*m,2*n)=k(2,4);K=zeros(20,20);%用來(lái)存儲(chǔ)整體剛度矩陣用來(lái)存儲(chǔ)整體剛度矩陣集成總剛的非對(duì)角線元素集成總剛的非對(duì)角線元素(這里的元素指這里的元素指2*2的小矩陣）的小矩陣）在下面的集成中，將在下面的集成中，將總剛看成總剛看成10*10的矩陣，每個(gè)元素為的矩陣，每個(gè)元素為2*2的小矩陣的小矩陣 m=Element(ie,2); %ie單元的末節(jié)點(diǎn)號(hào)單元的末

7、節(jié)點(diǎn)號(hào) n=Element(ie,1); %ie單元的首節(jié)點(diǎn)號(hào)單元的首節(jié)點(diǎn)號(hào) K(2*m-1,2*n-1)=k(3,1); K(2*m-1,2*n)=k(3,2); K(2*m,2*n-1)=k(4,1); K(2*m,2*n)=k(4,2);end 集成總剛的對(duì)角線元素（這里的元素指集成總剛的對(duì)角線元素（這里的元素指2*2的小矩陣）的小矩陣）for i=1:1:10 %按節(jié)點(diǎn)的順序循環(huán)按節(jié)點(diǎn)的順序循環(huán) for j=1:1:21 %對(duì)于每個(gè)節(jié)點(diǎn)，再按單元的順序循環(huán)對(duì)于每個(gè)節(jié)點(diǎn)，再按單元的順序循環(huán) k=PlaneTrussElementStiffness(j); if Element(j,1)=

8、I %如果如果i節(jié)點(diǎn)為節(jié)點(diǎn)為j單元的首節(jié)點(diǎn)單元的首節(jié)點(diǎn) K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1); K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2); K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1); K(2*i,2*i)=K(2*i,2*i)+k(2,2); end if Element(j,2)=i %如果如果i節(jié)點(diǎn)為節(jié)點(diǎn)為j單元的末節(jié)點(diǎn)單元的末節(jié)點(diǎn) K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(3,3); K(2*i-1,2*i)=K(2*i-1,2*i)+k(3,4); K(2*i,2*i-1)=K(2

9、*i,2*i-1)+k(4,3); K(2*i,2*i)=K(2*i,2*i)+k(4,4); end endend求解位移：求解位移：u=zeros(20);根據(jù)約束情況修改總剛，采用對(duì)角元素置根據(jù)約束情況修改總剛，采用對(duì)角元素置1法法for i=1:1:20 K(1,i)=0; K(2,i)=0; K(18,i)=0; K(i,1)=0;K(i,2)=0; K(i,18)=0; end %自由度自由度1、2、18被約束了被約束了,所在的行和列的其他元素都改為所在的行和列的其他元素都改為0K=K*1e15;%乘以一個(gè)大數(shù)，減小計(jì)算誤差乘以一個(gè)大數(shù)，減小計(jì)算誤差f=f*1e15;u=Kf;求解

10、求解K(1,1)=1;%對(duì)角線元素置對(duì)角線元素置1K(2,2)=1;K(18,18)=1;求解軸力：求解軸力：獲取單元兩端的節(jié)點(diǎn)號(hào)獲取單元兩端的節(jié)點(diǎn)號(hào)i = Element( ie, 1 ) ;%ie為單元號(hào)為單元號(hào)j = Element( ie, 2 ) ;獲取單元兩端的節(jié)點(diǎn)位移獲取單元兩端的節(jié)點(diǎn)位移uElement = zeros( 4, 1 ) ;uElement( 1:2 ) = u( (i-1)*2+1:(i-1)*2+2 ) ;uElement( 3:4 ) = u( (j-1)*2+1:(j-1)*2+2 ) ;k = PlaneTrussElementStiffness( ie

11、 ) ;nodef = k *uElement ;%整體坐標(biāo)下的節(jié)點(diǎn)力整體坐標(biāo)下的節(jié)點(diǎn)力T = TransformMatrix( ie ) ;%計(jì)算坐標(biāo)轉(zhuǎn)換矩陣計(jì)算坐標(biāo)轉(zhuǎn)換矩陣nodef = transpose( T ) * nodef ;%從整體坐標(biāo)轉(zhuǎn)換到局部坐標(biāo)從整體坐標(biāo)轉(zhuǎn)換到局部坐標(biāo)計(jì)算單元的節(jié)點(diǎn)力計(jì)算單元的節(jié)點(diǎn)力輸出求解結(jié)果：輸出求解結(jié)果：輸出位移輸出位移fprintf( 節(jié)點(diǎn)位移n ) ;for i=1:1:10 disp(節(jié)點(diǎn)號(hào),num2str(i), x方向位移：,num2str(u(2*i-1,1), y方向位移：,num2str(u(2*i,1);end輸出節(jié)點(diǎn)力輸出節(jié)點(diǎn)力f

12、printf( nn節(jié)點(diǎn)力n ) ; for ie=1:1:21 nodef=NodeForce( ie ); disp( 單元號(hào): ,num2str(ie), 節(jié)點(diǎn)號(hào)：,num2str(Element(ie,1), 節(jié)點(diǎn)號(hào)：,num2str(Element(ie,2), 軸力:,num2str(nodef(1) ) ; end例二：網(wǎng)架例二：網(wǎng)架思路分析思路分析網(wǎng)架是由多根桿件按照一定的網(wǎng)格形式通過(guò)節(jié)點(diǎn)連結(jié)而成的空間結(jié)構(gòu)。構(gòu)成網(wǎng)架的基本單元有三角錐，三棱體，正方體，截頭四角錐等。鑒于網(wǎng)架的形鑒于網(wǎng)架的形式較多，本程序提供一種通用的網(wǎng)架輸入方法，但錄入較為繁瑣，同時(shí)提供式較多，本程序提供一種通

13、用的網(wǎng)架輸入方法，但錄入較為繁瑣，同時(shí)提供一種正放四角錐網(wǎng)架的簡(jiǎn)易輸入方法作為典型。一種正放四角錐網(wǎng)架的簡(jiǎn)易輸入方法作為典型。考慮幾何非線性考慮幾何非線性。本程序采用荷載分級(jí)加載的方式考慮網(wǎng)架的幾何非線性。將總荷載分成1000份分步施加，求解各荷載步下的節(jié)點(diǎn)位移，修改網(wǎng)架相應(yīng)節(jié)點(diǎn)坐標(biāo)以及剛度矩陣，依次迭代求出網(wǎng)架的總位移。本程序的網(wǎng)架位移求解函數(shù)附在主程序后面，主程序運(yùn)行時(shí)調(diào)用該函數(shù)。幾點(diǎn)說(shuō)明幾點(diǎn)說(shuō)明用戶自定義輸入用戶自定義輸入幾何建模幾何建模正放四角錐網(wǎng)正放四角錐網(wǎng)架簡(jiǎn)易輸入架簡(jiǎn)易輸入定義荷載定義荷載定義邊界定義邊界條件條件網(wǎng)架分析網(wǎng)架分析位移位移應(yīng)變應(yīng)變應(yīng)力應(yīng)力位移求解函數(shù)位移求解函數(shù)單剛

14、矩陣單剛矩陣荷載矩陣荷載矩陣約束矩陣約束矩陣總剛矩陣總剛矩陣求解位移求解位移&分級(jí)加載，分級(jí)加載，通過(guò)修改通過(guò)修改節(jié)點(diǎn)坐標(biāo)，節(jié)點(diǎn)坐標(biāo)，迭代求解迭代求解幾何非線性幾何非線性e=input(選擇網(wǎng)架類型，0代表自由定義網(wǎng)架，1代表四角錐網(wǎng)架) %網(wǎng)架類型的選擇網(wǎng)架類型的選擇網(wǎng)架類型的選擇網(wǎng)架類型的選擇用戶自定義網(wǎng)架（網(wǎng)架信息的錄入，包括節(jié)點(diǎn)、單元、截面、彈性模量等）用戶自定義網(wǎng)架（網(wǎng)架信息的錄入，包括節(jié)點(diǎn)、單元、截面、彈性模量等）if e=0%選擇自定義網(wǎng)架選擇自定義網(wǎng)架Node=input(定義節(jié)點(diǎn)編號(hào)及對(duì)應(yīng)坐標(biāo)，按1 x1 y1 z1;2 x2 y2 z2;.輸入);%形成節(jié)點(diǎn)儲(chǔ)存矩陣形成節(jié)點(diǎn)

15、儲(chǔ)存矩陣 Men=input(定義單元與節(jié)點(diǎn)的關(guān)系，按1 node1 node2;2 node3 node4;.輸入,node1node2,依次類推);%形成單元儲(chǔ)存矩陣形成單元儲(chǔ)存矩陣 Msum=length(Men);%查找網(wǎng)架錄入的單元數(shù)查找網(wǎng)架錄入的單元數(shù) Cont1=input(定義單元實(shí)常數(shù)，若所有桿件截面面積和彈性模量不變，則輸入0，否則輸入1);定義單元屬性的輸入方式定義單元屬性的輸入方式 if Cont1=0 AE1=input(請(qǐng)輸入統(tǒng)一的截面面積與彈性模量，按A E輸入); AE=zeros(Msum,3); AE(:,1)=1:Msum;AE(:,2)=AE1(1,1)

16、;AE(:,3)=AE1(1,2); else AE=input(請(qǐng)輸入相應(yīng)單元的截面面積與彈性模量，按1,A1 E1;2,A2 E2;.輸入); end P=input(定義節(jié)點(diǎn)荷載，按node1 P1;node2 P2;.輸入); %網(wǎng)架荷載輸入網(wǎng)架荷載輸入 BC=input(定義邊界約束，按node1 Conx Cony Conz;node2 Conx Cony Conz);.輸入,Con代表x、y、z方向約束，取0為約束，取1無(wú)約束); %網(wǎng)架邊界條件網(wǎng)架邊界條件end單元屬性相同單元屬性相同單元屬性不同單元屬性不同荷載及邊界條件荷載及邊界條件正放四角錐網(wǎng)架定義正放四角錐網(wǎng)架定義if

17、e=1hu=input(輸入網(wǎng)架上層節(jié)點(diǎn)行數(shù)); %定義網(wǎng)架上層節(jié)點(diǎn)的行數(shù)定義網(wǎng)架上層節(jié)點(diǎn)的行數(shù)lu=input(輸入網(wǎng)架上層節(jié)點(diǎn)列數(shù)); %定義網(wǎng)架上層節(jié)點(diǎn)的列數(shù)定義網(wǎng)架上層節(jié)點(diǎn)的列數(shù)dis_xu=input(輸入網(wǎng)架上層節(jié)點(diǎn)列間距); %定義網(wǎng)架上層的行間距定義網(wǎng)架上層的行間距dis_yu=input(輸入網(wǎng)架上層節(jié)點(diǎn)行間距); %定義網(wǎng)架上層的列間距定義網(wǎng)架上層的列間距hd=hu-1; %網(wǎng)架下層節(jié)點(diǎn)的行數(shù)網(wǎng)架下層節(jié)點(diǎn)的行數(shù)ld=lu-1; %網(wǎng)架下層節(jié)點(diǎn)的列數(shù)網(wǎng)架下層節(jié)點(diǎn)的列數(shù)dis_xd=dis_xu; %網(wǎng)架下層的行間距網(wǎng)架下層的行間距dis_yd=dis_yu; %網(wǎng)架下層的行間

18、距網(wǎng)架下層的行間距dis_z=input(輸入網(wǎng)架上下層間距); %網(wǎng)架上下層間距網(wǎng)架上下層間距定義網(wǎng)架上層節(jié)點(diǎn)定義網(wǎng)架上層節(jié)點(diǎn)定義網(wǎng)架下層節(jié)點(diǎn)定義網(wǎng)架下層節(jié)點(diǎn)定義網(wǎng)架高度定義網(wǎng)架高度f(wàn)or i=1:hu for j=1:lu Node(i-1)*lu+j,2)=(j-1)*dis_xu; Node(i-1)*lu+j,3)=(i-1)*dis_yu; Node(i-1)*lu+j,4)=dis_z; endendfor i=1:hd for j=1:ld Node(i-1)*ld+j+hu*lu,2)=(j-1+0.5)*dis_xd; Node(i-1)*ld+j+hu*lu,3)=(i-

19、1+0.5)*dis_yd; Node(i-1)*ld+j+hu*lu,4)=0; endend網(wǎng)架上層節(jié)點(diǎn)編號(hào)與對(duì)應(yīng)坐標(biāo)網(wǎng)架上層節(jié)點(diǎn)編號(hào)與對(duì)應(yīng)坐標(biāo)網(wǎng)架下層節(jié)點(diǎn)編號(hào)與對(duì)應(yīng)坐標(biāo)網(wǎng)架下層節(jié)點(diǎn)編號(hào)與對(duì)應(yīng)坐標(biāo)Nsum=length(Node); %查詢網(wǎng)架的節(jié)點(diǎn)數(shù)查詢網(wǎng)架的節(jié)點(diǎn)數(shù)for i=1:Nsum %將節(jié)點(diǎn)編號(hào)錄入節(jié)點(diǎn)矩陣將節(jié)點(diǎn)編號(hào)錄入節(jié)點(diǎn)矩陣 Node(i,1)=i;endfor i=1:hu for j=1:lu-1 Men(i-1)*(lu-1)+j,2)=(i-1)*lu+j; Men(i-1)*(lu-1)+j,3)=(i-1)*lu+j+1; endendfor i=1:lu for

20、 j=1:hu-1 Men(i-1)*(hu-1)+j+(lu-1)*hu,2)=(j-1)*lu+i; Men(i-1)*(hu-1)+j+(lu-1)*hu,3)=j*lu+i; endend節(jié)點(diǎn)編號(hào)的錄入節(jié)點(diǎn)編號(hào)的錄入網(wǎng)架上層橫向單元的拓?fù)渚W(wǎng)架上層橫向單元的拓?fù)渚W(wǎng)架上層縱向單元的拓?fù)渚W(wǎng)架上層縱向單元的拓?fù)?for i=1:hd for j=1:ld-1 Men(i-1)*(ld-1)+(lu-1)*hu+(hu-1)*lu+j,2)=(i-1)*ld+j+hu*lu; Men(i-1)*(ld-1)+(lu-1)*hu+(hu-1)*lu+j,3)=(i-1)*ld+j+hu*lu+1;

21、 endendfor i=1:ld for j=1:hd-1 Men(i-1)*(hd-1)+(ld-1)*hd+(lu-1)*hu+(hu-1)*lu+j,2)=(j-1)*ld+i+hu*lu; Men(i-1)*(hd-1)+(ld-1)*hd+(lu-1)*hu+(hu-1)*lu+j,3)=j*ld+i+hu*lu; endend網(wǎng)架下層縱向單元的拓?fù)渚W(wǎng)架下層縱向單元的拓?fù)渚W(wǎng)架下層橫向單元的拓?fù)渚W(wǎng)架下層橫向單元的拓?fù)渚W(wǎng)架腹桿單元的拓?fù)渚W(wǎng)架腹桿單元的拓?fù)鋐or i=1:hd for j=1:ld Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)

22、*ld+(ld-1)*hd+1,2)=(i-1)*lu+j; Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+1,3)=(i-1)*ld+hu*lu+j;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+2,2)=(i-1)*lu+j+1;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+2,3)=(i-1)*ld+j+hu*lu;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(

23、lu-1)*hu+(hd-1)*ld+(ld-1)*hd+3,2)=i*lu+j;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+3,3)=(i-1)*ld+j+hu*lu;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+4,2)=i*lu+j+1;Men(i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+4,3)=(i-1)*ld+j+hu*lu; endend腹桿腹桿N腹桿腹桿N+1腹桿腹桿N

24、+2腹桿腹桿N+3單元編號(hào)錄入單元儲(chǔ)存矩陣單元編號(hào)錄入單元儲(chǔ)存矩陣 Msum=length(Men); %查詢網(wǎng)架單元數(shù) for i=1:Msum %將單元編號(hào)錄入單元將單元編號(hào)錄入單元儲(chǔ)存儲(chǔ)存矩陣矩陣 Men(i,1)=i;end定義截面屬性定義截面屬性E=2.1e11;%默認(rèn)材料為steelA1=input(請(qǐng)輸入網(wǎng)架上層單元的截面面積); %默認(rèn)網(wǎng)架上層單元截面尺寸相同默認(rèn)網(wǎng)架上層單元截面尺寸相同A2=input(請(qǐng)輸入網(wǎng)架下層單元的截面面積); %默認(rèn)網(wǎng)架下層單元截面尺寸相同默認(rèn)網(wǎng)架下層單元截面尺寸相同A3=input(請(qǐng)輸入網(wǎng)架腹桿單元的截面面積); %默認(rèn)網(wǎng)架腹桿單元截面尺寸相同默

25、認(rèn)網(wǎng)架腹桿單元截面尺寸相同AE=zeros(Msum,3); %定義單元屬性矩陣定義單元屬性矩陣m1=(hu-1)*lu+(lu-1)*hu; %上層單元截止編號(hào)上層單元截止編號(hào)m2=m1+(hd-1)*ld+(ld-1)*hd; %下層單元截止編號(hào)下層單元截止編號(hào)AE(1:m1,2)=A1; %將上層單元尺寸錄入將上層單元尺寸錄入AE矩陣矩陣AE(m1+1):m2,2)=A2; %將下層單元尺寸錄入將下層單元尺寸錄入AE矩陣矩陣AE(m2+1):Msum,2)=A3; %將腹桿單元尺寸錄入將腹桿單元尺寸錄入AE矩陣矩陣AE(:,1)=1:Msum; %將單元編號(hào)錄入將單元編號(hào)錄入AE矩陣矩陣

26、AE(:,3)=E; %將材料彈性模量錄入將材料彈性模量錄入AE矩陣矩陣定義荷載定義荷載cont2=input(定義節(jié)點(diǎn)荷載，若網(wǎng)架上層節(jié)點(diǎn)力與下層節(jié)點(diǎn)力均布，則輸入0，否則輸入1);if cont2=0 P1=input(請(qǐng)輸入網(wǎng)架上層節(jié)點(diǎn)荷載); P2=input(請(qǐng)輸入網(wǎng)架下層節(jié)點(diǎn)荷載); m3=hu*lu; P(1:Nsum,1)=1:Nsum; P(1:m3,2)=P1;P(m3+1):Nsum,2)=P2;else P=input(定義節(jié)點(diǎn)荷載，按node1 P1;node2 P2;.輸入); end定義邊界條件定義邊界條件cont3=input(定義邊界約束，若網(wǎng)架上層周邊節(jié)點(diǎn)全

27、約束，則輸入0，若下層周邊節(jié)點(diǎn)全約束，輸入1，否則輸入2);if cont3=0 n1=2*(hu+lu-2); BC=zeros(n1,4); BC(1:lu-2,1)=2:lu-1; BC(lu-1):(2*lu-4),1)=lu*(hu-1)+2:lu*hu-1; BC(2*lu-3):(2*lu-4+hu),1)=1:lu:lu*(hu-1)+1; BC(2*lu-3+hu):n1,1)=lu:lu:hu*lu;elseif cont3=1 n1=2*(hd+ld-2); BC=zeros(n1,4); BC(1:ld-2,1)=2:ld-1; BC(ld-1):(2*ld-4),1)

28、=ld*(hd-1)+2:ld*hd-1; BC(2*ld-3):(2*ld-4+hd),1)=1:ld:ld*(hd-1)+1; BC(2*ld-3+hd):n1,1)=ld:ld:hd*ld; for i=1:n1 BC(i,1)=BC(i,1)+hu*lu; endelse BC=input(定義邊界約束，按node1 Conx Cony Conz;node2 Conx Cony Conz);.輸入,Con代表x、y、z方向約束，取0為約束，取1無(wú)約束);endendNsum=length(Node);Msum=length(Men);Psum=length(P);BCsum=lengt

29、h(BC); %提取各矩陣的行數(shù)提取各矩陣的行數(shù)考慮幾何非線性分析網(wǎng)架考慮幾何非線性分析網(wǎng)架for i=1:Psum %將力分為將力分為1000份份 P(i,2)=P(i,2)/1000;endU=zeros(3*Nsum,1); %總位移矩陣總位移矩陣for i=1:1000 u,L1,Kz = grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum); for j=1:Nsum %根據(jù)節(jié)點(diǎn)位移修改網(wǎng)架的節(jié)點(diǎn)坐標(biāo)根據(jù)節(jié)點(diǎn)位移修改網(wǎng)架的節(jié)點(diǎn)坐標(biāo) Node(j,2)=Node(j,2)+u(3*j-2,1); Node(j,3)=Node(j,3)+u(3*j-1,

30、1); Node(j,4)=Node(j,4)+u(3*j,1); end U=U+u; %每次迭代位移的疊加每次迭代位移的疊加end迭代法修正剛度矩陣和網(wǎng)架位移迭代法修正剛度矩陣和網(wǎng)架位移求解網(wǎng)架桿件的應(yīng)力應(yīng)變求解網(wǎng)架桿件的應(yīng)力應(yīng)變L0=zeros(Msum,1); %所有根桿的最初長(zhǎng)度所有根桿的最初長(zhǎng)度f(wàn)or i=1:Msum %單元兩端的節(jié)點(diǎn)編號(hào)單元兩端的節(jié)點(diǎn)編號(hào) p=Men(i,2); q=Men(i,3); X1=Node(p,2); %單元端節(jié)點(diǎn)的坐標(biāo)單元端節(jié)點(diǎn)的坐標(biāo) Y1=Node(p,3); Z1=Node(p,4); X2=Node(q,2); Y2=Node(q,3); Z

31、2=Node(q,4); L0(i,1)=sqrt(X2-X1)2+(Y2-Y1)2+(Z2-Z1)2); %網(wǎng)架桿件的初始長(zhǎng)度網(wǎng)架桿件的初始長(zhǎng)度endLt=L1-L0; %所有桿件長(zhǎng)度的增量所有桿件長(zhǎng)度的增量strain=zeros(Msum,1); %定義應(yīng)變矩陣定義應(yīng)變矩陣stress=zeros(Msum,1); %定義應(yīng)力矩陣定義應(yīng)力矩陣for i=1:Msum E=AE(i,3); strain(i,1)=Lt(i,1)/L0(i,1); %第第i根桿件應(yīng)變根桿件應(yīng)變 stress(i,1)=E*strain(i,1); %第第i根桿件應(yīng)力根桿件應(yīng)力enddisp(U); %輸出網(wǎng)

32、架節(jié)點(diǎn)位移輸出網(wǎng)架節(jié)點(diǎn)位移disp(stress); %輸出網(wǎng)架桿件應(yīng)力輸出網(wǎng)架桿件應(yīng)力網(wǎng)架節(jié)點(diǎn)位移求解函數(shù)網(wǎng)架節(jié)點(diǎn)位移求解函數(shù)function u,L1,Kz = grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum)Kz=zeros(3*Nsum,3*Nsum); %定義剛度矩陣定義剛度矩陣L=zeros(Msum,1);for i=1:Msum %單元兩端的節(jié)點(diǎn)編號(hào)單元兩端的節(jié)點(diǎn)編號(hào) p=Men(i,2); q=Men(i,3); A=AE(i,2); E=AE(i,3);% 單元兩端節(jié)點(diǎn)坐標(biāo)單元兩端節(jié)點(diǎn)坐標(biāo) X1=Node(p,2);Y1=Node(p,3);Z1=Node(p,4);X2=Node(q,2);Y2=Node(