*DECK HEAT23
      SUBROUTINE INITAL
C...
C...  PI CONTROL OF A THREE-PASS SHELL AND TUBE HEAT EXCHANGER
C...
C...  THE FEASIBILITY OF CONTROLLING THE OUTLET TUBE TEMPERATURE OF A
C...  THREE-PASS SHELL AND TUBE HEAT EXCHANGER BY MANIUPLATING THE
C...  SHELL SIDE FLUID FLOW RATE IS TO BE INVESTIGATED WITH THE FOLLOW-
C...  ING PROGRAM.  THE PARTIAL DIFFERENTIAL EQUATIONS (PDES) WHICH
C...  MODEL THE EXCHANGER ARE
C...
C...   T1  + VT*T1  = CT*(TS - T1)
C...     T        Z
C...
C...   T2  - VT*T2  = CT*(TS - T2)
C...     T        Z
C...
C...   T3  + VT*T3  = CT*(TS - T3)
C...     T        Z
C...
C...   TS  - VS*TS  = CS*((T1 - TS) + (T2 - TS) + (T3 - TS))
C...     T        Z
C...
C...  WHERE
C...
C...   T1,T2,T3   TEMPERATURE OF THE FLUID IN TUBE PASSES 1, 2 AND 3,
C...              RESPECTIVELY
C...
C...   TS         TEMPERATURE OF THE FLUID IN THE SHELL
C...
C...   Z          POSITION ALONG THE EXCHANGER
C...
C...   T          TIME
C...
C...   VT,VS      AVERAGE VELOCITY OF THE TUBE AND SHELL SIDE FLUIDS,
C...              RESPECTIVELY
C...
C...   CT,CS      CONSTANTS DEFINING THE HEAT TRANSFER RATES TO THE
C...              TUBE SIDE FLUID AND SHELL SIDE FLUID, RESPECTIVELY
C...
C...  THE EXCHANGER IS TO BE SIMULATED SUBJECT TO THE FOLLOWING CON-
C...  DITIONS AND PARAMETERS:
C...
C...   T1(Z,0) = T2(Z,0) = T3(Z,0) = TS(Z,0) = 25
C...
C...   T1(0,T) = 100, T2(ZL,T) = T1(ZL,T), T3(0,T) = T2(0,T)
C...
C...   TS(ZL,T) = 220
C...
C...   QT = 1, AT = 1, AS = 1, CT = CS = 0.1, ZL = 10, N = 11
C...
C...  THE OUTLET TUBE TEMPERATURE, T3(ZL,T) WILL BE SENSED AND USED
C...  TO ADJUST THE SHELL SIDE FLOW RATE, QS, ACCORDING TO THE CON-
C...  TROLLER EQUATION
C...
C...                                T
C...     QS = QSS + KC*(E + (1/TI)*INT E*DT)
C...                                0
C...
C...  WHERE
C...
C...     QT       TUBE SIDE VOLUMETRIC FLOW RATE
C...
C...     AT       TOTAL TUBE CROSS SECTIONAL AREA IN EACH PASS, I.E.,
C...              QT = VT*AT
C...
C...     QS       SHELL SIDE VOLUMETRIC FLOW RATE
C...
C...     AS       SHELL CROSS SECTIONAL AREA, I.E., QS = AS*VS
C...
C...     QSS      SHELL SIDE VOLUMETRIC FLOW RATE AT STEADY STATE
C...
C...     KC       CONTROLLER GAIN
C...
C...     TI       CONTROLLER INTEGRAL TIME
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    T1(11),    T2(11),    T3(11),     TS(11),        EI
     2      /F/   T1T(11),   T2T(11),   T3T(11),    TST(11),         E
     3      /S/   T1Z(11),   T2Z(11),   T3Z(11),    TSZ(11)
     4      /C/        T0,        ZL,        VT,         QS,
     5                 CT,        CS,       TS0,        QSS,
     6                 KC,        TI,     T3SET,     ZP(11),
     7                  N,        IP,       ISP
      REAL             KC
C...
C...  SET THE MODEL PARAMETERS
      T0=25.
      ZL=10.
      VT=1.0
      QSS=0.5
      CT=0.1
      CS=0.1
      TS0=220.
      T3SET=170.
      IF(NORUN.EQ.1)THEN
         KC=0.
         TI=1.0E+10
      ELSE IF(NORUN.EQ.2)THEN
         KC=0.1
         TI=30.0
      ELSE
      END IF
      N=11
      IP=0
C...
C...  MODEL INITIAL CONDITIONS
      DO 1 I=1,N
      T1(I)=T0
      T2(I)=T0
      T3(I)=T0
      TS(I)=T0
1     CONTINUE
      EI=0.
C...
C...  T0 NOW BECOMES THE ENTERING TUBE TEMPERATUTE (T1(0,T) = T0) IN
C...  SUBROUTINE DERV
      T0=100.
C...
C...  INITIALIZE THE DERIVATIVES
      CALL DERV
C...
C...  INITIALIZE THE SPATIAL GRID USED IN SUBROUTINE PRINT
      ISP=0
      DO 2 I=1,N
      ZP(I)=ZL*FLOAT(I-1)/FLOAT(N-1)
2     CONTINUE
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    T1(11),    T2(11),    T3(11),     TS(11),        EI
     2      /F/   T1T(11),   T2T(11),   T3T(11),    TST(11),         E
     3      /S/   T1Z(11),   T2Z(11),   T3Z(11),    TSZ(11)
     4      /C/        T0,        ZL,        VT,         QS,
     5                 CT,        CS,       TS0,        QSS,
     6                 KC,        TI,     T3SET,     ZP(11),
     7                  N,        IP,       ISP
      REAL             KC
