平面弹性力学有限元源程序（FORTRAN）

hao1982

1. $DEBUG
   
2.                 PROGRAM PLANE
   
3.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
   
4.                 ALLOCATABLE::IJK(:,:),XY(:,:),BCA(:,:),SK(:,:),STR(:,:),MB(:,:),ZB(:),B(:)
   
5.                 ALLOCATABLE::DELD(:,:,:),TOD(:,:),DELST(:,:,:),TOST(:,:),DELSUP(:,:),TOTSUP(:)
   
6.                 DIMENSION EK(6,6)
   
7.         CHARACTER PN*40,FN*12
   
8. 
   
9.         WRITE(*,'(A)') ' 本程序为计算平面问题的有限元程序'
   
10.         WRITE(*,'(A)') ' 特点：(1)采用三结点三角形单元；'
    
11.         WRITE(*,'(A)') '       (2)采用等带宽存贮技术；'
    
12.         WRITE(*,'(A)') '       (3)采用高斯消元法解线性方程组。'
    
13.         WRITE(*,'(/A)') ' 输入计算问题名(PN):'
    
14.         READ(*,'(A)') PN
    
15.         CALL FNAME(PN,'.DAT',FN)
    
16.         WRITE(*,'(2A)') '  输入数据文件名为：',FN
    
17.         OPEN(5,FILE=FN,STATUS='OLD')
    
18.         CALL FNAME(PN,'.OUT',FN)
    
19.         WRITE(*,'(/2A)') ' 结果输出数据文件名为： ',FN
    
20.         OPEN(6,FILE=FN,STATUS='UNKNOWN')
    
21.         CALL FNAME(PN,'.OU1',FN)
    
22.         WRITE(*,'(/2A)') ' 参数输出数据文件名为： ',FN
    
23.         OPEN(7,FILE=FN,STATUS='UNKNOWN')
    
24. 
    
25.                 READ(5,*) NG,NE,MC,NX,NB,EO,VO,DENSITY,T
    
26.                 WRITE(6,120) NG,NE,MC,NX,NB
    
27.                 WRITE(6,130) EO,VO,DENSITY,T
    
28.                 READ(5,*) NWA,NWE,NWK,NWP
    
29.                 NT=2*NG
    
30.                 ALLOCATE (IJK(3,NE),XY(2,NG),BCA(7,NE),STR(3,NE),MB(2,NB),ZB(NB),B(NT))
    
31.                 ALLOCATE (DELD(2,NG,NX),TOD(2,NG),DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB))
    
32.                 CALL CLEAR(2,NG,TOD)
    
33.                 CALL CLEAR(3,NE,TOST)
    
34.                 CALL CLEAR1(NB,TOTSUP)
    
35.                 IF (NG.EQ.0) THEN
    
36.                   STOP 000
    
37.                 ENDIF
    
38.                 CALL INPUT(NG,NE,NB,IJK,XY,MB,ZB)
    
39.                 CALL CALND(NG,NE,IJK,ND)
    
40.                 ALLOCATE (SK(NT,ND))
    
41.                 IF(MC.EQ.0) THEN
    
42.                 E=EO
    
43.                 V=VO
    
44.                 ELSE
    
45.                 E=EO/(1-VO*VO)
    
46.                 V=VO/(1-VO)
    
47.                 ENDIF
    
48.                 IP=0
    
49.                 NX1=NX
    
50.                 A1=E/(1-V*V)/4.0
    
51.                 A2=0.5*(1-V)
    
52.                 CALL ABC(NG,NE,IJK,XY,BCA,NWA)
    
53.                 CALL CLEAR(NT,ND,SK)
    
54.                 DO 100 K=1,NE
    
55.                 CALL KE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE)
    
56.                 CALL SUMK(K,NE,NT,ND,IJK,EK,SK)
    
57. 100                CONTINUE
    
58.                 CALL CHECK(NT,ND,SK)
    
59. 110                CONTINUE
    
60.                 IP=IP+1
    
61.                 WRITE(6,'(///1X,A,I4)') "荷载工况=",IP
    
62.                 READ(5,*) NF,NP,NM
    
63.                 DO I=1,NT
    
64.                 B(I)=0.0
    
65.                 ENDDO
    
66.                 IF(NF.GT.0) THEN
    
67.                 CALL PF(NF,NT,B)
    
68.                 ENDIF
    
69.                 IF(NP.GT.0) THEN
    
70.                 CALL PP(NP,NG,NT,T,XY,B)
    
