结点三角形单元有限元程序MATLAB语言

1、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,ecode,NN,node,RN,RC

2、,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负方向上-1/2N的力作用在matlab环境下，

3、输入则程序运行结果如下：该程序求解的结点位移结果和结点应力结果与课本给出的结果一致。附录：%-平面三角形单元有限元法-%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 error(

4、9;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),:); %各三角形单元的节点坐标 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 r>k 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

13、 A(i,j)=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)