結(jié)點(diǎn)三角形單元有限元程序MATLAB語(yǔ)言

1、精選優(yōu)質(zhì)文檔-傾情為你奉上3結(jié)點(diǎn)三角形單元有限元程序（MATLAB語(yǔ)言）學(xué)號(hào)： 吳晴晴該程序包括以下6個(gè)部分：1 主程序tri_fem：用于數(shù)據(jù)的錄入和其他程序的調(diào)用；2 總剛程序Kf：計(jì)算結(jié)構(gòu)的總體剛度；3 各結(jié)點(diǎn)位移求解程序xf：求解各結(jié)點(diǎn)的位移；4 線性方程組求解程序Jordan：Gauss-Jordan法求解非約束結(jié)點(diǎn)的位移；5 應(yīng)力應(yīng)變程序ss：由各結(jié)點(diǎn)位移求解各單元內(nèi)的三個(gè)結(jié)點(diǎn)的應(yīng)力stress和應(yīng)變strain；6 數(shù)據(jù)錄入程序input：錄入材料、幾何尺寸、單元編號(hào)和結(jié)點(diǎn)編號(hào)、位移約束和已知載荷等。以課本P25頁(yè)例2.2為例，其input程序?yàn)閒unction E,v,t,EN

2、,ecode,NN,node,RN,RC,PN,PC=input()E=2.1e11; v=1/3; t=1; %楊氏模量Pa，泊松比，厚度EN=2; %單元數(shù)ecode=3 1 2; %單元編號(hào) 單元1 3-1-2;單元2 1-3-4 1 3 4; NN=4; %結(jié)點(diǎn)數(shù)node=0 0; %各結(jié)點(diǎn)坐標(biāo) 2 0; 2 1; 0 1; RN=2; %被約束的位移數(shù)RC=1 4; %被約束的結(jié)點(diǎn)PN=2; %有載荷的結(jié)點(diǎn)數(shù)%PC（1）表示有載荷的結(jié)點(diǎn)，PC(2)表示各結(jié)點(diǎn)的力，PC（3）表示載荷方向，0為x方向，1為y方向PC=2 3; -1/2 -1/2; 1 1; %結(jié)點(diǎn)2、3分別有y負(fù)方向上

3、-1/2N的力作用在matlab環(huán)境下，輸入則程序運(yùn)行結(jié)果如下：該程序求解的結(jié)點(diǎn)位移結(jié)果和結(jié)點(diǎn)應(yīng)力結(jié)果與課本給出的結(jié)果一致。附錄：%-平面三角形單元有限元法-%function x strain stress=tri_fem()E,v,t,EN,ecode,NN,node,RN,RC,PN,PC=input; %調(diào)入已定材料、幾何尺寸以及單元和結(jié)點(diǎn)編號(hào)及約束和載荷分布n m=size(ecode);if EN=n error(Wrong elementnumber EN or wrong elementcode ecode!); return;endn m=size(node);if NN=n

4、 error(Wrong nodenumber NN or wrong node-coordinate node!); return;ende=zeros(EN,6); A=zeros(EN,1); %面積for i=1:EN e(i,:)=node(ecode(i,1),:),node(ecode(i,2),:),node(ecode(i,3),:); %各三角形單元的節(jié)點(diǎn)坐標(biāo) D=1,node(ecode(i,1),:) 1,node(ecode(i,2),:) 1,node(ecode(i,3),:); A(i,1)=1/2*det(D);end% 形成單元參數(shù) b=zeros(EN,3

5、); c=zeros(EN,3); %各單元參數(shù)初始化for i=1:EN b(i,1)=e(i,4)-e(i,6); b(i,2)=e(i,6)-e(i,2); b(i,3)=e(i,2)-e(i,4); c(i,1)=-e(i,3)+e(i,5); c(i,2)=-e(i,5)+e(i,1); c(i,3)=-e(i,1)+e(i,3);end% 求得總剛，并引入約束和載荷求得各結(jié)點(diǎn)位移K=Kf(E,v,t,EN,ecode,NN,A,b,c); %調(diào)用函數(shù)Kf，求得結(jié)構(gòu)的總體剛度矩陣x=xf(NN,RN,RC,PN,PC,K); %調(diào)用函數(shù)xf,求得各結(jié)點(diǎn)位移 % 求解應(yīng)力應(yīng)變strai