71.                 ENDIF
    
72.                 IF(NM.GT.0) THEN
    
73.                 CALL PG(NM,NE,NT,DENSITY,T,IJK,BCA,B)
    
74.                 ENDIF
    
75.                 DO I=1,NB
    
76.                 I1=2*(MB(1,I)-1)+2-MB(2,I)
    
77.                 DELSUP(I,IP)=-B(I1)
    
78.                 ENDDO
    
79.                 IF(NF.NE.0.OR.NP.NE.0.OR.NM.NE.0) THEN
    
80.                 CALL DBC(NB,ND,NT,NX,NX1,MB,ZB,SK,B)
    
81.                 CALL GAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP)
    
82.                 CALL STRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR)
    
83.                 CALL TRES(IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST)
    
84.                 CALL SUPT(IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP,TOTSUP)
    
85.                 CALL OUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP)
    
86.                 ELSE
    
87.                 WRITE(*,'(2x,A,I4,A)') "荷载工况=",IP,"没有荷载!"
    
88.                 WRITE(6,'(2x,A,I4,A)') "荷载工况=",IP,"没有荷载!"
    
89.                 ENDIF
    
90.                 NX1=NX1-1
    
91.                 IF (NX1.GT.0) GOTO 110
    
92. 120                FORMAT(1X, '结点总数=', I3, 1X, '单元总数=', I3, 1X, '问题类型=', &
    
93.       & I1, 1X, '荷载工况数=', I2, 1X, '支承位移数=', I2)
    
94. 130                FORMAT(/1X, '弹性模量=', E10.3, 1X, '泊松比=', F5.3, 1X, '密度=', E10.3, &
    
95.       & 1X, '厚度=', F6.3)
    
96.                 END
    
97. 
    
98.                                
    
99.                 SUBROUTINE GAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP)
    
100.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
101.                 DIMENSION SK(NT,ND),A(NT,NT),B(NT)
     
102.                 CALL CLEAR(NT,NT,A)
     
103.                 DO I=1,NT
     
104.                  DO J=1,ND
     
105.                   IF((I+J-1).LE.NT) THEN
     
106.                   A(I,I+J-1)=SK(I,J)
     
107.                   ENDIF
     
108.                  ENDDO
     
109.                 ENDDO
     
110.                 IF(NWK.EQ.1.AND.NX1.EQ.NX) THEN
     
111.             WRITE(7,'(A)') "结构总刚(列出右上三角元素)"
     
112.                 DO I=1,NT
     
113.                 WRITE(7,100) (A(I,J),J=1,NT)
     
114.                 ENDDO
     
115.                 ENDIF
     
116.                 IF (NWP.EQ.1) THEN
     
117.                 WRITE(7,'(1X,A,I3,A)') "第",NX-NX1+1,"荷载工况的荷载列阵"
     
118.                 DO I=1,NT
     
119.                 WRITE(7,'(1x,E11.4)') B(I)
     
120.                 ENDDO
     
121.                 ENDIF
     
122.                 DO 50 K=1,NT-1
     
123.                  DO 60 I=K+1,NT
     
124.                   C1=A(K,I)/A(K,K)
     
125.                   DO 70 J=I,NT
     
126.                   A(I,J)=A(I,J)-C1*A(K,J)
     
127. 70                  CONTINUE
     
128.                   B(I)=B(I)-C1*B(K)
     
129. 60                 CONTINUE
     
130. 50                CONTINUE
     
131.                 B(NT)=B(NT)/A(NT,NT)
     
132.                 DO 80 I=NT-1,1,-1
     
133.                  DO 90 J=I+1,NT
     
134.                  B(I)=B(I)-A(I,J)*B(J)
     
135. 90                 CONTINUE
     
136.                  B(I)=B(I)/A(I,I)
     
137. 80                CONTINUE
     
138. 100                FORMAT(1x, 4000E10.3)
     
139.                 END
     
140. 
     
141. 
     
142.                 SUBROUTINE CHECK(NT,ND,SK)
     
143.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
144.                 DIMENSION SK(NT,ND)
     
145.                 M=0
     
146.                 DO I=1,NT
     
147.                  IF(SK(I,1).LE.0.0) THEN
     
148.                  M=M+1
     
149.                  ENDIF
     
150.                 ENDDO
     
151.                 IF (M.GT.0) THEN
     
152.                 WRITE(*,*) "总刚矩阵的对角元素为负值!!"
     
153.                 STOP 2222
     
154.                 ENDIF
     
