結點三角形單元有限元程序MATLAB語言

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

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

3、-1/2N的力作用在matlab環(huán)境下，輸入則程序運行結果如下：該程序求解的結點位移結果和結點應力結果與課本給出的結果一致。附錄：%-平面三角形單元有限元法-%function x strain stress=tri_fem()E,v,t,EN,ecode,NN,node,RN,RC,PN,PC=input; %調入已定材料、幾何尺寸以及單元和結點編號及約束和載荷分布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é)點坐標 D=1,node(ecode(i,1),:) 1,node(ecode(i,2),:) 1,node(ecode(i,3),:); A(i,1)=1/2*det(D);end% 形成單元參數 b=zeros(EN,3

5、); c=zeros(EN,3); %各單元參數初始化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% 求得總剛，并引入約束和載荷求得各結點位移K=Kf(E,v,t,EN,ecode,NN,A,b,c); %調用函數Kf，求得結構的總體剛度矩陣x=xf(NN,RN,RC,PN,PC,K); %調用函數xf,求得各結點位移 % 求解應力應變strai

6、n stress=ss(E,v,EN,ecode,A,b,c,x);% 單元剛度矩陣與結構剛度矩陣function K=Kf(E,v,t,EN,ecode,NN,A,b,c)Ke=zeros(6,6); %單元的剛度矩陣，初始為6*6階零矩陣K=zeros(NN*2,NN*2); %結構的總體剛度矩陣，初始為零矩陣for m=1:1:EN %m為單元號 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); %各單元剛度矩陣整體編號,并疊加 end end K=K+cell2mat(Ka);end %分塊矩陣K合成一個矩陣K% 引入位移向量和右端項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 %被約束結點位移為0%-引入已知結點載荷-%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%-引入已知結點載荷-%-求解非約束結點的位移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法求解非約束結點的位移X%-求解非約束結點的位移X-%-由X可得所有結點位移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消去法 將系數矩陣化成對角矩陣求解方程組的數值解function x=Jordan(A,b)%開始計算

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 %消元計算 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 % 求解應力應變function strain stress=ss(E,v,EN,ecode,A,b,c,x) strain=zeros(3,EN); %應變初始為0矩陣 stress=zeros(3,EN); %應力初始為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