平面刚架的计算源程序（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

1. c　　平面刚架的振动
   
2. IMPLICIT INTEGER*2 (I-N)
   
3. C   CHARACTER CH(80)*1,DFALE*12,OFALE*12,YN*1
   
4.         INTEGER*2 LI(1000)
   
5.         REAL*4 A(10000)
   
6.   DIMENSION RH(10000)
   
7.    OPEN(11,FILE='DFRAM.OUT')
   
8.     OPEN(7,FILE='DFRAM.IN')
   
9.   READ(7,*)NN,MS,NF,M1,NC
   
10.   WRITE(11,'(/2X,A//2X,A)')'输出数据如下：','一.输入数据部分：'
    
11.         WRITE(11,'(/2X,5A12/I10,4I12)')'结点数','杆件数',
    
12.      1   '支撑约束数','计算振型数','集中质量数',NN,MS,NF,M1,NC
    
13.   NJ=NN
    
14.   N3=NN*3
    
15.   NP=N3+1
    
16.   N=N3-NF
    
17.   ia=1
    
18.   ixx=ia+n3*n3
    
19.   iam=ixx+n3*np
    
20.   ivc=iam+m1
    
21.   ivtl=ivc+n*m1
    
22.   ix=ivtl+n3
    
23.   ive=ix+n3
    
24.   iz=ive+n
    
25.   iy=iz+n
    
26.   i1ma=iy+nn
    
27.   i2ma=i1ma+nc
    
28.   ia9=i2ma+nc
    
29.   is9=ia9+ms
    
30.   mal=is9+ms
    
31.    jns=1
    
32.    jig=jns+nf
    
33.    jjod=jig+n
    
34.    ji9=jjod+nc
    
35.    jj9=ji9+ms
    
36.    nal=jj9+ms
    
37.   if (mal.gt.10000-80.or.nal.gt.1000-80)
    
38.      1 stop '数组越界'
    
39.   READ (7,*)(A(IX+I),A(IY+I),I=0,NN-1)
    
40.   READ (7,*)(LI(I),I=1,NF),E
    
41. WRITE (11,'(/2X,A/A8,2A16)')'结点坐标：','结点号',
    
42.      1  'x方向坐标','y方向坐标'
    
43. WRITE (11,'(/2X,A,I3,A,F14.3,F16.3)')
    
44.      1   ('(',I+1,')',A(IX+I),A(IY+I),I=0,NN-1)
    
45. WRITE (11,'(1X,2A)')'约束自由度号(X方向位移号＝3X结点号－2；',
    
46.      1    'y方向位移号＝3X结点号－1;转角位移号＝3X结点号)'
    
47. WRITE (11,'(I4,9I6)')(li(i),i=1,NF)
    
48. WRITE(11,'(/2X,A/A10)')'材料常数：','弹性模量'
    
49. WRITE (11,'(E10.4)')E
    
50. C
    
51.      READ(7,*)(LI(JI9+I),LI(JJ9+I),A(IS9+I),A(IA9+I),RH(I+1),I=0,MS-1)
    
52. IF (NC.GT.0)READ (7,*)(LI(JJOD+I),A(I1MA+I),A(I2MA+I),I=0,NC-1)
    
53. WRITE(11,'(/2X,A/A8,6A12)')'各杆件单元信息：',
    
54.      1 '单元号','左端结点号','右端结点号','惯性','截面积','密度'
    
55. WRITE(11,'(2X,A,I3,A,I9,I12,E15.5,E12.5,E16.5)')('(',I+1,')',
    
56.      1    LI(JI9+I),LI(JJ9+I),A(IS9+I),A(IA9+I),RH(I+1),I=0,MS-1)
    
57. IF (NC.GT.0)THEN
    
58.     WRITE(11,'(/2X,A)')'集中质量和转动惯量：'
    
59.           WRITE(11,'(A16,2A12)')'质量所在结点号','质量','转动惯量'
    
60.     WRITE(11,'(I10,F17.0,F12.3)')
    
61.      1    (LI(JJOD+I),A(I1MA+I),A(I2MA+I),I=0,NC-1)
    
62.   ENDIF
    
63. C
    
64. 
    
