平面三角形单元Matlab程序

%变量说明%NPOIN NELEM NVFIX NFORCE NNODE%总结点数，单元数，受约束边界点数， 结点力数, 单元结点数%COORD LNODS YOUNG POISS THICK%结构结点坐标数组，单元定义数组，弹性模量， 泊松比， 厚度%FORCE FIXED BMATX DMATX SMATX %节点力数组，约束信息数组，单元应变矩阵，单元弹性矩阵， 单元应力矩阵%AREA HK ASLOD ASDISP FP1%单元面积，总体刚度矩阵，总体荷载向量， 结点位移向量，数据文件指针format short eclearFP1=fopen(C:UsersAdministratorDesktopinput.txt,rt);NPION=fscanf(FP1,%d,1); %结点个数（结点编码总数）NELEM=fscanf(FP1,%d,1); %单元个数（单元编码总数）NFORCE=fscanf(FP1,%d,1); %结点荷载个数NVFIX=fscanf(FP1,%d,1); %受约束边界点数YOUNG=fscanf(FP1,%e,1); %弹性模量POISS=fscanf(FP1,%f,1); %泊松比 THICK=fscanf(FP1,%d,1); %厚度LNODS=fscanf(FP1,%d,3,NELEM); %单元定义数组（单元结点号）COORD=fscanf(FP1,%f,2,NPION); %结点坐标数组FORCE=fscanf(FP1,%f,3,NFORCE); %结点力数组FIXED=fscanf(FP1,%d,3,inf); %约束信息数组%引用所需的全局变量%global NPION NELEM COORD LNODS YOUNG POISS%global BMATX DMATX SMATX AREA%生成弹性矩阵Da=YOUNG/(1-POISS2); DMATX(1,1)=1*a;DMATX(1,2)=POISS*a;DMATX(2,1)=POISS*a;DMATX(2,2)=1*a;DMATX(3,3)=(1-POISS)*a/2; for i=1:NELEM; %i为当前所计算的单元号 %计算当前单元的面积AREA=det(1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);. 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);. 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);)/2; end %生成应变矩阵Bfor j=0:2 b(j+1)=COORD(LNODS(i,(rem(j+1),3)+1),2)-COORD(LNODS(i,(rem(j+2),3)+1),2); c(j+1)=-COORD(LNODS(i,(rem(j+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1);endBMATX=b(1) 0 b(2) 0 b(3) 0;. 0 c(1) 0 c(2) 0 c(3);. c(1) b(1) c(2) b(2) c(3) b(3)/(2*AREA); SMATX=DMATX*BMATX; %求应力矩阵S=D*B% 所引用的全局变量：%global NPION NELEM NVFIX LNODS ASTIF THICK%global BMATX SMATX AREA FIXEDHK=seros(2*NPION,2*NPION); %张成特定大小总刚矩阵并置0 for i=1:NELEM EK=BMATX*SMATX*THICK*AREA; %求解单元刚度矩阵 a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号 for j=1:3 for k=1:3HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+EK(j*2-1:j*2,k*2-1:k*2); %跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中 end end end %将约束信息加入总刚（置0置1法）NUM=1; %计数器(当前已分析的节点数)i=0; %计数器(当前已处理的约束数)tmp(NVFIX)=0; %临时存被约束的序号while iNVFIX for j=-1:0 if FIXED(NUM,j+3)=1 %若发现约束 i=i+1; %计数器自增 tmp(i)=FIXED(NUM)*2+j; %求约束序号 end end NUM=NUM+1;endfor i=1:NVFIX HK(1:2*NPION,tmp(i)=0; %将一约束序号处的总刚列向量清0 HK(tmp(i),1:2*NPION)=0; %.将一约束序号处的总刚行向量清0 HK(tmp(i),tmp(i)=1; %将行列交叉处的元素置为1end %所需引用的全局变量%global ASLOD NPION NFORCE FORCEASLOD(1:2*NPION)=0; %张成特定大小的0向量for i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end %将有约束处的荷载置为0NUM=1; i=0; tmp(NVFIX)=0; while iNVFIX for j=-1:0 if FIXED(NUM,j+3)=1 i=i+1; tmp(i)=FIXED(NUM)*2+j; end end NUM=NUM+1; endfor i=1:NVFIX ASLOD(tmp(i)=0; endASDISP=HKASLOD; %计算结点位移向量 %所引用的全局变量 %global NELEM NPION SMATX ASDISP LNODS ELEDISP(1:6)=