平面刚架的计算源程序（wilson法）
下面给出求平面刚架的计算源程序，程序中的变量及数组说明如下：NN――结点总数
NE――单元总数
NC――支撑约束数
M1――要计算的振型数
NM――集中质量数
N3――结点位移总数，N3＝3＊NN
N――自由位移数，N＝N3－NF
X――一维数组，存放结点的X坐标值
Y――一维数组，存放结点的Y坐标值
NS――一维数组，存放位移约束信息
E――材料弹性模量
RH――材料密度
JOD――有集中质量的结点号
MAS1――集中质量值
MAS2――集中转动惯量值
I9――一维数组，存放单元始端结点号
J9――一维数组，存放单元终端结点号
S9――一维数组，存放单元截面惯性矩
A9――一维数组，存放单元截面面积
A――二维数组，存放整体刚度矩阵
SSM――二维数组，存放单元质量矩阵
XX――二维数组，存放整体质量矩阵
AM――一维数组，存放频率值
VC――二维数组，存放特征向量
DC――二维数组，存放坐标转换矩阵
ST――二维数组，存放单元刚度矩阵
ND――体系自由度数
THETA――θ值
DT――积分时间步长Δt
TMAX――最大积分时间
NEQ（L）――定义动力荷载的关键点数目
TC（I）――动力作用时刻
P（I）――对应于TC（I）时刻的动力值
UD（I）――位移
VA（I）――加速度
程序名：DFRAM.FOR
c 平面刚架的振动
IMPLICIT INTEGER*2 (I-N)
C CHARACTER CH(80)*1,DFALE*12,OFALE*12,YN*1
INTEGER*2 LI(1000)
REAL*4 A(10000)
DIMENSION RH(10000)
OPEN(11,FILE='DFRAM.OUT')
OPEN(7,FILE='DFRAM.IN')
READ(7,*)NN,MS,NF,M1,NC
WRITE(11,'(/2X,A//2X,A)')'输出数据如下：','一.输入数据部分：'
WRITE(11,'(/2X,5A12/I10,4I12)')'结点数','杆件数',
1 '支撑约束数','计算振型数','集中质量数',NN,MS,NF,M1,NC
NJ=NN
N3=NN*3
NP=N3+1
N=N3-NF
ia=1
ixx=ia+n3*n3
iam=ixx+n3*np
ivc=iam+m1
ivtl=ivc+n*m1
ix=ivtl+n3
ive=ix+n3
iz=ive+n
iy=iz+n
i1ma=iy+nn
i2ma=i1ma+nc
ia9=i2ma+nc
is9=ia9+ms
mal=is9+ms
jns=1
jig=jns+nf
jjod=jig+n
ji9=jjod+nc
jj9=ji9+ms
nal=jj9+ms
if (mal.gt.10000-80.or.nal.gt.1000-80)
1 stop '数组越界'
READ (7,*)(A(IX+I),A(IY+I),I=0,NN-1)
READ (7,*)(LI(I),I=1,NF),E
WRITE (11,'(/2X,A/A8,2A16)')'结点坐标：','结点号',
1'x方向坐标','y方向坐标'
WRITE (11,'(/2X,A,I3,A,F14.3,F16.3)')
1 ('(',I+1,')',A(IX+I),A(IY+I),I=0,NN-1)
WRITE (11,'(1X,2A)')'约束自由度号(X方向位移号＝3X结点号－2；',
1 'y方向位移号＝3X结点号－1;转角位移号＝3X结点号)'
WRITE (11,'(I4,9I6)')(li(i),i=1,NF)
WRITE(11,'(/2X,A/A10)')'材料常数：','弹性模量'
WRITE (11,'(E10.4)')E
C
READ(7,*)(LI(JI9+I),LI(JJ9+I),A(IS9+I),A(IA9+I),RH(I+1),I=0,MS-1)
IF (NC.GT.0)READ (7,*)(LI(JJOD+I),A(I1MA+I),A(I2MA+I),I=0,NC-1)
WRITE(11,'(/2X,A/A8,6A12)')'各杆件单元信息：',
1 '单元号','左端结点号','右端结点号','惯性','截面积','密度'
WRITE(11,'(2X,A,I3,A,I9,I12,E15.5,E12.5,E16.5)')('(',I+1,')',
1 LI(JI9+I),LI(JJ9+I),A(IS9+I),A(IA9+I),RH(I+1),I=0,MS-1)
IF (NC.GT.0)THEN
WRITE(11,'(/2X,A)')'集中质量和转动惯量：'
WRITE(11,'(A16,2A12)')'质量所在结点号','质量','转动惯量'
WRITE(11,'(I10,F17.0,F12.3)')
1 (LI(JJOD+I),A(I1MA+I),A(I2MA+I),I=0,NC-1)
ENDIF
C
CALL DFRAM(A(IA),A(IXX),A(IAM),A(IVC),A(IVTL),A(IX),A(IVE),
1A(IZ),A(IY),A(I1MA),A(I2MA),A(IA9),A(IS9),E,RH,LI(JNS),
2LI(JIG),LI(JJOD),LI(JI9),LI(JJ9),NN,MS,NF,M1,NC,NJ,N3,NP,N)
CLOSE(7)
CLOSE(11)
STOP
END
C
SUBROUTINE DFRAM(A,XX,AM,VC,VTL,X,VE,Z,Y,MAS1,MAS2,A9,S9,
1E,RH,NS,IG,JOD,I9,J9,NN,MS,NF,M1,NC,NJ,N3,NP,N)
IMPLICIT INTEGER*2 (I-N)
INTEGER*2 NS(NF),IG(N)
INTEGER*2 JOD(NC),I9(MS),J9(MS),PO
REAL*4 A(N3,N3),XX(N3,NP),AM(M1),VC(N,M1),VTL(N3),X(N3),
1 DC(6,6),SSM(6,6),S2(6,6),VE(N),Z(N),Y(NN),SC(N3,NP)
REAL*4 MAS1(NC),MAS2(NC),A9(MS),S9(MS),L
N2=N3
C
DO 3530 ME=1,MS
I=I9(ME)
J=J9(ME)
SA=S9(ME)
CA=A9(ME)
L=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)
C
C=(X(J)-X(I))/L
S=(Y(J)-Y(I))/L
DO 2370 II=1,6
DO 2370 JJ=1,6
DC(II,JJ)=0.0
2370 SSM(II,JJ)=0.0
DC(1,1)=C
DC(2,2)=C
DC(4,4)=C
DC(5,5)=C
DC(1,2)=S
DC(4,5)=S
DC(2,1)=-S
DC(5,4)=-S
DC(3,3)=1.0
DC(6,6)=1.0
C
I3=3*I
I2=I3-1
I1=I3-2
J3=J*3
J2=J3-1
J1=J3-2
C1=12.0*E*SA*S*S/L**3+C*C*CA*E/L
C2=12.0*E*SA*C*S/L**3-C*S*CA*E/L
C3=12.0*E*SA*C*C/L**3+S*S*CA*E/L
C4=6.0*E*SA*S/L**2
C5=6.0*E*SA*C/L**2
C6=4.0*E*SA/L
A(I1,I1)=A(I1,I1)+C1
A(I2,I1)=A(I2,I1)-C2
A(I1,I2)=A(I1,I2)-C2
A(I3,I1)=A(I3,I1)+C4
A(I1,I3)=A(I1,I3)+C4
A(J1,I1)=A(J1,I1)-C1
A(I1,J1)=A(I1,J1)-C1
A(J2,I1)=A(J2,I1)+C2
A(I1,J2)=A(I1,J2)+C2
A(J3,I1)=A(J3,I1)+C4
A(I1,J3)=A(I1,J3)+C4
A(I2,I2)=A(I2,I2)+C3
A(I3,I2)=A(I3,I2)-C5
A(I2,I3)=A(I2,I3)-C5
A(J1,I2)=A(J1,I2)+C2
A(I2,J1)=A(I2,J1)+C2
A(J2,I2)=A(J2,I2)-C3
A(I2,J2)=A(I2,J2)-C3
A(J3,I2)=A(J3,I2)-C5
A(I2,J3)=A(I2,J3)-C5
A(I3,I3)=A(I3,I3)+C6
A(J1,I3)=A(J1,I3)-C4
A(I3,J1)=A(I3,J1)-C4
A(J2,I3)=A(J2,I3)+C5
A(I3,J2)=A(I3,J2)+C5
A(J3,I3)=A(J3,I3)+.5*C6
A(I3,J3)=A(I3,J3)+.5*C6
A(J1,J1)=A(J1,J1)+C1
A(J2,J1)=A(J2,J1)-C2
A(J1,J2)=A(J1,J2)-C2
A(J3,J1)=A(J3,J1)-C4
A(J1,J3)=A(J1,J3)-C4
A(J2,J2)=A(J2,J2)+C3
A(J3,J2)=A(J3,J2)+C5
A(J2,J3)=A(J2,J3)+C5
A(J3,J3)=A(J3,J3)+C6
C
C1=RH*CA*L/6
C2=RH*CA*L/420
SSM(1,1)=2.0*C1
SSM(4,1)=C1
SSM(2,2)=156.0*C2
SSM(3,2)=-22.0*L*C2
SSM(5,2)=54.0*C2
SSM(6,2)=13.0*L*C2
SSM(2,3)=-22.0*L*C2
SSM(3,3)=4.0*L*L*C2
SSM(5,3)=-13.0*L*C2
SSM(6,3)=-3.0*L**2*C2
SSM(1,4)=C1
SSM(4,4)=2.0*C1
SSM(2,5)=54.0*C2
SSM(3,5)=-13.0*L*C2
SSM(5,5)=156.0*C2
SSM(6,5)=22.0*L*C2
SSM(2,6)=13.0*L*C2
SSM(3,6)=-3.0*L**2*C2
SSM(5,6)=22.0*L*C2
SSM(6,6)=4.0*L**2*C2
C
DO 3190 II=1,6
DO 3190 JJ=1,6
3190 S2(II,JJ)=0.0
DO 3230 II=1,6
DO 3230 JJ=1,6
DO 3230 KK=1,6
3230S2(II,JJ)=S2(II,JJ)
1 +SSM(II,KK)*DC(KK,JJ)
DO 3290 II=1,6
DO 3290 JJ=1,6
3290SSM(II,JJ)=0.0
DO 3350 II=1,6
DO 3350 JJ=1,6
DO 3350 KK=1,6
3350 SSM(II,JJ)=SSM(II,JJ)+DC(KK,II)*S2(KK,JJ)
C
I1=3*I-3
J1=3*J-3
DO 3520 II=1,3
DO 3520 JJ=1,3
MG=I1+II
ITR=I1+JJ
IAC=J1+II
IXK=J1+JJ
XX(MG,ITR)=XX(MG,ITR)+SSM(II,JJ)
XX(IAC,ITR)=XX(IAC,ITR)+SSM(II+3,JJ)
XX(MG,IXK)=XX(MG,IXK)+SSM(II,JJ+3)
XX(IAC,IXK)=XX(IAC,IXK)+SSM(II+3,JJ+3)
3520 CONTINUE
3530CONTINUE
C
IF(NC.EQ.0) GOTO 3660
DO 3650 II=1,NC
PO=JOD(II)
GMS=MAS1(II)
MI=MAS2(II)
K1=3*PO
J1=K1-1
I1=K1-2
XX(I1,I1)=XX(I1,I1)+GMS
XX(J1,J1)=XX(J1,J1)+GMS
XX(K1,K1)=XX(K1,K1)+MI
3650CONTINUE
C
3660MM=0
DO 3880 I=1,NF
NA=NS(I)
MM=MM+1
NA=NA-MM+1
NB=N3-MM
IF (NB.LT.NA) GOTO 3790
DO 3760 II=NA,NB
DO 3760 JJ=NA,NB
A(II,JJ)=A(II+1,JJ+1)
3760XX(II,JJ)=XX(II+1,JJ+1)
3790IF (NA.EQ.1) GOTO 3880
DO 3850 II=1,NA-1
DO 3850 JJ=NA,NB
A(II,JJ)=A(II,JJ+1)
XX(II,JJ)=XX(II,JJ+1)
A(JJ,II)=A(JJ+1,II)
3850 XX(JJ,II)=A(JJ+1,II)
3880CONTINUE
C
1640 N9=N-1
DO 1830 IX=1,N
DI=A(IX,1)
DII=0.0
IF (ABS(DI).GT.1.0E-16)DII=1.0/DI
IF (ABS(DI).LT.1.0E-16)WRITE(11,'(A)')
1 'MATRIX IS SINGULAR'
DO 1700 IY=1,N9
IY9=IY+1
1700 A(IX,IY)=A(IX,IY9)/DI
A(IX,N)=1.0/DI
DO 1820 IZ=1,N
IF (IZ.EQ.IX)GO TO 1820
O=A(IZ,1)
DO 1790 IY=1,N9
IY9=IY+1
1790 A(IZ,IY)=A(IZ,IY9)-A(IX,IY)*O
A(IZ,N)=-A(IX,N)*O
1820CONTINUE
1830 CONTINUE
C
3900 DO 4000 I=1,N
DO 3960 J=1,N
VE(J)=0.0
DO 3940 II=1,N
3940 VE(J)=VE(J)+A(I,II)*XX(II,J)
3960 CONTINUE
DO 3980 J=1,N
3980 A(I,J)= VE(J)
4000 CONTINUE
C
500 MN=N
NN=N
CALL EIG1(A,X,Z,NN,N2,0.0001,XM,NKK)
530 M=1
DO 560 I=1,NN
VC(I,M)=X(I)
560 XX(I,M)=VC(I,M)
AM(M)=XM
C
IF (M1.LT.2) GOTO 220
DO 1410 M=2,M1
DO 630 I=1,NN
GK4=ABS(XX(I,M-1)-1.)
IF (GK4.LT.0.001)IR=I
630 CONTINUE
IG(M-1)=IR
DO 680 I=1,NN
680 XX(MN-I+1,MN-M+3)=A(IR,I)
DO 730 I=1,NN
DO 730 J=1,NN
J1=MN-J+1
J2=MN-M+3
730 A(I,J)=A(I,J)-XX(I,M-1)*XX(J1,J2)
DO 910 I =1,NN
IF (I.EQ.IR)GOTO 910
IF (I.GT.IR)K1=I-1
IF (I.LE.IR)K1=I
DO 900 J=1,NN
IF (J.EQ.IR)GOTO 900
IF (J.GT.IR)K2=J-1
IF (J.LE.IR)K2=J
A(K1,K2)=A(I,J)
900 CONTINUE
910 CONTINUE
NN=NN-1
M3=NN
IF (M.NE.MN) GOTO 980
XM=A(1,1)
X(1)=1
970 GOTO 990
980 L1500=2
CALL EIG1(A,X,Z,NN,N2,0.0001,XM,NKK)
990 DO 1000 I=1,NN
1000 XX(I,M)=X(I)
AM(M)=XM
M4=M-1
M5=1000-M4
DO 1400 M8=M5,999
M6=M3+1
M2=1000-M8
M7=IG(M2)+1
IF(M6.LT.M7)GOTO 1150
N9=1000-M7
N8=1000-M6
DO 1130 I3=N8,N9
II=1000-I3
1130 X(II)=X(II-1)
1150 JJ=IG(M2)
X(JJ)=0.0
GSM=0.0
DO 1210 I=1,M6
J3=MN-I+1
J4=MN-M2+2
1210 GSM=GSM+XX(J3,J4)*X(I)
XK=(AM(M2)-XM)/GSM
DO 1250 I=1,M6
1250 X(I)=XX(I,M2)-XK*X(I)
GSM=0.0
DO 1300 I=1,M6
IF (ABS(GSM).LT.ABS(X(I)))GSM=X(I)
1300 CONTINUE
DO 1330 I=1,M6
1330 X(I)=X(I)/GSM
M3=M3+1
IF (M2.NE.1)GOTO 1400
DO 1380 I=1,M3
1380 VC(I,M)=X(I)
1400 CONTINUE
1410 CONTINUE
C
220 WRITE(11,'(/2X,A)')'计算结果部分：'
WRITE(11,'(/2X,A)')'输出各振型频率和振型'
DO 360 I=1,M1
TT=0.0
IF (AM(I).GT.1E-8)TT=SQRT(1.0/AM(I))/(2.0*3.14159)
WRITE(11,'(/2(2X,A,I2,A,F9.4))')'第',I,'振型频率：',TT,'振型：'
J=1
J0=1
DO 330 IT=1,N3
IF (J0.GT.NF)GOTO 320
DO 300 JS=J0,NF
IF(IT.EQ.NS(JS))GOTO 325
300 CONTINUE
320 VTL(IT)=VC(J,I)
J=J+1
GO TO 330
325 VTL(IT)=0.0
J0=J0+1
330 CONTINUE
WRITE(11,'(A8,A18,2A18/(I6,2E18.6,E20.6))')
1 '结点号','水平方向位移','垂直方向位移','转角',
2 (J,VTL(3*J-2),VTL(3*J-1),VTL(3*J),J=1,NJ)
360 CONTINUE
C
ND=(N3-NF)/3
DO 1500 ME=1,MS
I=I9(ME)
J=J9(ME)
IF(Y(I).EQ.Y(J))ND=ND-1
1500 CONTINUE
WRITE(11,*)'ND=',ND
DO 1510 I=1,ND
DO 1510 J=1,ND
AX=0.0
BX=0.0
SC(I,J)=AX*A(I,J)+BX*XX(I,J)
1510 CONTINUE
CALL STEPM(A,XX,SC,ND,N3)
RETURN
END
C SUBROUTINE EIG1(A,X,Z,NN,N2,D,XM,NKK)
IMPLICIT INTEGER*2 (I-N)
DIMENSION A(N2,N2),X(NN),Z(NN)
Y1=100000.0
NKK=1
DO 1320 I=1,NN
1320 X(I)=1.0
XM=-100000.0
1340 DO 1400 I=1,NN
SG=0.0
DO 1380 J=1,NN
1380 SG=SG+A(I,J)*X(J)
1400 Z(I)=SG
XM=0.0
DO 1450 I=1,NN
IF (ABS(XM).LT.ABS(Z(I)))XM=Z(I)
1450 CONTINUE
DO 1500 I=1,NN
1500 X(I)=Z(I)/XM
IF (ABS(Y1-XM)/XM.GT.D)GOTO 1530
GOTO 1550
1530 Y1=XM
NKK=NKK+1
IF (NKK.GT.200)GOTO 1623
GOTO 1340
1550 X3=0.0
DO 1580 I=1,NN
IF (ABS(X3).LT.ABS(X(I)))X3=X(I)
1580 CONTINUE
DO 1610 I=1,NN
1610 X(I)=X(I)/X3
RETURN
1623 WRITE(*)'NUMBER OF ITERATION NKK=',NKK
1630 FORMAT(A30,I2)
STOP
END
C
SUBROUTINE STEPM(SK,SM,SC,ND,N3)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION SK(N3,N3),SM(N3,N3),SC(N3,N3),F(N3,N3),X(N3,N3),
1 DUA(N3),UD(N3),UV(N3),UA(N3),TC(N3),P(N3),NEQ(N3)
ND1=ND+1
DO 10 I=1,N3
UD(I)=0.0
UV(I)=0.0
DO 10 J=1,N3
10 F(I,J)=0.0
WRITE(11,'(/1X,A)')'TIME PARA FOR ITEGRATION'
READ(7,*)THETA,DT,TMAX
WRITE(11,'(4X,3A6/4X,3F7.4)')
1 'THETA','DT','TMAX',THETA,DT,TMAX
WRITE(11,'(/1X,A)')'NO.OF LOAD DATA POINTS'
READ(7,*)(NEQ(L),L=1,ND)
WRITE(11,'(4X,A6,I2,A3,I5)')('NEQ(',L,')=',NEQ(L),L=1,ND)
ANN=0.0
II=1
DO 60 ID=1,ND
NE=NEQ(ID)
IF (NE.EQ.0) GOTO 60
IF(NE.GT.TMAX/DT)NE=TMAX/DT+0.5
READ(7,*)(TC(J),P(J),J=1,NE)
WRITE(11,'(/1X,A)')'DYNAMIC LOAD ON STRUCTURE'
WRITE(11,'(4X,A12,A15)')'TC(I)','P(I)'
WRITE(11,'(4X,F12.5,E15.6)')(TC(J),P(J),J=1,NE)
NT=TC(NE)/DT+0.5
NT1=NT+1
NT2=NT+2
F(ID,1)=P(1)
ANN=0.0
II=1
DO 50 I=2,NT2
AI=I-1
T=AI*DT
IF(T.GT.TC(NE))GOTO 60
IF(T.LE.TC(II+1))GOTO 40
ANN=-TC(II+1)+T-DT
II=II+1
40 ANN=ANN+DT
F(ID,I)=P(II)+(P(II+1)-P(II))*ANN/(TC(II+1)-TC(II))
50 CONTINUE
60 CONTINUE
NT=TMAX/DT+0.5
NT1=NT+1
DO 70 I=1,ND
X(I,ND1)=F(I,J)
DO 70 J=1,ND
70 X(I,J)=SM(I,J)
DO 80 I=1,ND
DO 80 J=1,ND
80 X(I,ND1)=X(I,ND1)-SC(I,J)*UV(J)-SK(I,J)*UD(J)
CALL SOLVE(ND,X,N3)
DO 90 I=1,ND
90 UA(I)=X(I,ND1)
WRITE(11,'(/1X,A)')'RESPONSE OF STRUCTURE'
WRITE(11,'(/3X,2A5,3A11)')'COR. ','TIME','DISPL.','VEL.','ACC.'
TU=THETA*DT
A1=3.0/TU
A2=6.0/TU
A3=TU/2.0
A4=A2/TU
DO 190 L=1,NT1
DO 110 I=1,ND
DO 110 J=1,ND
110 X(I,J)=SK(I,J)+A4*SM(I,J)+A1*SC(I,J)
AL=L
T=AL*DT
DO 130 I=1,ND
X(I,ND1)=(F(I,L+1)+F(I,L+2)-F(I,L+1))*(THETA-1.0)-F(I,L)
DO 120 J=1,ND
120 X(I,ND1)=X(I,ND1)+(SM(I,J)*A2+SC(I,J)*3.0)*UV(J)+(SM(I,J)*3.0
2 +A3*SC(I,J))*UA(J)
130 CONTINUE
CALL SOLVE(ND,X,N3)
DO 140 I=1,ND
DUA(I)=A4*X(I,ND1)-A2*UV(I)-3.0*UA(I)
DUA(I)=DUA(I)/THETA
DUV=DT*UA(I)+DT*DUA(I)/2.0
UD(I)=UD(I)+DT*UV(I)+DT*DT*UA(I)/2.0+DT*DT*DUA(I)/6.0
UV(I)=UV(I)+DUV
140 CONTINUE
DO 160 I=1,ND
X(I,ND1)=F(I,L+1)
DO 150 J=1,ND
X(I,ND1)=X(I,ND1)-SK(I,J)*UD(J)-SC(I,J)*UV(J)
150 X(I,J)=SM(I,J)
160 CONTINUE
CALL SOLVE(ND,X,N3)
DO 170 I=1,ND
UA(I)=X(I,ND1)
170 WRITE(11,'(1X,I3,1X,F8.6,3E12.4)')I,T,UD(I),
1 UV(I),UA(I)
190 CONTINUE
RETURN
END
C
SUBROUTINE SOLVE(N,A,N3)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(N3,N3+1)
M=1
EPS=1.0E-10
NPLUSM=N+M
DET=1.0
DO 40 K=1,N
DET=DET*A(K,K)
IF(DABS(A(K,K)).GT.EPS)GOTO 10
GOTO 60
10 KP1=K+1
DO 20 J=KP1,NPLUSM
20 A(K,J)=A(K,J)/A(K,K)
A(K,K)=1.0
DO 40 I=1,N
IF (I.EQ.K.OR.A(I,K).EQ.0)GOTO 40
DO 30 J=KP1,NPLUSM
30 A(I,J)=A(I,J)-A(I,K)*A(K,J)
A(I,K)=0.0
40 CONTINUE
60 RETURN
END
页:
[1]