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

3结点三角形单元有限元程序（MATLAB语言）学号：2011100290吴晴晴该程序包括以下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,PN,PC]=input()E=2.1e11;v=1/3;t=1;%杨氏模量Pa，泊松比，厚度EN=2;%单元数ecode=[312;%单元编号单元13-1-2;单元21-3-4134];NN=4;%结点数node=[00;%各结点坐标20;21;01];RN=2;%被约束的位移数RC=[14];%被约束的结点PN=2;%有载荷的结点数%PC（1）表示有载荷的结点，PC(2)表示各结点的力，PC（3）表示载荷方向，0为x方向，1为y方向PC=[23;-1/2-1/2;11];%结点2、3分别有y负方向上-1/2N的力作用在matlab环境下，输入则程序运行结果如下：该程序求解的结点位移结果和结点应力结果与课本给出的结果一致。附录：%%-------------平面三角形单元有限元法---------------------%%function[xstrainstress]=tri_fem()[E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input;%调入已定材料、几何尺寸以及单元和结点编号及约束和载荷分布[nm]=size(ecode);ifEN~=nerror('WrongelementnumberENorwrongelementcodeecode!');return;end[nm]=size(node);ifNN~=nerror('WrongnodenumberNNorwrongnode-coordinatenode!');return;ende=zeros(EN,6);A=zeros(EN,1);%面积fori=1:ENe(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);c=zeros(EN,3);%各单元参数初始化fori=1:ENb(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,求得各结点位移%%求解应力应变[strainstress]=ss(E,v,EN,ecode,A,b,c,x);%%单元刚度矩阵与结构刚度矩阵functionK=Kf(E,v,t,EN,ecode,NN,A,b,c)Ke=zeros(6,6);%单元的刚度矩阵，初始为6*6阶零矩阵K=zeros(NN*2,NN*2);%结构的总体刚度矩阵，初始为零矩阵form=1:1:EN%m为单元号fori=1:1:3forj=1:1:3Ke(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(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;endendKe=E*t/(4*(1-v^2)*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块fori=1:1:3forj=1:1:3Ka(ecode(m,i),ecode(m,j))=Kb(i,j);%各单元刚度矩阵整体编号,并叠加endendK=K+cell2mat(Ka);end%分块矩阵K合成一个矩阵K%%引入位移向量和右端项functionx=xf(NN,RN,RC,PN,PC,K)x=ones(NN*2,1);%位移初始为0向量fori=1:RNx(RC(i)*2-1)=0;x(RC(i)*2)=0;end%被约束结点位移为0%%----------------引入已知结点载荷-------------%px=zeros(NN*2,1);%载荷初始为0向量fori=1:PNifPC(3,i)==1px(PC(1,i)*2)=PC(2,i);elseifPC(3,i)==0px(PC(1,i)*2-1)=PC(2,i);endendend%%----------------引入已知结点载荷-------------%%%----------------求解非约束结点的位移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=zeros(2*RN,2*(NN-RN));BNa=mat2cell(BN,ones(1,RN)*2,ones(1,NN-RN)*2);m=1;fori=1:1:NNifx(2*i)==1M(m)=i;m=m+1;endendfori=1:1:NN-RNforj=1:1:NN-RNANa(i,j)=Ka(M(i),M(j));bna(i,1)=pxa(M(i),1);endendfori=1:RNforj=1:NN-RNBNa(i,j)=Ka(RC(i),M(j));endendAN=cell2mat(ANa);bn=cell2mat(bna);BN=cell2mat(BNa);X=Jordan(AN,bn);%利用Gauss-Jordan法求解非约束结点的位移X%%----------------求解非约束结点的位移X-------------%%----------------由X可得所有结点位移x-------------%BN=BN*X;m=1;n=1;fori=1:1:NNifx(2*i)==1x(2*i-1)=X(m);x(2*i)=X(m+1);m=m+2;elseifx(2*i)==0px(2*i-1)=BN(n);px(2*i)=BN(n+1);n=n+2;endendend%%列主元Jordan消去法将系数矩阵化成对角矩阵求解方程组的数值解functionx=Jordan(A,b)%开始计算，赋初值[n,m]=size(A);x=zeros(n,1);fork=1:n%选主元max1=0;fori=k:nifabs(A(i,k))max1max1=abs(A(i,k));r=i;endend%交换两行ifrkforj=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%消元计算b(k)=b(k)/A(k,k);forj=k+1:nA(k,j)=A(k,j)/A(k,k);endfori=1:nifi~=kforj=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend%输出xfori=1:nx(i)=b(i);end%%求解应力应变function[strainstress]=ss(E,v,EN,ecode,A,b,c,x)strain=zeros(3,EN);%应变初始为0矩阵stress=zeros(3,EN);%应力初始为0矩阵D=E/(1-v^2)*[1v0;v10;00(1-v)/2];form=1:1:ENB=[b(m,1)0b(m,2)0b(m,3)0;0c(m,1)0c(m,2)0c(m,3);c(m,1)b(m,1)c(m,2)b(m,2)c(m,3)b(m,3)];B=B/2/A(m,1);%应变矩阵S=D*B;%应力矩阵X=[x(2*ecode(m,1)-1),x(2*ecode(m,1)),x(2*ecode(m,2)-1),x(2*ecode(m,2)),x(2*ecode(m,3)-1),x(2*ecode(m,3))]';strain(:,m)=B*X;stress(:,m)=S*X;end