6、n stress=ss(E,v,EN,ecode,A,b,c,x);% 單元?jiǎng)偠染仃嚺c結(jié)構(gòu)剛度矩陣function K=Kf(E,v,t,EN,ecode,NN,A,b,c)Ke=zeros(6,6); %單元的剛度矩陣，初始為6*6階零矩陣K=zeros(NN*2,NN*2); %結(jié)構(gòu)的總體剛度矩陣，初始為零矩陣for m=1:1:EN %m為單元號(hào) for i=1:1:3 for j=1:1:3 Ke(2*i-1,2*j-1)=b(m,i)*b(m,j)+(1-v)*c(m,i)*c(m,j)/2; Ke(2*i-1,2*j)=v*c(m,i)*b(m,j)+(1-v)*b(m,i)*c(

7、m,j)/2; Ke(2*i,2*j-1)=v*b(m,i)*c(m,j)+(1-v)*c(m,i)*b(m,j)/2; Ke(2*i,2*j)=c(m,i)*c(m,j)+(1-v)*b(m,i)*b(m,j)/2; end end Ke=E*t/(4*(1-v2)*A(EN).*Ke; %獲得單元m的剛度矩陣 Kb=mat2cell(Ke,ones(1,3)*2,ones(1,3)*2); %將單元矩陣Ke分為3*3塊 set1=ones(1,NN)*2; Ka=mat2cell(zeros(NN*2,NN*2),set1,set1); %將總剛K分為NN*NN塊 for i=1:1:3

8、for j=1:1:3 Ka(ecode(m,i),ecode(m,j)=Kb(i,j); %各單元?jiǎng)偠染仃囌w編號(hào),并疊加 end end K=K+cell2mat(Ka);end %分塊矩陣K合成一個(gè)矩陣K% 引入位移向量和右端項(xiàng)function x=xf(NN,RN,RC,PN,PC,K)x=ones(NN*2,1); %位移初始為0向量for i=1:RN x(RC(i)*2-1)=0; x(RC(i)*2)=0;end %被約束結(jié)點(diǎn)位移為0%-引入已知結(jié)點(diǎn)載荷-%px=zeros(NN*2,1); %載荷初始為0向量for i=1:PN if PC(3,i)=1 px(PC(1,i)

9、*2)=PC(2,i); else if PC(3,i)=0 px(PC(1,i)*2-1)=PC(2,i); end endend%-引入已知結(jié)點(diǎn)載荷-%-求解非約束結(jié)點(diǎn)的位移X-%set1=ones(1,NN)*2;Ka=mat2cell(K,set1,set1); pxa=mat2cell(px,set1,1);AN=zeros(2*(NN-RN),2*(NN-RN);ANa=mat2cell(AN,ones(1,NN-RN)*2,ones(1,NN-RN)*2);bn=zeros(2*(NN-RN),1);bna=mat2cell(bn,ones(1,NN-RN)*2,1);BN=ze

10、ros(2*RN,2*(NN-RN);BNa=mat2cell(BN,ones(1,RN)*2,ones(1,NN-RN)*2);m=1; for i=1:1:NN if x(2*i)=1 M(m)=i; m=m+1; endend for i=1:1:NN-RN for j=1:1:NN-RN ANa(i,j)=Ka(M(i),M(j); bna(i,1)=pxa(M(i),1); end end for i=1:RN for j=1:NN-RN BNa(i,j)=Ka(RC(i),M(j); end end AN=cell2mat(ANa); bn=cell2mat(bna); BN=ce

11、ll2mat(BNa); X=Jordan(AN,bn); %利用Gauss-Jordan法求解非約束結(jié)點(diǎn)的位移X%-求解非約束結(jié)點(diǎn)的位移X-%-由X可得所有結(jié)點(diǎn)位移x-% BN=BN*X; m=1; n=1; for i=1:1:NN if x(2*i)=1 x(2*i-1)=X(m); x(2*i)=X(m+1); m=m+2; else if x(2*i)=0 px(2*i-1)=BN(n); px(2*i)=BN(n+1); n=n+2; end end end % 列主元Jordan消去法 將系數(shù)矩陣化成對(duì)角矩陣求解方程組的數(shù)值解function x=Jordan(A,b)%開(kāi)始計(jì)算

12、，賦初值n,m=size(A);x=zeros(n,1);for k=1:n %選主元 max1=0; for i=k:n if abs(A(i,k)max1 max1=abs(A(i,k); r=i; end end %交換兩行 if rk for j=k:n z=A(k,j); A(k,j)=A(r,j); A(r,j)=z; end z=b(k); b(k)=b(r); b(r)=z; end %消元計(jì)算 b(k)=b(k)/A(k,k); for j=k+1:n A(k,j)=A(k,j)/A(k,k); end for i=1:n if i=k for j=k+1:n A(i,j)=

13、A(i,j)-A(i,k)*A(k,j); end b(i)=b(i)-A(i,k)*b(k); end endend%輸出xfor i=1:n x(i)=b(i);end % 求解應(yīng)力應(yīng)變function strain stress=ss(E,v,EN,ecode,A,b,c,x) strain=zeros(3,EN); %應(yīng)變初始為0矩陣 stress=zeros(3,EN); %應(yīng)力初始為0矩陣 D=E/(1-v2)*1 v 0;v 1 0;0 0 (1-v)/2; for m=1:1:EN B=b(m,1) 0 b(m,2) 0 b(m,3) 0; 0 c(m,1) 0 c(m,2) 0 c(m,3); c(m,1) b(m,1) c(m