65.       CALL DFRAM(A(IA),A(IXX),A(IAM),A(IVC),A(IVTL),A(IX),A(IVE),
    
66.      1  A(IZ),A(IY),A(I1MA),A(I2MA),A(IA9),A(IS9),E,RH,LI(JNS),
    
67. 2  LI(JIG),LI(JJOD),LI(JI9),LI(JJ9),NN,MS,NF,M1,NC,NJ,N3,NP,N)
    
68. CLOSE(7)
    
69. CLOSE(11)
    
70. STOP
    
71. END
    
72. C
    
73.       SUBROUTINE DFRAM(A,XX,AM,VC,VTL,X,VE,Z,Y,MAS1,MAS2,A9,S9,
    
74.      1  E,RH,NS,IG,JOD,I9,J9,NN,MS,NF,M1,NC,NJ,N3,NP,N)
    
75. IMPLICIT INTEGER*2 (I-N)
    
76. INTEGER*2 NS(NF),IG(N)
    
77. INTEGER*2 JOD(NC),I9(MS),J9(MS),PO
    
78. REAL*4 A(N3,N3),XX(N3,NP),AM(M1),VC(N,M1),VTL(N3),X(N3),
    
79.      1 DC(6,6),SSM(6,6),S2(6,6),VE(N),Z(N),Y(NN),SC(N3,NP)
    
80. REAL*4 MAS1(NC),MAS2(NC),A9(MS),S9(MS),L
    
81. N2=N3
    
82. 
    
83. C
    
84.       DO 3530 ME=1,MS
    
85.    I=I9(ME)
    
86.    J=J9(ME)
    
87.    SA=S9(ME)
    
88.    CA=A9(ME)
    
89. L=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)
    
90. 
    
91. C
    
92.          C=(X(J)-X(I))/L
    
93.    S=(Y(J)-Y(I))/L
    
94.       DO 2370 II=1,6
    
95.          DO 2370 JJ=1,6
    
96.                 DC(II,JJ)=0.0
    
97. 2370            SSM(II,JJ)=0.0
    
98.            DC(1,1)=C
    
99.      DC(2,2)=C
    
100. DC(4,4)=C
     
101. DC(5,5)=C
     
102. DC(1,2)=S
     
103. DC(4,5)=S
     
104. DC(2,1)=-S
     
105. DC(5,4)=-S
     
106. DC(3,3)=1.0
     
107. DC(6,6)=1.0
     
108. C
     
109.       I3=3*I
     
110. I2=I3-1
     
111. I1=I3-2
     
112. J3=J*3
     
113. J2=J3-1
     
114. J1=J3-2
     
115. C1=12.0*E*SA*S*S/L**3+C*C*CA*E/L
     
116. C2=12.0*E*SA*C*S/L**3-C*S*CA*E/L
     
117. C3=12.0*E*SA*C*C/L**3+S*S*CA*E/L
     
118. C4=6.0*E*SA*S/L**2
     
119. C5=6.0*E*SA*C/L**2
     
120. C6=4.0*E*SA/L
     
121. A(I1,I1)=A(I1,I1)+C1
     
122. A(I2,I1)=A(I2,I1)-C2
     
123. A(I1,I2)=A(I1,I2)-C2
     
124. A(I3,I1)=A(I3,I1)+C4
     
125. A(I1,I3)=A(I1,I3)+C4
     
126. A(J1,I1)=A(J1,I1)-C1
     
127. A(I1,J1)=A(I1,J1)-C1
     
128. A(J2,I1)=A(J2,I1)+C2
     
129. A(I1,J2)=A(I1,J2)+C2
     
130. A(J3,I1)=A(J3,I1)+C4
     
131. A(I1,J3)=A(I1,J3)+C4
     
132. A(I2,I2)=A(I2,I2)+C3
     
133. A(I3,I2)=A(I3,I2)-C5
     
134. A(I2,I3)=A(I2,I3)-C5
     
135. A(J1,I2)=A(J1,I2)+C2
     
136. A(I2,J1)=A(I2,J1)+C2
     
137. A(J2,I2)=A(J2,I2)-C3
     
138. A(I2,J2)=A(I2,J2)-C3
     
139. A(J3,I2)=A(J3,I2)-C5
     
140. A(I2,J3)=A(I2,J3)-C5
     
141. A(I3,I3)=A(I3,I3)+C6
     
