*DECK HEAT18
      SUBROUTINE INITAL
C...
C...  DYNAMICS OF A REGENERATOR WITH PLOTTING OF THE SPATIAL PROFILES
C...
C...  THE PARTIAL DIFFERENTIAL EQUATIONS DESCRIBING THE DYNAMICS OF A
C...  DISTRIBUTED REGENERATOR ARE
C...
C...     TF  = -Q/(S*F)TF  - (AV*H)/(F*RHOF*CPF)(TF - TS)            (1)
C...       T             Z
C...
C...     TS  = (AV*H)/((1-F)*RHOS*CPS)(TF - TS)                      (2)
C...       T
C...
C...  WITH AUXILIARY CONDITIONS
C...
C...     TF(0,T) = F(T), TF(Z,0) = TF0, TS(Z,0) = TS0          (3)(4)(5)
C...
C...  WHERE
C...
C...     TF        FLUID TEMPERATURE
C...
C...     TS        SOLID TEMPERATURE
C...
C...     Z         AXIAL POSITION
C...
C...     T         TIME
C...
C...     Q         FLUID VOLUMETRIC FLOW RATE
C...
C...     S         REGENERATOR CROSS SECTIONAL AREA
C...
C...     F         SOLID VOID FRACTION
C...
C...     AV        HEAT TRANSFER AREA BETWEEN THE FLUID AND SOLID PER
C...               UNIT VOLUME OF REGENERATOR
C...
C...     H         HEAT TRANSFER COEFFICIENT
C...
C...     RHOF      FLUID DENSITY
C...
C...     RHOS      SOLID DENSITY
C...
C...     CPF       FLUID HEAT CAPACITY
C...
C...     CPS       SOLID HEAT CAPACITY
C...
C...     TF0       INITIAL FLUID TEMPERATURE
C...
C...     TS0       INITIAL SOLID TEMPERATURE
C...
C...     F(T)      ENTERING FLUID TEMPERATURE AS A FUNCTION OF TIME (FOR
C...               THE FOLLOWING ANALYSIS, F(T) = TFE, A CONSTANT INLET
C...               TEMPERATURE)
C...
C...  TO SIMPLIFY THE ANALYSIS, THE FOLLOWING CONSTANTS ARE DEFINED
C...
C...     A = Q/(S*F*DZ), B = (AV*H)/(F*RHOF*CPF)
C...
C...     C = (AV*H)/((1-F)*RHOS*CPS)
C...
C...  NOTE THAT A*DZ IS THE FLUID VELOCITY, AND 1/B AND 1/C CAN BE 
C...  CONSIDERED HEAT TRANSFER TIME CONSTANTS FOR THE FLUID AND SOLID
C...  RESPECTIVELY.  ALSO, TYPICALLY, B/C IS LARGE, I.E., THE SOLID
C...  PDE IS MUCH SLOWER THAN THE GAS PDE OR IN OTHER WORDS, THE
C...  MODEL IS A STIFF SET OF PDES.
C...
C...  THESE CONCEPTS ARE ILLUSTRATED WITH TWO RUNS OF THE MODEL: (1)
C...  HAV = 0 FOR WHICH ONLY GAS FLOW TAKES PLACE (A FAST PROCESS)
C..   AND (2) HAV NE 0 FOR WHICH HEAT TRANSFER TO THE SOLID TAKES
C...  PLACE (A SLOW PROCESS).  FOR THE LATTER, THE TIME SCALE IS 25000
C...  TIMES LONGER (NOTE THE TIME SCALES DEFINED IN THE DATA).  ALSO,
C...  BECAUSE OF THE STIFFNESS, THE TIME INTEGRATION SHOULD BE DONE
C...  WITH AN IMPLICIT INTEGRATOR (IT IS INSTRUCTIVE TO RUN THIS PRO-
C...  BLEM WITH AN EXPLICIT INTEGRATOR AND COMPARE THE RUN TIME WITH
C...  THAT FOR AN IMPLICIT INTEGRATOR - THE LATTER IS MUCH SHORTER).
C...
C...  THE FLUID AND SOLID TEMPERATURES, TF(0,T), TS(0,T), TF(ZL,T) AND
C...  TS(ZL,T), ARE PLOTTED AS A FUNCTION OF TIME IN SUBROUTINE PRINT.
C...  ALSO, THE SPATIAL PROFILES, TF(Z,T) AND TS(Z,T) FOR FOUR SELECTED
C...  TIMES IN SUBROUTINE PRINT.
C...
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/    TF(21),    TS(21)
     2      /F/ DTFDT(21), DTSDT(21)
     3      /S/ DTFDZ(21)
     4      /C/         A,         B,         C,       TF0,       TS0,
     5                TFE,         N,        DZ,        ZL,        IP
C...
C...  THE REGENERATOR IS SIMULATED BY THE NUMERICAL METHOD OF LINES
C...  IN WHICH THE SPATIAL DERIVATIVE TF  IS REPLACED WITH A FIRST-
C...                                    Z
C...  ORDER, TWO-POINT DIFFERENCE.  THIS IS IN EFFECT SIMULATING
C...  THE REGENERATOR WITH 20 PERFECTLY MIXED SECTIONS.
C...
C...  SET THE MODEL PARAMETERS
      N=21
      ZL=2.
      DZ=0.1
      TF0=20.
      TS0=20.
      TFE=500.
      Q=1.0
      S=0.5
      F=0.3
      RHOF=0.15
      RHOS=3000.
      CPF=5000.
      CPS=5000.
