*DECK HEAT30
      SUBROUTINE INITAL
C...
C...  DYNAMICS OF A CROSS FLOW HEAT EXCHANGER
C...
C...  CONSIDER A CROSS FLOW HEAT EXCHANGER WITH THE TUBES RUNNING
C...  HORIZONTALLY IN 28 PASSES (IN THE Y DIRECTION) AND THE SHELL
C...  RUNNING IN ONE PASS (IN THE Z DIRECTION)
C...
C                                      + TH(4)
C                                      .
C            ......................... . ............................
C            .                         .                            .
C            .                         .                            .
C    TC(28)  .                                                      .
C       +.....................................................      .
C            .                      PASS 28                   .     .
C            .      ..........................................      .
C            .     .                PASS 27                         .
C            .      ..........................................      .
C            .                      PASS 26                   .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                      PASS I+1                  .     .
C            .      ..........................................      .
C            .     .                 PASS I                         .
C            .      ..........................................      .
C            .                      PASS I-1                  .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                                                .
C            .      ..........................................      .
C            .                                                .     .
C            .      ..........................................      .
C            .     .                PASS 3                          .
C            .      ..........................................      .
C     TC(0)  .                      PASS 2                    .     .
C       ....+.................................................      .
C            .                      PASS 1                          .
C        +   .                                                      .
C        .   .                         +                            .
C        .   .                         .                            .
C        .   ......................... . ............................
C        Z   Y ---+                    .
C                                      . TH(0)
C...
C...  HOT FLUID ENTERS THE SHELL AT Z = 0 WITH TEMPERATURE TH(0) AND
C...  LEAVES THE SHELL AT Z = ZL WITH TEMPERATURE TH(4).  COLD FLUID
C...  ENTERS THE FIRST TUBE PASS AT Y = 0 WITH TEMPERATURE TC(0) AND
C...  LEAVES THE 28TH TUBE PASS AT Y = YL (MEASURED ALONG THE TUBE
C...  PASSES) WITH TEMPERATURE TC(28).
C...
C...  IN ORDER TO ACCOUNT FOR THE DISTRIBUTED CHARACTERISTIC OF THE
C...  EXCHANGER (THE SPATIAL DISTRIBUTION OF THE FLUID TEMPERATURES),
C...  THE SHELL IS DIVIDED INTO FOUR PERFECTLY MIXED SECTIONS, AND
C...  THE TUBE PASSES ARE DIVIDED INTO 28 PERFECTLY MIXED SECTIONS.
C...  THUS THE SUBSCRIPTING OF THE EXITING TEMPERATURES, TH(4) AND
C...  TC(28).  THE NUMBERS OF SECTIONS HAVE BEEN CHOSEN ARBITRARILY.
C...  THE MODEL SHOULD THEREFORE BE RUN WITH DIFFERING NUMBERS OF
C...  SECTIONS TO DETERMINE HOW SENSITIVE THE NUMERICAL SOLUTION IS TO
C...  THE NUMBERS OF SECTIONS USED.
C...
C...  THUS, THERE ARE SEVEN TUBE SECTIONS FOR EACH SHELL SECTION.  THE
C...  TEMPERATURES WITHIN THE INDIVIDUAL SHELL AND TUBE SECTIONS CAN
C...  NOW BE COMPUTED BY APPLYING AN ENERGY CONSERVATION EQUATION TO
C...  EACH SECTION.
C...
C...  A SECTION OF TUBE OF LENGTH DY IS DEPICTED BELOW
C...
C...                         +.... DY ...+
C...             ......................................
C...                         .           .
C...                         .    Q(I)   .
C...          TC(I-1) ......+.     +     .......+ TC(I)
C...                         .     .     .
C...                         .     .     .
C...              .....................................
C...                               .
C...                             TH(J)
C...
C...  AN ENERGY BALANCE ON THE TUBE SIDE FLUID IN THE SECTION OF LENGTH
C...  DY, ASSUMING PERFECT MIXING IN THIS SECTION, GIVES
C...
C...     DY*AT*RHOM*CPM*DTC(I)/DT =
C...                                                            (1)
C...        VM*AT*RHOM*CPM*(TC(I-1) - TC(I)) + Q(I)
C...
C...  WHERE
C...
C...     Q(I) = PI*DT*DY*U*(TH(J) - TC(I))                      (2)
C...
C...  THE NOTATION FOR EQUATIONS (1) AND (2) IS DEFINED BELOW.
C...
C...  A SECTION OF SHELL OF LENGTH DZ IS DEPICTED BELOW
C...
C...                         .           .
C...                         .           .
C...                         .   TH(J)   .
C...                         .     +     .
C...                         .     .     .
C...                         .     .     .
C...                         .     .     .
C...                       + .............
C...                       . .           .
C...                       . .           .
C...                       . .           .    7
C...                      DZ .      +....... SUM Q(I)
C...                       . .           .   I=1
C...                       . .           .
C...                       . .           .
C...                       + .............
C...                         .     +     .
C...                         .     .     .
C...                         .     .     .
C...                         .     .     .
C...                         .  TH(J-1)  .
C...                         .           .
C...                         .           .
C...
C...  AN ENERGY BALANCE ON A SECTION OF THE SHELL OF LENGTH DZ, ASSUMING
C...  PERFECT MIXING IN THE SECTION GIVES
C...
C...     DZ*AS*RHOH*CPH*DTH(J)/DT =
C...                                                            (3)
C...                                            7
C...        VH*AS*RHOH*CPH*(TH(J-1) - TH(J)) + SUM Q(I)
C...                                           I=1
C...
C...  EQUATIONS (1) AND (3) HAVE THE INITIAL CONDITIONS
C...
C...     TC(I) = 52, I = 0, 1, 2, ..., 28                       (4)
C...
C...     TH(J) = 94, J = 0, 1, 2, ..., 4                        (5)
C...
C...  NOTE THAT I AND J ARE THE INDICES FOR THE TUBE AND SHELL SECTIONS,
C...  RESPECTIVELY, AND THERE ARE 28 AND 4 SECTIONS OF THE TUBE AND
C...  SHELL, RESPECTIVELY (THESE NUMBERS OF SECTIONS CAN EASILY BE
C...  CHANGED IN THE SUBSEQUENT PROGRAMMING).
C...
C...  A PI CONTROLLER FOR THE EXITING TUBE TEMPERATURE IS ALSO INCLUDED
C...  IN THE MODEL THROUGH THE EQUATION
C...
C...     VH = VHSS + AKC2*(E2 + (1/TI2)*E2I)                    (6)
C...
C...  THERE E2 IS THE ERROR IN THE EXITING TUBE TEMPERATURE
C...
C...     E2 = TCSET - TC(28)                                    (7)
C...
C...  THUS, THE MODEL CONSISTS OF 29 + 5 + 1 = 35 ORDINARY DIFFERENTIAL
C...  EQUATIONS.
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1     /Y/   TC(0:28),   TH(0:4),       E2I
     2     /F/  TCT(0:28),  THT(0:4),        E2
     3     /C/         NC,        NH,        VM,        PI,        AS,
     4                 AT,      RHOH,       CPH,      RHOM,       CPM,
     5                 ZL,        DZ,        YL,        DY,      AKC2,
     6                TI2,     Q(28),         U,        DT,        VH,
     7             CON(5),     TCSET,        IP