155.                 END
     
156. 
     
157.                 SUBROUTINE STRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR)
     
158.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
159.                 DIMENSION D1(6),D2(3),B(NT),S(3,6),IJK(3,NE),BCA(7,NE),STR(3,NE)
     
160.                 CALL CLEAR(3,NE,STR)
     
161.                 DO 10 K=1,NE
     
162.                  DO I=1,3
     
163.                  II=IJK(I,K)
     
164.                  D1(2*(I-1)+1)=B(2*(II-1)+1)
     
165.                  D1(2*(I-1)+2)=B(2*(II-1)+2)
     
166.                  ENDDO
     
167.                  CALL CLEAR(3,6,S)
     
168.                  C1=2*A1/BCA(7,K)
     
169.                  DO 20 I=1,3
     
170.                   S(1,2*(I-1)+1)=C1*BCA(I,K)
     
171.                   S(1,2*(I-1)+2)=C1*V*BCA(I+3,K)
     
172.                   S(2,2*(I-1)+1)=C1*V*BCA(I,K)
     
173.                   S(2,2*(I-1)+2)=C1*BCA(I+3,K)
     
174.                   S(3,2*(I-1)+1)=C1*A2*BCA(I+3,K)
     
175.                   S(3,2*(I-1)+2)=C1*A2*BCA(I,K)
     
176. 20                 CONTINUE
     
177.                  CALL MATMUL(3,6,1,S,D1,D2)
     
178.                  DO I=1,3
     
179.                  STR(I,K)=D2(I)
     
180.                  ENDDO
     
181. 10                CONTINUE
     
182.                 END
     
183. 
     
184.                 SUBROUTINE OUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP)
     
185.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
186.                 DIMENSION MB(2,NB),DELD(2,NG,NX),TOD(2,NG)
     
187.                 DIMENSION DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB)
     
188.                 CHARACTER VE*6
     
189.                 WRITE(6,20) IP
     
190.                 WRITE(6,30)
     
191.                 DO I=1,NG
     
192.                 WRITE(6,40) I,DELD(1,I,IP),TOD(1,I),DELD(2,I,IP),TOD(2,I)
     
193.                 ENDDO
     
194.                 WRITE(6,50)
     
195.                 DO J=1,NE
     
196.                 WRITE(6,60) J,(DELST(L,J,IP),TOST(L,J),L=1,3)
     
197.                 ENDDO
     
198.                 WRITE(6,70)
     
199.                 DO J=1,NB
     
200.                 IF (MB(2,J).EQ.1) THEN
     
201.                 VE='x方向'
     
202.                 ELSE
     
203.                 VE='Y方向'
     
204.                 ENDIF
     
205.                 WRITE(6,80) MB(1,J),VE,DELSUP(J,IP),TOTSUP(J)
     
206.                 ENDDO
     
207. 20                FORMAT(/1X, '荷载工况=', I3, '的计算结果')
     
208. 30                FORMAT(/, 1X, '结点号', 5X, 'X位移增量', 5X, 'X累计位移', 5x, &
     
209.       & 'Y位移增量', 5X, 'Y累计位移')
     
210. 40                FORMAT(1X, I3, 6X, F10.4, 4x, F10.4, 4X, F10.4, 4x, F10.4)
     
211. 50                FORMAT(/, 1X, '单元号', 6x, 'X应力增量', 6x, 'X累计应力', 6x, &
     
212.       & 'Y应力增量', 6x, 'Y累计应力', 6x, '剪应力增量', 6x, &
     
213.       & '累计剪应力')
     
214. 60                FORMAT(2X, I3, 6X, E10.3, 5X, E10.3, 5X, E10.3, 5X, E10.3, 6X, E10.3, &
     
215.       & 6X, E10.3)
     
216. 70                FORMAT(/, 1X, '支座结点号', 6x, '支反力方向', 6x, '支反力增量', &
     
217.       & 6x, '支反力累计量')
     
218. 80                FORMAT(5X, I3, 12X, A, 6x, E10.3, 8X, E10.3)
     
219.             END                                                           
     
220. 
     
221. 
     
222.                 SUBROUTINE MATMUL(M,N,L,A,B,C)
     
223.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
224.         DIMENSION A(M,N),B(N,L),C(M,L)
     
225.         DO 100 I=1,M
     
226.          DO 100 J=1,L
     
227.           C(I,J)=0.0
     
228.           DO 100 K=1,N
     
229. 100     C(I,J)=C(I,J)+A(I,K)*B(K,J)
     
