一平面桁架静力分析源程序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
升级会员 