      SUBROUTINE INITAL
C...
C...  CONVECTIVE COOLING OF A MOVING POLYMER SHEET
C...
C...  CONSIDER THE CONVECTIVE COOLING OF A MOVING POLYMER SHEET ILLUS-
C...  TRATED BELOW
C                             +----- DZ -----+
C        ...........................................................
C            +                .              .
C            .                .              .
C            D    V, TP  ----+.    TP(Z,T)   . ----+ V, TP
C            .         AT Z   .              .            AT Z+DZ
C            +                .       .      .
C        ...........................................................
C      Z = 0                  Z       .     Z+DZ                 Z = ZL
C                                     .
C                                     +
C
C                                    TA
C...
C...  AN ENERGY BALANCE ON THE INCREMENTAL SECTION OF THE POLYMER SHEET
C...  GIVES
C...
C...     (DZ*W*D*RHO*CP*TP ) = DZ*W*V*RHO*CP*TP
C...                      T                    AT Z
C...
C...                         - DZ*W*V*RHO*CP*TP
C...                                           AT Z+DZ
C...
C...                         + 2*W*DZ*U*(TA - TP)
C...
C...  THE 2 IN THE HEAT TRANSFER TERM ACCOUNTS FOR COOLING ON BOTH
C...  FLAT SURFACES OF THE POLYMER.
C...
C...  DIVISION BY D*W*DZ*RHO*CP FOLLOWED BY DZ ---+ 0 GIVES
C...
C...     TP  = -V*TP  + 2*U*/(D*CP*RHO)*(TA - TP)                  (1)
C...       T        Z
C...
C...  THE BOUNDARY AND INITIAL CONDITIONS FOR EQUATION (1) ARE
C...
C...     TP(0,T) = T0 + (TE - T0)*(1 - EXP(-T/2.0))                (2)
C...
C...     TP(Z,0) = T0                                              (3)
C...
C...  WHERE
C...
C...     TP        POLYMER TEMPERATURE
C...
C...     TA        AMBIENT TEMPERATURE
C...
C...     T         TIME
C...
C...     Z         POSITION ALONG THE POLYMER
C...
C...     V         VELOCITY OF THE POLYMER SHEET
C...
C...     D         THICKNESS OF THE POLYMER SHEET
C...
C...     W         WIDTH OF THE POLYMER SHEET
C...
C...     U         HEAT TRANSFER COEFFICIENT BETWEEN THE FLAT SURFACES
C...               OF THE POLYMER SHEET AND THE SURROUNDING COOLING
C...               MEDIUM
C...
C...     CP        HEAT CAPACITY OF THE POLYMER
C...
C...     RHO       DENSITY OF THE POLYMER
C...
C...  THE PROBLEM, THEN, IS TO COMPUTE TP(Z,T) TO DETERMINE IF THE
C...  POLYMER IS COOLED SUFFICIENTLY AT Z = ZL WHERE ZL IS THE LENGTH
C...  OF THE COOLING SECTION.
C...
C...  THE NUMERICAL VALUES OF THE MODEL PARAMETERS ARE
C...
C...     V = 10 CM/SEC
C...
C...     D = 0.5 CM
C...
C...     CP = 0.8 CAL/GM-C
C...
C...     RHO = 1.2 GM/CM**3
C...
C...     U = 0.1 CAL/SEC-CM**2-C
C...
C...     TE = 400 C
C...
C...     T0 = 25 C
C...
C...     TA = 25
C...
C...     ZL = 100 CM
C...
C...  THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS ARE COMPARED
C...  IN SUBROUTINE PRINT.  THE STEADY STATE ANALYTICAL SOLUTION TO
C...  EQUATIONS (1) TO (3) CAN BE OBTAINED IN THE FOLLOWING WAY.
C...
C...  AT STEADY STATE, EQUATION (1) REDUCES TO
C...
C...     V*DTP/DZ = 2*U/(D*RHO*CP)*(TA - TP)
C...
C...  WHICH CAN BE REARRANGED TO
C...
C...     DTP/(TA - TP) = E*DZ, E = 2*U*/(D*RHO*CP*V)
C...
C...  INTEGRATION OF BOTH SIDES GIVES
C...
C...       TP                     Z
C...      INT DTP/(TP - TA) = -E*INT DX
C...       TE                     0
C...
C...  EVALUATION OF THE INTEGRALS AT THE LIMITS GIVES
C...
C...    LN((TP - TA)/(TE - TA)) = -E*Z
C...
C...  OR
C...
C...     (TP - TA)/(TE - TA) = EXP(-E*Z)
C...
C...  OR
C...
C...     TP = TA + (TE - TA)*EXP(-E*Z)
C...
C...  THIS STEADY STATE ANALYTICAL SOLUTION CAN BE COMPARED WITH THE
C...  STEADY STATE NUMERICAL SOLUTION.
C...
C...  FOR THE SPECIAL CASE U = 0, E = 0 IN THE PRECEDING SOLUTION WHICH
C...  REDUCES TO
C...
C...     TP = TE
C...
C...  AS EXPECTED.
C...
C...  THE NUMERICAL SOLUTION IS PROGRAMMED IN SUBROUTINE DERV FOR TWO
C...  APPROXIMATIONS OF THE SPATIAL DERIVATIVE, TP , IN EQUATION (1)
C...                                              Z
C...     (1)  A FIVE-POINT CENTERED APPROXIMATION USED IN SUBROUTINE
C...          DSS004.
C...
C...     (2)  A FIVE-POINT BIASED UPWIND APPROXIMATION USED IN SUBROU-
C...          TINE DSS020.
C...
C...  AN INTERESTING EXPERIMENT WITH THIS PROGRAM IS TO SET THE HEAT
C...  TRANSFER COEFFICIENT, U, TO ZERO.  FOR THIS CASE, EQUATION (1)
C...  REDUCES TO THE ADVECTION EQUATION, WHICH, SUPERFICIALLY, IS
C...  QUITE SIMPLE, BUT ACTUALLY IS DIFFICULT TO SOLVE NUMERICALLY
C...  BECAUSE IT PROPAGATES DISCONTINUITIES.  THESE DISCONTINUITIES
C...  PRODUCE DIFFUSION AND OSCILLATION WHICH ARE EVIDENT IN THE
C...  NUMERICAL SOLUTIONS.
C...
C...  IF FUNCTIONS OF T LESS STRINGENT THAN DISCONTINUITIES ARE APPLIED
C...  AS BOUNDARY CONDITIONS (AT Z = 0), THE FIVE-POINT BIASED UPWIND
C...  APPROXIMATION WILL GENERALLY GIVE BETTER SOLUTIONS THAN THE FIVE-
C...  POINT CENTERED APPROXIMATION.  THIS POINT HAS BEEN STUDIED EXTEN-
C...  SIVELY, AND A TUTORIAL PAPER IS AVAILABLE FROM W. E. SCHIESSER.
C...
C...  ONE EASY WAY TO REMOVE DISCONTINUITIES FROM THE PROBLEM (EQUATIONS
C...  (1) TO (3)) IS TO SET T0 = TE, FOR WHICH THE FOLLOWING PROGRAM
C...  PRODUCES A SMOOTH SOLUTION EVEN FOR U = 0.  OTHERWISE (T0 NE TE),
C...
C...  FOR U = 0, THE PROBLEM IS ESSENTIALLY IMPOSSIBLE SINCE TP  IN EQ-
C...                                                           Z
C...  UATION (1) IS INFINITE.  HOWEVER, EVEN FOR THIS CASE, A NUMERICAL
C...  SOLUTION IS COMPUTED BY THE PROGRAM WHICH CAN BE COMPARED WITH THE
C...  EXACT SOLUTION, I.E., FOR U = 0, EQUATION (1) BECOMES THE ADVEC-
C...  TION EQUATION,
C...
C...     TP  + V*TP  = 0
C...       T       Z
C...
C...  WHICH, FOR THE INITIAL AND BOUNDARY CONDITIONS
C...
C...     TP(Z,0) = T0, TP(0,T) = T0 + (TE - T0)*H(T)
C...
C...  HAS THE EXACT SOLUTION
C...
C...     TP(Z,T) = T0 + (TE - T0)*H(T - Z/V)
C...
C...  WHERE H(T) IS THE HEAVISIDE UNIT STEP FUNCTION
C...
C...            0, T LT 0
C...
C...     H(T) =
C...
C...            1, T GT 0
C...
C...  THIS EXACT SOLUTION FOR U = 0 (A STEP FUNCTION TRAVELING AT
C...  THE VELOCITY V) CAN BE COMPARED WITH THE NUMERICAL SOLUTION.
C...  IN PARTICULAR, AT Z = ZL, THE EXACT SOLUTION IS JUST A STEP AT
C...  T = ZL/V FROM T0 TO TE, AND, OF COURSE, THE NUMERICAL SOLUTION
C...  WILL ONLY APPROXIMATE THIS EXACT SOLUTION.
C...
C...  ONE OTHER IMPORTANT PARAMETER (BESIDE THE HEAT TRANSFER COEFFI-
C...  CIENT, U) CAN BE STUDIED WITH THE FOLLOWING PROGRAM, THE NUMBER
C...  OF GRID POINTS N.  IN THE PRESENT CASE, 41 WAS SELECTED RATHER
C...  ARBITRARILY.
C...
      COMMON/T/          T,     NSTOP,     NORUN
     1      /Y/     TP(41)
     2      /F/    TPT(41)
     3      /S/    TPZ(41)
     4      /C/          V,         D,        CP,       RHO,         U,
     5                  TE,        T0,        TA,        ZL,         E,
     6                   N,        IP
