有限元计算结构力学fortran程序

1、圏丄箝支辺兴爭计算结构力学程序 计算结构力学编程大作业 时间： 2007年6月 1* !关于程序的说明 1* !一、功能： !1、可计算包括节点力，一般非节点力，支座沉降、温度荷载作用、制造误差的平 !面桁架、梁、刚架及其组合结构的节点位移与杆端力； !2、可同时计算多种工况下的节点位移与杆端力。 1* 1* 1 !二、变量说明： 1NE单元数； 1N结构中自由度数； !NJ节点数； 1NS一一特殊节点数，包扌舌支座节点、主从节点(1节点不做主节点)、连接桁架的较节点(没有转角)； ! NAI一一结构的单元截面类型数： ! MT一一单元截面类型号； 1NL一一荷载工况数； ! H截面高度； !

2、 E弹性模量； ! JC单元定位向量数组； !X(NJ),Y(NJ)节点的X, Y坐标值； ! JE(NE,2)一一单元两端节点码数组； !AI(NAI,2)按单元类型顺序存放A与I, Al (I, 1)第I类单元的截面积，Al (I, 2)第I类单元的 !惯性矩； ! MT(NE)单元所属单元类型号； !JS(NS,4)特殊节点信息，JS(I,1)结点码；JS(I,2), JS(I,3), JS(I,4)U, V, CETA 约束信息， !有约束为1,没有约束为0；从节点某位移同主节点时位移时，该位移约束信息填主节点码； 1 !PJ(NP,3)节点荷载信息数组；PJ (I, 1)节点力所在节

3、点号；PJ (L 2)节点力作用坐标方向： !坐标方向U, V, M分别为1, 2, 3； PJ (L 3)节点力的大小(含正负号)：U, V方向集中力时, !与坐标轴正向同向为正，M按右手法则为正；本程序推导过程取y轴向下为正。 1 ! PF(NF,4)一一非节点荷载数组，并给出以卞类型说明： !前6类型数据输法(梯形等可以用叠加法计算)： ! PF(I,1)-单元码；PF(I,2)-类型；PF(I,3)-荷载大小；PF(I,4)-c 值； !1一一垂直于单元的均布力，人小为q,以坐标轴正向为正,c为荷载末端距1节点距离； ! 2一一非节点集中力P, c为荷载距1节点距离； -! 3一非节点

4、集中力距M, c为荷载距1节点距离，右手法则判正负; ! 4一一三角形荷载，c为荷载距1节点距离，1端为0,距离1端c时力为q; !j端为0的三角形，可按叠加法处理。 ! 5一一沿杆轴向均布力,大小为q,c为荷载末端距1节点距离； ! 6一沿杆轴向集中力,大小为q,c为荷载末端距1节点距离； 1 !从第7到第9类型(支座沉降)数据输法:PF(I,1)-单元码;PF(I,2)-类型；PF(I,3)-位移大小(含正负)，坐 标轴正向位为正，转角按右手法则：PF(I,4)-沉降所在的单元位移分量，1端为1-3, j端为4-6； 1 ! 7一一沿轴向支座沉降； ! 8一垂直于轴向支座沉降； ! 9一一

5、支座转动； !10一一制造误差,PF(I,1)一制造误差所在单元，PF(I,2)-类型；PF(I,3)-误差人小(含正负)，正负取决于消除 !误差时端点的运动方向，PF(I,4)一误差所在坐标号； !11一一温度荷载，PF(I,1)荷载所在单元，数据形式为:ElementNo.l，如2单元上有温度荷载，则PF(I,1)=2.1; !PF(I,2)一温度变化值tl, PF(I,3)温度变化值t2, PF(I,4)一材料线膨胀系数； t !TK(NN)一采用一维存储结构刚度矩阵，上半带元素(每列第一个非零元素到对角元)； ! KD一一主元地址数组，表示结构刚度矩阵的主元在TK中的序号,KD中最后一

6、个数是TK中元素的总个数: ! H结构刚度矩阵上半带的非对角元素在TK中的地址，JI=KD(J)-J+I： ! JN(NJ,3)一一结点位移分量编号数组，用于存放结点三个位移的位移分量号码， !JN(I,1)，JN(I,2), JN(I,2)-分别为结点I的U, V, CETA分量的位移分量(坐标)号码； 1 ! P(N)一节点荷载列阵；在回代求位移时存放位移量； ! F(N)一求得的杆端力列阵； ! FO(6)一一等效节点荷载列阵； m* *平面结构分析源程序内容* m* PROGRAM PFF DIMENSION X(50),Y(50),JE(50,2)MT(30),AI(l 0、2) J

