基于Matlab语言的按平面三角形单元划分的结构有限元程序设计

format short e clear

FP1=fopen('LinearTriangleElement of Thin ','rt'); NELEM=fscanf(FP1,'%d',1) NPION=fscanf(FP1,'%d',1) NVFIX=fscanf(FP1,'%d',1) NFORCE=fscanf(FP1,'%d',1) YOUNG=fscanf(FP1,'%e',1) POISS=fscanf(FP1,'%f',1) THICK=fscanf(FP1,'%f',1)

LNODS=fscanf(FP1,'%d',[ NELEM,3]) COORD=fscanf(FP1,'%f',[ NPION,2]) FORCE=fscanf(FP1,'%f',[ NFORCE,3]) FIXED=fscanf(FP1,'%d',[ NVFIX,3])

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 for i=1:NELEM %生成弹性矩阵D

D= YOUNG/(1-POISS^2)*[1 POISS 0; POISS 1 0;

0 0 (1-POISS)/2] %计算当前单元的面积A

A=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 %生成应变矩阵B for 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); end

B=[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*A); B1( :,:,i)=B; %求应力矩阵S=D*B

S=D*B;

ESTIF=B'*S*THICK*A; a=LNODS(i,:); for j=1:3 for k=1:3

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…

=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2); end end end

%将约束信息加入总体刚度矩阵 for i=1:NVFIX if FIXED(i,2)==1

ASTIF(:,(FIXED(i,1)*2-1))=0; ASTIF((FIXED(i,1)*2-1),:)=0;

ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; end

if FIXED(i,3)==1

ASTIF( :,FIXED(i,1)*2)=0; ASTIF(FIXED(i,1)*2,:)=0;

ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; end end

%生成荷载向量

ASLOD(1:2*NPION)=0; for i=1:NFORCE

ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3) end

%求解内力

ASDISP=ASTIF\\ASLOD' ELEDISP(1:6)=0; for i=1:NELEM for j=1:3

ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); end i

STRESS=D*B1(:, :, i)*ELEDISP' end

fclose(FP1);