TALK=T;RUN(1,1) idealised shell-and-tube heat exchanger

   PHOTON USE
   p;;;;
   
   vi y
   msg(vectors and contours at symmetry plane
   vec y 1 sh
   pause
   msg(shell-side velocity vectors
   vec off
   msg(shell-side pressure contours
   con ps y 1 fi;.001
   pause
   msg(shell-side temperature contours
   con tems y 1 fi;.001
   pause
   msg(tube-side temperature contours
   con temt y 1 fi;.001
   pause
   
   ENDUSE
   DISPLAY
    A 3D shell-and-tube heat-exchanger, with one tube-side pass,
    and two shell-side baffles, in steady flow.
 
    The shell-side fluid is saturated water; the tube-side fluid
    is SAE_5W-30_engine_oil.

    In-Form is used for many purposes.
    Note by dbs on 24.08.08.
    This file, dated  20.05.05, was probably that which was later 
    used for the 2005 ASME paper. The latter however was developed
    further, not least by the use of MUSES to enable the stresses in
    the tubes to be calculated.
    The diagram has today been changed because it still referred to
    a five-pass exchanger!
    Also the nomenclature has been streamlined so thats rather than 
    1 refers to the shell side and t to the tube side. 
   ENDDIS
    GROUP 1. Run title
TEXT(3D Shell and (1-pass) Tube Heat Exchanger
TITLE
    User-defined variables
    TSIN = inlet temperature of shell fluid
    TTIN = inlet temperature of tube fluid
    FLOS = mass-flow rate of shell fluid
    FLOt = mass-flow rate of tube fluid
    cps  = specific heat of shell fluid
    cpt  = specific heat of tube fluid
    RESCO = flow-resistance coefficient of shell-side fluid
           caused by the presence of the tube bank
 
    The heat exchanger is a rectangular box, 1m high (x),
    1m wide (y) and 4m long (z). A uniform cartesian grid is used.
 
    Only one half of the exchanger is included in the
    calculation domain, because of the symmetry of the
    geometry about the central y=constant plane.
    
                                       ^ shell outlet
                                     ! ! !
       //////////////////////////////!   !///
       ------------------------!------   ----
    >>>>>>>>>>>>>>>>>>>>>>>>>>>!>>>>>>>>>>> !    1 tube-side pass
    ^  ----   ------!------------------------
    !  ///! ^ !//////////////////////////////
    X     ! ! !
    !     shell inlet
    O------ Z ------>
  SAVE1BEGIN        
REAL(TSIN,TTIN,FLOS,FLOT,DIAM,PITCH,areadvol)
diam=0.01; pitch=0.0125
areadvol=3.1416*diam/pitch**2
real(cpscon,cptcon)
cpscon=4138.    ! average specific heats
cptcon=1920.
TSIN=10; TTIN=80.0 ! degrees Celsius
FLOS=100; FLOT=100 ! kg/s mass flows of the fluids
flot=50
REAL(TUBVEL,TUBDEN)
TUBDEN=842
TUBVEL=(FLOT/(2*TUBDEN*XULAST*YVLAST))*(PITCH/DIAM)**2
TUBVEL

REAL(RESCO);RESCO=1.E2 ! this is a guessed constant value
                       ! it ought to be computed point by point
 
xulast=1.0          ! the height
yvlast=0.5          ! one half of the width
zwlast=4.0          ! the length   
nx=30; ny=5; nz=30  ! the computational grid
steady=t
#unigrid
  SAVE1END        
 
    GROUP 7. Variables stored, solved , named
 
    The shell-side fluid is a single-phase one, for which
    five variables must be solved: 
    
    p1, pressure
    u1, x-direction velocity
    v1, y-direction velocity
    w1, z-direction velocity
    h1; enthalpy    
    for the tube-side fluid, only the enthalpy needs be 
    computed. Its velocity along the tube is supposed known. 
    Nonuniformities within the header regions are neglected.
 
SOLVE(P1,U1,V1,W1,H1,ht) ! 1 is shell-side, 2 is tube side
                         ! enthalpies are solved, not temperatures

name(h1)=hs              ! name change, for consistency
name(p1)=ps              ! name change, for consistency

solutn(ps,y,y,y,p,p,p)
  The space occupied by the tubes, and therefore denied to the
  shell-side fluid, is reflected by way of 'porosity' factors.
STORE(tems,temt,VABS,reys,nuss,prns,coes,reyt,nust,prnt,coet,coeu)
  ! tems for shell-side fluid because it is stored, not solved
STORE(TEMM,EPOR,NPOR,HPOR)
STORE(enul,rho1,cps,COND)   ! shell-side properties 
  ! rho1 retained rather than rhos because it is reserved for
  ! phase-1 density in In-Form property statement
STORE(enup,rhot,cpt,cont)   ! tube-side properties
  ! enup is used rather than enut because the latter has a 
  ! reserved meaning
  SAVE7BEGIN 
(stored var tems is hs/cpscon) 
(stored var temt is ht/cptcon)
   
(stored var temm is (coes*temt+coet*tems)/(coes+coet))
(stored var reys is diam*vabs/enul)
(stored var prns is cps*rho1*enul/COND)
(stored var nuss is 0.2*reys^0.6*prns^0.33)
  
(stored var coes is areadvol*nuss*COND/diam)
(stored var reyt is diam*tubvel/enup)
(stored var prnt is cpt*rhot*enup/cont)
(stored var nust is max(2.0,0.328*(reyt*prnt)^0.33))
  (stored var coet is areadvol*nust*cont/diam)
  or, to investigate the importance of the viscosity-ratio effect
  use the following. see below for enum
  
(stored var coet is (enum/enup)^0.14*areadvol*nust*cont/diam)
   (stored var coeU is 1/(1/coes+1/coet))
(stored var coeU is coes*coet/(coes+coet)) ! as above but simpler
                           ! differences in the results were tiny
  
  to compute the volumetric heat flux
(stored var flux is coeu*(temt-tems))
  SAVE7END 
    GROUP 8. Terms (in differential equations) & devices
 
    The "diffusion" terms are here cut out for all variables, but
    could be restored were plausible formulae suggested.
    The built-in sources (kinetic heating) for the enthalpies are
    cut out as .
    Also, the convection terms are cut out for the tube-side 
    fluid because the neighbour-cell sources do what is necessary.
TERMS(U1,Y,Y,N,Y,Y,Y); TERMS(V1,Y,Y,N,Y,Y,Y)
TERMS(W1,Y,Y,N,Y,Y,Y); TERMS(hs,N,Y,N,Y,Y,Y)
TERMS(ht,N,N,N,Y,N,N)
    
    GROUP 9 Material properties
  SAVE9BEGIN  
char(rho_expression)    
char(emu_expression)    
char(enu_expression)    
char(cp_expression)    
char(COND_expression)
char(vn)
vn=tems
   
  get these from core-library case 089
  shell-side fluid is saturated water    
rho_expression=POL5(:VN:,2446.,-20.6741,.11576,-3.12895e-4 ,4.0505$
E-7,-2.054E-10)
 cp_expression=exp((8.29041-.012557*:VN:)/(1.-1.52373e-3*:VN:))
 emu_expression=1.e-7*exp((1.12646-.039638*:VN:) /(1.-7.29769E-3*:V$
N:))

(stored var COND is .001*POL5(:VN:,62.282,-1.768417,.03499,-1.15706E-4,$
  1.53599E-7,-7.7477E-11))
  mesg(rho_expression = :rho_expression:
  mesg(cp_expression  = :cp_expression:
  mesg(emu_expression = :emu_expression:
  mesg(COND_expression= :COND_expression:
  
(property rho1 is :rho_expression:)
(property enul is :emu_expression:/(rho1))
(stored var cps is :cp_expression:)


vn=temt
  tube-side fluid is SAE_5W-30_engine_oil
 rho_expression=1052.3-0.6420*:VN:
 cp_expression=753.7+3.65*:VN:
 emu_expression=10.^(POL4(:VN:,58.2987,-.53817,1.92827e-3, -3.16448$
E-6,1.97922E-9)-2)
 COND_expression=0.1447-2.3073E-5*:VN:
   rho_expression
   cp_expression
   emu_expression
   COND_expression
                               
(stored var rhot is :rho_expression:)
(stored var enup is :emu_expression:/(rhot))
(stored var cpt is :cp_expression:)
(stored var cont is :COND_expression:)
vn=temm
 emu_expression=10.^(POL4(:VN:,58.2987,-.53817,1.92827e-3, -3.16448$
E-6,1.97922E-9)-2)
(stored var enum is :emu_expression:/(rhot))
  SAVE9END  
    GROUP 11. Initialization of variable or porosity fields
  
  specification of porosities, dependent on tube spacing
FIINIT(EPOR)=0.5; FIINIT(NPOR)=0.5; FIINIT(HPOR)=0.5
  note that no calculation of porosities based upon the bundle
  geometry is actually provided here.
fiinit(tems) = TSIN+273 
fiinit(temt) = TTIN+273.0 
fiinit(hs)   = cpscon*fiinit(tems)
fiinit(ht)   = cptcon*fiinit(temt)
    GROUP 13. Boundary CONDitions and special sources
 
integer(izz)
izz = nz/10
izz
    West boundary; shell fluid inlet ;  cells 3 to 5 at low-x
PATCH(INLET1,CELL,1,1,3,5,1,izz,1,lstep)
COVAL(INLET1,ps,FIXFLU,FLOS*5/(3*2.0*NY*izz))  ! divide by 2 because only half 
                                ! exchanger is considered
COVAL(INLET1,hs,ONLYMS,(TSIN+273)*cpscon)    ! shell-fluid specific heat
 
    East boundary; shell fluid outlet; 3 out of 5 cells at high-x
PATCH(OUTLET1,EAST,NX,NX,3,5,NZ-izz+1,NZ,1,lstep)
COVAL(OUTLET1,ps,FIXP,0.0)
 
  FLOT/(2*NY*NX) is the mass flow rate of tube fluid per cell 
               in cross-section
  
    High boundary, tube fluid inlet; NY cells in high-z wall
PATCH(INLET2,CELL,1,NX,1,NY,NZ,NZ,1,1000)
COVAL(INLET2,ht,FLOT/(2*NY*NX),(TTIN+273)*cptcon)  ! tube-fluid specific heat
 
    Note how the giving of special names to patches,
    beginning NE (for neighbour), coupled with LOCNE in the
    "value" location, produces sources which simulate along-
    the-tube convection fluid-to-metal heat transfer etc.
    by activating special calls to built-in open-source coding, 
    in file gxnepat for.
    
    Flow of tube fluid in first pass
PATCH(NEhs,CELL,1,nx,1,NY,1,NZ-1,1,lstep)
COVAL(NEhs,ht,FLOT/(2*NY*nx),LOCNE)   
  ! (NY*nx) = cell number
 
  
patch(heatexch,VOLUME,1,NX,1,NY,1,NZ,1,1000)
  SAVE13BEGIN
  (source of hs at heatexch is coeu*(temt-tems) with line)
  (source of ht at heatexch is coeu*(tems-temt) with line)
  ! with line has no effect because the depedent variables do not 
  ! appear in the source expression
(source of hs at heatexch is coeu*(temt-tems))
(source of ht at heatexch is coeu*(tems-temt))
  SAVE13END
    
    
    Baffle 1 at NZ/3
PATCH(BAFFLE1,HIGH,1,3*NX/5,1,NY,NZ/3,NZ/3,1,1000)
COVAL(BAFFLE1,W1,FIXVAL,0.0)       ! no leakage is allowed for
 
    Baffle 2 at 2*NZ/3
PATCH(BAFFLE2,HIGH,2*NX/5+1,NX,1,NY,2*NZ/3,2*NZ/3,1,1000)
COVAL(BAFFLE2,W1,FIXVAL,0.0)
 
    Resistance to flow exerted by tubes, throughout the shell.
PATCH(RESIST,PHASEM,1,NX,1,NY,1,NZ,1,1000)   ! 2,NZ  ??
COVAL(RESIST,U1,RESCO,0.0); COVAL(RESIST,V1,RESCO,0.0)
COVAL(RESIST,W1,0.5*RESCO,0.0)
 
    GROUP 15. Termination of sweeps
LSWEEP=200
    GROUP 16. Termination of iterations
LITER(ps)=100
    GROUP 17. Under-relaxation devices
RELAX(U1,FALSDT,0.1);RELAX(V1,FALSDT,0.1)
RELAX(W1,FALSDT,0.1)
relax(ps,linrlx,0.5)

varmin(rho1)= 500
varmin(rhot)= 500
varmax(rho1)= 5000
varmax(rhot)= 5000
    
    GROUP 21. Print-out of variables
    Print-out of porosities is suppressed.
OUTPUT(EPOR,N,N,N,N,N,N);OUTPUT(NPOR,N,N,N,N,N,N)
OUTPUT(HPOR,N,N,N,N,N,N)

  SAVE22BEGIN
(longname of hs print as shell-side_fluid_enthalpy)  
(longname of tems print as shell-side_fluid_temperature)  
(longname of rho1 print as shell-side_fluid_density)  
(longname of cps print as shell-side_fluid_specific_heat)  
(longname of cond print as shell-side_fluid_thermal_conductivity)  
(longname of enul print as shell-side_fluid_kinematic_viscosity)  
(longname of ht print as tube-side_fluid_enthalpy)  
(longname of temt print as tube-side_fluid_temperature)  
(longname of rhot print as tube-side_fluid_density)  
(longname of cpt print as tube-side_fluid_specific_heat)  
(longname of cont print as tube-side_fluid_thermal_conductivity)  
(longname of enup print as tube-side_fluid_kinematic_viscosity)  
  SAVE22END
 
    GROUP 21. Spot-value print-out
IXMON=5;IYMON=5
    GROUP 23. Field print-out and plot control
IPLTL=LSWEEP;IPROF=1;ORSIZ=0.4;XZPR=T 
 
PATCH(TABMAP,CONTUR,1,NX,1,NY,NZ/2,NZ/2,1,1000)
PLOT(TABMAP,hs,0.0,10.0);PLOT(TABMAP,ht,0.0,10.0)
 
PATCH(TABYEQ3,CONTUR,1,NX,3,3,1,NZ,1,1000)
PLOT(TABYEQ3,ps,0.0,10.0);PLOT(TABYEQ3,hs,0.0,10.0)
PLOT(TABYEQ3,ht,0.0,10.0)
LIBREF=800
tstswp=-1
#endpause
  #conprom
#maxabs
    GROUP 24. Dumps for restarts