平面桁架计算程序(共9页)

1、精选优质文档-倾情为你奉上! PTA说明: 1,数据输入文件为DATAIN.TXT中,数据依次为(NE,NJ,NR,NP,NP); (X,Y)*NJ ; ! (IJ1,IJ2,A,EI)*NE;(JR1,JR2,JR3,JR4)*NR;(PJ1,PJ2,PJ3)*NP;(PF1,PF2,PF3,PF4)*NF ! 2,NE单元总数；NJ节点总数；NR约束总数；NP节点荷载总数；NF非节点荷载数；X,Y节点坐标； ! IJ(NE,1),IJ(NE,2) 单元期终点；A 单元面积 ；ZI 截面惯性矩；JR(NR,1)约束结点号; ! JR(NE,2:4)横向竖向转动约束(1); PJ(NP,1:3

2、)分别为结点号、荷载类型、荷载值； ! PF(NF,1:4)分别为单元号、荷载类型，荷载值、荷载据起点距离 !主程序 PROGRAM Plane_Truss_Analysis implicit none integer NE,NJ,NP,NF,NR,N real E integer,allocatable,DIMENSION(:,:):IJ,JR real,allocatable,DIMENSION(:):A,P,X,Y !动态数组定义 real,allocatable,DIMENSION(:,:):PJ,PF,TK OPEN(1,FILE='datain.TXT',STATU

3、S='OLD') !打开文件大datain.txt并存到标号1中 OPEN(2,FILE='dataout.TXT',STATUS='NEW') !建立文件dataout存储计算结果 READ(1,*)NE,NJ,NR,NP,NF,E N=NJ*2 !位移总数 allocate(X(1:NJ),Y(1:NJ),IJ(1:NE,2),JR(1:NR,4),A(1:NE),PJ(1:NP,3),PF(1:NF,4),TK(1:N,1:N),P(1:N) WRITE(2,10)NE,NJ,NR,NP,NF,E !打印表头10 FORMAT(/1X,&#

4、39; *平面桁架力计算PTA*'/4X,'单元数NE=',I2,12X,'结点数NJ=',I2,14X,'支座数NR=',I2,/4X,'结点荷载数NP=',I2,8X,'非节点荷载数NF=',I2,8X,'弹性模量E=',E12.4) CALL INPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,JR,PJ,PF) !数据输入 CALL TSM(NE,NJ,E,X,Y,IJ,A,TK,N) !形成结构原始刚度矩阵 CALL JLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF

5、,P,N) !形成结构综合节点荷载阵列 CALL ISC(NR,JR,TK,P,N) !引入结构约束条件 CALL GAUSS(TK,P,N) !高斯消去法计算结构坐标下的节点位移 CALL MVN(NE,NJ,NF,E,X,Y,IJ,A,PF,P,N) !计算单元杆端内力 CLOSE(1) CLOSE(2) deallocate(X,Y,IJ,JR,A,P,PJ,PF,TK) END PROGRAM Plane_Truss_Analysis !原始数据输入 SUBROUTINE INPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,JR,PJ,PF) DIMENSION X(NJ),Y

