一平面桁架静力分析源程序PTSAPFOR.docx
《一平面桁架静力分析源程序PTSAPFOR.docx》由会员分享,可在线阅读,更多相关《一平面桁架静力分析源程序PTSAPFOR.docx(55页珍藏版)》请在冰豆网上搜索。
一平面桁架静力分析源程序PTSAPFOR
一、平面桁架静力分析源程序——PTSAP.FOR
C
C
CPLANETRUSSSTRUCTURALANALYSISPROGRAM
C
CMAINPROGRAM
DIMENSIONJE(2,100),JS(3,50),X(100),Y(100),E(100),A(100),
&PJ(3,50),P(200),F(100)
REALKB(200,50)
OPEN(5,FILE='PTR.DAT',STATUS='OLD')
OPEN(6,FILE='PTW.DAT',STATUS='NEW')
CINPUTCONTROLPARAMETERSOFSTRUCTURE
READ(5,*)NJ,NE,NS,NPJ
CALLDAT(NJ,NE,NS,NPJ,NW,NEQ,X,Y,JE,JS,PJ,A,E)
CALLSSM(NEQ,NE,NJ,NW,JE,X,Y,E,A,KB)
CALLFPJ(NEQ,NPJ,PJ,P)
CALLCSE(NEQ,NS,NW,JS,KB,P)
CALLBGS(NEQ,NJ,NW,KB,P)
CALLTRS(NEQ,NE,NJ,NS,NPJ,JE,JS,X,Y,E,A,PJ,P,F)
CLOSE(5)
CLOSE(6)
END
C
CSUBPROGRAM-DATAINPUTANDPRINT
SUBROUTINEDAT(NJ,NE,NS,NPJ,NW,NEQ,X,Y,JE,JS,PJ,A,E)
DIMENSIONJE(2,NE),JS(3,NS),X(NJ),Y(NJ),E(NE),A(NE),PJ(3,NPJ)
READ(5,*)(X(I),Y(I),I=1,NJ)
READ(5,*)((JE(I,J),I=1,2),J=1,NE)
READ(5,*)((JS(I,J),I=1,3),J=1,NS)
READ(5,*)((PJ(I,J),I=1,3),J=1,NPJ)
READ(5,*)(A(I),E(I),I=1,NE)
NW=0
DO70I=1,NE
IW=IABS(JE(1,I)-JE(2,I))
IF(IW.GT.NW)NW=IW
70CONTINUE
NW=2*(NW+1)
NEQ=2*NJ
WRITE(6,75)
75FORMAT(1X,'PLANETRUSSSTRUCTUREANALYSIS'/)
COUTPUTCONTROLPARAMETERSOFSTRUCTURE
WRITE(6,130)
WRITE(6,140)NE,NJ,NS,NPJ,NEQ,NW
COUTPUTDATAOFELEMENTSANDJOINTS
WRITE(6,150)
WRITE(6,78)
78FORMAT(/1X,'NODALCOODINATES')
WRITE(6,80)(I,X(I),Y(I),I=1,NJ)
WRITE(6,90)(J,(JE(I,J),I=1,2),A(J),E(J),J=1,NE)
WRITE(6,100)((JS(I,J),I=1,3),J=1,NS)
WRITE(6,110)((PJ(I,J),I=1,3),J=1,NPJ)
80FORMAT(7X,'NODE',7X,'X',14X,'Y',/(1X,I10,2E15.6))
90FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',7X,'A',14X,'E'
&/(1X,3I10,2E15.6))
100FORMAT(/1X,'CONSTRAINTEDNODES'/7X,'NODE',
&8X,'XX',8X,'YY'/(1X,3I10))
110FORMAT(/1X,'JOINTLOADS'/7X,'NODE',8X,'PX',11X,'PY'
&/(1X,F10.0,2E15.6))
130FORMAT(1X,'STRUCTURALPARAMETERS'//1X,'TOTAL',1X,'NO.OF')
140FORMAT(1X,'ELEMENTS=',I5/1X,'JOINTS=',I5/1X,'CONSTRAINTEDJOINTS=
&',I5/1X,'JOINTLOADS=',I5/1X,'DEGREESOFFREEDOM=',I5/
&1X,'BANDWIDTH=',I5)
150FORMAT(/1X,'DATAOFELEMENTS&JOINTS')
RETURN
END
CSUBPROGRAM-1WHOLESTIFFNESSMATRIX
SUBROUTINESSM(NEQ,NE,NJ,NW,JE,X,Y,E,A,KB)
DIMENSIONJE(2,NE),X(NJ),Y(NJ),E(NE),A(NE),KB(NEQ,NW),KE(4,4)
REALKB,KE
DO10I=1,NEQ
DO10J=1,NW
KB(I,J)=0.0
10CONTINUE
DO20M=1,NE
CALLESM(M,NE,NJ,JE,X,Y,E,A,KE)
DO20I=1,2
DO20II=1,2
IM=2*(I-1)+II
IW=2*(JE(I,M)-1)+II
DO20J=1,2
DO20JJ=1,2
JM=2*(J-1)+JJ
JW=2*(JE(J,M)-1)+JJ
JL=JW-IW+1
IF(JL.GT.0)THEN
KB(IW,JL)=KB(IW,JL)+KE(IM,JM)
ENDIF
20CONTINUE
RETURN
END
CSUBPROGRAM-2JOINTLOADVECTOR
SUBROUTINEFPJ(NEQ,NPJ,PJ,P)
DIMENSIONPJ(3,NPJ),P(NEQ)
DO10I=1,NEQ
10P(I)=0.0
DO20I=1,NPJ
J=PJ(1,I)
P(2*J-1)=PJ(2,I)
P(2*J)=PJ(3,I)
20CONTINUE
RETURN
END
CSUBPROGRAM-3INTRODUCESUPPORTCONDITION
SUBROUTINECSE(NEQ,NS,NW,JS,KB,P)
DIMENSIONJS(3,NS),P(NEQ)
REALKB(NEQ,NW)
DO10I=1,NS
J=JS(1,I)
DO20K=1,2
IF(JS(K+1,I).NE.0)THEN
L=2*(J-1)+K
KB(L,1)=1.0
DO30JJ=2,NW
30KB(L,JJ)=0.0
IF(L.GE.NW)JM=NW
IF(L.LT.NW)JM=L
DO40JJ=2,JM
II=L-JJ+1
40KB(II,JJ)=0.0
P(L)=0.0
ENDIF
20CONTINUE
10CONTINUE
RETURN
END
CSUBPROGRAM-4SOLVEEQUATIONS
SUBROUTINEBGS(NEQ,NJ,NW,KB,P)
DIMENSIONKB(NEQ,NW),P(NEQ)
REALKB
N1=NEQ-1
DO10K=1,N1
IM=K+NW-1
IF(NEQ.LT.IM)IM=NEQ
I1=K+1
DO20I=I1,IM
L=I-K+1
C=KB(K,L)/KB(K,1)
JM=NW-L+1
DO30J=1,JM
M=J+I-K
30KB(I,J)=KB(I,J)-C*KB(K,M)
20P(I)=P(I)-C*P(K)
10CONTINUE
P(NEQ)=P(NEQ)/KB(NEQ,1)
DO40K=1,N1
I=NEQ-K
JM=K+1
IF(NW.LT.JM)JM=NW
DO50J=2,JM
L=J+I-1
50P(I)=P(I)-KB(I,J)*P(L)
40P(I)=P(I)/KB(I,1)
WRITE(6,70)
70FORMAT(/5X,'NODALDISPLACEMENTS')
WRITE(6,60)(I,P(2*I-1),P(2*I),I=1,NJ)
60FORMAT(/7X,'NODE',10X,'U',14X,'V'/(1X,I10,2E15.6))
RETURN
END
C
CSUBPROGRAM-5INTERNALFORCES&SUPPORTFORCES
SUBROUTINETRS(NEQ,NE,NJ,NS,NPJ,JE,JS,X,Y,E,A,PJ,P,F)
DIMENSIONJE(2,NE),P(NEQ),F(NE),JS(3,NS),X(NJ),Y(NJ),E(NE)
&,A(NE),PJ(3,NPJ)
WRITE(6,10)
10FORMAT(/5X,'ELEMENTINTERNALFORCE'//7X,'ELEMENT',6X,'N',11X,
&'STRESS')
DO20M=1,NE
CALLSCL(M,NE,NJ,JE,X,Y,BL,SI,CO)
I=JE(1,M)
J=JE(2,M)
DU=P(2*J-1)-P(2*I-1)
DV=P(2*J)-P(2*I)
EAL=E(M)*A(M)/BL
F(M)=EAL*(DU*CO+DV*SI)
SG=F(M)/A(M)
WRITE(6,30)M,F(M),SG
30FORMAT(1X,I10,2E15.6)
20CONTINUE
WRITE(6,40)
40FORMAT(/5X,'SUPPORTFORCES'//7X,'NODE',7X,'HR',13X,'VR')
DO50II=1,NS
HR=0.0
VR=0.0
I=JS(1,II)
DO60J=1,NPJ
L=PJ(1,J)
IF(I.EQ.L)THEN
HR=HR-PJ(2,J)
VR=VR-PJ(3,J)
ENDIF
60CONTINUE
DO70M=1,NE
IF(I.EQ.JE(1,M).OR.I.EQ.JE(2,M))THEN
B=-1.0
IF(I.EQ.JE(2,M))B=1.0
CALLSCL(M,NE,NJ,JE,X,Y,BL,SI,CO)
HR=HR+B*F(M)*CO
VR=VR+B*F(M)*SI
ENDIF
70CONTINUE
WRITE(6,80)I,HR,VR
50CONTINUE
80FORMAT(1X,I10,2E15.6)
RETURN
END
c
CSUBPROGRAM-01ELEMENTSTIFFNESSMATRIX
SUBROUTINEESM(M,NE,NJ,JE,X,Y,E,A,KE)
DIMENSIONJE(2,NE),X(NJ),Y(NJ),E(NE),A(NE)
REALKE(4,4)
DO10I=1,4
DO10J=1,4
10KE(I,J)=0.0
CALLSCL(M,NE,NJ,JE,X,Y,BL,SI,CO)
EAL=E(M)*A(M)/BL
S1=EAL*CO*CO
S2=EAL*SI*CO
S3=EAL*SI*SI
KE(1,1)=S1
KE(1,2)=S2
KE(2,1)=S2
KE(2,2)=S3
DO20I=1,2
DO20J=1,2
KE(I,J+2)=-KE(I,J)
KE(I+2,J)=-KE(I,J)
20KE(I+2,J+2)=KE(I,J)
RETURN
END
CSUBPROGRAM-001ELEMENTCONSTANTS
SUBROUTINESCL(M,NE,NJ,JE,X,Y,BL,SI,CO)
DIMENSIONJE(2,NE),X(NJ),Y(NJ)
I=JE(1,M)
J=JE(2,M)
DX=X(J)-X(I)
DY=Y(J)-Y(I)
BL=SQRT(DX*DX+DY*DY)
SI=DY/BL
CO=DX/BL
END
C
C
二、空间桁架静力分析源程序——sTSAP.FOR
c
cspactrussstructuralanalysisprogram
c
DIMENSIONJE(2,100),JS(4,50),X(100),Y(100),Z(100),E(100),A(100),
&PJ(4,50)
REAL*8SKB(200,50),P(200),F(100)
OPEN(5,FILE='STR.DAT',STATUS='OLD')
OPEN(6,FILE='STW.DAT',STATUS='NEW')
Cinputcontrolparametersofstructure
READ(5,*)NJ,NE,NS,NPJ
CALLDAT(NJ,NE,NS,NPJ,NW,NEQ,X,Y,Z,JE,JS,PJ,A,E)
CALLSSM(NEQ,NE,NJ,NW,JE,X,Y,Z,E,A,SKB)
CALLFPJ(NEQ,NPJ,PJ,P)
CALLCSE(NEQ,NS,NW,JS,SKB,P)
CALLBGS(NEQ,NJ,NW,SKB,P)
CALLTRS(NEQ,NE,NJ,NS,NPJ,JE,JS,X,Y,Z,E,A,PJ,P,F)
CLOSE(5)
CLOSE(6)
END
c
cSUBROUTINE---DATINPUTandPRINT
c
SUBROUTINEDAT(NJ,NE,NS,NPJ,NW,NEQ,X,Y,Z,JE,JS,PJ,A,E)
DIMENSIONJE(2,ME),JS(4,NS),X(NJ),Y(NJ),Z(NJ),E(NE),A(NE),
&PJ(4,NPJ)
READ(5,*)(X(I),Y(I),Z(I),I=1,NJ)
READ(5,*)((JE(I,J),I=1,2),J=1,NE)
READ(5,*)((JS(I,J),I=1,4),J=1,NS)
READ(5,*)((PJ(I,J),I=1,4),J=1,NPJ)
READ(5,*)(A(I),E(I),I=1,NE)
NW=0
DO70I=1,NE
IW=IABS(JE(1,I)-JE(2,I))
IF(IW.GT.NW)NW=IW
70CONTINUE
NW=3*(NW+1)
NEQ=3*NJ
WRITE(6,75)
75FORMAT(1X,'SPACETRUSSSTRUCTUREANALYSIS'/)
Coutputcontrolparametersofstructure
WRITE(6,130)
WRITE(6,140)NE,NJ,NS,NPJ,NEQ,NW
Coutputdataofelementsandjoints
WRITE(6,150)
WRITE(6,78)
78FORMAT(/1X,'NODALCOODINATES')
WRITE(6,80)(I,X(I),Y(I),Z(I),I=1,NJ)
WRITE(6,90)(J,(JE(I,J),I=1,2),A(J),E(J),J=1,NE)
WRITE(6,100)((JS(I,J),I=1,4),J=1,NS)
WRITE(6,110)((PJ(I,J),I=1,4),J=1,NPJ)
80FORMAT(7X,'NODE',7X,'X',14X,'Y',14X,'Z'/(1X,I10,3E15.6))
90FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',7X,'A',14X,'E'
&/(1X,3I10,2E15.6))
100FORMAT(/1X,'CONSTRAINTEDNODES'/7X,'NODE',8X,'XX',8X,'YY',8X,'ZZ'
&/(1X,4I10))
110FORMAT(/1X,'JOINTLOADS'/7X,'NODE',8X,'PX',11X,'PY',13X,'PZ'
&/(1X,F10.0,3E15.6))
130FORMAT(1X,'STRUCTURALPARAMETERS'//1X,'TOTAL',1X,'NO.OF')
140FORMAT(1X,'ELEMENTS=',I5/1X,'JOINTS=',I5/1X,'CONSTRAINTEDJOINTS=
&',I5/1X,'JOINTLOADS=',I5/1X,'DEGREESOFFREEDOM=',I5/
&1X,'BANDWIDTH=',I5)
150FORMAT(/1X,'DATAOFELEMENTS&JOINTS')
RETURN
END
c
Csubprogram-1wholestiffnessmatrix
c
SUBROUTINESSM(NEQ,NE,NJ,NW,JE,X,Y,Z,E,A,SKB)
DIMENSIONJE(2,NE),X(NJ),Y(NJ),Z(NJ),E(NE),A(NE)
REAL*8SKB(NEQ,NW),SKE(6,6)
DO10I=1,NEQ
DO10J=1,NW
SKB(I,J)=0.0
10CONTINUE
DO20M=1,NE
CALLESM(M,NE,NJ,JE,X,Y,Z,E,A,SKE)
DO20I=1,2
DO20II=1,3
IM=3*(I-1)+II
IW=3*(JE(I,M)-1)+II
DO20J=1,2
DO20JJ=1,3
JM=3*(J-1)+JJ
JW=3*(JE(J,M)-1)+JJ
JL=JW-IW+1
IF(JL.GT.0)THEN
SKB(IW,JL)=SKB(IW,JL)+SKE(IM,JM)
ENDIF
20CONTINUE
RETURN
END
c
Csubprogram-01elementstiffnessmatrix
c
SUBROUTINEESM(M,NE,NJ,JE,X,Y,Z,E,A,SKE)
DIMENSIONJE(2,NE),X(NJ),Y(NJ),Z(NJ),E(NE),A(NE)
REAL*8SKE(6,6),S1,S2,S3,S4,S5,S6,BL,CA,CB,CR,EAL
DO10I=1,6
DO10J=1,6
10SKE(I,J)=0.0
CALLSCL(M,NE,NJ,JE,X,Y,Z,BL,CA,CB,CR)
EAL=E(M)*A(M)/BL
S1=EAL*CA*CA
S2=EAL*CA*CB
S3=EAL*CA*CR
S4=EAL*CB*CB
S5=EAL*CB*CR
S6=EAL*CR*CR
SKE(1,1)=S1
SKE(1,2)=S2
SKE(1,3)=S3
SKE(2,1)=S2
SKE(2,2)=S4
SKE(2,3)=S5
SKE(3,1)=S3
SKE(3,2)=S5
SKE(3,3)=S6
DO20I=1,3
DO20J=1,3
SKE(I,J+3)=-SKE(I,J)
SKE(I+3,J)=-SKE(I,J)
20SKE(I+3,J+3)=SKE(I,J)
RETURN
END
c
Csubprogram-001elementconstants
c
SUBROUTINESCL(M,NE,NJ,JE,X,Y,Z,BL,CA,CB,CR)
DIMENSIONJE(2,NE),X(NJ),Y(NJ),Z(NJ)
REAL*8DX,DY,DZ,BL,CA,CB,CR
I=JE(1,M)
J=JE(2,M)
DX=X(J)-X(I)
DY=Y(J)-Y(I)
DZ=Z(J)-Z(I)
BL=SQRT(DX*DX+DY*DY+DZ*DZ)
CA=DX/BL
CB=DY/BL
CR=DZ/BL
END
Csubprogram-2jointloadvector
SUBROUTINEFPJ(NEQ,NPJ,PJ,P)
DIMENSIONPJ(4,NPJ)
REAL*8P(NEQ)
DO10I=1,NEQ
10P(I)=0.0
DO20I=1,NPJ
J=PJ(1,I)
P(3*J-2)=PJ(2,I)
P(3*J-1)=PJ(3,I)
P(3*J)=PJ(4,I)
20CONTINUE
RETURN
END
C
Csubprogram-3introducesupportcondition
c
SUBROUTINECSE(NEQ,NS,NW,JS,SKB,P)
DIMENSIONJS(4,NS)
REAL*8SKB(NEQ,NW),P(NEQ)
DO10I=1,NS
J=JS(1,I)
DO20K=1,3
IF(JS(K+1,I).NE.0)THEN
L=3*(J-1)+K
SKB(L,1)=1.0
DO30JJ=2,NW
30SKB(L,JJ)=0.0
IF(L.GE.NW)JM=NW
IF(L.LT.NW)JM=L
DO40JJ=2,JM
II=L-JJ+1
40SKB(II,JJ)=0.0
P(L)=0.0
ENDIF
20CONTINUE
10CONTINUE
RETURN
END
c
Csubprogram-4solveequations
c
SUBROUTINEBGS(NEQ,NJ,NW,SKB,P)
REAL*8SKB(NEQ,NW),P(NEQ),C
N1=NEQ-1
DO10K=1,N1
IM=K+NW-1
IF(NEQ.LT.IM)IM=NEQ
I1=K+1
DO20I=I1,IM
L=I-K+1
C=SKB(K,L)/SKB(K,1)
JM=NW-L+1
DO30J=1,JM
M=J