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

《有限元计算结构力学fortran程序》由会员分享，可在线阅读，更多相关《有限元计算结构力学fortran程序（21页珍藏版）》请在装配图网上搜索。

1、计算结构力学编程大作业时间：2007年6月$ 力、mt，+ + + + + + + + + + + + +力力力壬*关于程序的说明esse $ 力、mt，+ + + + + + + + + + + + +力力力壬*!一、功能：!1、可计算包括节点力，一般非节点力，支座沉降、温度荷载作用、制造误差的平!面桁架、梁、刚架及其组合结构的节点位移与杆端力；!2、可同时计算多种工况下的节点位移与杆端力。，力力力力力力力力力力 $ & 力力力力力 $ 力、mmmmx “$ & 力力“!*力$力力$力、+ + + + + + + + + + + + +!*!二、变量说明：! NE单元数；! N结构中自由度数

2、；! NJ节点数；! NS特殊节点数，包括支座节点、主从节点(1节点不做主节点)、连接桁架的铰节点(没有转角)；! NAI结构的单元截面类型数；! MT单元截面类型号；! NL荷载工况数；! H截面高度；! E弹性模量；! JC单元定位向量数组；! X(NJ),Y(NJ)节点的X，Y坐标值；! JE(NE,2)单元两端节点码数组；! AI(NAI,2)按单元类型顺序存放A与I，AI(I，1)第I类单元的截面积，AI(I，2)第I类单元的!惯性矩；! MT(NE)单元所属单元类型号；! JS(NS,4)特殊节点信息，JS(I,1)一结点码；JS(I,2)，JS(I,3)，JS(I,4)U，V，C

3、ETA约束信息，!有约束为1，没有约束为0；从节点某位移同主节点时位移时，该位移约束信息填主节点码；! PJ(NP,3)节点荷载信息数组；PJ(I，1)节点力所在节点号；PJ(I，2)节点力作用坐标方向：!坐标方向U，V，M分别为1，2, 3； PJ(I，3)一节点力的大小(含正负号)；U，V方向集中力时，!与坐标轴正向同向为正，M按右手法则为正；本程序推导过程取y轴向下为正。! PF(NF,4)非节点荷载数组，并给出以下类型说明：!前6类型数据输法(梯形等可以用叠加法计算)：! PF(I,1)-单元码；PF(I,2)-类型;PF(I,3)-荷载大小；PF(I,4)-c 值；! 1垂直于单元的

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

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

6、的主元在TK中的序号，KD中最后一个数是TK中元素的总个数;! JI结构刚度矩阵上半带的非对角元素在TK中的地址，JI=KD(J)-J+I；! JN(NJ,3)结点位移分量编号数组，用于存放结点三个位移的位移分量号码，!JN(I,1)，JN(I,2)，JN(I,2)-分别为结点I的U，V，CETA分量的位移分量(坐标)号码；! P(N)节点荷载列阵；在回代求位移时存放位移量；! F(N)求得的杆端力列阵；! FO(6)等效节点荷载列阵；个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个

7、个个个个个个个个e力力力力力力力力力力力力力力力力“*平面结构分析源程序内容*PROGRAM PFFDIMENSION X(50),Y(50),JE(50,2),MT(30),AI(10,2),JS(20,4),PJ(50,3),PF(50,4),JN(50,3),& KD(150),TK(1000),P(150),F(6),H(50)DOUBLE PRECISION TK,P,FCHARACTER *200 TLOPEN(1,FILE=INDAT.DAT,STATUS=OLD)OPEN(2,FILE=OUTDAT.DAT,STATUS=NEW)READ(1,70) TLREAD(1,70)

8、TLREAD(1,*)NE,NJ,NS,NAI,NL,E10304050602070WRITE(2,10)NE,NJ,NS,NAI,NL,EFORMAT(5X,PLANE FRAME STRUCTURE ANALYSIS/5X,*/2X,CONTROL PARAMETERS &OF STRUCTURE/5X,NE=,I2,8X,NJ=,I2,8X,NS=,I2,8X,NAI=,I2,/5X,NL=,I2,8X,E=,E12.4)CALL INPUT(NE,NJ,NS,NAI,X,YJE,MT,AI,JS,H) !读入数据文件CALL DJN(NJ,NS,JS,JN,N)!计算结构自由度数N,形