142. A(J1,I3)=A(J1,I3)-C4
     
143. A(I3,J1)=A(I3,J1)-C4
     
144. A(J2,I3)=A(J2,I3)+C5
     
145. A(I3,J2)=A(I3,J2)+C5
     
146. A(J3,I3)=A(J3,I3)+.5*C6
     
147. A(I3,J3)=A(I3,J3)+.5*C6
     
148. A(J1,J1)=A(J1,J1)+C1
     
149. A(J2,J1)=A(J2,J1)-C2
     
150. A(J1,J2)=A(J1,J2)-C2
     
151. A(J3,J1)=A(J3,J1)-C4
     
152. A(J1,J3)=A(J1,J3)-C4
     
153. A(J2,J2)=A(J2,J2)+C3
     
154. A(J3,J2)=A(J3,J2)+C5
     
155. A(J2,J3)=A(J2,J3)+C5
     
156. A(J3,J3)=A(J3,J3)+C6
     
157. C　　　
     
158.       C1=RH*CA*L/6
     
159. C2=RH*CA*L/420
     
160. SSM(1,1)=2.0*C1
     
161. SSM(4,1)=C1
     
162. SSM(2,2)=156.0*C2
     
163. SSM(3,2)=-22.0*L*C2
     
164. SSM(5,2)=54.0*C2
     
165. SSM(6,2)=13.0*L*C2
     
166. SSM(2,3)=-22.0*L*C2
     
167. SSM(3,3)=4.0*L*L*C2
     
168. SSM(5,3)=-13.0*L*C2
     
169. SSM(6,3)=-3.0*L**2*C2
     
170. SSM(1,4)=C1
     
171. SSM(4,4)=2.0*C1
     
172. SSM(2,5)=54.0*C2
     
173. SSM(3,5)=-13.0*L*C2
     
174. SSM(5,5)=156.0*C2
     
175. SSM(6,5)=22.0*L*C2
     
176. SSM(2,6)=13.0*L*C2
     
177. SSM(3,6)=-3.0*L**2*C2
     
178. SSM(5,6)=22.0*L*C2
     
179. SSM(6,6)=4.0*L**2*C2
     
180. C
     
181.       DO 3190 II=1,6
     
182.    DO 3190 JJ=1,6
     
183. 3190        S2(II,JJ)=0.0
     
184.       DO 3230 II=1,6
     
185.    DO 3230 JJ=1,6
     
186.       DO 3230 KK=1,6
     
187. 3230  S2(II,JJ)=S2(II,JJ)
     
188.      1     +SSM(II,KK)*DC(KK,JJ)
     
189. DO 3290 II=1,6
     
190.     DO 3290 JJ=1,6
     
191. 3290  SSM(II,JJ)=0.0
     
192.       DO 3350 II=1,6
     
193.     DO 3350 JJ=1,6
     
194.         DO 3350 KK=1,6
     
195. 3350    SSM(II,JJ)=SSM(II,JJ)+DC(KK,II)*S2(KK,JJ)
     
196. C
     
197.       I1=3*I-3
     
198. J1=3*J-3
     
199.   DO 3520 II=1,3
     
200.        DO 3520 JJ=1,3
     
201.        MG=I1+II
     
202.    ITR=I1+JJ
     
203.    IAC=J1+II
     
204.    IXK=J1+JJ
     
205. XX(MG,ITR)=XX(MG,ITR)+SSM(II,JJ)
     
206.       XX(IAC,ITR)=XX(IAC,ITR)+SSM(II+3,JJ)
     
207. XX(MG,IXK)=XX(MG,IXK)+SSM(II,JJ+3)
     
208. XX(IAC,IXK)=XX(IAC,IXK)+SSM(II+3,JJ+3)
     
209. 3520      CONTINUE
     
210. 3530  CONTINUE
     
211. C
     
212.       IF(NC.EQ.0) GOTO 3660
     
213. DO 3650 II=1,NC
     
214.    PO=JOD(II)
     
215. GMS=MAS1(II)
     
216. MI=MAS2(II)
     
217. K1=3*PO
     
218. J1=K1-1
     
219. I1=K1-2
     
220. XX(I1,I1)=XX(I1,I1)+GMS
     
