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

2. Matlab语言程序代码

%******************************************************************************

%初始化及数据调用

format short e %设定输出类型 clear %清除内存变量

FP1=fopen('','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',[3,NELEM])' %单元定义数组（单元结点号）

%相应为单元结点号（编码）、按逆时针顺序输入

COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组

%坐标：x，y坐标（共NPOIN组）

FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组（受力结点编号, x方向,y方向） FIXED=fscanf(FP1,'%d',[3,NVFIX])' %约束信息（约束点，x约束，y约束） %有约束为1，无约束为0

%*****************************************************************************

%生成单元刚度矩阵并组成总体刚度矩阵

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0

%***************************************************************************** for i=1:NELEM

%生成弹性矩阵D

D= [1 POISS 0; POISS 1 0;

0 0 (1-POISS)/2]*YOUNG/(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; %对角元素为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; %对角元素为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); %关闭数据文件

五、 程序运行结果

NELEM = 4 NPION = 6 NVFIX = 2 NFORCE = 1 YOUNG = 0 POISS =

THICK =

LNODS =

1 2 6 2 3 4 2 4 5 2 5 6 COORD =

0 0 1 0 2 0 2 1 1 1 0 1 FORCE =

4 0 -150 FIXED =

1 1 1 6 1 1 D =

+008 +007 0 +007 +008 0

0 0 +007 A =

D =

+008 +007 0 +007 +008 0

0 0 +007 A = D =

+008 +007 0 +007 +008 0

0 0 +007 A = D =

+008 +007 0 +007 +008 0

0 0 +007 A =

ASDISP =

0 0

0 0 （说明：以上12个ASDISP输出结果数据表示节点位移向量，即各节点分别在X，Y方向上产生的位移。） i = 1 STRESS = +005 +004 +005 i = 2 STRESS = +004