230.         END
     
231. 
     
232.                 SUBROUTINE CLEAR(M,N,A)
     
233.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
234.                 DIMENSION A(M,N)
     
235.                 DO I=1,M
     
236.                 DO J=1,N
     
237.                 A(I,J)=0.0
     
238.                 ENDDO
     
239.                 ENDDO
     
240.                 END
     
241. 
     
242.                 SUBROUTINE ABC(NG,NE,IJK,XY,BCA,NWA)
     
243.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
244.                 DIMENSION X(2,5),XY(2,NG),IJK(3,NE),B(7),BCA(7,NE)
     
245.                 DO 10 I=1,NE
     
246.                  DO 20 K=1,3
     
247.                  DO 20 J=1,2
     
248.                  X(J,K)=XY(J,IJK(K,I))
     
249. 20                 CONTINUE
     
250.                  DO 30 J=1,2
     
251.                  X(J,4)=X(J,1)
     
252.                  X(J,5)=X(J,2)
     
253. 30                 CONTINUE
     
254.                  DO 40 K=1,3
     
255.                  B(K)=X(2,K+1)-X(2,K+2)
     
256.                  B(K+3)=X(1,K+2)-X(1,K+1)
     
257. 40                 CONTINUE
     
258.                  B(7)=(B(1)*B(5)-B(4)*B(2))*0.5
     
259.                  IF (B(7).LE.0.0) THEN
     
260.                  WRITE(*,*) "单元号=",I
     
261.                  WRITE(*,*) "单元的结点I,J AND K编号：",(IJK(K,I),K=1,3)
     
262.                  STOP 1111
     
263.                  ELSE
     
264.                  DO J=1,7
     
265.                  BCA(J,I)=B(J)
     
266.                  ENDDO
     
267.                  ENDIF
     
268. 10                CONTINUE
     
269.                 IF(NWA.EQ.1) THEN
     
270.                 WRITE(7,100)
     
271.                 DO I=1,NE
     
272.                 WRITE(7,110) I,(BCA(J,I),J=1,7)
     
273.                 ENDDO
     
274.                 ENDIF
     
275. 100                 FORMAT(1X, '单元号', 3x, 'bi', 7x, 'bj', 7x, 'bk', 7x, 'ci', 7x, &
     
276.       & 'cj', 7x, 'ck', 7x, '面积A')
     
277. 110                 FORMAT(2X, I3, 3X, F6.3, 3X, F6.3, 3X, F6.3, 3X, F6.3, 3X, F6.3, &
     
278.       & 3X, F6.3, 3X, F6.3)
     
279.              END
     
280. 
     
281.                  SUBROUTINE SUMK(K,NE,NT,ND,IJK,EK,SK)
     
282.              IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
283.                  DIMENSION        IJ(3),EK(6,6),IJK(3,NE),SK(NT,ND)
     
284.                  DO I=1,3
     
285.                  IJ(I)=IJK(I,K)
     
286.                  ENDDO
     
287.                  DO 10 I=1,3
     
288.                   DO 20 J=1,3
     
289.                   IF (IJ(I).GT.IJ(J)) GOTO 20
     
290.                   M=2*(IJ(I)-1)+1
     
291.                   N=2*(IJ(J)-1)+1-(2*(IJ(I)-1)+1)+1
     
292.                   MO=2*(I-1)+1
     
293.                   NO=2*(J-1)+1
     
294.                   SK(M,N)=SK(M,N)+EK(MO,NO)
     
295.                   SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1)
     
296.                   SK(M+1,N)=SK(M+1,N)+EK(MO+1,NO+1)
     
297.                   IF (IJ(I).EQ.IJ(J)) GOTO 20
     
298.                   SK(M+1,N-1)=SK(M+1,N-1)+EK(MO+1,NO)
     
299. 20                  CONTINUE
     
300. 10             CONTINUE
     
301.                  END
     
302. 
     
303.                 SUBROUTINE KE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE)
     
304.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
305.                 DIMENSION  IJ(3),IJK(3,NE),BCA(7,NE),EK(6,6)
     
306.                 CALL CLEAR(6,6,EK)
     
307.                 DO I=1,3
     
308.                 IJ(I)=IJK(I,K)
     
309.                 ENDDO
     
310.                 DO 10 I=1,3
     
311.                  DO 20 J=I,3
     
312.                  IO=2*(I-1)
     
313.                  JO=2*(J-1)
     