221. XX(J1,J1)=XX(J1,J1)+GMS
     
222. XX(K1,K1)=XX(K1,K1)+MI
     
223. 3650  CONTINUE
     
224. C
     
225. 3660  MM=0
     
226.       DO 3880 I=1,NF
     
227. NA=NS(I)
     
228. MM=MM+1
     
229. NA=NA-MM+1
     
230. NB=N3-MM
     
231. IF (NB.LT.NA) GOTO 3790
     
232.   DO 3760 II=NA,NB
     
233.      DO 3760 JJ=NA,NB
     
234.       A(II,JJ)=A(II+1,JJ+1)
     
235. 3760  XX(II,JJ)=XX(II+1,JJ+1)
     
236. 3790  IF (NA.EQ.1) GOTO 3880
     
237.          DO 3850 II=1,NA-1
     
238.      DO 3850 JJ=NA,NB
     
239.             A(II,JJ)=A(II,JJ+1)
     
240. XX(II,JJ)=XX(II,JJ+1)
     
241. A(JJ,II)=A(JJ+1,II)
     
242. 3850 XX(JJ,II)=A(JJ+1,II)
     
243. 3880  CONTINUE
     
244. C
     
245. 1640 N9=N-1
     
246. DO 1830 IX=1,N
     
247. DI=A(IX,1)
     
248. DII=0.0
     
249. IF (ABS(DI).GT.1.0E-16)DII=1.0/DI
     
250. IF (ABS(DI).LT.1.0E-16)WRITE(11,'(A)')
     
251.      1   'MATRIX IS SINGULAR'
     
252. DO 1700 IY=1,N9
     
253. IY9=IY+1
     
254. 1700 A(IX,IY)=A(IX,IY9)/DI
     
255.       A(IX,N)=1.0/DI
     
256. DO 1820 IZ=1,N
     
257. IF (IZ.EQ.IX)GO TO 1820
     
258. O=A(IZ,1)
     
259. DO 1790 IY=1,N9
     
260. IY9=IY+1
     
261. 1790 A(IZ,IY)=A(IZ,IY9)-A(IX,IY)*O
     
262. A(IZ,N)=-A(IX,N)*O
     
263. 1820  CONTINUE
     
264. 1830 CONTINUE
     
265. C
     
266. 3900 DO 4000 I=1,N
     
267. DO 3960 J=1,N
     
268. VE(J)=0.0
     
269.   DO 3940 II=1,N
     
270. 3940   VE(J)=VE(J)+A(I,II)*XX(II,J)
     
271. 3960 CONTINUE
     
272. DO 3980 J=1,N
     
273. 3980 A(I,J)= VE(J)
     
274. 4000 CONTINUE
     
275. C
     
276. 500 MN=N
     
277. NN=N
     
278. CALL EIG1(A,X,Z,NN,N2,0.0001,XM,NKK)
     
279. 530 M=1
     
280. DO 560 I=1,NN
     
281. VC(I,M)=X(I)
     
282. 560 XX(I,M)=VC(I,M)
     
283. AM(M)=XM
     
284. C
     
285. IF (M1.LT.2) GOTO 220
     
286. DO 1410 M=2,M1
     
287. DO 630 I=1,NN
     
288. GK4=ABS(XX(I,M-1)-1.)
     
289. IF (GK4.LT.0.001)IR=I
     
290. 630 CONTINUE
     
291. IG(M-1)=IR
     
292. DO 680 I=1,NN
     
293. 680 XX(MN-I+1,MN-M+3)=A(IR,I)
     
294. DO 730 I=1,NN
     
295. DO 730 J=1,NN
     
296. J1=MN-J+1
     
297. J2=MN-M+3
     
298. 730 A(I,J)=A(I,J)-XX(I,M-1)*XX(J1,J2)
     
299. DO 910 I =1,NN
     
300. IF (I.EQ.IR)GOTO 910
     
301. IF (I.GT.IR)K1=I-1
     
302. IF (I.LE.IR)K1=I
     
303. DO 900 J=1,NN
     
304. IF (J.EQ.IR)GOTO 900
     
305. IF (J.GT.IR)K2=J-1
     
306. IF (J.LE.IR)K2=J
     
307. A(K1,K2)=A(I,J)
     
308. 900 CONTINUE
     
