PROGRAM ONEDEE
C=======================================================================
C ONE DIMENSIONAL TIME DEPENDENT CONVECTION-DIFFUSION EQUATION
C METHOD: FINITE ELEMENT METHOD
C ELEMENT: LINEAR
C EIJI FUKUMORI SEPTEMBER 1985
C=======================================================================
PARAMETER ( IR = 1, IW = 2, MXE=2000, MXN=MXE+1 )
IMPLICIT REAL * 8 ( A-H, O-Z )
DIMENSION X(MXN), RHO(MXN), V(MXN), DB(MXN), S(MXN), C(MXN),
* A(MXN,3), G(MXN,3), H(MXN,3)
C=======================================================================
C---------- READ IN DATA ---------
OPEN ( IR, FILE ='ONEDEE.DAT', STATUS='OLD')
READ (IR,*) NNODE
READ (IR,*) DT
READ (IR,*) TMAX
READ (IR,*) IB1, BV1
READ (IR,*) IBN, BVN
C
IF ( TMAX .LT. DT ) STOP'ERROR TMAX < DT'
IF ( DT .LE. 0. ) STOP'ERROR DT <= 0.'
DO I = 1 , NNODE
READ (IR,*) X(I),RHO(I),V(I),DB(I)
END DO
READ (IR,*) ( S(I) , I = 1 , NNODE )
CLOSE (IR)
C=======================================================================
C-------- RESET MAIRICES [G] AND [H]
DO J = 1 , 3
DO I = 1 , NNODE
G(I,J) = 0.D0
H(I,J) = 0.D0
END DO
END DO
C=======================================================================
C-------- EVALUATION OF FINITE ELEMENT EQUATION
NE = NNODE -1
DO I = 1 , NE
J = I + 1
RL = X(I+1) - X(I)
DR = RHO(J) - RHO(I)
DDB = DB(J) - DB(I)
AVERHO = 0.5 * ( RHO(I) + RHO(J) )
AVEDB = 0.5 * ( DB(I) + DB(J) )
AVEV = 0.5 * ( V(I) + V(I) )
C------ ADDITIVE VISCOSITY
FKX = 1.00 * AVEV * RL / 2.
GAMMAX = FKX / AVEDB
ADDVX = OPT(GAMMAX) * FKX
C
E1 = 1./RL * ( ADDVX + RHO(I)*DB(I) + 0.5 * ( DB(1) * DR +
* RHO(I) * DDB ) + 1./3. * DR * DDB )
E2 = 0.5 * ( AVEDB * DR )
E3 = 0.5 * AVERHO * AVEV
E4 = (1./DT) * ( RL/2.) * AVERHO
C
B11 = - E3 + E1 - E2
B12 = E3 - E1 - E2
B21 = - E3 - E1 + E2
B22 = E3 + E1 + E2
C
G11 = E4 + 0.5 * B11
G12 = 0.5 * B12
G21 = 0.5 * B21
G22 = E4 + 0.5 * B22
C
H11 = E4 - 0.5 * B11
H12 = - 0.5 * B12
H21 = - 0.5 * B21
H22 = E4 - 0.5 * B22
C
C----- ASSEMBLY OF MATRICES, [G] AND [H]
C [G]{S} = [H](S) : {S} = {S} AT T+DT
C (S) = {S} AT T
G(I,2) = G(I,2) + G11
G(I,3) = G(I,3) + G12
G(J,1) = G(J,1) + G21
G(J,2) = G(J,2) + G22
H(I,2) = H(I,2) + H11
H(I,3) = H(I,3) + H12
H(J,1) = H(J,1) + H21
H(J,2) = H(J,2) + H22
END DO
C=======================================================================
OPEN ( IW, FILE ='SOLUTION.FEM', STATUS='UNKNOWN')
C=======================================================================
ITIME = 0
CALL PRINT ( DT, IW, MXN, NNODE, ITIME, X , S )
C
C-------- TIME MARCHING PROCESS STARTS HERE.
NTIME = TMAX/DT + 0.1
DO ITIME = 1 , NTIME
TIME = DT * ITIME
C------ RESET RIGHT HAND SIDE OF SYSTEM EQUATION.
DO I = 1 , NNODE
C (I) = 0.
END DO
C------ COPY MATRIX [G] TO [A]
DO J = 1 , 3
DO I = 1 , NNODE
A(I,J) = G(I,J)
END DO
END DO
C------ COMPUTATION OF RIGHT RAND SIDE
C(1) = H(1,2) * S(1) + H(1,3) * S(2)
C(NNODE) = H(NNODE,2)* S(NNODE) + H(NNODE,1) * S(NNODE-1)
DO I = 2 , NNODE-1
C(I) = S(I-1) * H(I,1) + S(I) * H(I,2) + S(I+1) * H(I,3)
END DO
C================= IMPLEMENTATION OF BOUNDARY CONDITIONS ===============
C IF IBX = 1, FIRST KIND
C IF IBX = 2, SECOND KIND
C
C--------- TOP BOUNDARY
IF ( IB1 .EQ. 1 ) THEN
A(1,1) = 0.D0
A(1,2) = 1.D0
A(1,3) = 0.D0
C(1) = BV1
ENDIF
IF ( IB1 .EQ. 2 ) THEN
C(1) = C(1) + BV1
ENDIF
C--------- BOTTOM
IF ( IBN .EQ. 1 ) THEN
A(NNODE,1) = 0.D0
A(NNODE,2) = 1.D0
A(NNODE,3) = 0.D0
C(NNODE) = BVN
ENDIF
IF ( IBN .EQ. 2 ) THEN
C(NNODE) = C(NNODE) - BVN
ENDIF
C====================== SOLVING SYSTEM EQUATION ========================
CALL SYSTEM ( MXN, NNODE, A, C )
DO I = 1 , NNODE
S(I) = C(I)
END DO
CALL PRINT ( DT, IW, MXN, NNODE, ITIME, X , S )
END DO
C=======================================================================
CLOSE (IW)
C=======================================================================
STOP
END
C
C
SUBROUTINE PRINT ( DT, IW, MXN, NNODE, ITIME, X , S )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN) , S(MXN)
DATA ISKIP / 40 /
IF ( ITIME .EQ. 0 ) THEN
WRITE (IW,*) NNODE
WRITE (IW,200) ( X(I) , I = 1 , NNODE )
END IF
IF ( ITIME/ISKIP*ISKIP - ITIME .NE. 0 ) RETURN
TIME = DT * ITIME
WRITE (IW,200) TIME
WRITE (IW,200) ( S(I) , I = 1 , NNODE )
200 FORMAT ( F10.3 )
RETURN
END
C
C
FUNCTION OPT ( GAMMA )
IMPLICIT REAL*8 ( A-H , O-Z )
IF ( GAMMA .LT. 0.01 ) THEN
OPT = GAMMA / 3.
RETURN
END IF
IF ( GAMMA .GT. 10. ) THEN
OPT = 1.
RETURN
END IF
E = DEXP ( GAMMA )
COTH = ( E + 1./E ) / ( E - 1./E )
OPT = COTH - 1./ GAMMA
RETURN
END
C
SUBROUTINE SYSTEM ( MXN, NNODE, A , B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,3) , B(MXN)
B(1) = B(1) / A(1,2)
A(1,2) = A(1,3) / A(1,2)
DO I = 2 , NNODE
P = A(I,2) - A(I,1) * A(I-1,2)
A(I,2) = A(I,3) / P
B(I) = ( B(I) - A(I,1)*B(I-1) ) / P
END DO
C------ BACK SUBSTITUTION ----
DO I = NNODE-1, 1, -1
B(I) = B(I) - A(I,2) * B(I+1)
END DO
RETURN
END