6、(NJ),IJ(NE,2),A(NE),JR(NR,4),PJ(NP,3),PF(NF,4) READ(1,*) (X(I),Y(I),I=1,NJ) !从1中顺序读取节点坐标 READ(1,*) (IJ(I,1),IJ(I,2),A(I),I=1,NE) !从1中顺序读取单元期终编号，截面积，惯性矩 READ(1,*) (JR(I,J),J=1,3),I=1,NR) !从1中顺序读取支座信息 IF(NP.GT.0)READ(1,*)(PJ(I,J),J=1,3),I=1,NP) !判断并从1中顺序读取结点荷载编码，方向，大小 IF(NF.GT.0)READ(1,*)(PF(I,J),J=1,

7、4),I=1,NF) !判断并从1中顺序读取非结点荷载编码，荷载类型，大小，位子参数 WRITE(2,10) (I,X(I),Y(I),I=1,NJ) !将结点号、结点坐标，按10格式写入2中 WRITE(2,20) (I,IJ(I,1),IJ(I,2),A(I),I=1,NE) !将单元号、单元起终号、单元截面性质，按20格式写入2中 WRITE(2,30) (JR(I,J),J=1,3),I=1,NR) !将结点支座信息按格式30写入2中 IF(NP.GT.0)WRITE(2,40)(PJ(I,J),J=1,3),I=1,NP) IF(NF.GT.0)WRITE(2,50)(PF(I,J)

8、,J=1,4),I=1,NF) 10 FORMAT(/4X,'1>. 结点坐标'/8X,'结点号',12X,'X',12X,'Y'/(6X,I4,5X,2F12.4) 20 FORMAT(/4X,'2>. 单元信息'/8X,'单元号',6X,'单元起点I',6X,'单元终点J',10X,'A'/(8X,I4,10X,I4,10X,I4,6x,F12.4) 30 FORMAT(/4X,'3>. 结点约束信息'/8X,&#

9、39;约束结点',5X,'水平XR',5X,'竖直YR'/(4X,3I10) 40 FORMAT(/4X,'4>. 节点荷载'/8X,'结点号',10X,'方向',10X,'荷载值'/(9X,F5.0,9X,F5.0,10X,F12.4) 50 FORMAT(/4X,'5>. 非结点荷载'/8X,'单元号',6X,'荷载类型',6X,'荷载值',6X,'荷载位置'/(6X,F6.0,6X,F6.0,4

10、X,2F12.4) END SUBROUTINE INPUT !计算单元常数 SUBROUTINE LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO) DIMENSION X(NJ),Y(NJ),IJ(NE,2) I=IJ(M,1) J=IJ(M,2) DX=X(J)-X(I) DY=Y(J)-Y(I) BL=SQRT(DX*DX+DY*DY) SI=DY/BL CO=DX/BL END SUBROUTINE LSC !计算局部坐标系下单刚 SUBROUTINE ESM(M,NE,E,A,BL,EK) DIMENSION A(NE),EK(4,4) EK=0.0 EK(1,1)=E*A(

11、M)/BL EK(1,3)=-EK(1,1) EK(3,1)=-EK(1,1) EK(3,3)=EK(1,1) END SUBROUTINE ESM!坐标转换矩阵 SUBROUTINE CTM(SI,CO,T)DIMENSION T(4,4)T=0.0T(1,1:2)=(/CO,SI/)T(2,1:2)=(/-SI,CO/)T(3:4,3:4)=T(1:2,1:2)returnend SUBROUTINE CTM !计算结构坐标系下单刚 SUBROUTINE TTKT(EK,T) DIMENSION EK(4,4),T(4,4),TE(4,4) DO I=1,4 DO J=1,4 TE(I,J

12、)=0.0 DO K=1,4 TE(I,J)=TE(I,J)+T(K,I)*EK(K,J) ENDDO ENDDO ENDDO DO I=1,4 DO J=1,4 EK(I,J)=0.0 DO K=1,4 EK(I,J)=EK(I,J)+TE(I,K)*T(K,J) ENDDO ENDDO ENDDO END SUBROUTINE TTKT !形成结构的原始刚度矩阵SUBROUTINE TSM(NE,NJ,E,X,Y,IJ,A,TK,N)DIMENSION X(NJ),Y(NJ),IJ(NE,2),A(NE),TK(N,N),EK(4,4),T(4,4),LV(4)TK=0.0 !N*N阶原始

13、刚度矩阵元素全赋值0 DO M=1,NE !以单元号循环，逐个将元素放入结构的原始刚度矩阵中 CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO) !计算第M个单元的常数 CALL ESM(M,NE,E,A,BL,EK) !计算第M个单元在局部坐标系下的单元刚度矩阵 CALL CTM(SI,CO,T) !计算第M个单元的坐标转换矩阵 CALL TTKT(EK,T) !计算第M个单元在结构坐标系下的单元刚度矩阵 DO K=1,2 !求单元位移分量在位移列阵中的编号 LV(K)=2*(IJ(M,1)-1)+K !单元定位向量 LV(2+K)=2*(IJ(M,2)-1)+K ENDDO

