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

+004 +004 i = 3 STRESS = +004 +004 +004 i = 4 STRESS = +005 +004 +004 （说明：以上四组STRESS输出结果数据表示单元应力SX,SY,SXY，i为单元号。）

六、 ANSYS建模比较分析

利用ANSYS有限元分析软件，完全按照前面Matlab程序设计的各项条件参数以及单元划分方式建立ANSYS模型，其求解结果与以上程序计算结果比较，验证程序是否正确。

ANSYS模型节点单元如下图所示：

ANSYS模型求解变形图如下所示：

ANSYS求解节点位移结果列表显示如下：（单位：m）

PRINT U NODAL SOLUTION PER NODE

***** POST1 NODAL DEGREE OF FREEDOM LISTING *****

THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEM

NODE UX UY UZ USUM 1 2 3 4 5

6 MAXIMUM ABSOLUTE VALUES

NODE 4 4 0 4 VALUE

ANSYS求解单元应力结果列表显示如下：（单位：KN/m2）

PRINT S ELEMENT SOLUTION PER ELEMENT

***** POST1 ELEMENT NODAL STRESS LISTING *****

THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATES ELEMENT= 1 PLANE182

NODE SX SY SZ SXY SYZ 1 +06 -33586. +06 2 +06 -33586. +06 6 +06 -33586. +06 ELEMENT= 2 PLANE182

NODE SX SY SZ SXY SYZ 2 -55503. -55503. -55503. 3 -55503. -55503. -55503. 4 -55503. -55503. -55503. ELEMENT= 3 PLANE182

NODE SX SY SZ SXY SYZ 2 55503. -49526. -94497. 4 55503. -49526. -94497. 5 55503. -49526. -94497. ELEMENT= 4 PLANE182

NODE SX SY SZ SXY SYZ 2 +06 -27040. -17932. 5 +06 -27040. -17932. 6 +06 -27040. -17932.

结论

通过比较Matlab语言设计程序运行结果和ANSYS建模分析结果可知，两种方式算出的结果完全一致，说程序设计正确。所以本程序适用于按三角形单元划分的平面结构有限元分析。

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',[3,NELEM]) COORD=fscanf(FP1,'%f',[2,NPION]) FORCE=fscanf(FP1,'%f',[3,NFORCE]) FIXED=fscanf(FP1,'%d',[3,NVFIX])

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);