C...
C...  MODEL PARAMETERS
C...
C...  NUMBER OF TUBE AND SHELL SECTIONS
      NC=28
      NH=4
C...
C...  TUBE FLUID VELOCITY
      VM=135.0
C...
C...  THE CONSTANT PI
      PI=3.1415927
C...
C...  SHELL AND TUBE CROSS SECTIONAL AREAS
      AS=929.0
      AT=11.43
C...
C...  SHELL FLUID DENSITY AND HEAT CAPACITY
      RHOH=1.0
      CPH=1.0
C...
C...  TUBE FLUID DENSITY AND HEAT CAPACITY
      RHOM=1.2
      CPM=0.55
C...
C...  SHELL LENGTH AND SECTION LENGTH
      ZL=305.0
      DZ=ZL/FLOAT(NH)
C...
C...  TUBE LENGTH AND SECTION LENGTH
      YL=30480.0
      DY=YL/FLOAT(NC)
C...
C...  CONTROLLER SET POINT
      TCSET=75.0
C...
C...  CONTROLLER GAIN AND INTEGRAL TIME
C...
C...     NO CONTROL
         IF(NORUN.EQ.1)THEN
            AKC2=0.
            TI2=1.0E+10
C...
C...     PI CONTROL  
C...     THIS CONTROLLER NEEDS SOME ADDITIONAL TUNING TO CLOSELY 
C...     ACHIEVE THE SET POINT TCSET = 75.  ALSO, THE CONTROLLER
C...     IS APPARENTLY INOPERABLE WITH TCSET = 82.6 UNLESS THE 
C...     TEMPERATURE OF THE SHELL FLUID IS RAISED ABOVE 94.  THE
C...     LIMITS ON VH PROGRAMMED AT THE BEGINNING OF SUBROUTINE
C...     DERV PROBABLY HAVE SOME BEARING ON THESE CONCLUSIONS
         ELSE IF(NORUN.EQ.2)THEN
            AKC2=0.2
            TI2=1.0E-02
         END IF
