有限元程序设计.doc

1、.平面问题有限元程序设计理学院 学号 兔兔设计人 完成日期 一、 程序功能说明本程序适用于节点荷载作用下的桁架分析问题，当有节间荷载存在时可按照静力等效原理将其转化为节点荷载。可求解平面桁架在静力荷载作用下的内力和位移。二、 框图的设计开 始 输 入 数 据 数 组 定 义计算各杆截面面积和半带宽 调用形成单刚矩阵UNIT调用形成半带宽存贮的结构原始刚度矩阵TOTAL有节点荷载否输入节点荷载值，并将其送入相应的荷载列阵P（N）中考虑结构是否自重将杆自重引起的等效荷载叠加到P（N）中支座处理、解方程，并输出U（N）、V（N）调用UNIT，求各单元杆端内力 结 束单元循环没有有否是 三、 程序的标

2、识符及数组说明NPOIN 最大节点数NELEM 最大单元数 NLOAD 节点的荷载总数NZERO 节点的约束位移总数WT 结构的自重EE 材料的弹性模量LL 一维数组，用于存放单元杆件的长度AA 一维数组，用于存放单元杆件的面积COORD 节点坐标数组LNODE 单元节点数组BH 二维数组，用于存放单元截面尺寸NRES 二维数组，用于存放约束的位移值JP 二维数组，用于存放节点的荷载值ESTIF 四维数组，用于存放整体坐标系下的单元刚度矩阵ASTIF 二十维数组，用于存放半带宽结构原始刚度矩阵P 用于存放节点的荷载列阵U 用于存放节点x方向的位移值V 用于存放节点y方向的位移值四、 源程序IN

3、TEGER E,NELEM,Z,H REAL LL,ESTIF,ASTIF,JP DIMENSION COORD(3,2),LNODE(3,2),AA(200),BH(3,2),RES(3,2), &LL(200),ESTIF(4,4),ASTIF(400,20),JP(1,2),P(400),U(200), &V(200) OPEN(2,FILE='D:NMXJIA.DAT',STATUS='NEW') C 输入已知数据DATA NPOIN,NELEM,NJP,NRES,EE,WT/3,3,1,3,21000,0/DATA COORD/0,6,

4、0,0,0,6/DATA LNODE/1,2,1,2,3,3/DATA BH/3*2,3*10/DATA RES/3*0,1,2,4/C 计算各单元面积DO 200 E=1,NELEMAA(E)=BH(E,1)*BH(E,2) CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY)200 CONTINUE C计算半带宽 L2=2*NPOIN NHBW=0 DO 210 E=1,NELEM M=ABS(LNODE(E,1)-LNODE(E,2) IF(NHBW.LT.M) NHBW=M210 CONTINUEWRITE(2,*) '半带宽' N

5、HBW=2*(NHBW+1) WRITE(2,220) NHBW220 FORMAT(1X,'NHBW=', I2)C单元循环 DO 300 I1=1,L2 DO 300 J1=1,NHBW300 ASTIF(I1,J1)=0.0 DO 400 E=1,NELEM CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) CALL TOTAL(E,LNODE,ESTIF,ASTIF) 400 CONTINUE DO 560 N=1,L2560 P(N)=0.0 IF(NJP.EQ.0) GOTO 650DATA JP/10,5/ DO 630

