PROGRAM HEAT1DIM
C======================================================================
C ONE DIMENSIONAL THERMAL DIFFUSION ONLY HEAT EQUATION
C VARIABLE DEFINITION:
C FK = K/(RHO*CP); TL = DOMAIN LENGTH; DX = LENGTH OF ELEMENT
C IALPHA = 1: EXPLICIT
C ALPHA = 0: IMPLICIT
C ALPHA = 0.5: CRANK-NICOLSON
C DEC 16, 1995 WRITTEN BY EIJI FUKUMORI
C======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=2, MXE=10, MXN=MXE+1,NBW=ND )
DIMENSION NODEX(MXE,ND),X(MXN),A(MXN,NBW),RHS(MXN),
* IBTYPE(2), BV(2), T(MXN)
C======================================================================
C------- PRESET VARIABLES
TL = 4.D0*DATAN(1.D0)
WRITE (*,*) ' DOAMIN LENGTH(TL) = ', TL
NE = 10
NE = NE/2*2
IF ( NE .LT. 2 ) NE = 2
NNODE = NE + 1
DX = TL / NE
WRITE (*,*) ' NUMBER OF ELEMENTS(NE) =', NE
C
WRITE (*,*) ' ONE DIMENSIONAL THERMAL DIFFUSION PROBLEM'
C======================================================================
C-------- TYPE-IN DATA FROM KEYBOARD
C----------- RECOMMENDED FK=1., ALPHA = 0.5, DT = 0.1D0
WRITE(*,100)
100 FORMAT (1X,' K(TRY 1. IF UNDECIDED) = ' $ )
READ (*,*) FK
IF ( FK .LE. 0.D0 ) FK = 1.
C======================================================================
WRITE (*,*) ' TYPE IN ALPHA'
WRITE (*,*) ' IF ALPHA = 1. THEN EXPLICIT'
WRITE (*,*) ' IF ALPHA = 0. THEN IMPLICIT'
WRITE (*,*) ' IF ALPHA = 0.5 THEN CRANK-NICOLSON'
WRITE(*,110)
110 FORMAT (1X,' ALPHA = ' $ )
READ (*,*) ALPHA
IF ( ALPHA .LT. 0.) ALPHA = 0.
IF ( ALPHA .GT. 1.) ALPHA = 1.
C======================================================================
WRITE(*,120)
120 FORMAT (1X,' DT(TRY 0.001) = ' $ )
READ (*,*) DT
IF ( DT .LT. 0. ) DT = 0.0001
C======================================================================
WRITE (*,130)
130 FORMAT ( 1X, ' NUMBER OF INCREMENTS IN TIME(NSTEP) = ' $ )
READ (*,*) NSTEP
IF ( NSTEP .LT. 10 ) NSTEP = 10
C======================================================================
C-------- ELEMENT
DO I = 1 , NE
NODEX(I,1) = I
NODEX(I,2) = I + 1
END DO
C-------- COORDINATE
X(1) = 0.
DO I = 2 , NNODE
X(I) = (I-1)*DX
END DO
C-------- BOUNDARY CONDITION
IBTYPE(1) = 1
BV(1) = 0.
IBTYPE(2) = 1
BV(2) = 0.
C-------- INITIAL CONDITION
TIME = 0.D0
DO I = 1 , NNODE
T(I) = DSIN(X(I))
END DO
C-------- PRINTING INITIAL VALUE AT A SELECTED COORDINATE
C------- NODESLT = SELECTED NODE FOR COMPARISON BETWEEN THEM
NODESLT = NNODE/2+1
EXACT = DEXP(-FK*TIME)*DSIN(X(NODESLT))
WRITE (*,210)
210 FORMAT (1X,3X,'CURRENT TIME',1X,'T(X=TL/2,TIME)',10X,'EXACT',
* 5X,'DIFFERENCE' )
WRITE (*,200) TIME, T(NODESLT), EXACT
200 FORMAT ( 1X, 4F15.10 )
C======================================================================
C
C --------TIME MARCING START HERE--------
C
C======================================================================
DO ITIME = 1 , NSTEP
CALL MATRIX (MXE,MXN,ND,NBW,NE,NNODE,NODEX,X,A,RHS,DT,FK,ALPHA,T)
CALL FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
CALL SYSTEM ( MXN, NBW, NNODE, A, RHS )
DO I = 1 , NNODE
T(I) = RHS(I)
END DO
TIME = TIME + DT
NODESLT = NNODE/2+1
C------- NODESLT = SELECTED NODE FOR COMPARISON BETWEEN THEM
EXACT = DEXP(-FK*TIME)*DSIN(X(NODESLT))
WRITE (*,200) TIME, T(NODESLT), EXACT, RHS(NODESLT)-EXACT
END DO
STOP' NORMAL TERMINATION'
END
C
C
SUBROUTINE INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),IBTYPE(2), BV(2)
OPEN ( 1,FILE='BUCKLE.DAT', STATUS='OLD')
READ(1,*) P
READ(1,*) NE
DO I = 1 , NE
READ(1,*) IEL, (NODEX(IEL,J),J=1,ND),EI(IEL)
END DO
NNODE = NE + 1
DO I = 1 , NNODE
READ(1,*) NODE, X(NODE)
END DO
READ(1,*) IBTYPE(1), BV(1)
READ(1,*) IBTYPE(2), BV(2)
CLOSE (1)
RETURN
END
C
C
SUBROUTINE MATRIX (MXE,MXN,ND,NBW,NE,NNODE,NODEX,X,A,RHS,
* DT, FK, ALPHA, T )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),X(MXN),A(MXN,NBW),RHS(MXN),T(MXN)
DO I = 1 , NNODE
A(I,1) = 0.
A(I,2) = 0.
RHS(I) = 0.
END DO
DO IEL = 1 , NE
I = NODEX(IEL,1)
J = NODEX(IEL,2)
RL = X(J) - X(I)
A(I,1) = A(I,1) + RL/(2.*DT)+(1.-ALPHA)*FK/RL
A(I,2) = A(I,2) - (1.-ALPHA)*FK/RL
A(J,1) = A(J,1) + RL/(2.*DT)+(1.-ALPHA)*FK/RL
RHS(I)=RHS(I)+(RL/(2.*DT)-ALPHA*FK/RL)*T(I)+ALPHA*FK/RL*T(J)
RHS(J)=RHS(J)+(RL/(2.*DT)-ALPHA*FK/RL)*T(J)+ALPHA*FK/RL*T(I)
END DO
RETURN
END
C
C
SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION IBTYPE(2),BV(2),A(MXN,NBW), RHS(MXN)
IF ( IBTYPE(1) .EQ. 1 ) THEN
A(1,1) = 1.
RHS(1) = BV(1)
RHS(2) = RHS(2) - BV(1)*A(1,2)
A(1,2) = 0.
ELSE
RHS(1) = RHS(1) - BV(1)
END IF
IF ( IBTYPE(2) .EQ. 1 ) THEN
A(NNODE,1) = 1.
RHS(NNODE) = BV(2)
RHS(NNODE-1) = RHS(NNODE-1) - BV(2)*A(NNODE-1,2)
A(NNODE-1,2) = 0.
ELSE
RHS(NNODE) = RHS(NNODE) - BV(2)
END IF
RETURN
END
C
C
SUBROUTINE SYSTEM ( MXN, NBW, NNODE, A , B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,NBW) , B(MXN)
B(1) = B(1) / A(1,1)
A(1,1) = A(1,2) / A(1,1)
DO I = 2 , NNODE
P = A(I,1) - A(I-1,2) * A(I-1,1)
A(I,1) = A(I,2) / P
B(I) = ( B(I) - A(I-1,2)*B(I-1) ) / P
END DO
C------ BACK SUBSTITUTION ----
DO I = NNODE-1, 1,-1
B(I) = B(I) - A(I,1) * B(I+1)
END DO
RETURN
END