C...
C...  SET THE MODEL PARAMETERS
      V=10.
      D=0.5
      CP=0.8
      RHO=1.2
      U=0.10
      TE=400.
      T0=25.
      TA=25.
      ZL=100.
      N=41
C...
C...  COMPUTE THE SPATIAL EIGENVALUE E
      E=2.*U/(D*RHO*CP*V)
C...
C...  MODEL INITIAL CONDITIONS
      DO 1 I=1,N
      TP(I)=T0
1     CONTINUE
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/          T,     NSTOP,     NORUN
     1      /Y/     TP(41)
     2      /F/    TPT(41)
     3      /S/    TPZ(41)
     4      /C/          V,         D,        CP,       RHO,         U,
     5                  TE,        T0,        TA,        ZL,         E,
     6                   N,        IP
C...
C...  BOUNDARY CONDITION
      TP(1)=T0+(TE-T0)*(1.0-EXP(-T/2.0))
      TPT(1)=(TE-T0)*(1.0/2.0)*EXP(-T/2.0)
C...
C...  COMPUTE THE SPATIAL DERIVATIVE BY FIVE-POINT CENTERED DIFFERENCES
      IF(NORUN.EQ.1)CALL DSS004(0.,ZL,N,TP,TPZ)
C...
C...  COMPUTE THE SPATIAL DERIVATIVE BY FIVE-POINT BIASED UPWIND DIFFER-
C...  ENCES
      IF(NORUN.EQ.2)CALL DSS020(0.,ZL,N,TP,TPZ,V)