309. 910 CONTINUE
     
310. NN=NN-1
     
311. M3=NN
     
312. IF (M.NE.MN) GOTO 980
     
313. XM=A(1,1)
     
314. X(1)=1
     
315. 970 GOTO 990
     
316. 980 L1500=2
     
317. CALL EIG1(A,X,Z,NN,N2,0.0001,XM,NKK)
     
318. 990 DO 1000 I=1,NN
     
319. 1000 XX(I,M)=X(I)
     
320. AM(M)=XM
     
321. M4=M-1
     
322. M5=1000-M4
     
323. DO 1400 M8=M5,999
     
324. M6=M3+1
     
325. M2=1000-M8
     
326. M7=IG(M2)+1
     
327. IF(M6.LT.M7)GOTO 1150
     
328. N9=1000-M7
     
329. N8=1000-M6
     
330. DO 1130 I3=N8,N9
     
331. II=1000-I3
     
332. 1130 X(II)=X(II-1)
     
333. 1150 JJ=IG(M2)
     
334. X(JJ)=0.0
     
335. GSM=0.0
     
336. DO 1210 I=1,M6
     
337. J3=MN-I+1
     
338. J4=MN-M2+2
     
339. 1210 GSM=GSM+XX(J3,J4)*X(I)
     
340. XK=(AM(M2)-XM)/GSM
     
341. DO 1250 I=1,M6
     
342. 1250 X(I)=XX(I,M2)-XK*X(I)
     
343. GSM=0.0
     
344. DO 1300 I=1,M6
     
345. IF (ABS(GSM).LT.ABS(X(I)))GSM=X(I)
     
346. 1300 CONTINUE
     
347. DO 1330 I=1,M6
     
348. 1330 X(I)=X(I)/GSM
     
349. M3=M3+1
     
350. IF (M2.NE.1)GOTO 1400
     
351. DO 1380 I=1,M3
     
352. 1380 VC(I,M)=X(I)
     
353. 1400 CONTINUE
     
354. 1410 CONTINUE
     
355. C
     
356. 220 WRITE(11,'(/2X,A)')'计算结果部分：'
     
357. WRITE(11,'(/2X,A)')'输出各振型频率和振型'
     
358. DO 360 I=1,M1
     
359. TT=0.0
     
360. IF (AM(I).GT.1E-8)TT=SQRT(1.0/AM(I))/(2.0*3.14159)
     
361. WRITE(11,'(/2(2X,A,I2,A,F9.4))')'第',I,'振型频率：',TT,'振型：'
     
362. J=1
     
363. J0=1
     
364. DO 330 IT=1,N3
     
365. IF (J0.GT.NF)GOTO 320
     
366. DO 300 JS=J0,NF
     
367. IF(IT.EQ.NS(JS))GOTO 325
     
368. 300 CONTINUE
     
369. 320 VTL(IT)=VC(J,I)
     
370. J=J+1
     
371. GO TO 330
     
372. 325 VTL(IT)=0.0
     
373. J0=J0+1
     
374. 330 CONTINUE
     
375. WRITE(11,'(A8,A18,2A18/(I6,2E18.6,E20.6))')
     
376.      1 '结点号','水平方向位移','垂直方向位移','转角',
     
377.      2 (J,VTL(3*J-2),VTL(3*J-1),VTL(3*J),J=1,NJ)
     
378. 360 CONTINUE
     
379. C
     
380. ND=(N3-NF)/3
     
381. DO 1500 ME=1,MS
     
382. I=I9(ME)
     
383. J=J9(ME)
     
384. IF(Y(I).EQ.Y(J))ND=ND-1
     
385. 1500 CONTINUE
     
386. WRITE(11,*)'ND=',ND
     
387. DO 1510 I=1,ND
     
388. DO 1510 J=1,ND
     
389. AX=0.0
     
390. BX=0.0
     
391. SC(I,J)=AX*A(I,J)+BX*XX(I,J)
     
392. 1510 CONTINUE
     
393. CALL STEPM(A,XX,SC,ND,N3)
     
394. RETURN
     
395. END
     
396. C        SUBROUTINE EIG1(A,X,Z,NN,N2,D,XM,NKK)
     
397. IMPLICIT INTEGER*2 (I-N)
     