14、 DO L=1,4 !单刚按行循环 I=LV(L) !单刚第L行放入总刚第I行 DO K=1,4 !单刚按列循环 J=LV(K) !单刚第K列放入总刚第J列 TK(I,J)=TK(I,J)+EK(L,K) !对号入座 ENDDO ENDDO ENDDO END SUBROUTINE TSM !计算单元固端力的子程序 SUBROUTINE EFF(L,PF,NF,BL,FO) DIMENSION PF(NF,4),FO(6) NO=INT(PF(L,2) Q=PF(L,3) ! 非节点荷载类型号 C=PF(L,4) B=BL-C C1=C/BL C2=C1*C1 C3=C1*C2 DO I=1,

15、4 FO(I)=0.0 end doGOTO(10,20),NO10 FO(1)=-Q*BL/2 !垂直于单元的均布荷载 FO(3)=FO(1) RETURN20 FO(1)=-Q*B/BL !垂直于单元的集中荷载 FO(3)=-Q*C1 END !形成结构综合节点荷载列阵SUBROUTINE JLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N)DIMENSION X(NJ),Y(NJ),IJ(NE,2),PJ(NP,3),PF(NF,4),P(N),FO(4),PE(4)P=0.0IF (NP.GT.0)THEN DO I=1,NP J=INT(PJ(I,1) !节点号 L=2

16、*(J-1)+INT(PJ(I,2) P(L)=PJ(I,3) ENDDO ENDIF IF (NF.GT.0)THEN DO L=1,NF M=INT(PF(L,1) CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO) CALL EFF(L,PF,NF,BL,FO) PE(1)=-FO(1)*CO+FO(2)*SI PE(2)=-FO(1)*SI-FO(2)*CO PE(3)=-FO(3)*CO+FO(4)*SI PE(4)=-FO(3)*SI-FO(4)*CO I=IJ(M,1) J=IJ(M,2) P(2*I-1)=P(2*I-1)+PE(1) P(2*I)=P(2*I)

17、+PE(2) P(2*J-1)=P(2*J-1)+PE(3) P(2*J)=P(2*J)+PE(4) ENDDO ENDIF END SUBROUTINE JLP !主1副零法引入结构的支承条件 SUBROUTINE ISC(NR,JR,TK,P,N) DIMENSION TK(N,N),P(N),JR(NR,4) DO I=1,NR J=JR(I,1) DO K=1,2 IF(JR(I,K+1).NE.0)THEN L=2*(J-1)+K DO JJ=1,N TK(L,JJ)=0.0 TK(JJ,L)=0.0 ENDDO TK(L,L)=1.0 P(L)=0.0 ENDIF ENDDO EN

18、DDO END SUBROUTINE ISC !输出结点位移 计算单元杆端力 SUBROUTINE MVN(NE,NJ,NF,E,X,Y,IJ,A,PF,P,N) DIMENSION X(NJ),Y(NJ),IJ(NE,2),A(NE),P(N),PF(NF,4),FO(4),F(4),D(4),TD(4),T(4,4),EK(4,4) WRITE(2,10)10 FORMAT(/4X,'6>. 结点位移值'/8X,'结点号',10X,'水平U',10X,'竖直V') DO I=1,NJ WRITE(2,15)I,P(2*I

19、-1),P(2*I)15 FORMAT(6X,I6,6X,2F12.4) ENDDO WRITE(2,25)25 FORMAT(/4X,'7>. 单元杆端力'/8X,'单元号',9X,'轴力N') DO M=1,NE CALL LSC(M,NE,NJ,X,Y,IJ,BL,SI,CO) CALL ESM(M,NE,E,A,BL,EK) CALL CTM(SI,CO,T) I=IJ(M,1) J=IJ(M,2) DO K=1,2 D(K)=P(2*(I-1)+K) D(K+2)=P(2*(J-1)+K) ENDDO DO I=1,4 TD(I)

20、=0.0 DO J=1,4 TD(I)=TD(I)+T(I,J)*D(J) ENDDO ENDDO DO I=1,4 F(I)=0.0 DO J=1,4 F(I)=F(I)+EK(I,J)*TD(J) ENDDO ENDDO IF(NF.GT.0)THEN DO L=1,NF I=INT(PF(L,1) IF(M.EQ.I)THEN CALL EFF(L,PF,NF,BL,FO) DO J=1,4 F(J)=F(J)+FO(J) ENDDO ENDIF ENDDO ENDIF WRITE(2,80)M,F(1),F(3)80 FORMAT(2X,I8,4X,'NX1=',F12

21、.4,3X,'NX2=',F12.4) ENDDO END SUBROUTINE MVN !GAUSS解法 SUBROUTINE GAUSS(A,B,N) DIMENSION A(N,N),B(N) DO K=1,N-1 DO I=K+1,N A1=A(K,I)/A(K,K) DO J=K+1,N A(I,J)=A(I,J)-A1*A(K,J) ENDDO B(I)=B(I)-A1*B(K) ENDDO ENDDO B(N)=B(N)/A(N,N) DO I=N-1,1,-1 DO J=I+1,N B(I)=B(I)-A(I,J)*B(J) ENDDO B(I)=B(I)/A(I,I) ENDDO END SUBROUTINE GAUSS 例子：1：数据文件8,6,3,2,0,10,0,1,1,2,0,3,1,4,0,2,-0.1,2,11,3,12,3,12,4,13,4,13,5,14,5,13,6,1,1,16,1,15,0,12,2,-204,2,-20 2：运行结果*平面桁架力计算PTA* 单元数NE= 8 结点数NJ= 6 支座数NR= 3 结点荷载数NP= 2 非节点荷载数NF= 0 弹性模量E= .10