7、S(20,4)fJ(50,3)fF(50、4)，JN(50,3)， & KD(150),TK( 1000/P(l 50),F(6人H(50) DOUBLE PRECISION TK.RF CHARACTER *200 TL OPEN(hFILE=,INDAr.DATSTATUS=,OLD,) OPEN(2.FILE=,OUTDAT.DArSTArUS=TSIEW,) READ(l,70) TL READ(IJO) TL READ(1，*)NE,NJ，NS，NAI,NL、E 3 计算结构力学程序 10 30 40 50 60 20 70 WRITE(2 J 0)NE,NJ,NS,NAI,NL、E

8、 FORNIAT(5X/PLANE FRAME STRUCTORE ANALYSIS75X； * * * *7/2X/CONTROL PARAMETERS &OF STRUCTURE75X,?.Q*C*C*(0.5-2.0*C1/3.0+().25*C2) FO(5)=-Q*C*C2*(1.0-0.5*C1) FO(6)=Q*C*C*C1*(1.0/3.0-0.25*C1) RETURN 20FO(2)=.Q*B*B*(1.0+2.0*C 1)/BL/BL FO(3)=-Q*C *B *B/BL/BL FO(5)=-Q*C2*(1.0+2.0*B/BL) FO(6)=Q*C2*B RETURN

9、 30FO(2)=6.0*Q*C1 *B/BLBL FO(3)=Q*B*(2.03.0*B/BL)ZBL FO(5)=-6.0*Q*C1 *B/BL/BL FO(6 Q*C1*(2.O.3.O*C1) RETURN 40FO(2)=-Q*C*(0.5-0.75*C2+0.4*C3) FO(3)=.Q*C*C*(L0/3.0-0.5*Cl+0.2*C2) FO(5)=-Q*C*C2*(0.75.0.4*C1) FO(6 =Q*C*C*C 1 *(0.25-0.2*C 1) RETURN 50FO(1)=-Q*C*(1.0-0.5*C1) FO(4)=-0.5*Q*C*C1 RETURN 60FO

10、(1)=-Q*B/BL FO(4)=.Q*C1 70 L=INT(C) K=MT(M) S=E*AI(K1)*Q,BL FO(L)=S IF(L.EQ.1)FO(4)=-S IF(L.EQ.4)FO(1)=-S RETURN 80 L=INT(C) K=MT(M) FO(L)=12.0*E*AI(K,2)*Q/BL/BL/BL IF(L.EQ.2)FO(5)=-FO(2) IF(L.EQ.5)FO(2)=-FO(5) FO(3)=0.5*BL*FO(2) FO(6)=FO(3) RETURN 90 L=INT(C) K=MT(M) S=2.0*E*AI(K,2)*Q/BL FO(L)=20*S

11、 IF(LEQ3)FO(6)=S IF(LEQ6)FO(3)=S FO(2)=3.0*S/BL FO(5)=-FO(2) RETURN 100 L=INT(C) K=MT(M) S=E*AI(K,1)*Q,BL FO(L)=S IF(L.EQ.1)FO(4)=-S IF(L.EQ.4)FO(1)=-S RETURN 110 K=MT(M) FO(1 尸05*E*AI(KJ)*PF(I,4)*(PF(I,2)+PF(I,3) F0(2)=00 FO(3)=-E*AI(K7)*PF(L4)*(PF(I,2)+PF(I,3)/H(M) FO(4)=-FO(1) F0(5)=00 FO(6)=-FO(

12、3) END * 算例： 如图所示平面刚架，尺寸见图。其中梁的截面积人= 0.2MJ厶= 0.01AT；梁的截面积A =0.25A/2, /=0.06M4, = 3x107A/M2o支座沉降当做单独的工况计算。 输入数据文件(INDAT.DAT)： COMPOSITED PLANE FRAME STRUCTURE ANALYSIS NE NJ NS NAI NL E 673223.0E7 * COORDINATES OF JOINT JOINT x(l) Y(I) 1 0.0 0.0 2 0.0 0.0 3 0.0 5.0 4 8.0 0.0 存 b 8.0 5.0 6 0.0 10.0 7

13、8.0 10.0 * INFORMATION OF ELEMENTS ELEMENT JOINT-I JOINT-J AI-TYPE H 1 2 4 1 0 2 1 3 2 0 3 4 2 0 4 3 存 1 0 存 3 6 2 0 657 * INFORMATION OF A & OF ELEMENT TYPE-NO.AI 10.20.01 20.250.006 * INFORMATION OF SPECIAL JOINT JOINT U V CETA 1 2 2 0 6 1 1 1 7 1 1 1 * FIRST LOAD CASE INFORMATION NPNF 21 * JOINT

14、FORCE NO. JOINT COORDINATE LOAD 1 1.0 1.0 18.0 23.01.0 24.0 * NON-JOINT FORCE NO. ELEMENT TYPE LOAD 1 1.0 1.0 20.0 8.0 * SECOND LOAD CASE INFORMATION NPNF 02 * JOINT FORCE NO. JOINT COORDINATE LOAD * NON-JOINT FORCE NO. ELEMENT TYPE LOAD 16.07.0 0.03 4.0 2 6.0 8.0 -0.02 5.0 输出计算结果(OUTDAT.DAT)： PLANE

15、 FRAME STRUCTURE ANALYSIS * CONTROL PARAMETERS OF STRUCTURE NE=6 NJ= 7NS= 3 NAI=2 NL=2 E= .3000E 十 08 COORDINATES OF JOINTS JOINT X Y 1 .0000 .0000 2 .0000 .0000 3 .0000 5.0000 4 &0000 .0000 5 &0000 5.0000 6 .0000 10.0000 7 &0000 10.0000 INFORMATION OF ELEMENTS ELEMENT JOINT-I JOINT-J TYPE 1 2 4 1 2

16、 1 3 2 3 4 存 5 2 4 3 5 1 5 3 6 2 6 5 7 2 TYPE A I 1 .200000 .010000 2 .250000 .006000 INFORMATION OF SPECIAL JOINTS JOINT uvceta 1220 6 111 7111 LOAD DATA LOAD CASE= 1 NP= 2NF= 1 JOINT LOAD JOINTXYMLOAD 1.1.18.0000 3.1.24.0000 NON-JOINT LOAD ELEMENT TYPE LOAD 1. 1. 20.0000 &0000 JOINT DISPLACEMENTS

17、JOINTU Vceta 1 198620D-02 2 .198620D-02 3 175088D-02 .846572D-04 .846572D-04 .383719D-04 -.485791D-05 .105183D-02 150910D-03 4.195921D-02 .128676D-03-.664935D-03 5 175352D-02 6 .000000D 十 00 7 000000D+00 .682948D-04 .000000D-K)0 000000D+00 278622D03 .000000D+00 .000000D+00 MEMBER-END FORCES OF ELEME

18、NTS ELEMENT N V M 1 Nl= .2024D-02 Vl= -.6943D 十 02 Ml= -.4336D-08 N2= -.2024D-02 V2= -.9057D+02 M2= .8458D 十 02 2 Nl= .6943D-02 Vl= .2243D+01 Ml= .8882D-14 N2= -.6943D-02 V2= -.2243D+01 M2= .1122D+02 3 Nl= .9057D+02 Vl= -.2024D 十 02 Ml= -.8458D+02 N2= -.9057D+02 V2= .2024D+02 M2= -.1664D+02 4 Nl= J9

19、79D十01 Vl= .1187D 十 02 Ml= .4269D+02 N2= .1979D十01 V2= -.1187D 十 02 M2= .5227D 十 02 Nl= .5756D-02 Vl= -.2374D 十 02 Ml= -.5391D+02 N2= -.5756D+02 V2= .2374D+02 M2= -.6477D 十 02 6 Nl= .1024D-03 Vl= -.1826D 十 02 Ml= -.3563D+02 N2= -.1024D-03 V2= .1826D+02 M2= -.5569D 十 02 LOAD DATA LOAD CASE= 2 NP= 0NF

20、= 2 NON-JOINT LOAD ELEMENTTYPELOADC 6.7.03004.0000 MEMBER-END FORCES OF ELEMENTS ELEMENT N V M 1 Nl= .2843D 十 02 Vl= .4692D+00 Ml= -.2611D-08 N2= -.2843D 十 02 V2= -.4692D+00 M2= .3754D 十 01 2 Nl= -.4692D 十 00 Vl= .2843D+02 Ml= .1705D-12 N2= .4692D 十 00 V2= -.2843D 十 02 M2= .1422D 十 03 3 Nl= .4692D-00 Vl= -.2844D 十 0