C...
C...  CONTROLLER EQUATIONS
      E=T3SET-T3(N)
      QS=QSS+KC*(E+(1./TI)*EI)
      IF(QS.LT.0.)QS=0.
C...
C...  SINCE AS = 1, VS = QS/AS = QS
      VS=QS
C...
C...  BOUNDARY CONDITIONS
      T1(1)=T0
      TS(N)=TS0
      T2(N)=T1(N)
      T3(1)=T2(1)
C...
C...  MODEL SPATIAL DERIVATIVES
      CALL DSS020(0.,ZL,N,T1,T1Z, 1.)
      CALL DSS020(0.,ZL,N,T2,T2Z,-1.)
      CALL DSS020(0.,ZL,N,T3,T3Z, 1.)
      CALL DSS020(0.,ZL,N,TS,TSZ,-1.)
C...
C...  MODEL PDES
      DO 1 I=1,N
      T1T(I)=-VT*T1Z(I)+CT*(TS(I)-T1(I))
      T2T(I)= VT*T2Z(I)+CT*(TS(I)-T2(I))
      T3T(I)=-VT*T3Z(I)+CT*(TS(I)-T3(I))
      TST(I)= VS*TSZ(I)+CS*((T1(I)-TS(I))+(T2(I)-TS(I))+(T3(I)-TS(I)))
1     CONTINUE
C...
C...  ZERO THE TEMPORAL DERIVATIVES AT THE BOUNDARIES
      T1T(1)=0.
      T2T(N)=0.
      T3T(1)=0.
      TST(N)=0.
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    T1(11),    T2(11),    T3(11),     TS(11),        EI
     2      /F/   T1T(11),   T2T(11),   T3T(11),    TST(11),         E
     3      /S/   T1Z(11),   T2Z(11),   T3Z(11),    TSZ(11)
     4      /C/        T0,        ZL,        VT,         QS,
     5                 CT,        CS,       TS0,        QSS,
     6                 KC,        TI,     T3SET,     ZP(11),
     7                  N,        IP,       ISP
      REAL             KC
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION   T3P(101),  TSP(101),  QSP(101),   TP(101),
     1          T1SP(7,11),T2SP(7,11),T3SP(7,11),TSSP(7,11)
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)WRITE(NO,1)
1     FORMAT(9X,1HT,4X,8H T1(0,T),4X,8HTS(ZL,T),
     14X,8HT3(ZL,T),4X,8H TS(0,T),4X,8H   QS(T))
C...
C...  PRINT THE SOLUTION
      IP=IP+1
      IF((IP-1)/10*10.NE.(IP-1))GO TO 3
      WRITE(NO,2)T,T1(1),TS(N),T3(N),TS(1),QS
2     FORMAT(F10.1,4F12.1,F12.3)
C...
C...  STORE THE SOLUTION FOR PLOTTING
3     TP(IP)=T
      T3P(IP)=T3(N)
      TSP(IP)=TS(1)
      QSP(IP)=QS
C...
C...  THE SPATIAL PROFILES ARE STORED EVERY 20TH CALL TO PRINT
      IF((IP-1)/20*20.NE.(IP-1))GO TO 4
      ISP=ISP+1
      DO 5 I=1,N
      T1SP(ISP,I)=T1(I)
      T2SP(ISP,I)=T2(I)
      T3SP(ISP,I)=T3(I)
      TSSP(ISP,I)=TS(I)
5     CONTINUE
C...
C...  PLOT THE SOLUTION AT THE END OF THE RUN.  THE ENTERING SHELL
C...  FLUID TEMPERATURE IS FIRST STORED TO SCALE THE TEMPERATURE PLOTS
4     IF(IP.LT.101)RETURN
      DO 6 I=1,N
      T1SP(7,I)=TS0
      T2SP(7,I)=TS0
      T3SP(7,I)=TS0
      TSSP(7,I)=TS0
6     CONTINUE
      CALL TPLOTS(1,IP,TP,T3P)
      WRITE(NO,7)
7     FORMAT(1H ,//,31H T3(ZL,T) VS T                 )
      CALL TPLOTS(1,IP,TP,TSP)
      WRITE(NO,8)
8     FORMAT(1H ,//,31H TS(0,T) VS T                  )
      CALL TPLOTS(1,IP,TP,QSP)
      WRITE(NO,9)
9     FORMAT(1H ,//,31H QS VS T                       )
      CALL TPLOTS(7,N,ZP,T1SP)
      WRITE(NO,12)
12    FORMAT(1H ,//,31H T1(Z,T) VS Z, T = 0, 40,.. 200)
      CALL TPLOTS(7,N,ZP,T2SP)
      WRITE(NO,13)
13    FORMAT(1H ,//,31H T2(Z,T) VS Z, T = 0, 40,.. 200)
      CALL TPLOTS(7,N,ZP,T3SP)
      WRITE(NO,14)
14    FORMAT(1H ,//,31H T3(Z,T) VS Z, T = 0, 40,.. 200)
      CALL TPLOTS(7,N,ZP,TSSP)
      WRITE(NO,15)
15    FORMAT(1H ,//,31H TS(Z,T) VS Z, T = 0, 40,.. 200)
      RETURN
      END