314.                  EK(IO+1,JO+1)=BCA(I,K)*BCA(J,K)+A2*BCA(I+3,K)*BCA(J+3,K)
     
315.                  EK(IO+1,JO+2)=V*BCA(I,K)*BCA(J+3,K)+A2*BCA(I+3,K)*BCA(J,K)
     
316.                  EK(IO+2,JO+2)=BCA(I+3,K)*BCA(J+3,K)+A2*BCA(I,K)*BCA(J,K)
     
317.                  IF (I.EQ.J) GOTO 20
     
318.                  EK(IO+2,JO+1)=V*BCA(I+3,K)*BCA(J,K)+A2*BCA(I,K)*BCA(J+3,K)
     
319. 20                 CONTINUE
     
320. 10                CONTINUE
     
321.                 DO I=2,6
     
322.                  DO J=1,I-1
     
323.                  EK(I,J)=EK(J,I)
     
324.                  ENDDO
     
325.                 ENDDO
     
326.                 DO I=1,6
     
327.                  DO J=1,6
     
328.                  EK(I,J)=EK(I,J)*A1*T/BCA(7,K)
     
329.                  ENDDO
     
330.                 ENDDO
     
331.                 IF(NWE.EQ.1) THEN
     
332.              WRITE(7,40)        K, "号单元的单刚"
     
333.                 DO I=1,6
     
334.                  WRITE(7,50) (EK(I,J),J=1,6)
     
335.                 ENDDO
     
336.                 ENDIF
     
337. 40                FORMAT(1X, I3, A)
     
338. 50                FORMAT(1X, 6(E10.4, 4X))
     
339.                 END
     
340. 
     
341.                 SUBROUTINE PF(NF,NT,B)
     
342.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
343.                 DIMENSION MF(2,NF),ZF(NF),B(NT)
     
344.                 DO I=1,NF
     
345.                 READ(5,*) K,MF(1,I),MF(2,I),ZF(I)
     
346.                 ENDDO
     
347.                 WRITE(6,'(/A,I3)') "集中荷载个数NF=",NF
     
348.                 WRITE(6,20)
     
349.                 DO I=1,NF
     
350.                 WRITE(6,30) MF(1,I),MF(2,I),ZF(I)
     
351.                 ENDDO
     
352.                 DO 10 I=1,NF
     
353.                  II=2*(MF(1,I)-1)-MF(2,I)+2
     
354.                  B(II)=B(II)+ZF(I)
     
355. 10                CONTINUE
     
356. 20                FORMAT(1X, '结点号', 4x, '方向', 4x, '荷载值')
     
357. 30                FORMAT(2X, I3, 5X, I3, 4X, E10.4)
     
358.                 END               
     
359. 
     
360.         SUBROUTINE FNAME(PN,FN2,FN)
     
361.         CHARACTER PN*40,FN2*4,FN*12
     
362. !       去掉PN中前面的空格
     
363.         DO 10 I=1,40
     
364.          IF(PN(I:I).EQ.' ') GOTO 10
     
365.          IP=I
     
366.          GOTO 20
     
367. 10      CONTINUE
     
368. 20      CONTINUE
     
369.         FN(1:8)=PN(IP:IP+7)
     
370. !       去掉FN中后面的空格
     
371.         DO 30 I=8,1,-1
     
372.          IF(FN(I:I).EQ.' ') GOTO 30
     
373.          IL=I
     
374.          GOTO 40
     
375. 30      CONTINUE
     
376. 40      CONTINUE
     
377. !       生成文件名FN=PN+FN2
     
378.         FN(IL+1:IL+4)=FN2(1:4)
     
379.         END
     
380. 
     
381.                 SUBROUTINE PP(NP,NG,NT,T,XY,B)
     
382.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
383.                 DIMENSION MP(2,NP),ZP(NP),XY(2,NG),B(NT)
     
384.                 DO I=1,NP
     
385.                 READ(5,*) K,MP(1,I),MP(2,I),ZP(I)
     
386.                 ENDDO
     
387.                 WRITE(6,'(/A,I3)') "分布荷载数NP=",NP
     
388.                 WRITE(6,20)
     
389.                 DO I=1,NP
     
390.                 WRITE(6,30) I,MP(1,I),MP(2,I),ZP(I)
     
391.                 ENDDO
     
392.                 DO 10 I=1,NP
     
393.                  II=2*(MP(1,I)-1)
     
394.                  JJ=2*(MP(2,I)-1)
     
395.                  PX=0.5*ZP(I)*T*(XY(2,MP(1,I))-XY(2,MP(2,I)))
     
