! under construction (but running) ! ED WILSON'S 1964 TRANSIENT HEAT CONDUCTION & CONVECTION PROGRAM ! FOR PLANE AND AXISYMMETRIC STRUCTURES ! (THE FIRST NON-STRUCTURAL FINITE ELEMENT PROGRAM) ! CONVERTED TO F90 FROM F66 BY ED AKIN, JAN 2000 ! COMMON B (500), X (500), Y (500), T (500), D (500), KODE (500),& HED (20), LM (5), IX (3), E (3, 3), KX (4), P (5), S (5, 5), & DD (5), XCOND (20), SPHT (20), DENS (20), QX (20) INTEGER OK COMMON / SYMARG / NUMNP, MBAND, A (500, 27), Q (500) !****************************************************************** ! READ AND PRINT OF CONTROL INFORMATION !****************************************************************** 50 READ (5, 1000, IOSTAT = OK) HED, NUMNP, NUMEL, NUMCBC, KAT, & NUMMAT, NDT, INTER, DT 1000 FORMAT (20A4/7I5,E15.6) IF ( OK < 0 ) CALL EXIT (0) ! END OF FILE IF (KAT) 54, 56, 54 54 WRITE (6, 2010) 2010 FORMAT ('AXISYMMETRIC SOLID BODY') GOTO 58 56 WRITE (6, 2011) 2011 FORMAT ('TWO DIMENSIONAL PLANE BODY') 58 WRITE (6, 2000) HED, NUMNP, NUMEL, NUMCBC, NUMMAT, NDT, INTER, DT 2000 FORMAT (20A4,//, 'NUMBER OF NODAL POINTS--',I4,/,& & 'NUMBER OF ELEMENTS------',I4,/,'NUMBER OF CONVECTION BC-',I4,/,& & 'NUMBER OF MATERIALS-----',I4,/,'NUMBER OF INCREMENTS----',I4,/,& & 'OUTPUT INTERVAL---------',I4,/,'TIME INTERVAL------',E10.3) ! READ (5, 1003) (M, XCOND (M), SPHT (M), DENS (M), QX (M), & N = 1, NUMMAT) 1003 FORMAT (I10,4F10.0) WRITE (6, 2009) (M, XCOND (M), SPHT (M), DENS (M), QX (M), & M = 1, NUMMAT) 2009 FORMAT (' M',14X, 'K',14X, 'C',14X, 'D',14X, 'Q',/, & (I4,4E15.6)) !****************************************************************** ! READ OR GENERATE NODAL POINT INFORMATION !****************************************************************** WRITE (6, 2001) 2001 FORMAT ('N.P. NO. CODE',14X, 'X',14X, 'Y',14X, 'T') L = 1 ! 60 READ (5, 1001) N, KODE (N), X (N), Y (N), T (N) 1001 FORMAT (2I5,3F10.0) DIFF = N + 1 - L IF (N - L) 65, 80, 70 65 WRITE (6, 2020) N 2020 FORMAT ('CARD NO. ',I4, ' OUT OF ORDER') GOTO 60 70 DX = (X (N) - X (L - 1) ) / DIFF DY = (Y (N) - Y (L - 1) ) / DIFF 75 KODE (L) = 0 X (L) = X (L - 1) + DX Y (L) = Y (L - 1) + DY T (L) = 0.0 80 WRITE (6, 2002) L, KODE (L), X (L), Y (L), T (L) 2002 FORMAT (2I10,3E15.6) L = L + 1 IF (N - L) 90, 80, 75 90 IF (NUMNP + 1 - L) 100, 100, 60 100 CONTINUE !****************************************************************** ! FORM CONDUCTIVITY MATRIX FOR COMPLETE BODY !****************************************************************** DO 111 I = 1, NUMNP D (I) = 0.0 B (I) = 0.0 DO 110 J = 1, 20 A (I, J) = 0.0 110 END DO 111 END DO MBAND = 0 NUM = 0 WRITE (6, 2003) 2003 FORMAT (' N I J K L MATERIAL' ) ! DO 200 N = 1, NUMEL ! ! 1. READ OR GENERATE ELEMENT PROPERTIES ! IF (NUM - N) 120, 121, 121 120 READ (5, 1002) NUM, KX (1), KX (2), KX (3), KX (4), MTYPE 1002 FORMAT (6I5) ! 121 DO 122 I = 1, 4 LM (I) = LM (I) + 1 122 END DO ! IF (NUM - N) 123, 124, 126 123 WRITE (6, 2021) NUM 2021 FORMAT ('BAD CARD NO. ',I4) GOTO 120 ! 124 DO 125 I = 1, 4 LM (I) = KX (I) 125 END DO COND = XCOND (MTYPE) ! 126 WRITE (6, 2004) N, LM (1), LM (2), LM (3), LM (4), MTYPE 2004 FORMAT (5I5,I10) ! ! 2. FORM ELEMENT CONDUCTIVITY MATRIX ! DO 149 I = 1, 5 DD (I) = 0.0 P (I) = 0.0 DO 150 J = 1, 5 S (I, J) = 0.0 150 END DO 149 END DO ! I = LM (1) J = LM (2) K = LM (3) L = LM (4) LM (5) = I ! XX = (X (I) + X (J) + X (K) + X (L) ) / 4. YY = (Y (I) + Y (J) + Y (K) + Y (L) ) / 4. ! DO 152 K = 1, 4 ! I = LM (K) J = LM (K + 1) IF (I - J) 135, 152, 135 135 AJ = X (J) - X (I) AK = XX - X (I) BJ = Y (J) - Y (I) BK = YY - Y (I) C = BJ - BK DX = AK - AJ ! XMUL = 1.0 IF (KAT) 136, 137, 136 136 XMUL = XMUL * (X (I) + X (J) + XX) / 3.0 ! 137 XLAM = AJ * BK - AK * BJ COMM = .5 * COND * XMUL / XLAM QQ = XLAM * XMUL * QX (MTYPE) / 4.0 QSTORE = XLAM * XMUL * SPHT (MTYPE) * DENS (MTYPE) / 4.0 ! E (1, 1) = C**2 + DX**2 E (1, 2) = BK * C - AK * DX E (1, 3) = - BJ * C + AJ * DX E (2, 1) = E (1, 2) E (2, 2) = BK**2 + AK**2 E (2, 3) = - BJ * BK - AJ * AK E (3, 1) = E (1, 3) E (3, 2) = E (2, 3) E (3, 3) = BJ**2 + AJ**2 ! IX (1) = K IX (2) = K + 1 IF (K - 4) 145, 140, 145 140 IX (2) = 1 145 IX (3) = 5 ! DO 148 I = 1, 3 II = IX (I) P (II) = P (II) + QQ DD (II) = DD (II) + QSTORE DO 151 J = 1, 3 JJ = IX (J) S (II, JJ) = S (II, JJ) + E (I, J) * COMM 151 END DO 148 END DO ! 152 END DO ! DO 154 I = 1, 4 DO 155 J = 1, 4 S (I, J) = S (I, J) - S (I, 5) * S (J, 5) / S (5, 5) 155 END DO 154 END DO ! ! 3. ADD ELEMENT CONDUCTIVITY TO COMPLETE CONDUCTIVITY MATRIX ! DO 174 L = 1, 4 I = LM (L) D (I) = D (I) + DD (L) B (I) = B (I) + P (L) DO 175 M = 1, 4 J = LM (M) - I + 1 IF (27 - J) 123, 158, 158 158 IF (MBAND-J) 160, 165, 165 160 MBAND = J 165 IF (J) 175, 175, 170 170 A (I, J) = A (I, J) + S (L, M) 175 END DO 174 END DO ! 200 END DO !****************************************************************** ! BOUNDARY CONDITIONS !****************************************************************** ! ! 1. CONVECTION BOUNDARY CONDITIONS ! IF (NUMCBC) 220, 220, 205 205 WRITE (6, 2006) 2006 FORMAT (' I J H TEMPERATURE') DO 215 N = 1, NUMCBC WRITE (6, 2004) NUMNP, NUMCBC, N READ (5, 1007) I, J, H, TEMP 1007 FORMAT (2I5,2F10.0) WRITE (6, 2007) I, J, H, TEMP 2007 FORMAT (2I5,2E15.6) XL = SQRT ( (X (J) - X (I) ) **2 + (Y (J) - Y (I) ) **2) IF (KAT) 206, 207, 206 206 XL = XL * (X (I) + X (J) ) / 2. 207 TEMP = H * XL * TEMP / 2. H = H * XL / 4. B (I) = B (I) + TEMP B (J) = B (J) + TEMP A (I, 1) = A (I, 1) + H A (J, 1) = A (J, 1) + H K = J - I + 1 IF (K) 212, 212, 210 210 A (I, K) = A (I, K) + H GOTO 215 212 K = I - J + 1 A (J, K) = A (J, K) + H 215 END DO 220 CONTINUE ! ! 2. TEMPERATURE BOUNDARY CONDITIONS ! DO 300 N = 1, NUMNP B (N) = B (N) + T (N) IF (KODE (N) ) 225, 300, 225 225 DO 250 M = 2, MBAND K = N - M + 1 IF (K) 235, 235, 230 230 B (K) = B (K) - A (K, M) * T (N) A (K, M) = 0.0 235 L = N + M - 1 IF (NUMNP - L) 245, 240, 240 240 B (L) = B (L) - A (N, M) * T (N) 245 A (N, M) = 0.0 250 END DO A (N, 1) = 1.0 B (N) = T (N) 300 END DO !****************************************************************** ! SOLVE FOR NODAL POINT TEMPERATURES !****************************************************************** ! FORM EFFECTIVE CONDUCTIVITY MATRIX FOR TIME INCREMENT ! IF (DT) 700, 700, 304 304 DT2 = 1.0 / DT DO 320 N = 1, NUMNP T (N) = 0.0 IF (KODE (N) ) 320, 305, 320 305 IF (D (N) ) 310, 320, 310 310 D (N) = DT2 * D (N) A (N, 1) = A (N, 1) + D (N) 320 END DO CALL SYMSOL (1) ! ! CALCULATE TEMPERATURE AT THE END OF EACH TIME INCREMENT ! TIME = 0.0 LL = 0 ! DO 600 L = 1, NDT ! ! 1. CALCULATE EFFECTIVE LOAD MATRIX ! DO 400 I = 1, NUMNP Q (I) = B (I) IF (KODE (I) ) 400, 395, 400 395 Q (I) = B (I) + D (I) * T (I) 400 END DO ! ! 2. SOLVE FOR TEMPERATURES ! CALL SYMSOL (2) ! DO 500 I = 1, NUMNP T (I) = Q (I) 500 END DO ! TIME = TIME+DT LL = LL + 1 IF (LL - INTER) 600, 550, 550 550 WRITE (6, 2005) TIME, (N, T (N), N = 1, NUMNP) 2005 FORMAT ('TIME=',E12.5/ (6(I6,E14.6)) ) LL = 0 ! 600 END DO GOTO 50 ! ! STEADY STATE SOLUTION ! 700 DO 800 I = 1, NUMNP Q (I) = B (I) 800 END DO CALL SYMSOL (1) CALL SYMSOL (2) WRITE (6, 2005) TIME, (N, Q (N), N = 1, NUMNP) GOTO 50 END SUBROUTINE SYMSOL (KKK) ! COMMON / SYMARG / NN, MM, A (500, 27), B (500) ! GOTO (1000, 2000), KKK ! ! REDUCE MATRIX ! 1000 DO 280 N = 1, NN DO 261 L = 2, MM C = A (N, L) / A (N, 1) I = N + L - 1 IF (NN - I) 260, 240, 240 240 J = 0 DO 250 K = L, MM J = J + 1 A (I, J) = A (I, J) - C * A (N, K) 250 END DO 260 A (N, L) = C 261 END DO 280 END DO GOTO 500 ! ! REDUCE VECTOR ! 2000 DO 291 N = 1, NN DO 286 L = 2, MM I = N + L - 1 IF (NN - I) 290, 285, 285 285 B (I) = B (I) - A (N, L) * B (N) 286 END DO 290 B (N) = B (N) / A (N, 1) 291 END DO ! ! BACK SUBSTITUTION ! N = NN 300 N = N - 1 IF (N) 350, 500, 350 350 DO 400 K = 2, MM L = N + K - 1 IF (NN - L) 400, 370, 370 370 B (N) = B (N) - A (N, K) * B (L) 400 END DO GOTO 300 ! 500 RETURN ! END SUBROUTINE SYMSOL