9、成结点位移分量数组JNCALL ADE(NE,NJ,N,JE,JN,KD,NN)!形成主元地址数组 KD(N)CALL SSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK)!形成总刚，一维存储数组 TK(NN)CALL UTDU3(TK,NN,KD,N)!对总刚进行UTDU分解，以用于解方程组DO 20 LC=1,NL!对工况循环READ(1,70)TLREAD(1,70)TLREAD(1,70)TLREAD(1,*)NP,NF!读入工况信息WRITE(2,30)LC,NP,NFFORMAT(/2X,LOAD DATA/10X,LOAD CASE=,I3/10

10、X,NP=,I3,8X,NF=,I3)CALL NLV(NE,NJ,NAI,E,N,NP,NF,X,YJE,JN,PJ,PF,MT,AI,P,H)!形成总荷载列阵 P(N)CALL BACK3(TK,NN,P,N,KD,JN,NJ) !解方程组并输出结点位移，存放在数组P(N)中WRITE(2,40)FORMAT(/4X,MEMBER-END FORCES OF ELEMENTS/4X,ELEMENT,13X,N,17X,V,17X,M)DO 60 M=1,NECALL MQN(M,NE,NJ,NAI,N,NF,E,X,YJE,MT,AI,JN,PF,P,F,H) !计算单元杆端力，存放在数组

11、 F(6)中WRITE(2,50)M,(F(I),I=1,6)!输出杆端力FORMAT(/1X,I10,3X,N1=,D12.4,3X,V1=,D12.4,3X,M1=,D12.4/14X,N2=,&D12.4,3X,V2=,D12.4,3X,M2=,D12.4)CONTINUECONTINUEFORMAT(A)CLOSE(1)CLOSE(2)ENDSUBROUTINE INPUT(NE,NJ,NS,NAI,X,YJE,MT,AI,JS,H)!读入数据文件DIMENSION X(NJ),Y(NJ),JE(NE,2),MT(NE),AI(NAI,2),JS(NS,4),H(NE)INTEGER