C...
C...  HEAT TRANSFER COEFFICIENT AND TUBE DIAMETER
      U=0.00264
      DT=3.81
C...
C...  SHELL FLUID VELOCITY
      VH=0.958
C...
C...  PRECOMPUTE SOME CONSTANTS
      CON(1)=PI*DT*DY*U
      CON(2)=DY*AT*RHOM*CPM
      CON(3)=VM*AT*RHOM*CPM
      CON(4)=DZ*AS*RHOH*CPH
      CON(5)=VH*AS*RHOH*CPH
C...
C...  INITIAL CONDITIONS
      DO 1 I=0,NC
      TC(I)=52.0
1     CONTINUE
      DO 2 I=0,NH
      TH(I)=94.0
2     CONTINUE
      E2I=0.
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1     /Y/   TC(0:28),   TH(0:4),       E2I
     2     /F/  TCT(0:28),  THT(0:4),        E2
     3     /C/         NC,        NH,        VM,        PI,        AS,
     4                 AT,      RHOH,       CPH,      RHOM,       CPM,
     5                 ZL,        DZ,        YL,        DY,      AKC2,
     6                TI2,     Q(28),         U,        DT,        VH,
     7             CON(5),     TCSET,        IP
C...
C...  MODEL ALGEBRA
C...
C...     ERROR IN EXITING TUBE SIDE TEMPERATURE
         E2=TCSET-TC(NC)
C...
C...     SHELL SIDE VELOCITY
         VH=0.958+AKC2*(E2+(1./TI2)*E2I)
C...
C...     LIMIT THE SHELL SIDE VELOCITY FROM THE CONTROLLER EQUATION
C...     (THESE LIMITS REFLECT THE PHYSICAL REALITY OF WHAT A CONTROL
C...     SYSTEM COULD PROVIDE)
         IF(VH.LT.0.)VH=0.
         IF(VH.GT.3.0)VH=3.0
C...
C...     PUT THE SHELL SIDE VELOCITY IN THE MODEL EQUATIONS
         CON(5)=VH*AS*RHOH*CPH
C...
C...  BOUNDARY CONDITIONS
C...     TUBE SIDE
         TC(0)=52.0
C...
C...     SHELL SIDE
         TH(0)=94.0
C...
C...  DIFFERENTIAL EQUATIONS
C...
C...  STEP THROUGH THE PASTEURIZER
C...
C...  INITIALIZE THE SHELL SIDE SECTION NUMBER, HEAT TRANSFER RATE
         J=1
         QZ=0.
C...
C...  STEP ALONG THE NC TUBE SECTIONS
         DO 10 I=1,NC
            Q(I)=CON(1)*(TH(J)-TC(I))
            TCT(I)=(1./CON(2))*(CON(3)*(TC(I-1)-TC(I))+Q(I))
            QZ=QZ+Q(I)
C...
C...  CHECK THE SHELL SIDE SECTION NUMBER, AND USE THE CORRESPONDING
C...  ENERGY BALANCE.  FOR A NEW SHELL SIDE SECTION, INCREMENT THE
C...  SHELL SIDE NUMBER AND INITIALIZE THE HEAT TRANSFER RATE FOR THE
C...  NEW SHELL SIDE SECTION
         IF((I.EQ. 7).OR.
     1      (I.EQ.14).OR.
     2      (I.EQ.21).OR.
     3      (I.EQ.28))THEN
            THT(J)=(1./CON(4))*(CON(5)*(TH(J-1)-TH(J))-QZ)
            J=J+1
            QZ=0.