398. DIMENSION A(N2,N2),X(NN),Z(NN)
     
399. Y1=100000.0
     
400. NKK=1
     
401. DO 1320 I=1,NN
     
402. 1320 X(I)=1.0
     
403. XM=-100000.0
     
404. 1340 DO 1400 I=1,NN
     
405. SG=0.0
     
406. DO 1380 J=1,NN
     
407. 1380 SG=SG+A(I,J)*X(J)
     
408. 1400 Z(I)=SG
     
409. XM=0.0
     
410. DO 1450 I=1,NN
     
411. IF (ABS(XM).LT.ABS(Z(I)))XM=Z(I)
     
412. 1450 CONTINUE
     
413. DO 1500 I=1,NN
     
414. 1500 X(I)=Z(I)/XM
     
415. IF (ABS(Y1-XM)/XM.GT.D)GOTO 1530
     
416. GOTO 1550
     
417. 1530 Y1=XM
     
418. NKK=NKK+1
     
419. IF (NKK.GT.200)GOTO 1623
     
420. GOTO 1340
     
421. 1550 X3=0.0
     
422. DO 1580 I=1,NN
     
423. IF (ABS(X3).LT.ABS(X(I)))X3=X(I)
     
424. 1580 CONTINUE
     
425. DO 1610 I=1,NN
     
426. 1610 X(I)=X(I)/X3
     
427. RETURN
     
428. 1623 WRITE(*)'NUMBER OF ITERATION NKK=',NKK
     
429. 1630 FORMAT(A30,I2)
     
430. STOP
     
431. END
     
432. C
     
433. SUBROUTINE STEPM(SK,SM,SC,ND,N3)
     
434. IMPLICIT REAL*8(A-H,O-Z)
     
435. DIMENSION SK(N3,N3),SM(N3,N3),SC(N3,N3),F(N3,N3),X(N3,N3),
     
436.      1 DUA(N3),UD(N3),UV(N3),UA(N3),TC(N3),P(N3),NEQ(N3)
     
437. ND1=ND+1
     
438. DO 10 I=1,N3
     
439. UD(I)=0.0
     
440. UV(I)=0.0
     
441. DO 10 J=1,N3
     
442. 10 F(I,J)=0.0
     
443. WRITE(11,'(/1X,A)')'TIME PARA FOR ITEGRATION'
     
444. READ(7,*)THETA,DT,TMAX
     
445. WRITE(11,'(4X,3A6/4X,3F7.4)')
     
446.      1 'THETA','DT','TMAX',THETA,DT,TMAX
     
447. WRITE(11,'(/1X,A)')'NO.OF LOAD DATA POINTS'
     
448. READ(7,*)(NEQ(L),L=1,ND)
     
449. WRITE(11,'(4X,A6,I2,A3,I5)')('NEQ(',L,')=',NEQ(L),L=1,ND)
     
450. ANN=0.0
     
451. II=1
     
452. DO 60 ID=1,ND
     
453. NE=NEQ(ID)
     
454. IF (NE.EQ.0) GOTO 60
     
455. IF(NE.GT.TMAX/DT)NE=TMAX/DT+0.5
     
456. READ(7,*)(TC(J),P(J),J=1,NE)
     
457. WRITE(11,'(/1X,A)')'DYNAMIC LOAD ON STRUCTURE'
     
458. WRITE(11,'(4X,A12,A15)')'TC(I)','P(I)'
     
459. WRITE(11,'(4X,F12.5,E15.6)')(TC(J),P(J),J=1,NE)
     
460. NT=TC(NE)/DT+0.5
     
461. NT1=NT+1
     
462. NT2=NT+2
     
463. F(ID,1)=P(1)
     
464. ANN=0.0
     
465. II=1
     
466. DO 50 I=2,NT2
     
467. AI=I-1
     
468. T=AI*DT
     
469. IF(T.GT.TC(NE))GOTO 60
     
470. IF(T.LE.TC(II+1))GOTO 40
     
471. ANN=-TC(II+1)+T-DT
     
472. II=II+1
     
473. 40 ANN=ANN+DT
     
474. F(ID,I)=P(II)+(P(II+1)-P(II))*ANN/(TC(II+1)-TC(II))
     
475. 50 CONTINUE
     
