有限元程序设计004.doc

有限元程序 设 计 报 告 程序功能说明 本程序包含四个子程序，其主要功能分别是： 1.MODPS（DMATX，YOUNG，POISS） 该子程序是形成弹性矩阵[D]的，[D]阵的平面应力问题表达式为： (4-56) 如果要计算平面应变问题，可将上式中弹性模量 、泊松比 值分别用 、 值代替。 2、BMATPS（INELE，COORD，LNODE，BMATX，AREA，NPOIN，NELEM） 该子程序用来计算应变矩阵[B]，[B]的表达式为： (4-57) 3、DBE（DMATX，BMATX，DMATX） 该子程序用来形成应力矩阵[S]，[S]阵的表达式为 (4-58) 4、GUASS（NZERO，LZERO，ASTIF，NHBW，ALOAD，NP2，NPOIT） 该子程序是用来进行支座处理和解方程的，最后输出节点处位移分量。支座处理的基本原理是如果位移受到约束，则将其对应刚度对角元素变为1，该对角元素所对应的行和列的其余元素变为0，并且将右端项约束位移对应元素充为0。解方程应用带消去法，解出位移后，将节点位移分量存放在原来存放右端项的标识符ALOAD中，最后将ALOAD（节点位移）输出。 程序框图设计 程序标识符及数组的说明 数组及变量说明 改正COORD 节点坐标数组 LNODE 单元节点数组 ALOAD 荷载（右端项）数组，解出位移后，将节点处的位移充于其中 ESTIF 单元刚度矩阵，为6×6阶方阵 ASTIF 总刚度矩阵，带状存贮 DBMAT 应力矩阵[S]阵，为3×6阶 DMATX 弹性矩阵[D]阵，为3×3阶 BMATX 应变矩阵[S]阵，为3×6阶 STRES 单元应力矩阵，为3×1阶 JJS 受载节点矩阵 ZX Ｘ方向已知载荷向量 ZY Y方向已知载荷向量 DISP 单元节点位移矩阵 LZERO 约束位移矩阵 变量的说明 NPOIT 最大节点数 NELEM 最大单元数 NNODE 单元节点数（三角元中为3） NLOAD 最大受载节点数 NZERO 最大约束位移个数 YOUNG 弹性模量 POISS 泊松比 TE 板厚（假设t=1） NHBW 半带宽 NP2 位移总数 AREA 单元面积 4. 源程序 FOR FINITE ELEMENT DIMENSION COORD(6,2),LNODE(4,3),ALOAD(12),ESTIF(6,6), 1ASTIF(12,8),DBMAT(3,6),STRES(3),BMATX(3,6),DMATX(3,3), 1jjS(1),ZX(1),ZY(1),DISP(6),LZERO(6) c 输入已知数据 OPEN(1,FILE=D:\NMX\DATA.DAT) DATA NPOIN,NELEM,NNODE,NLOAD,NZERO/6,4,3,1,6/ DATA YOUNG,POISS,TE/210000.0,0.3,1.0/ DATA ZX/0.0/ DATA ZY/-1.0/ DATA JJS/1/ DATA LNODE/1,4,2,6,2,5,5,3,3,2,3,5/ DATA COORD/0.0,0.0,0.5,0.0,0.5,1.0,1.0,0.5,0.5,0.0,0.0,0.0/ DATA LZERO/1,3,7,8,10,12/ C 计算半带宽 NHBW=0 DO 11 INELE=1,NELEM DO 11 I=1,NNODE DO 11 J=1,NNODE LN=IABS(LNODE(INELE,I)-LNODE(INELE,J)) IF (LN.GT.NHBW) NHBW=LN 11 CONTINUE NHBW=(NHBW+1)*2 WRITE(1,*)半带宽 WRITE(1,12) NHBW 12 FORMAT(1x,NHBW=,I3) NP2