!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 6.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Please Note: ! ! ! ! (1) This computer program is written by Tao Pang in conjunction with ! ! his book, "An Introduction to Computational Physics," published ! ! by Cambridge University Press in 1997. ! ! ! ! (2) No warranties, express or implied, are made for this program. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! PROGRAM N_T_FIELD ! ! This program solves the time dependent temperature field ! around a nuclear waste rod in a two-dimensional model. ! Copyright (c) Tao Pang 1997. ! IMPLICIT NONE INTEGER, PARAMETER :: M1=1001,M=2001,N=19 INTEGER :: I,J REAL :: HT,TC,XK,RA,RB,H,H2,T0,S0,CS,A0,B0,C0,R,T REAL, DIMENSION (N) :: A,B,C,Y,G,W,V,U REAL, DIMENSION (M,N) :: S REAL, DIMENSION (M,N+1) :: X ! HT = 1.0/(M1-1) TC = 1.0 XK = 3153.6 RA = 25.0 RB = 100.0 H = RB/(N+1) H2 = H*H T0 = 10.0 S0 = HT*XK*T0/RA**2 CS = HT*XK/H2 ! DO I = 1, M DO J = 1, N+1 X(I,J) = 0.0 END DO END DO ! DO I = 1, N R = H*I A(I) = 2.0*(1.0+CS)*R B(I) = -(1.0+0.5/I)*CS*R C(I) = -(1.0-0.5/I)*CS*R IF (R.LT.RA) THEN S(1,I) = S0*R ELSE S(1,I) = 0.0 END IF END DO ! ! Assign the source of the radiation heat ! DO I = 2, M T = HT*(I-1) DO J = 1, N R = H*J IF (R.LT.RA) THEN S(I,J) = R*S0/EXP(T/TC) ELSE S(I,J) = 0.0 END IF END DO END DO ! DO I = 2, M ! ! Find the values from the last time step ! A0 = 2.0*(1.0-CS)*H B0 = (1.0+0.5)*CS*H C0 = (1.0-0.5)*CS*H G(1) = A0*X(I-1,1)+B0*X(I-1,2) & +C0*(4.0*X(I-1,1)-X(I-1,2))/3.0+S(I-1,1)+S(I,1) DO J = 2, N R = J*H A0 = 2.0*(1.0-CS)*R B0 = (1.0+0.5/J)*CS*R C0 = (1.0-0.5/J)*CS*R G(J) = A0*X(I-1,J)+B0*X(I-1,J+1)+C0*X(I-1,J-1)+(S(I-1,J)+S(I,J)) END DO ! ! Find the elements in L and U ! W(1) = A(1) V(1) = C(1) U(1) = B(1)/W(1) DO J = 2, N V(J) = C(J) W(J) = A(J)-V(J-1)*U(J-1) U(J) = B(J)/W(J) END DO ! ! Find the solution of the temperature ! Y(1) = G(1)/W(1) DO J = 2, N Y(J) = (G(J)-V(J-1)*Y(J-1))/W(J) END DO ! X(I,N) = Y(N) DO J = N-1,1,-1 X(I,J) = Y(J)-U(J)*X(I,J+1) END DO END DO ! I = M DO J = 1, N+1 R = J*H WRITE (6,"(2F16.8)") R,X(I,J) END DO END PROGRAM N_T_FIELD