+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);
相关推荐：
loading