396.                  PY=0.5*ZP(I)*T*(XY(1,MP(2,I))-XY(1,MP(1,I)))
     
397.                  B(II+1)=B(II+1)+PX
     
398.                  B(II+2)=B(II+2)+PY
     
399.                  B(JJ+1)=B(JJ+1)+PX
     
400.                  B(JJ+2)=B(JJ+2)+PY
     
401. 10                CONTINUE
     
402. 20                FORMAT(1X, '序号', 4x, '首结点号', 4x, '末结点号', 7x, '荷载值')
     
403. 30                FORMAT(1X, I3, 8x, I3, 9X, I3, 7X, F10.3)
     
404.                 END               
     
405. 
     
406.                 SUBROUTINE PG(NM,NE,NT,DENSITY,T,IJK,BCA,B)
     
407.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
408.                 DIMENSION IJK(3,NE),BCA(7,NE),MG(NM),B(NT)
     
409.                 IF (NM.GE.NE) THEN
     
410.                 NM=NE
     
411.                 DO I=1,NM
     
412.                 MG(I)=I
     
413.                 ENDDO
     
414.                 ELSE
     
415.                 READ(5,*) (MG(I),I=1,NM)
     
416.                 ENDIF
     
417.                 WRITE(6,'(/A,I3)') "计自重单元数NM=",NM
     
418.                 WRITE(6,'(A)') "计自重的单元号为："
     
419.                 K1=0
     
420.                 K2=0
     
421.                 DO WHILE((NM-K2).GE.10)
     
422.                 WRITE(6,100) (MG(K2+J),J=1,10)
     
423.                 K1=K1+1
     
424.                 K2=10*K1
     
425.                 ENDDO
     
426.                 WRITE(6,100) (MG(K2+J),J=1,NM-K2)
     
427.                 DO 10 I=1,NM
     
428.                  II=2*(IJK(1,MG(I))-1)
     
429.                  JJ=2*(IJK(2,MG(I))-1)
     
430.                  KK=2*(IJK(3,MG(I))-1)
     
431.                  PE1=-DENSITY*T*BCA(7,MG(I))/3.0
     
432.                  B(II+2)=B(II+2)+PE1
     
433.                  B(JJ+2)=B(JJ+2)+PE1
     
434.                  B(KK+2)=B(KK+2)+PE1
     
435. 10                CONTINUE
     
436. 100                FORMAT(10(1X,I3))
     
437.                 END               
     
438. 
     
439.                 SUBROUTINE INPUT(NG,NE,NB,IJK,XY,MB,ZB)
     
440.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
441.                 DIMENSION IJK(3,NE),XY(2,NG),MB(2,NB),ZB(NB)
     
442.                 DO I=1,NE
     
443.                 READ(5,*) K,(IJK(J,I),J=1,3)
     
444.                 ENDDO
     
445.                 DO I=1,NG
     
446.                 READ(5,*) K,XY(1,I),XY(2,I)
     
447.                 ENDDO
     
448.                 DO I=1,NB
     
449.                  READ(5,*) K,MB(1,I),MB(2,I),ZB(I)
     
450.                 ENDDO
     
451.                 WRITE(6,15)
     
452.                 DO K=1,NE
     
453.                 WRITE(6,20) K,(IJK(J,K),J=1,3)
     
454.                 ENDDO
     
455.                 WRITE(6,30)
     
456.                 DO K=1,NG
     
457.                 WRITE(6,35) K,XY(1,K),XY(2,K)
     
458.                 ENDDO
     
459.                 WRITE(6,40)
     
460.                 DO K=1,NB
     
461.                  WRITE(6,45) K,MB(1,K),MB(2,K),ZB(K)
     
462.                 ENDDO
     
463. 15                FORMAT(/1X, '单元号', 4X, 'I结点号', 4x, 'J结点号', 4x, 'K结点号')
     
464. 20                FORMAT(2X, I3, 7X, I3, 8X, I3, 8X, I3)
     
465. 30                FORMAT(/1X, '结点号', 8X, 'X坐标', 8X, 'Y坐标')
     
466. 35                FORMAT(2X, I3, 6X, F8.3, 5X, F8.3)
     
467. 40                FORMAT(/1X, '序号', 4x, '支承结点号', 6x, '支承方向', 6x, &
     
468.       & '支承位移')
     
469. 45                FORMAT(2X, I2, 8X, I4, 11X, I2, 7X, F10.3)
     
470.                 END
     
