*DECK HEAT24
      SUBROUTINE INITAL
C...
C...  PARAMETER ESTIMATION IN HEAT TRANSFER SYSTEMS
C...
C...  THE FOLLOWING PROGRAMMING ILLUSTRATES A METHOD BY WHICH PARAMETERS
C...  CAN BE ESTIMATED IN HEAT TRANSFER SYSTEMS.  SPECIFICALLY, THE HEAT
C...  TRANSFER COEFFICIENT OF A THREE-PASS SHELL AND TUBE HEAT EXCHANGER
C...  IS TO BE ESTIMATED FROM THE STEADY STATE INLET AND OUTLET TEMPERA-
C...  TURES.
C...
C...  THE PARTIAL DIFFERENTIAL EQUATIONS WHICH MODEL THE HEAT EXCHANGER
C...  ARE
C...
C...   T1  + VT*T1  = UT*(TS - T1)
C...     T        Z
C...
C...   T2  - VT*T2  = UT*(TS - T2)
C...     T        Z
C...
C...   T3  + VT*T3  = UT*(TS - T3)
C...     T        Z
C...
C...   TS  - VS*TS  = US*((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...   UT,US      CONSTANTS DEFINING THE HEAT TRANSFER RATES TO THE
C...              TUBE SIDE FLUID AND SHELL SIDE FLUID, RESPECTIVELY
C...
C...  NOTE THAT THE ESTIMATION PROCEDURE IS BASED ON THE USE OF A
C...  DYNAMIC MODEL FOR THE HEAT EXCHANGER WHICH WILL CONVERGE TO THE
C...  GIVEN STEADY STATE INLET AND OUTLET TEMPERATURES.
C...
C...  UT AND US ARE A GROUP OF CONSTANTS FOR THE TUBE AND SHELL SIDES,
C...  RESPECTIVELY, WHICH CONTAIN THE OVERALL HEAT TRANSFER COEFFI-
C...  CIENT.  UT AND US ARE TO BE ESTIMATED FROM THE STEADY STATE INLET
C...  AND OUTLET TEMPERATURES OF THE HEAT EXCHANGER, INDICATED BELOW
C...  BELOW, ALONG WITH THE OTHER PARAMETERS OF THE EXCHANGER
C...
C...     INLET TEMPERATURES
C...
C...        T1(0,T) = 100, TS(ZL,T) = 220
C...
C...     OUTLET TEMPERATURES
C...
C...        T3(ZL,T) = 170
C...
C...        AT STEADY STATE, THE OTHER OUTLET TEMPERATURE, TS(O,T),
C...        IS SET BY THE OVERALL HEAT BALANCE FOR THE EXCHANGER
C...
C...        VT*(T3(ZL,T) - T1(0,T)) = VS*(TS(ZL,T) - TS(0,T))
C...
C...        TWO CASES ARE PROGRAMMED
C...
C...        (1)  VS = 1
C...
C...             FOR VT = 1, FROM THE OVERALL ENERGY BALANCE,
C...
C...             TS(0,T) = 150
C...
C...        (2)  VS = 0.5
C...
C...             FOR VT = 1, FROM THE OVERALL ENERGY BALANCE,
C...
C...             TS(0,T) = 80
C...
C...     OTHER PARAMETERS
C...
C...        ZL = 10, N = 11
C...
C...  FOR CASE (1), VS = 1, THE OUTLET TEMPERATURE TS(0,T) = 150 IS
C...  PHYSICALLY REALIZABLE AND THEREFORE THE ESTIMATION OF THE HEAT
C...  TRANSFER COEFFICIENT IS POSSIBLE.  FOR CASE (2), VS = 0.5, THE
C...  OUTLET TEMPERATURE TS(0,T) = 80 IS PHYSICALLY UNREALIZABLE SINCE
C...  IT IS BELOW BOTH ENTERING TEMPERATURES.  THE PARAMETER ESTIMATION
C...  WILL STILL EXECUTE, BUT THE FINAL HEAT TRANSFER COEFFICIENT GIVES
C...  AN OUTLET TEMPERATURE, TS(0,T), WHICH DOES NOT SATISFY THE OVER-
C...  ALL ENERGY BALANCE.  CASE (2) WAS INCLUDED TO DEMONSTRATE A PIT-
C...  FALL THAT CAN OCCUR IN PARAMETER ESTIMATION, I.E., THE ESTIMATED
C...  PARAMETER MIGHT APPEAR REASONABLE, BUT ONE OR MORE BASIC REQUIRE-
C...  MENTS OF THE RESULTING MODEL HAS BEEN VIOLATED, IN THIS CASE, THE
C...  OVERALL ENERGY BALANCE.
C...
C...  THE BOUNDARY CONDITIONS FOR THE PDES ARE (IN ADDITION TO THE INLET
C...  TEMPERATURE CONDITIONS
C...
C...        T2(ZL,T) = T1(ZL,T), T3(0,T) = T2(0,T)
C...
C...  THE INITIAL CONDITIONS FOR THE DYNAMIC MODEL ARE TAKEN AS
C...
C...        T1(Z,0) = T2(Z,0) = T3(Z,0) = TS(Z,0) = 25
C...
C...  THESE INITIAL CONDITIONS ARE SELECTED ARBITRARILY, AND HAVE NO
C...  EFFECT ON THE FINAL STEADY STATE OF THE MODEL.
C...
C...  THE OUTLET TUBE TEMPERATURE, T3(ZL,T), IS COMPARED WITH THE DESIR-
C...  ED VALUE, 170.  THE DIFFERENCE OR ERROR, E = 170 - T3(ZL,T), IS
C...  USED TO ADJUST UT ACCORDING TO THE EQUATION
C...
C...                                   T
C...        UT = UTE + KC*(E + (1/TI)*INT E*DT)
C...                                   0
C...
C...  WHERE
C...
C...     UT       VALUE OF UT FROM THE PARAMETER ESTIMATION ALGORITHM
C...
C...     UTE      INITIAL ESTIMATE OF UT
C...
C...     E        ERROR IN THE OUTLET TUBE TEMPERATURE = 170 - T3(ZL,T)
C...
C...     KC       GAIN OF THE PARAMETER ESTIMATION EQUATION
C...
C...     TI       INTEGRAL TIME OF THE PARAMETER ESTIMATION EQUATION
C...
C...  THE EQUATION FOR UT WILL ADJUST UT UNTIL THE ERROR, E, IS REDUCED
C...  TO ZERO CORRESPONDING TO THE CORRECT OUTLET TUBE TEMPERATURE.
C...  AT THIS POINT, THE OUTLET SHELL TEMPERATURE, TS(0,T), SHOULD HAVE
C...  THE CORRECT VALUE (FOR CASE (1), WHILE FOR CASE (2), THE OVERALL
C...  HEAT BALANCE WILL NOT BE SATISFIED).
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,       TS0,        VT,         VS,
     5                 UT,        US,       UTE,         KC,
     6                 TI,     T3SET,        ZL,     ZP(11),
     7                  N,        IP,       ISP
      REAL             KC
C...
C...  SET THE MODEL PARAMETERS
      T0=25.
      TS0=220.
      VT=1.0
      IF(NORUN.EQ.1)VS=1.0
      IF(NORUN.EQ.2)VS=0.5
      ZL=10.
      N=11
      IP=0
C...
C...  SET THE INITIAL ESTIMATE OF UT
      UTE=0.5
C...
C...  DEFINE THE REMAINDER OF THE PARAMETER ESTIMATION ALGORITHM
      KC=0.05
      TI=100.
      T3SET=170.
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,       TS0,        VT,         VS,
     5                 UT,        US,       UTE,         KC,
     6                 TI,     T3SET,        ZL,     ZP(11),
     7                  N,        IP,       ISP
      REAL             KC
C...
C...  PARAMETER ADJUSTMENT
      E=T3SET-T3(N)
      UT=UTE+KC*(E+(1./TI)*EI)
      IF(UT.LT.0.)UT=0.
C...
C...  US AND UT ARE RELATED BY A SINGLE PROPORTIONALITY CONSTANT, A,
C...  I.E., US = A*UT, WHICH INCLUDES THE RATIO OF THE TOTAL TUBE CROSS
C...  SECTIONAL AREA TO THE TOTAL SHELL CROSS SECTIONAL AREA (FOR FLOW
C...  IN BOTH CASES).  FOR A GIVEN HEAT EXCHANGER, THIS PROPORTIONALITY
C...  CONSTANT IS DEFINED BY THE EXCHANGER DIMENSIONS AND GEOMETRY.  IN
C...  THE PRESENT CASE, FOR SIMPLICITY, WE WILL TAKE A = 1
      US=UT
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)+UT*(TS(I)-T1(I))
      T2T(I)= VT*T2Z(I)+UT*(TS(I)-T2(I))
      T3T(I)=-VT*T3Z(I)+UT*(TS(I)-T3(I))
      TST(I)= VS*TSZ(I)+US*((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,       TS0,        VT,         VS,
     5                 UT,        US,       UTE,         KC,
     6                 TI,     T3SET,        ZL,     ZP(11),
     7                  N,        IP,       ISP
      REAL             KC
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION   T3P(101),  TSP(101),  UTP(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   UT(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),UT
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)
      UTP(IP)=UT
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,UTP)
      WRITE(NO,9)
9     FORMAT(1H ,//,31H UT VS T                       )
      CALL TPLOTS(7,N,ZP,T1SP)
      WRITE(NO,12)
12    FORMAT(1H ,//,31H T1(Z,T) VS Z, T = 0, 60,.. 300)
      CALL TPLOTS(7,N,ZP,T2SP)
      WRITE(NO,13)
13    FORMAT(1H ,//,31H T2(Z,T) VS Z, T = 0, 60,.. 300)
      CALL TPLOTS(7,N,ZP,T3SP)
      WRITE(NO,14)
14    FORMAT(1H ,//,31H T3(Z,T) VS Z, T = 0, 60,.. 300)
      CALL TPLOTS(7,N,ZP,TSSP)
      WRITE(NO,15)
15    FORMAT(1H ,//,31H TS(Z,T) VS Z, T = 0, 60,.. 300)
C...
C...  ******************************************************************
C...
C...  THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE-DEPENDENT
C...  AND WILL HAVE TO BE CHANGED FOR OTHER COMPUTERS, OR CONVERTED TO
C...  COMMENTS
C...
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(TP,T3P,IP,3H*T*,10H*T3(ZL,T)*,2H**,1)
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(TP,TSP,IP,3H*T*,10H*TS(0,T) *,2H**,1)
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(TP,UTP,IP,3H*T*,10H* UT(T)  *,2H**,1)
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(ZP,T1SP,N,3H*Z*,10H*T1(Z,T) *,2H**,7)
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(ZP,T2SP,N,3H*Z*,10H*T2(Z,T) *,2H**,7)
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(ZP,T3SP,N,3H*Z*,10H*T3(Z,T) *,2H**,7)
C...  CALL QIKSAZ(3,3)
C...  CALL QIKPLT(ZP,TSSP,N,3H*Z*,10H*TS(Z,T) *,2H**,7)
C...
C...  ******************************************************************
C...
      RETURN
      END