476. 60 CONTINUE
     
477. NT=TMAX/DT+0.5
     
478. NT1=NT+1
     
479. DO 70 I=1,ND
     
480. X(I,ND1)=F(I,J)
     
481. DO 70 J=1,ND
     
482. 70 X(I,J)=SM(I,J)
     
483. DO 80 I=1,ND
     
484. DO 80 J=1,ND
     
485. 80 X(I,ND1)=X(I,ND1)-SC(I,J)*UV(J)-SK(I,J)*UD(J)
     
486. CALL SOLVE(ND,X,N3)
     
487. DO 90 I=1,ND
     
488. 90 UA(I)=X(I,ND1)
     
489. WRITE(11,'(/1X,A)')'RESPONSE OF STRUCTURE'
     
490. WRITE(11,'(/3X,2A5,3A11)')'COR. ','TIME','DISPL.','VEL.','ACC.'
     
491. TU=THETA*DT
     
492. A1=3.0/TU
     
493. A2=6.0/TU
     
494. A3=TU/2.0
     
495. A4=A2/TU
     
496. DO 190 L=1,NT1
     
497. DO 110 I=1,ND
     
498. DO 110 J=1,ND
     
499. 110 X(I,J)=SK(I,J)+A4*SM(I,J)+A1*SC(I,J)
     
500. AL=L
     
501. T=AL*DT
     
502. DO 130 I=1,ND
     
503. X(I,ND1)=(F(I,L+1)+F(I,L+2)-F(I,L+1))*(THETA-1.0)-F(I,L)
     
504. DO 120 J=1,ND
     
505. 120 X(I,ND1)=X(I,ND1)+(SM(I,J)*A2+SC(I,J)*3.0)*UV(J)+(SM(I,J)*3.0
     
506.      2 +A3*SC(I,J))*UA(J)
     
507. 130 CONTINUE
     
508. CALL SOLVE(ND,X,N3)
     
509. DO 140 I=1,ND
     
510. DUA(I)=A4*X(I,ND1)-A2*UV(I)-3.0*UA(I)
     
511. DUA(I)=DUA(I)/THETA
     
512. DUV=DT*UA(I)+DT*DUA(I)/2.0
     
513. UD(I)=UD(I)+DT*UV(I)+DT*DT*UA(I)/2.0+DT*DT*DUA(I)/6.0
     
514. UV(I)=UV(I)+DUV
     
515. 140 CONTINUE
     
516. DO 160 I=1,ND
     
517. X(I,ND1)=F(I,L+1)
     
518. DO 150 J=1,ND
     
519. X(I,ND1)=X(I,ND1)-SK(I,J)*UD(J)-SC(I,J)*UV(J)
     
520. 150 X(I,J)=SM(I,J)
     
521. 160 CONTINUE
     
522. CALL SOLVE(ND,X,N3)
     
523. DO 170 I=1,ND
     
524. UA(I)=X(I,ND1)
     
525. 170 WRITE(11,'(1X,I3,1X,F8.6,3E12.4)')I,T,UD(I),
     
526.      1 UV(I),UA(I)
     
527. 190 CONTINUE
     
528. RETURN
     
529. END
     
530. C
     
531. SUBROUTINE SOLVE(N,A,N3)
     
532. IMPLICIT REAL*8(A-H,O-Z)
     
533. DIMENSION A(N3,N3+1)
     
534. M=1
     
535. EPS=1.0E-10
     
536. NPLUSM=N+M
     
537. DET=1.0
     
538. DO 40 K=1,N
     
539. DET=DET*A(K,K)
     
540. IF(DABS(A(K,K)).GT.EPS)GOTO 10
     
541. GOTO 60
     
542. 10 KP1=K+1
     
543. DO 20 J=KP1,NPLUSM
     
544. 20 A(K,J)=A(K,J)/A(K,K)
     
545. A(K,K)=1.0
     
546. DO 40 I=1,N
     
547. IF (I.EQ.K.OR.A(I,K).EQ.0)GOTO 40
     
548. DO 30 J=KP1,NPLUSM
     
549. 30 A(I,J)=A(I,J)-A(I,K)*A(K,J)
     
550. A(I,K)=0.0
     
551. 40 CONTINUE
     
552. 60 RETURN
     
553. END

复制代码