6、K1=1,NJP NN=JP(K1,2)+0.1 630 P(NN)=JP(K1,1) 650 IF(WT.LE.0.0) GOTO 750 DO 700 E=1,NELEM N1=LNODE(E,1) N2=LNODE(E,2) P(2*N1)=P(2*N1)-WT*AA(E)*LL(E)/2.0 P(2*N2)=P(2*N2)-WT*AA(E)*LL(E)/2.0700 CONTINUE WRITE(2,710)710 FORMAT(/4X,'荷载总数',8X,'水平荷载',8X,'铅垂荷载') DO 730 K=1,NO730 WRITE(

7、2,740) K,P(2*K-1),P(2*K)740 FORMAT(4X,I2,8X,F8.3,8X,F8.3) 750 DO 800 I1=1,NRES Z=RES(I1,2)+1E-5 ASTIF(Z,1)=ASTIF(Z,1)*1E8 P(Z)=ASTIF(Z,1)*RES(I1,1)800 CONTINUE DO 850 K1=1,L2-1 IF(L2.GT.(K1+NHBW-1) THEN IM=K1+NHBW-1 ELSE IM=L2 ENDIF DO 850 I1=K1+1,IM L1=I1-K1+1 C1=ASTIF(K1,L1)/ASTIF(K1,1) DO 830 J1=

8、1,NHBW-L1+1 MM=J1+I1-K1 ASTIF(I1,J1)=ASTIF(I1,J1)-C1*ASTIF(K1,MM)830 CONTINUE P(I1)=P(I1)-C1*P(K1)850 CONTINUE P(L2)=P(L2)/ASTIF(L2,1) DO 900 I1=L2-1,1,-1 IF(NHBW.GT.(L2-I1+1) THEN JM=L2-I1+1 ELSE JM=NHBW ENDIF DO 880 J1=2,JM H=J1+I1-1 P(I1)=P(I1)-ASTIF(I1,J1)*P(H)880 CONTINUE P(I1)=P(I1)/ASTIF(I1,1

9、)900 CONTINUE WRITE(2,910)910 FORMAT(/10X,'节点位移',10X,'水平位移',10X,'铅垂位移'/) DO 930 N=1,NO U(N)=P(2*N-1) V(N)=P(2*N)930 WRITE(2,950) N,U(N),V(N)950 FORMAT(15X,I2,6X,F12.7,6X,F12.7) WRITE(2,970)970 FORMAT(/4X,'单元号',8X,'节点号',8X,'N(KN)',8X,'Q(KN)'/) D

10、O 980 E=1,NELEM CALL UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) N1=LNODE(E,1) N2=LNODE(E,2) ULNODE=U(N1)-U(N2) VLNODE=V(N1)-V(N2) D1=ESTIF(1,1)*ULNODE+ESTIF(1,2)*VLNODE D2=ESTIF(1,2)*ULNODE+ESTIF(2,2)*VLNODE FI=CX*D1+CY*D2 FJ=-FI TI=-CY*D1+CX*D2 TJ=-TI WRITE(2,990) E,N1,FI,TI,N2,FJ,TJ990 FORMAT(4X,I2,

11、12X,I2,8X,F8.4,8X,F8.4/18X,I2,8X,F8.4,8X,F8.4)980 CONTINUE WRITE(2,1000)1000FORMAT(/28X,'结束',/15X,35('*')/) STOP END 子程序TOTAL形成总刚度矩阵SUBROUTINE TOTAL(E,LNODE,ESTIF,ASTIF) INTEGER E,DH,ZL,DL REAL ESTIF,ASTIF DIMENSION LNODE(3,2),ESTIF(4,4),ASTIF(400,20) DO 40 I1=1,2 DO 40 II=1,2 KH=2*(

12、I1-1)+II DH=2*(LNODE(E,I1)-1)+II DO 40 J1=1,2 DO 40 JJ=1,2 KL=2*(J1-1)+JJ ZL=2*(LNODE(E,J1)-1)+JJ DL=ZL-DH+1 IF(DL.GT.0) ASTIF(DH,DL)=ASTIF(DH,DL)+ESTIF(KH,KL)40 CONTINUE RETURN END 子程序UNIT形成单刚SUBROUTINE UNIT(E,EE,COORD,LNODE,AA,ESTIF,LL,CX,CY) INTEGER E REAL LL,ESTIF DIMENSION COORD(3,2),LNODE(3,2)

13、,AA(200),LL(200),ESTIF(4,4) N1=LNODE(E,1) N2=LNODE(E,2) CX=COORD(N2,1)-COORD(N1,1) CY=COORD(N2,2)-COORD(N1,2) LL(E)=SQRT(CX*CX+CY*CY) CX=CX/LL(E) CY=CY/LL(E) EAL=EE*AA(E)/LL(E) ESTIF(1,1)=EAL*CX*CX ESTIF(1,2)=EAL*CX*CY ESTIF(2,2)=EAL*CY*CY ESTIF(2,1)=ESTIF(1,2) DO 10 I=1,2 DO 10 J=1,2 ESTIF(I,J+2)=-

14、ESTIF(I,J) ESTIF(I+2,J)=-ESTIF(I,J)10 ESTIF(I+2,J+2)=ESTIF(I,J) RETURN END五、 算例如图1所示桁架，已知杆件材料的弹性模量,杆件截面高度H=10cm，截面的宽度为b=2cm，不计各杆的自重，求在荷载作用下，各杆的轴力。图1桁架例图计算模型的输入数据有以下各量：NPOINNRESNELEMNLOADYOGWT3331210000.0节点单元12LNODE(1,I)12LNODE(2,I)23LNODE(3,I)13编号数组123RES(I,1)0.00.00.0RES(I,2)124节点 I123COORD(I,1)0.06.00.0COORD(I，2)0.00.06.0JP(I,1)10.0JP(I,2)5单元号I123BH(I,1)101010BH(I,