471. 
     
472.                 SUBROUTINE DBC(NB,ND,NT,NX,NX1,MB,ZB,SK,B)
     
473.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
474.                 DIMENSION  MB(2,NB),ZB(NB),SK(NT,ND),B(NT)
     
475.                 DO 5 I=1,NB
     
476.                  N=2*MB(1,I)-MB(2,I)
     
477.                  Z=ZB(I)
     
478.                  IF (ABS(Z).LT.1.0E-10) THEN
     
479.                   IF (NX.NE.NX1) GOTO 30
     
480.                   SK(N,1)=1.0
     
481.                   DO J=2,ND
     
482.                   SK(N,J)=0.0
     
483.                   ENDDO
     
484.                   DO 10 K=2,ND
     
485.                   IF (N.LT.K) GOTO 20
     
486.                   M=N-K+1
     
487.                   SK(M,K)=0.0
     
488. 10                  CONTINUE
     
489. 20                  CONTINUE
     
490. 30                  B(N)=0.0
     
491.                  ELSE
     
492.                   IF (NX.NE.NX1) GOTO 40
     
493.                   SK(N,1)=SK(N,1)*1.0E+15
     
494. 40                  B(N)=SK(N,1)*Z
     
495.                  ENDIF
     
496. 5                CONTINUE
     
497.                 END
     
498. 
     
499.                 SUBROUTINE CALND(NG,NE,IJK,ND)
     
500.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
501.                 DIMENSION  IJK(3,NE)
     
502.                 ND=0
     
503.                 DO 10 I=1,NG
     
504.                  M1=0
     
505.                  DO 20 J=1,NE
     
506.                   DO 30 K=1,3
     
507.                   IF(IJK(K,J).EQ.I) THEN
     
508.                   M2=0
     
509.                    DO 40 L=1,3
     
510.                     M3=IJK(L,J)-I
     
511.                     IF(M3.GT.M2) THEN
     
512.                     M2=M3
     
513.                     ENDIF
     
514. 40                   CONTINUE
     
515.                   GOTO 25
     
516.                   ENDIF
     
517. 30                  CONTINUE
     
518. 25                 IF(M2.GT.M1) M1=M2
     
519. 20                CONTINUE
     
520.                 IF(M1.GT.ND)ND=M1
     
521. 10                CONTINUE
     
522.                 ND=2*(ND+1)
     
523.                 WRITE(7,100) ND
     
524. 100                FORMAT(1X,'总刚矩阵的半带宽ND为：',I4)
     
525.                 END
     
526. 
     
527.                 SUBROUTINE TRES(IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST)
     
528.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
529.                 DIMENSION B(NT),STR(3,NE),DELD(2,NG,NX),TOD(2,NG),DELST(3,NE,NX),TOST(3,NE)
     
530.                 DO I=1,NG
     
531.                 IO=2*(I-1)
     
532.                 DELD(1,I,IP)=B(IO+1)
     
533.                 DELD(2,I,IP)=B(IO+2)
     
534.                 TOD(1,I)=TOD(1,I)+DELD(1,I,IP)
     
535.                 TOD(2,I)=TOD(2,I)+DELD(2,I,IP)
     
536.                 ENDDO
     
537.                 DO I=1,NE
     
538.                 DO J=1,3
     
539.                 DELST(J,I,IP)=STR(J,I)
     
540.                 TOST(J,I)=TOST(J,I)+DELST(J,I,IP)
     
541.                 ENDDO
     
542.                 ENDDO
     
543.                 END
     
544. 
     
545. 
     
546.                 SUBROUTINE SUPT(IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP,TOTSUP)
     
547.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
548.                 DIMENSION IJK(3,NE),MB(2,NB),BCA(7,NE),DELD(2,NG,NX),EK(6,6)
     
549.                 DIMENSION D1(6),DELSUP(NB,NX),TOTSUP(NB)
     
550.                 DO 10 K=1,NB
     
551.                 I1=MB(1,K)
     
552.                 I2=MB(2,K)
     
553.                 RF=0.0
     
554.                  DO 20 I=1,NE
     
555.                   DO 30 J=1,3
     
556.                   I3=IJK(J,I)
     
557.                   IF (I1.NE.I3) THEN
     
558.                    GOTO 30
     
559.                   ELSE
     
560.                    I4=2*(J-1)+2-I2
     
561.                    CALL KE(I,NE,T,V,A1,A2,IJK,BCA,EK,0)
     
562.                    DO M=1,3
     