C...
C...  NO HEAT TRANSFER (JUST FLOW OF THE GAS THROUGH THE REGENERATOR)
      IF(NORUN.EQ.1)HAV=0.0E+05
C...
C...  HEAT TRANSFER TO THE SOLID INCLUDED
      IF(NORUN.EQ.2)HAV=1.0E+05
C...
C...  CONSTANTS A, B, C
      A=Q/(S*F*DZ)
      B=HAV/(F*RHOF*CPF)
      C=HAV/((1.0-F)*RHOS*CPS)
C...
C...  SET THE MODEL INITIAL CONDITIONS
      DO 1 I=1,N
      TF(I)=TF0
      TS(I)=TS0
1     CONTINUE
C...
C...  INITIALIZE THE REST OF THE MODEL
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/    TF(21),    TS(21)
     2      /F/ DTFDT(21), DTSDT(21)
     3      /S/ DTFDZ(21)
     4      /C/         A,         B,         C,       TF0,       TS0,
     5                TFE,         N,        DZ,        ZL,        IP
C...
C...  IMPLEMENT THE BOUNDARY CONDITION
      TF(1)=TFE
      DTFDT(1)=0.
C...
C...  THE NMOL ODES ARE PROGRAMMED EXPLICITLY
1     DO 11 I=1,N
      DTSDT(I)=C*(TF(I)-TS(I))
      IF(I.EQ.1)GO TO 11
      DTFDT(I)=A*(TF(I-1)-TF(I))+B*(TS(I)-TF(I))
11    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/    TF(21),    TS(21)
     2      /F/ DTFDT(21), DTSDT(21)
     3      /S/ DTFDZ(21)
     4      /C/         A,         B,         C,       TF0,       TS0,
     5                TFE,         N,        DZ,        ZL,        IP
C...
C...  DIMENSION THE ARRAYS FOR THE TIME PLOTS
      DIMENSION TP(4,101), TIME(101)
C...
C...  DIMENSION THE ARRAYS FOR THE SPATIAL PROFILES
      DIMENSION  TZ(8,21),    ZP(21)
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION AND THE VELOCITY OF
C...  THE GAS, VGAS, THE RECIPROCALS OF THE TIME CONSTANTS FOR THE
C...  GAS AND SOLID, B AND C
      IP=IP+1
      VGAS=A*DZ
      IF(IP.EQ.1)WRITE(NO,1)VGAS,B,C
1     FORMAT('     VGAS = ',E10.3,' M/SEC  ',
     +       '   1/TGAS = ',E10.3,' 1/SEC  ',
     +       ' 1/TSOLID = ',E10.3,' 1/SEC  ',//,
     +  9X,1HT,5X,7HTF(0,T),5X,7HTS(0,T),4X,8HTF(ZL,T),4X,8HTS(ZL,T))
C...
C...  PRINT THE NUMERICAL SOLUTION AT Z = 0, ZL EVERY TENTH CALL TO 
C...  SUBROUTINE PRINT
      IF(((IP-1)/10*10).EQ.(IP-1))THEN
         WRITE(NO,2)T,TF(1),TS(1),TF(N),TS(N)
2        FORMAT(F10.2,4F12.2)
      END IF
C...
C...  STORE THE NUMERICAL SOLUTION FOR THE TIME PLOTS
      TP(1,IP)=TF(1)
      TP(2,IP)=TS(1)
      TP(3,IP)=TF(N)
      TP(4,IP)=TS(N)
      TIME(IP)=T
C...
C...  STORE THE NUMERICAL SOLUTION FOR THE SPATIAL PROFILES
C...  AT T = 0, 1/4, 2/4 AND 3/4 OF THE TOTAL RUN TIMES
      IF((IP-1)*(IP-26)*(IP-51)*(IP-76).NE.0)GO TO 5
      IF(IP.EQ.  1)NSP=1
      IF(IP.EQ. 26)NSP=3
      IF(IP.EQ. 51)NSP=5
      IF(IP.EQ. 76)NSP=7
      DO 4 IZ=1,N
      TZ(NSP  ,IZ)=TF(IZ)
      TZ(NSP+1,IZ)=TS(IZ)
      ZP(IZ)=FLOAT(IZ-1)*DZ
4     CONTINUE
C...
C...  PLOT THE TIME RESPONSE
5     IF(IP.LT.101)RETURN
      CALL TPLOTS(4,IP,TIME,TP)
C...
C...  LABEL THE TIME PLOT
      WRITE(NO,3)
3     FORMAT(1H ,//,53H 1 - TF(0,T), 2 - TS(0,T), 3 - TF(ZL,T), 4 - TS(Z
     1L,T))
C...
C...  PLOT THE SPATIAL PROFILES
      CALL TPLOTS(8,N,ZP,TZ)
C...
C...  LABEL THE SPATIAL PROFILE PLOT
      WRITE(NO,6)
6     FORMAT(' 1 - TF(Z,  0)   2 - TS(Z,  0)',/,
     1       ' 3 - TF(Z,1/4)   4 - TS(Z,1/4)',/,
     2       ' 5 - TF(Z,2/4)   6 - TS(Z,2/4)',/,
     3       ' 7 - TF(Z,3/4)   8 - TS(Z,3/4)')
      RETURN
      END

