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

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 0.0000 0.0000 0.0000 0.0000 2 -0.80607E-03-0.15848E-02 0.0000 0.17780E-02 3 -0.10281E-02-0.44727E-02 0.0000 0.45893E-02 4 0.11937E-02-0.46947E-02 0.0000 0.48441E-02 5 0.86670E-03-0.18880E-02 0.0000 0.20774E-02 6 0.0000 0.0000 0.0000 0.0000 MAXIMUM ABSOLUTE VALUES

NODE 4 4 0 4 VALUE 0.11937E-02-0.46947E-02 0.0000 0.48441E-02

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 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 2 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182

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

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

NODE SX SY SZ SXY SYZ 2 0.16793E+06 -27040. 0.0000 -17932. 0.0000

5 0.16793E+06 -27040. 0.0000 -17932. 0.0000 6 0.16793E+06 -27040. 0.0000 -17932. 0.0000

结论

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

format short e clear

FP1=fopen('LinearTriangleElement of Thin plate.txt','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