563.                    IO=IJK(M,I)
     
564.                    D1(2*(M-1)+1)=DELD(1,IO,IP)
     
565.                    D1(2*(M-1)+2)=DELD(2,IO,IP)
     
566.                    ENDDO
     
567.                    DO L=1,6
     
568.                    RF=RF+EK(I4,L)*D1(L)
     
569.                    ENDDO
     
570.                    GOTO 20
     
571.                   ENDIF
     
572. 30                  CONTINUE
     
573. 20                 CONTINUE
     
574.                 DELSUP(K,IP)=DELSUP(K,IP)+RF
     
575.                 TOTSUP(K)= TOTSUP(K)+DELSUP(K,IP)
     
576. 10                CONTINUE
     
577.                 END               
     
578.                                    
     
579.                 SUBROUTINE CLEAR1(N,B)
     
580.                 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
     
581.                 DIMENSION B(N)
     
582.                 DO I=1,N
     
583.                 B(I)=0.0
     
584.                 ENDDO
     
585.                 END

复制代码

PFEM,EXE平面问题有限元程序的数据文件输入格式：

基本输入数据文件名：*.DAT
计算结果输出数据文件名：*.OUT
中间参数输出数据文件名：*.OU1
下面按在文件中的前后顺序进行说明：
1 基本控制参数信息：NG,NE,MC,NX,NB,EO,VO,DENSITY ,T(共计5个整形数，4个实型数)
  NG：结构的结点总数；
  NE：结构的单元总数；
  MC：平面问题的类型，MC=0，为平面应力，MC=1，为平面应变；
  NX：荷载工况数；
  NB：支承位移数；
  EO：材料弹性模量(Pa)；
  VO：材料泊松比；
  DENSITY ：容重(N/m3)
  T ：材料厚度(m)；
2 打印输出控制参数：NWA,NEW,NWK,NWP(4个整形数)
  等于1时，输出，否则不输出。
3 单元结点信息：(K,(IJK(I,K),I=1,3),K=1,NE) (每行4个整形数，共计NE行)
  K：单元号；
  IJK(1,K)：K单元I结点编号；
  IJK(2,K)：K单元J结点编号；
  IJK(3,K)：K单元K结点编号；
4  结点坐标信息：((K,XY(1,K),XY(2,K)),K=1,NG)(每行3个整形数，共计NG行)
  K：结点号
  XY(1,K)：K结点X坐标；
  XY(2,K)：K结点Y坐标；
5 支承信息：((K,MB(1,K),MB(2,K),ZB(K)),K=1,NB)(每行3个整形数，1个实型数，共计NB行)
  K：支承位移序号；
  MB(1,K)：第K个支承位移所在的结点号；
  MB(2,K)：第K个支承位移的坐标方向；
  ZB(K)：  第K个支承位移的数值；
6 按NX荷载工况数输入荷载信息：每一荷载工况如下
  (1) NF,NP,NM(3个整型数)
     NF：集中荷载个数；
     NP：分布荷载个数；
     NM：计自重单元数；
  (2) 若NF≠0，则输入下面数据
      K,MF(1,K),MF(2,K),ZF(K)(每行3个整形数，1个实型数，共计NF行)
       K：集中荷载序号；
       MF(1,K)：第K个集中荷载作用的结点号；
       MF(2,K)：第K个集中荷载的坐标方向；
       ZF(K)：  第K个集中荷载的数值；
  (3) 若NP≠0，则输入下面数据
      K,MP(1,K),MP(2,K),ZP(K)(每行3个整形数，1个实型数，共计NP行)
       K：分布荷载序号；
       MP(1,K)：第K个分布荷载作用的结点号；
       MP(2,K)：第K个分布荷载的坐标方向；
       ZP(K)：  第K个分布荷载的数值；
  (4) 若NM≠0，则输入下面数据
      若NM≥NE，则表示计所有单元的自重，不需输入计自重的单元号；
      若NM<NE，则需要输入计自重的单元号；

例题验证

1 (选自徐芝伦编《弹性力学简明教程》(第二版))设有对角受压的正方形薄板，荷载沿厚度均匀分布，为2N/m，由于xz和yz面均为该薄板的对称面，所以只需取1/4部分作为计算对象，图1(a,b)，将该对象划分为4个单元，共有6个结点。对称面上的结点没有垂直于对称面的位移分量，因此，在1、2、4三个结点上设置了水平链杆支座，在4、5、6三个结点设置了铅直链杆支座，如图1c所示。