《有限元计算程序》.doc

有 限 元 计 算 程 序 班级： 结构工程 学号： 094811109 姓名： 王 兰 指导老师： 周 文 伟 【算例一】如图所示的悬臂梁，受均布荷载q=1N/mm2作用。 QUOTE 按图划分为80个单元，55个节点，分布载荷已事先分配相应的节点，用有限元法计算此结构位移及应力。 解： 此为平面应力问题，每节点有x、y两方向位移，NF=2；三角单元有3个节点，NE=3；总计55个节点，NN=55.由图可见，任一单元中两节点编号之最大差值为6，因而最大半带宽MB=2（6+1）=14；根部有5个固定节点，每点有两项位移，则ND=10；共有80个单元，LE=80；平面内任一点有两个坐标，LS=2；弹性模量取E=210Gpa；泊松比UM=0.3；板厚T=10，对平面问题，T可为任意值；MS=110，ME=6，用以上数据可确定程序中数组体积。 由图，按节点顺序，依次给出其x、y坐标值，确定数组XY。分布荷载q已分配导1、6、11、16、21、26、31、36、41、46、51等节点，沿y轴向下如图，其余各节点荷载皆为零（固定点的约束力虽为未知，但在处理固定位移约束时，将由程序处理，此处可填任意值，现也给为零值），由此决定数组Q。 逐个单元按逆时针方向给出其3个节点号码（单元内第一节点可任意），确定数组NX。 最后5个节点被固定，相应的位移自由度序号为101—110，用以确定数组MC。输入数据文件INP.DAT中依次给出上述全部数据。 1、源程序 PROGRAM MAIN DIMENSION SK(300,30),EK(12,12),Q(300),MC(55),XY(3,100),XYE(3,4),QE(12),NX(4,100) OPEN(7,FILE=INP.DAT) REWIND 7 READ(7,*)NF,NE,NN,MB,ND,LE,LS READ(7,*)E,UM,T 10 FORMAT (7I5) 12 FORMAT(3F15.2) WRITE(*,600)NF,NE,NN,MB,ND,LE,LS,E,UM,T ME=NE*NF MS=NN*NF CALL INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND) WRITE(*,102)((XY(I,J),I=1,LS),J=1,NN) 102 FORMAT(10X,XY/,(2X,6F12.3)) WRITE(*,101)(Q(I),I=1,MS) 101 FORMAT(10X,Q/,(2X,6F12.3)) WRITE(*,500)((NX(I,J),I=1,NE),J=1,LE) WRITE(*,400)(MC(I),I=1,ND) 500 FORMAT(10X,NX/,(2X,12I6)) 600 FORMAT(10X,NF NE NN MB ND LE LS E UM T/7(2X,I4),3(2X,F8.4)) 400 FORMAT(10X,MC/,(2X,10I6)) CALL STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,UM,T) CALL SOLVE (SK,Q,MS,MB) OPEN (9,FILE=OUT.DAT) REWIND 9 WRITE(9,200) WRITE(9,250)(Q(I),I=1,MS) 200 FORMAT(5X,DISPLACEMENT) 250 FORMAT(2X,6E14.5) CALL STRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T) STOP 1000 END SUBROUTINE INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND) DIMENSION XY(LS,NN),Q(MS),NX(NE,LE),MC(ND) READ(7,*)XY READ(7,*)Q READ(7,*)NX READ(7,*)MC CLOSE(7) 10 FORMAT(6F11.2) 20 FORMAT(12I5) RETURN END SUBROUTINE STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,UM,T) DIMENSION SK(MS,MB),EK(ME,ME),Q(MS),NX(NE,LE),MC(ND),XY(LS,NN),XYE(LS,NE) DO 35 I=1,MS DO 35 J=1,MB 35 SK(I,J)=0. DO 200 L=1,LE DO 40 J=1,N