12、NOREAD(1,70) TLREAD(1,70) TLREAD(1,70) TLREAD(1,*)(NO,X(I),Y(I),I=1,NJ)READ(1,70) TLREAD(1,70) TLREAD(1,70) TLREAD(1,*)(NO,JE(I,1),JE(I,2),MT(I),H(I),I=1,NE)READ(1,70) TLREAD(1,70) TL1020304070!102040506030READ(1,70) TLREAD(1,*)(NO,(AI(I,J),J=1,2),I=1,NAI)READ(1,70) TLREAD(1,70) TLREAD(1,70) TLREAD(

13、1,*)(JS(I,J),J=1,4),I=1,NS)WRITE(2,10)(I,X(I),Y(I),I=1,NJ)WRITE(2,20)(I,JE(I,1),JE(I,2),MT(I),I=1,NE)WRITE(2,30)(I,(AI(I,J),J=1,2),I=1,NAI)WRITE(2,40)(JS(I,J),J=1,4),I=1,NS)FORMAT(/2X,COORDINATES OF JOINTS/6X,JOINT,11X,X,11X,Y/(6X,I4,5X,2F12.4)FORMAT(/2X,INFORMATION OF ELEMENTS/6X,ELEMENT,& 4X,JOINT

14、-I,4X,JOINT-J,5X,TYPE/(2X,4I10)FORMAT(/7X,TYPE,10X,A,12X,I/(8X,I2,5X,2F12.6)FORMAT(/2X,INFORMATION OF SPECIAL JOINTS/6X,JOINT,4X,u,4x,v,4x,ceta/(6X,4I5)FORMAT(A)ENDDetermine number of joint displacements.SUBROUTINE DJN(NJ,NS,JS,JN,N)!计算结构自由度数N，形成数组JN。JN为结点位移分量编号数组，用于存放结点三个位移的位移分量号码DIMENSION JS(NS,4)

15、,JN(NJ,3)DO 10 I=1,NJDO 10 J=1,3JN(I,J)=0DO 20 I=1,NSL=JS(I,1)DO 20 J=1,3JN(L,J)=JS(I,J+1)N=0ID=0DO 30 I=1,NJDO 30 J=1,3IF(JN(I,J)-1)40,50,60N=N+1JN(I,J)=NGOTO 30JN(I,J)=0GOTO 30ID=1CONTINUE!对JN置0!取特殊节点的节点码!将 JS 中位移约束信息 JS(I,2),JS(I,3),JS(I,4)附给 JN(L,1),JN(L,2),JN(L,3)!自由度置0!由条件的小于0，等于0，大于0来执行40，50，

16、60语句!小于0时，JN(I,J)=0，有自由度，自由度累加!等于0, JN(I,J)=1,没有自由度!大于0,结点I是从结点,自由度不增加，RETURNIF(ID.EQ.0)GOTO 80DO 70 I=1,NS!判断是否有主从结点L=JS(I,1)DO 70 J=1,3K=JS(I,J+1)IF(K.LE.1)GOTO 70!若K1,则L是从结点，K是主结点编号JN(L,J)=JN(K,J)!从结点L的第J个位移分量的编号应取主结点K的第J个位移分量的编号70 CONTINUE80 RETURNENDDetermine address of diagonal elements of tot

17、al stiffness matrix.SUBROUTINE ADE(NE,NJ,N,JE,JN,KD,NN)!形成主元地址数组 KD(N)DIMENSION JE(NE,2),JN(NJ,3),KD(N),JC(6)DO 10 I=1,N!主元地址数组KD(N)先置0!对每根杆件循环!形成单元定位向量JC!单元定位向量中的最小非零分量MIN，赋初值为N!对定位向量的六个分量循环!取定位向量的分量!排除零分量!求定位向量的最小非零分量!对定位向量的六个分量循环!取定位向量的分量!排除零分量!计算半带宽!第J列的真正半带宽10KD(I)=0DO 20 M=1,NECALL EJC(M,NE,NJ

18、,JE,JN,JC)MIN=NDO 30 I=1,6J=JC(I)IF(J.EQ.0) GOTO 30IF(J.LT.MIN) MIN=J30 CONTINUEDO 20 I=1,6J=JC(I)IF(J.EQ.0) GOTO 20NW=J-MIN+1IF(NW.GT.KD(J) KD(J)=NW20 CONTINUENN=1!NN置初值为1，NN为存放总刚度矩阵中非零元素的一维数组TK(NN)的元素个数DO 40 J=2,N!对总刚度矩阵K列循环NN=NN+KD(J)KD(J)=NN!累加各列半带宽，形成KD数组 40 CONTINUEEND! Assemble structure stif

19、fness matrix stored as a vector.SUBROUTINE SSM(NE,NJ,NAI,E,N,NN,X,YJE,MT,AI,JN,KD,TK) !形成存储总刚的一维数组 TK(NN)DIMENSION X(NJ),Y(NJ),JE(NE,2),JN(NJ,3),MT(NE),AI(NAI,2),KD(N),TK(NN),KE(6,6),JC(6)DOUBLE PRECISION TK,KE,BL,SI,CODO 10 I=1,NN10TK(I)=0.0!数组 TK(NN)S 0DO 20 M=1,NE !对每根杆件循环CALL SCL(M,NE,NJ,X,YJE,B

20、L,SI,CO)!求单元常数CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE)!形成 global 坐标下的单元刚度矩阵 KE(6,6)计算结构力学程序CALL EJC(M,NE,NJ,JE,JN,JC)!形成单元定位向量JCDO 30 L=1,6!对单元刚度矩阵的行循环I=JC(L)!取单刚的第L行在总刚中的行码IF(I.EQ.0) GOTO 30DO 40 K=1,6!对单元刚度矩阵的列循环J=JC(K)!取单刚的第K列在总刚中的列码IF(J.LT.I) GOTO 40!排除总刚的下三角元素JI=KD(J)-J+I!计算总刚第I行第J列元素在数组TK中的地址TK(

21、JI)=TK(JI)+KE(L,K)!叠加单刚元素到总刚中40CONTINUE30CONTINUE20CONTINUE!WRITE(2,50) TK!50FORMAT(TK,3X,2D15.4)ENDSUBROUTINE SCL (M,NE,NJ,X,YJE,BL,SI,CO) !计算单元常数DIMENSION X(NJ),Y(NJ),JE(NE,2)DOUBLE PRECISION BL,SI,COI=JE(M,1)!取M单元起始结点号J=JE(M,2)!取M单元终止结点号DX=X(J)-X(I)DY=Y(J)-Y(I)BL=SQRT(DX*DX+DY*DY) !计算单元杆长SI=DY/BL

22、!计算单元转角SIN值CO=DX/BL!计算单元转角COS值ENDDecompose total stiffness matrix using Crout method.SUBROUTINE UTDU3(A,NN,ID,N)!对总刚进行UTDU分解，用于解方程组，对应实参(TK,NN,KD,N)DIMENSION A(NN),ID(N)DOUBLE PRECISION A!对列循环!取主对角元素在一维数组A(NN)中的地址码,A(NN)对应TK(NN)DO 30 J=2,NJJ=ID(J)J1=J-1!第J列第一个非零元素的行号!排除下三角元素!对行循环MJ=J-JJ+ID(J1)+1IF(M

23、J.GT.J1) GOTO 30DO 20 I=MJ,JII=ID(I)I1=I-1MIJ=MJIF(I.EQ.J) GOTO 5MI=I-II+ID(I1)+1IF(MIJ.LT.MI) MIJ=MI5S=0.0IF(MIJ.GT.Il) GOTO 15DO 10 K=MIJ,I1KI=II-I+KKK=ID(K)KJ=JJ-J+K10 S=S+A(KI)*A(KK)*A(KJ)15 IF(I.EQ.J) THENA(JJ)=A(JJ)-S !求对角阵D的第J个元素，存放在一维数组A(NN)中ELSEJI=JJ-J+IA(JI)=(A(JI)-S)/A(II)!求三角阵U的第I行第J列元素，

24、存放在一维数组A(NN)中END IF20 CONTINUE30 CONTINUEEND ! Form total joint load vector.SUBROUTINE NLV(NE,NJ,NAI,E,N,NP,NF,X,YE,JN,PJ,PF,MT,AI,P,H) !形成总荷载列阵 P(N)DIMENSION X(NJ),Y(NJ),JE(NE,2),JN(NJ,3),PJ(NP,3),PF(NF,4),MT(NE),AI(NAI,2),& P(N),FO(6),FE(6),JC(6),H(NE)DOUBLE PRECISION P,FO,FE,BL,SI,COINTEGER NODO

25、10 I=1,NP(I)=0.0!总荷载列阵P(N)置010 CONTINUEREAD(1,70)TLREAD(1,70)TLREAD(1,70)TLIF(NP.GT.0) THENREAD(1,*)(NO,(PJ(I,J),J=1,3),I=1,NP)!读入结点荷载信息WRITE(2,20)(PJ(I,J),J=1,3),I=1,NP)20 FORMAT(/2X,JOINT LOAD/6X,JOINT,8X,XYM,12X,LOAD/(6X,F5.0,6X,F5.0,6X,F12.4)DO 30 I=1,NPJ=INT(PJ(I,1)!取结点荷载所在的结点号K=INT(PJ(I,2) !取结

26、点荷载在结点上的方向L=JN(J,K)!取结点荷载的位移编号IF(L.NE.0) P(L)=PJ(I,3)!把结点荷载叠加到总荷载列阵P(N)中30 CONTINUEEND IFREAD(1,70)TLREAD(1,70)TLREAD(1,70)TLIF(NF.GT.0) THENREAD(1,*)(NO,(PF(I,J),J=1,4),I=1,NF)!读入非结点荷载信息WRITE(2,40)(PF(I,J),J=1,4),I=1,NF)40 FORMAT(/2X,NON-JOINT LOAD/6X,ELEMENT,8X,TYPE,8X,LOAD,&12X,C/(6X,F6.0,6X,F6.0

27、,4X,2F12.4)DO 50 I=1,NF!对非结点荷载循环M=INT(PF(I,1)!取非结点荷载所在的单元码CALL SCL(M,NE,NJ,X,YJE,BL,SI,CO)!计算单元常数CALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)将非结点荷载转化为等效结点荷载，存放在FO(6 )中CALL EJC(M,NE,NJ,JE,JN,JC)!形成单元定位数组!将等效结点荷载转化到global坐标下,y轴向下的坐标系。!将等效结点荷载转化到global坐标下!将等效结点荷载转化到global坐标下!将等效结点荷载转化到global坐标下!将等效结点荷载转化

28、到global坐标下!将等效结点荷载转化到global坐标下FE(1)=-FO(1)*CO+FO(2)*SIFE(2)=-FO(1)*SI-FO(2)*COFE(3)=-FO(3)FE(4)=-FO(4)*CO+FO(5)*SIFE(5)=-FO(4)*SI-FO(5)*COFE(6)=-FO(6)DO 60 J=1,6L=JC(J)!将等效结点荷载叠加到总荷载列阵P(N)中IF(L.NE.0) P(L)=P(L)+FE(J)60CONTINUE 50 CONTINUEEND IF70 FORMAT(A)END ! Solve equations and print joint displac

29、ements.SUBROUTINE BACK3(A,NN,B,N,ID,JN,NJ) !解方程组并输出结点位移，存放在数组B(N)中。这里B(N) 数组对应主程序中的P(N)数组DIMENSION A(NN),B(N),ID(N),JN(NJ,3),D(3)DOUBLE PRECISION A,D,BDO 20 J=2,N!对列循环JJ=ID(J)J1=J-1MJ=J-JJ+ID(J1)+1!第J列第一个非零元素的行号IF(MJ.GT.J1) GOTO 20!若MJJ-1，则z的第J个值即为B的第J个值DO 10 I=MJ,J1JI=JJ-J+I10B(J)=B(J)-A(JI)*B(I) !

30、求z的第J个值，存放在B(N)数组。这里B(N)数组对应主程序中的P(N)数组20 CONTINUEDO 30 I=1,NII=ID(I)B(I)=B(I)/A(II)!求y的第J个值，存放在B(N)数组30 CONTINUE3540506070758010NW=1!半带宽置初值为1DO 35 I=2,NIW=ID(I)-ID(I-1)IF(IW.GT.NW) NW=IW!计算最大列半带宽CONTINUEN1=N-1DO 50 I=N1,1,-1!对行循环MIN=I+NW-1IF(N.LT.MIN) MIN=NDO 40 J=I+1,MIN!对列循环JJ=ID(J)MJ=J-JJ+ID(J-1

31、)+1!第J列第一个非零元素的行号IF(I.LT.MJ) GOTO 40 !判断是否为带内元素JI=JJ-J+IB(I)=B(I)-A(JI)*B(J) !求乂的第I个值，存放在B(N)数组。这里x的值即为结点位移CONTINUECONTINUEWRITE(2,60)FORMAT(/2X,JOINT DISPLACEMENTS/5X,JOINT,12X,U,14X,V,11X,ceta)DO 80 I=1,NJDO 70 J=1,3D(J)=0.0L=JN(I,J) !取位移号IF(L.NE.0) D(J)=B(L)!给位移赋值CONTINUEWRITE(2,75) I,D(1),D(2),D

32、(3) !输出第I号结点的位移值FORMAT(2X,I6,4X,3D15.6)CONTINUEENDSUBROUTINE MQN(M,NE,NJ,NAI,N,NF,E,X,YJE,MT,AI,JN,PF,P,F,H)计算单元杆端力，存放在数组 F(6)中DIMENSION X(NJ),Y(NJ),JE(NE,2),AI(NAI,2),P(N),JN(NJ,3),PF(NF,4),JC(6),FO(6),D(6),FD(6),&F(6),KE(6,6),MT(NE),H(NE)DOUBLE PRECISION P,FO,D,FD,F,KE,BL,SI,COCALL SCL(M,NE,NJ,X,Y

33、JE,BL,SI,CO)!求单元常数CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成 global 坐标下的单元刚度矩阵 KE(6,6)CALL EJC(M,NE,NJ,JE,JN,JC)!形成单元定位向量 JCDO 10 I=1,6L=JC(I)!取位移编号D(I)=0.0IF(L.NE.0) D(I)=P(L) !取global坐标下的杆端位移CONTINUEDO 20 I=1,6FD(I)=0.03020504010DO 30 J=1,6FD(I)=FD(I)+KE(I,J)*D(J) !计算杆端力CONTINUECONTINUEF(1)=FD(1)*

34、CO+FD(2)*SI!将杆端力转化到local坐标下F(2)=-FD(1)*SI+FD(2)*CO!将杆端力转化到local坐标下F(3)=FD(3)!将杆端力转化到local坐标下F(4)=FD(4)*CO+FD(5)*SI!将杆端力转化到local坐标下F(5)=-FD(4)*SI+FD(5)*CO!将杆端力转化到local坐标下F(6)=FD(6)!将杆端力转化到local坐标下IF(NF.GT.0)THENDO 40 I=1,NFL=INT(PF(I,1)IF(M.EQ.L)THENCALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)DO 50 J=1

35、,6F(J)=F(J)+FO(J)!如果有非结点荷载，杆端力中要加入等效结点荷载CONTINUEENDIFCONTINUEENDIFENDSUBROUTINE EJC(M,NE,NJ,JE,JN,JC)!形成单元定位向量JC，JC存放的是单元两端点的位移号DIMENSION JE(NE,2),JN(NJ,3),JC(6)I=JE(M,1) J=JE(M,2) DO 10 K=1,3 JC(K)=JN(I,K) JC(K+3)=JN(J,K)CONTINUE END!取M单元起始结点号!取M单元终止结点号!取单元位移号赋值给定位向量JC!取单元位移号赋值给定位向量JCSUBROUTINE ESM

36、(M,NE,NAI,E,MT,AI,BL,SI,CO,KE)!形成 global 坐标下的单元刚度矩阵 KE(6,6)DIMENSION MT(NE),AI(NAI,2),KE(6,6)DOUBLE PRECISION BL,SI,CO,KEI=MT(M)!取M单元截面类型号，MT单元截面类型号C1=E*AI(I,1)/BLC2=2.0*E*AI(I,2)/BLC3=3.0*C2/BLC4=2.0*C3/BLS1=C1*CO*CO+C4*SI*SIS2=(C1-C4)*SI*COS3=C3*SIS4=C1*SI*SI+C4*CO*COS5=C3*COS6=C2KE(1,1)=S1KE(1,2)

37、=S2KE(1,3)=-S3KE(1,4)=-S1KE(1,5)=-S2KE(1,6)=-S3KE(2,2)=S4KE(2,3)=S5KE(2,4)=-S2KE(2,5)=-S4KE(2,6)=S5KE(3,3)=2.0*S6KE(3,4)=S3KE(3,5)=-S5KE(3,6)=S6KE(4,4)=S1KE(4,5)=S2KE(4,6)=S3KE(5,5)=S4KE(5,6)=-S5KE(6,6)=2.0*S6DO 10 I=1,5DO 10 J=I+1,610 KE(J,I)=KE(I,J)! WRITE(2,20)(KE(I,J),J=1,6),I=1,6)!20FORMAT(KE,3

38、X,6D15.4)ENDSUBROUTINE EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)将非结点荷载转化为等效结点荷载存放在FO(6 )中DIMENSION PF(NF,4),MT(NE),AI(NAI,2),FO(6),H(NE)DOUBLE PRECISION BL,FOIF(PF(I,1)-M).GT.0.0)GOTO 110 !判断此非结点荷载是否为温度荷载NO=INT(PF(I,2)!取此非结点荷载的类型号Q=PF(I,3)C=PF(I,4)IF(NO.GT.6) GOTO 3!判断此非结点荷载是否为1 6号荷载B=BL-CC1=C/BLC2=C1*C

39、1C3=C1*C23DO 5 J=1,65FO(J)=0.0GOTO(10,20,30,40,50,60,70,80,90,100),NO !对不同的非结点荷载分别计算其等效结点荷载10FO(2)=-Q*C*(1.0-C2+C3/2.0)FO(3)=-Q*C*C*(0.5-2.0*C1/3.0+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)RETURN20FO(2)=-Q*B*B*(1.0+2.0*C1)/BL/BLFO(3)=-Q*C*B*B/BL/BLFO(5)=-Q*C2*(1.0+2.0*B/BL)F

40、O(6)=Q*C2*BRETURN30FO(2)=6.0*Q*C1*B/BL/BLFO(3)=Q*B*(2.0-3.0*B/BL)/BLFO(5)=-6.0*Q*C1*B/BL/BLFO(6)=Q*C1*(2.0-3.0*C1)RETURN40FO(2)=-Q*C*(0.5-0.75*C2+0.4*C3)FO(3)=-Q*C*C*(1.0/3.0-0.5*C1+0.2*C2)FO(5)=-Q*C*C2*(0.75-0.4*C1)FO(6)=Q*C*C*C1*(0.25-0.2*C1)RETURN50FO(1)=-Q*C*(1.0-0.5*C1)FO(4)=-0.5*Q*C*C1RETURN60

41、FO(1)=-Q*B/BLFO(4)=-Q*C170 L=INT(C)K=MT(M)S=E*AI(K,1)*Q/BLFO(L)=SIF(L.EQ.1)FO(4)=-SIF(L.EQ.4)FO(1)=-SRETURN80 L=INT(C)K=MT(M)FO(L)=12.0*E*AI(K,2)*Q/BL/BL/BLIF(L.EQ.2)FO(5)=-FO(2)IF(L.EQ.5)FO(2)=-FO(5)FO(3)=0.5*BL*FO(2)90100110FO(6)=FO(3)RETURNL=INT(C)K=MT(M)S=2.0*E*AI(K,2)*Q/BLFO(L)=2.0*SIF(L.EQ.3)F

42、O(6)=SIF(L.EQ.6)FO(3)=SFO(2)=3.0*S/BLFO(5)=-FO(2)RETURNL=INT(C)K=MT(M)S=E*AI(K,1)*Q/BLFO(L)=SIF(L.EQ.1)FO(4)=-SIF(L.EQ.4)FO(1)=-SRETURNK=MT(M)FO(1)=0.5*E*AI(K,1)*PF(I,4)*(PF(I,2)+PF(I,3)FO(2)=0.0FO(3)=-E*AI(K,2)*PF(I,4)*(PF(I,2)+PF(I,3)/H(M)FO(4)=-FO(1)FO(5)=0.0FO(6)=-FO(3)END*算例：如图所示平面刚架，尺寸见图。其中梁的截

43、面积气=0.2M2,/广0.01M4；梁的截面积A. =0.25M2 ,I = 0.06M4，E = 3x107N/M2。支座沉降当做单独的工况计算。输入数据文件（INDAT.DAD：COMPOSITED PLANE FRAME STRUCTURE ANALYSISNE NJ NS NAI NL E673223.0E7士 士 士 士 士 士 士 士 士 士 士 士 士 士 士 W W W W W 士 士 士 士 士 士 士 士 士 士 士 士v 士 士 士 士 士 士 士 士 士 士 士 士& 4-女个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个

44、个个个个个个个个个COORDINATES OF JOINTJOINTX(I)Y(I)10.00.020.00.030.05.048.00.058.05.060.010.078.010.0士 士 士 士 士 士 士 W W W v 士 士 士 士 士 士 士 士 士 士 士 士 士 士V 士 士 士 士 士 士 士 士 士 士 士 士 士 士v 士 士 士 士 士 士 士INFORMATION OF ELEMENTSELEMENTJOINT-IJOINT-JAI-TYPEH124102132034520435105362065720士 士 士 士 士 士 士 W W W 士 士 士 士 士 士

45、士 士 士 士 士 士 士 士 vV 士 士 士 士 士 士 士 士 士 士 士 士 士 士 vv 士 士 士 士 士 士 士个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个INFORMATION OF A & I OF ELEMENTTYPE-NO.AI1 0.20.012 0.250.006士 士 士 士 士 士 士 W W W v 士 士 士 士 士 士 士 士 士 士 士 士 士 士 v V 士 士 士 士 士 士 士 士 士 士 士 士 士 士 v v 士 士 士 士 士 士 士 个个个个个个个个个个个个个个个个个个个

46、个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个INFORMATION OF SPECIAL JOINTJOINTUVCETA122061117111士 士 士 士 士 士 士 W W W v 士 士 士 士 士 士 士 士 士 士 士 士 士 士 vV 士 士 士 士 士 士 士 士 士 士 士 士 士 士 vv 士 士 士 士 士 士 士个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个FIRST LOAD CASE INFORMATIONNPNF21士 士 士 士 士 W W W v 士 士 士

47、士 士个个个个个个个个个个个个个个个JOINT FORCE NO. 1 2JOINT 1.0 3.0COORDINATE1.01.0LOAD18.024.0必必必必必必e e e 必必必必必 个个个个个个个个个个个个个个个NON-JOINT FORCELOAD20.08.0NO. ELEMENT TYPE11.01.0必必必必必必必必e e e 必必必必必必必必必必必必必必必e e e 必必必必必必必必必必必必必必必e e e 必必必必必必必 个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个个SECOND LOAD CASE IN

48、FORMATIONNPNF02必必必必必必e e e 必必必必必*JOINT FORCENO. JOINT COORDINATE LOAD士 士 士 士 士 W W W v 士 士 士 士 士NON-JOINT FORCENO.ELEMENTTYPELOADC16.07.00.034.026.08.0-0.025.0输出计算结果(O UTDAT.DAT)PLANE FRAME STRUCTURE ANALYSIS士 士 W W W W W v 士 士 士 士 士 士 士 士 士 士 士 士 vv 士 士CONTROL PARAMETERS OF STRUCTURENE= 6 NJ= 7 NS=

49、 3 NAI= 2NL= 2 E= .3000E+08COORDINATES OF JOINTSJOINTXY1.0000.00002.0000.00003.00005.000048.0000.000058.00005.00006.000010.000078.000010.0000INFORMATION OF ELEMENTSELEMENT123456JOINT-I214335JOINT-J435567TYPE122122TYPEAI1.200000.0100002.250000.006000INFORMATION OF SPECIAL JOINTSJOINTuvceta12206111711

50、1LOAD DATALOAD CASE= 1NP=2NF=1JOINT LOADJOINTXYMLOAD1.1.18.00003.1.24.0000NON-JOINT LOADELEMENT1.TYPE1.LOAD 20.0000C8.0000JOINT DISPLACEMENTSJOINTU1.198620D-022.198620D-023.175088D-024.195921D-025.175352D-026.000000D+007.000000D+00V.846572D-04.846572D-04.383719D-04.128676D-03.682948D-04.000000D+00.0

51、00000D+00ceta-.485791D-05.105183D-02.150910D-03-.664935D-03.278622D-03.000000D+00.000000D+00MEMBER-END FORCES OF ELEMENTSELEMENTNVM1N1=.2024D+02V1=-.6943D+02M1=-.4336D-08N2=-.2024D+02V2=-.9057D+02M2=.8458D+022N1=.6943D+02V1=.2243D+01M1=.8882D-14N2=-.6943D+02V2=-.2243D+01M2=.1122D+023N1=.9057D+02V1=-

52、.2024D+02M1=-.8458D+02N2=-.9057D+02V2=.2024D+02M2=-.1664D+024N1=-.1979D+01V1=.1187D+02M1=.4269D+02N2=.1979D+01V2=-.1187D+02M2=.5227D+025N1=.5756D+02V1=-.2374D+02M1=-.5391D+02N2=-.5756D+02V2=.2374D+02M2=-.6477D+026N1=.1024D+03V1=-.1826D+02M1=-.3563D+02N2=-.1024D+03V2=.1826D+02M2=-.5569D+02LOAD DATALO

53、AD CASE= 2NP= 0 NF= 2NON-JOINT LOADELEMENTTYPE LOADC6.7.03004.0000MEMBER-END FORCES OF ELEMENTSELEMENTNVM1N1=.2843D+02V1=.4692D+00M1=-.2611D-08N2=-.2843D+02V2=-.4692D+00M2=.3754D+012N1=-.4692D+00V1=.2843D+02M1=.1705D-12N2=.4692D+00V2=-.2843D+02M2=.1422D+033N1=.4692D+00V1=-.2844D+02M1=-.3754D+01N2=-.

54、4692D+00V2=.2844D+02M2=-.1384D+034N1=-.1429D+03V1=-.2946D+02M1=-.1935D+02N2=.1429D+03V2=.2946D+02M2=-.2163D+035N1=.2899D+02V1=-.1144D+03M1=-.1228D+03N2=-.2899D+02V2=.1144D+03M2=-.4493D+036N1=-.2899D+02V1=.1144D+03M1=.3548D+03N2=.2899D+02V2=-.1144D+03M2=.2174D+036.8.-.02005.0000JOINT DISPLACEMENTSJOINTU1.340494D-012.340494D-013.179583D-014.340115D-015.181488D-016.000000D+007.000000D+00V.190150D-04.190150D-04.193278D-04.299810D-01.299807D-01.000000D+00.000000D+00ceta.256001D-02.372856D-02.453461D-02.377861D-02.190823D-02.000000D+00.000000D+00