C...
C...  PDE (1)
      DO 1 I=2,N
      TPT(I)=-V*TPZ(I)+2.*U/(D*RHO*CP)*(TA-TP(I))
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/          T,     NSTOP,     NORUN
     1      /Y/     TP(41)
     2      /F/    TPT(41)
     3      /S/    TPZ(41)
     4      /C/          V,         D,        CP,       RHO,         U,
     5                  TE,        T0,        TA,        ZL,         E,
     6                   N,        IP
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      REAL ZP(41), TPP(6,41)
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(1H ,//,9X,1HT,4X,8H TP(0,T),4X,8HTP(ZL,T))
C...
C...  PRINT THE SOLUTION (IN TIME)
      WRITE(NO,2)T,TP(1),TP(N)
2     FORMAT(F10.2,2F12.2)
C...
C...  STORE THE SOLUTION FOR PLOTTING (IN SPACE)
      DO 7 I=1,N
      ZP(I)=ZL*FLOAT(I-1)/FLOAT(N-1)
      TPP(IP,I)=TP(I)
7     CONTINUE
C...
C...  PLOT THE SOLUTION AS SPATIAL PROFILES
      IF(IP.EQ.6)THEN
      CALL SPLOTS(IP,N,ZP,TPP)
C...
C...  COMPUTE AND PRINT THE ANALYTICAL STEADY STATE SOLUTION ALONG WITH
C...  THE NUMERICAL SOLUTION FOR THE SECOND RUN
      WRITE(NO,5)
5     FORMAT(1H ,//,23H STEADY STATE SOLUTIONS,//,
     1 10H         Z,12H   TP (ANAL),12H    TP (NUM))
      DO 3 I=1,N
      Z=ZL*FLOAT(I-1)/FLOAT(N-1)
      TPA=TA+(TE-TA)*EXP(-E*Z)
      WRITE(NO,4)Z,TPA,TP(I)
4     FORMAT(F10.1,2F12.1)
3     CONTINUE
      END IF
      RETURN
      END

