C WILSON'S HEAT CONDUCTION PROGAM (1966) C FOR PLANE AND AXISYMMETRIC STRUCTURES C (THE FIRST NON-STRUCTURAL FINITE ELEMENT PROGRAM) C COMMON B(500),X(500),Y(500),T(500),D(500),KODE(500), 1 HED(20),LM(5),IX(3),E(3,3),KX(4),P(5),S(5,5),DD(5), 2 XCOND(20),SPHT(20),DENS(20),QX(20) COMMON /SYMARG/ NUMNP,MBAND,A(500,27),Q(500) C*********************************************************************** C READ AND PRINT OF CONTROL INFORMATION C*********************************************************************** 50 READ (5,1000) HED,NUMNP,NUMEL,NUMCBC,KAT,NUMMAT,NDT,INTER,DT IF (KAT) 54,56,54 54 WRITE (6,2010) GO TO 58 56 WRITE (6,2011) 58 WRITE (6,2000) HED,NUMNP,NUMEL,NUMCBC,NUMMAT,NDT,INTER,DT C READ (5,1003) (M,XCOND(M),SPHT(M),DENS(M),QX(M),N=1,NUMMAT) WRITE(6,2009) (M,XCOND(M),SPHT(M),DENS(M),QX(M),M=1,NUMMAT) C*********************************************************************** C READ OR GENERATE NODAL POINT INFORMATION C*********************************************************************** WRITE (6,2001) L=1 C 60 READ (5,1001) N,KODE(N),X(N),Y(N),T(N) DIFF=N+1-L IF (N-L) 65,80,70 65 WRITE (6,2020) N GO TO 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) L=L+1 IF (N-L) 90,80,75 90 IF (NUMNP+1-L) 100,100,60 100 CONTINUE C*********************************************************************** C FORM CONDUCTIVITY MATRIX FOR COMPLETE BODY C*********************************************************************** DO 111 I=1,NUMNP D(I)=0.0 B(I)=0.0 DO 110 J=1,20 A(I,J)=0.0 110 CONTINUE 111 CONTINUE MBAND=0 NUM=0 WRITE (6,2003) C DO 200 N=1,NUMEL C C 1. READ OR GENERATE ELEMENT PROPERTIES C IF (NUM-N) 120,121,121 120 READ (5,1002) NUM,KX(1),KX(2),KX(3),KX(4),MTYPE C 121 DO 122 I=1,4 LM(I)=LM(I)+1 122 CONTINUE C IF (NUM-N) 123,124,126 123 WRITE (6,2021) NUM GO TO 120 C 124 DO 125 I=1,4 LM(I)=KX(I) 125 CONTINUE COND=XCOND(MTYPE) C 126 WRITE (6,2004) N,LM(1),LM(2),LM(3),LM(4),MTYPE C C 2. FORM ELEMENT CONDUCTIVITY MATRIX C DO 149 I=1,5 DD(I)=0.0 P(I)=0.0 DO 150 J=1,5 S(I,J)=0.0 150 CONTINUE 149 CONTINUE C I=LM(1) J=LM(2) K=LM(3) L=LM(4) LM(5)=I C XX=(X(I)+X(J)+X(K)+X(L))/4. YY=(Y(I)+Y(J)+Y(K)+Y(L))/4. C DO 152 K=1,4 C 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 C XMUL=1.0 IF (KAT) 136,137,136 136 XMUL=XMUL*(X(I)+X(J)+XX)/3.0 C 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 C 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 C IX(1)=K IX(2)=K+1 IF (K-4) 145,140,145 140 IX(2)=1 145 IX(3)=5 C 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 CONTINUE 148 CONTINUE C 152 CONTINUE C 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 CONTINUE 154 CONTINUE C C 3. ADD ELEMENT CONDUCTIVITY TO COMPLETE CONDUCTIVITY MATRIX C 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 CONTINUE 174 CONTINUE C 200 CONTINUE C*********************************************************************** C BOUNDARY CONDITIONS C*********************************************************************** C C 1. CONVECTION BOUNDARY CONDITIONS C IF (NUMCBC) 220,220,205 205 WRITE (6,2006) DO 215 N=1,NUMCBC WRITE(6,2004) NUMNP,NUMCBC,N READ (5,1007) I,J,H,TEMP WRITE (6,2007) I,J,H,TEMP 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 GO TO 215 212 K=I-J+1 A(J,K)=A(J,K)+H 215 CONTINUE 220 CONTINUE C C 2. TEMPERATURE BOUNDARY CONDITIONS C 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 CONTINUE A(N,1)=1.0 B(N)=T(N) 300 CONTINUE C*********************************************************************** C SOLVE FOR NODAL POINT TEMPERATURES C*********************************************************************** C FORM EFFECTIVE CONDUCTIVITY MATRIX FOR TIME INCREMENT C 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 CONTINUE CALL SYMSOL(1) C C CALCULATE TEMPERATURE AT THE END OF EACH TIME INCREMENT C TIME=0.0 LL=0 C DO 600 L=1,NDT C C 1. CALCULATE EFFECTIVE LOAD MATRIX C 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 CONTINUE C C 2. SOLVE FOR TEMPERATURES C CALL SYMSOL(2) C DO 500 I=1,NUMNP T(I)=Q(I) 500 CONTINUE C TIME=TIME+DT LL=LL+1 IF(LL-INTER) 600,550,550 550 WRITE (6,2005) TIME,(N,T(N),N=1,NUMNP) LL=0 C 600 CONTINUE GO TO 50 C C STEADY STATE SOLUTION C 700 DO 800 I=1,NUMNP Q(I)=B(I) 800 CONTINUE CALL SYMSOL(1) CALL SYMSOL(2) WRITE (6,2005) TIME,(N,Q(N),N=1,NUMNP) GO TO 50 C C FORMAT STATEMENTS C 1000 FORMAT (20A4/7I5,E15.6) 1001 FORMAT (2I5,3F10.0) 1002 FORMAT (6I5) 1003 FORMAT (I10,4F10.0) 1007 FORMAT (2I5,2F10.0) C 2000 FORMAT (1H0 20A4// 25H0NUMBER OF NODAL POINTS-- I4/ 1 25H NUMBER OF ELEMENTS------ I4 / 25H NUMBER OF CONVECTION BC-I4/ 2 25H NUMBER OF MATERIALS----- I4 / 25H NUMBER OF INCREMENTS----I4/ 3 25H OUTPUT INTERVAL--------- I4 / 20H TIME INTERVAL------ E10.3) 2001 FORMAT (20H0 N.P. NO. CODE 14X,1HX,14X,1HY,14X,1HT) 2002 FORMAT (2I10,3E15.6) 2003 FORMAT (35H0 N I J K L MATERIAL ) 2004 FORMAT (5I5,I10) 2005 FORMAT (6H0TIME= E12.5/ (I6,E14.6,I6,E14.6,I6,E14.6,I6,E14.6, 1 I6,E14.6,I6,E14.6)) 2006 FORMAT (40H0 I J H TEMPERATURE ) 2007 FORMAT (2I5,2E15.6) 2009 FORMAT (4H0 M 14X 1HK 14X 1HC 14X 1HD 14X 1HQ/ (I4,4E15.6)) 2010 FORMAT (25H1AXISYMMETRIC SOLID BODY ) 2011 FORMAT (27H1TWO DIMENSIONAL PLANE BODY ) C 2020 FORMAT (10H0CARD NO. I4, 13H OUT OF ORDER ) 2021 FORMAT (13H0BAD CARD NO. I4) C END SUBROUTINE SYMSOL (KKK) C COMMON /SYMARG/ NN,MM,A(500,27),B(500) C GO TO (1000,2000),KKK C C REDUCE MATRIX C 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 CONTINUE 260 A(N,L)=C 261 CONTINUE 280 CONTINUE GO TO 500 C C REDUCE VECTOR C 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 CONTINUE 290 B(N)=B(N)/A(N,1) 291 CONTINUE C C BACK SUBSTITUTION C 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 CONTINUE GO TO 300 C 500 RETURN C END