平面桁架计算程序

!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)分别为结点号、荷载类型、荷载值；!PF(NF,1:4)分别为单元号、荷载类型，荷载值、荷载据起点距离!主程序PROGRAMPlane_Truss_AnalysisimplicitnoneintegerNE,NJ,NP,NF,NR,NrealEinteger,allocatable,DIMENSION(:,:)::IJ,JRreal,allocatable,DIMENSION(:)::A,P,X,Y!动态数组定义real,allocatable,DIMENSION(:,:)::PJ,PF,TKOPEN(1,FILE='datain.TXT',STATUS='OLD')!打开文件大datain.txt并存到标号1中OPEN(2,FILE='dataout.TXT',STATUS='NEW')!建立文件dataout存储计算结果READ(1,*)NE,NJ,NR,NP,NF,EN=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!打印表头10FORMAT(/1X,'***********平面桁架力计算PTA***********'//4X,'单元数NE=',I2,12X,'结点数NJ=',I2,14X,'支座数NR=',I2,/4X,'结点荷载数NP=',I2,8X,'非节点荷载数NF=',I2,8X,'弹性模量E=',E12.4)CALLINPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,JR,PJ,PF)!数据输入CALLTSM(NE,NJ,E,X,Y,IJ,A,TK,N)!形成结构原始刚度矩阵CALLJLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N)!形成结构综合节点荷载阵列CALLISC(NR,JR,TK,P,N)!引入结构约束条件CALLGAUSS(TK,P,N)!高斯消去法计算结构坐标下的节点位移CALLMVN(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)ENDPROGRAMPlane_Truss_Analysis!原始数据输入SUBROUTINEINPUT(NE,NJ,NR,NP,NF,X,Y,IJ,A,JR,PJ,PF)DIMENSIONX(NJ),Y(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,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),J=1,4),I=1,NF)10FORMAT(//4X,'1.结点坐标'/8X,'结点号',12X,'X',12X,'Y'/(6X,I4,5X,2F12.4))20FORMAT(//4X,'2.单元信息'/8X,'单元号',6X,'单元起点I',6X,'单元终点J',10X,'A'/(8X,I4,10X,I4,10X,I4,6x,F12.4))30FORMAT(//4X,'3.结点约束信息'/8X,'约束结点',5X,'水平XR',5X,'竖直YR'/(4X,3I10))40FORMAT(//4X,'4.节点荷载'/8X,'结点号',10X,'方向',10X,'荷载值'/(9X,F5.0,9X,F5.0,10X,F12.4))50FORMAT(//4X,'5.非结点荷载'/8X,'单元号',6X,'荷载类型',6X,'荷载值',6X,'荷载位置'/(6X,F6.0,6X,F6.0,4X,2F12.4))ENDSUBROUTINEINPUT!计算单元常数SUBROUTINELSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)DIMENSIONX(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/BLCO=DX/BLENDSUBROUTINELSC!计算局部坐标系下单刚SUBROUTINEESM(M,NE,E,A,BL,EK)DIMENSIONA(NE),EK(4,4)EK=0.0EK(1,1)=E*A(M)/BLEK(1,3)=-EK(1,1)EK(3,1)=-EK(1,1)EK(3,3)=EK(1,1)ENDSUBROUTINEESM!坐标转换矩阵SUBROUTINECTM(SI,CO,T)DIMENSIONT(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)returnendSUBROUTINECTM!计算结构坐标系下单刚SUBROUTINETTKT(EK,T)DIMENSIONEK(4,4),T(4,4),TE(4,4)DOI=1,4DOJ=1,4TE(I,J)=0.0DOK=1,4TE(I,J)=TE(I,J)+T(K,I)*EK(K,J)ENDDOENDDOENDDODOI=1,4DOJ=1,4EK(I,J)=0.0DOK=1,4EK(I,J)=EK(I,J)+TE(I,K)*T(K,J)ENDDOENDDOENDDOENDSUBROUTINETTKT!形成结构的原始刚度矩阵SUBROUTINETSM(NE,NJ,E,X,Y,IJ,A,TK,N)DIMENSIONX(NJ),Y(NJ),IJ(NE,2),A(NE),TK(N,N),EK(4,4),T(4,4),LV(4)TK=0.0!N*N阶原始刚度矩阵元素全赋值0DOM=1,NE!以单元号循环，逐个将元素放入结构的原始刚度矩阵中CALLLSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)!计算第M个单元的常数CALLESM(M,NE,E,A,BL,EK)!计算第M个单元在局部坐标系下的单元刚度矩阵CALLCTM(SI,CO,T)!计算第M个单元的坐标转换矩阵CALLTTKT(EK,T)!计算第M个单元在结构坐标系下的单元刚度矩阵DOK=1,2!求单元位移分量在位移列阵中的编号LV(K)=2*(IJ(M,1)-1)+K!单元定位向量LV(2+K)=2*(IJ(M,2)-1)+KENDDODOL=1,4!单刚按行循环I=LV(L)!单刚第L行放入总刚第I行DOK=1,4!单刚按列循环J=LV(K)!单刚第K列放入总刚第J列TK(I,J)=TK(I,J)+EK(L,K)!对号入座ENDDOENDDOENDDOENDSUBROUTINETSM!计算单元固端力的子程序SUBROUTINEEFF(L,PF,NF,BL,FO)DIMENSIONPF(NF,4),FO(6)NO=INT(PF(L,2))Q=PF(L,3)!非节点荷载类型号C=PF(L,4)B=BL-CC1=C/BLC2=C1*C1C3=C1*C2DOI=1,4FO(I)=0.0enddoGOTO(10,20),NO10FO(1)=-Q*BL/2!垂直于单元的均布荷载FO(3)=FO(1)RETURN20FO(1)=-Q*B/BL!垂直于单元的集中荷载FO(3)=-Q*C1END!形成结构综合节点荷载列阵SUBROUTINEJLP(NE,NJ,NP,NF,X,Y,IJ,PJ,PF,P,N)DIMENSIONX(NJ),Y(NJ),IJ(NE,2),PJ(NP,3),PF(NF,4),P(N),FO(4),PE(4)P=0.0IF(NP.GT.0)THENDOI=1,NPJ=INT(PJ(I,1))!节点号L=2*(J-1)+INT(PJ(I,2))P(L)=PJ(I,3)ENDDOENDIFIF(NF.GT.0)THENDOL=1,NFM=INT(PF(L,1))CALLLSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)CALLEFF(L,PF,NF,BL,FO)PE(1)=-FO(1)*CO+FO(2)*SIPE(2)=-FO(1)*SI-FO(2)*COPE(3)=-FO(3)*CO+FO(4)*SIPE(4)=-FO(3)*SI-FO(4)*COI=IJ(M,1)J=IJ(M,2)P(2*I-1)=P(2*I-1)+PE(1)P(2*I)=P(2*I)+PE(2)P(2*J-1)=P(2*J-1)+PE(3)P(2*J)=P(2*J)+PE(4)ENDDOENDIFENDSUBROUTINEJLP!主1副零法引入结构的支承条件SUBROUTINEISC(NR,JR,TK,P,N)DIMENSIONTK(N,N),P(N),JR(NR,4)DOI=1,NRJ=JR(I,1)DOK=1,2IF(JR(I,K+1).NE.0)THENL=2*(J-1)+KDOJJ=1,NTK(L,JJ)=0.0TK(JJ,L)=0.0ENDDOTK(L,L)=1.0P(L)=0.0ENDIFENDDOENDDOENDSUBROUTINEISC!输出结点位移计算单元杆端力SUBROUTINEMVN(NE,NJ,NF,E,X,Y,IJ,A,PF,P,N)DIMENSIONX(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)10FORMAT(//4X,'6.结点位移值'/8X,'结点号',10X,'水平U',10X,'竖直V')DOI=1,NJWRITE(2,15)I,P(2*I-1),P(2*I)15FORMAT(6X,I6,6X,2F12.4)ENDDOWRITE(2,25)25FORMAT(/4X,'7.单元杆端力'/8X,'单元号',9X,'轴力N')DOM=1,NECALLLSC(M,NE,NJ,X,Y,IJ,BL,SI,CO)CALLESM(M,NE,E,A,BL,EK)CALLCTM(SI,CO,T)I=IJ(M,1)J=IJ(M,2)DOK=1,2D(K)=P