C...
C...        CHECK FOR THE CALCULATION OF THE TIME DERIVATIVES FOR THE
C...        LAST TUBE AND SHELL SECTIONS
            IF(J.GT.NH)GO TO 10
         END IF
10       CONTINUE
C...
C...  BOUNDARY TEMPERATURES ARE INVARIANT
C...
C...     TUBE SIDE
         TCT(0)=0.
C...
C...     SHELL SIDE
         THT(0)=0.
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1     /Y/   TC(0:28),   TH(0:4),       E2I
     2     /F/  TCT(0:28),  THT(0:4),        E2
     3     /C/         NC,        NH,        VM,        PI,        AS,
     4                 AT,      RHOH,       CPH,      RHOM,       CPM,
     5                 ZL,        DZ,        YL,        DY,      AKC2,
     6                TI2,     Q(28),         U,        DT,        VH,
     7             CON(5),     TCSET,        IP
      REAL        TP(101),TEP(2,101),   EP(101),  EIP(101),  VHP(101),
     1             YP(29), TCP(6,29)
C...
C...  INITIALIZE A COUNTER FOR THE TUBE FLUID TEMPERATURE PROFILES
      IF(IP.EQ.0)IP1=0
C...
C...  PRINT THE SOLUTION EVERY 20TH CALL TO SUBROUTINE PRINT
      IP=IP+1
      IF(((IP-1)/20*20).EQ.(IP-1))THEN
      WRITE(NO,1)T
1     FORMAT(1H ,//,5H T = ,E10.3)
      WRITE(NO,2)( I, TC(I), TCT(I),I=0,NC)
2     FORMAT(1H ,//,10H         I,10H     TC(I),2X,10H    TCT(I),/,
     1       (I10,F10.1,E12.3))
      WRITE(NO,3)( J, TH(J), THT(J),J=0,NH)
3     FORMAT(1H ,//,10H         J,10H     TH(J),2X,10H    THT(J),/,
     1       (I10,F10.1,E12.3))
C...
C...  STORE THE SOLUTION FOR PLOTTING THE TEMPERATURE PROFILES
      IP1=IP1+1
      DO 4 I=0,NC
      IF(IP1.EQ.1)YP(I+1)=DY*FLOAT(I)
      TCP(IP1,I+1)=TC(I)
4     CONTINUE
      END IF
C...
C...  STORE THE SOLUTION FOR PLOTTING THE TIME RESPONSE
      TP(IP)=T
      TEP(1,IP)=TC(NC)
      TEP(2,IP)=TH(NH)
      EP(IP)=E2
      EIP(IP)=(1.0/TI2)*E2I
      VHP(IP)=VH
C...
C...  PLOT THE SOLUTION AT THE END OF THE RUN
      IF(IP.LT.101)RETURN
C...
C...  EXITING TEMPERATURES VS T
      CALL TPLOTS(2,IP,TP,TEP)
      WRITE(NO,10)
10    FORMAT(1H ,//,
     1 34H 1 - TC(NC), 2 - TH(NH) VS T (SEC))
C...
C...  ERROR VS T
      CALL TPLOTS(1,IP,TP,EP)
      WRITE(NO,11)
11    FORMAT(1H ,//,
     1 34H1 - E VS T (SEC)                  )
C...
C...  INTEGRAL TERM IN CONTROLLER EQUATION VS T
      CALL TPLOTS(1,IP,TP,EIP)
      WRITE(NO,12)
12    FORMAT(1H ,//,
     1 34H1 - (1/TI2) INT E DT VS T (SEC)   )
C...
C...  SHELL SIDE VELOCITY VS TIME
      CALL TPLOTS(1,IP,TP,VHP)
      WRITE(NO,13)
13    FORMAT(1H ,//,
     1 34H 1 - VH VS T (SEC)                )
C...
C...  PLOT THE TUBE TEMPERATURE PROFILES
      CALL TPLOTS(IP1,NC+1,YP,TCP)
      WRITE(NO,21)
21    FORMAT(1H ,//,
     1 56H TUBE FLUID TEMPERATURE VS Y, T = 0, 200, 400, ..., 1000)
      RETURN
      END

