PHOTON USE
  p;parphi
  5 1
  ;;;;
  
  up z
  msg pressures
  con p1 x 1 fi;0.001;pause;con off;red
  msg phase-1 vol fract
  con r1 x 1 fi;0.001;pause;con off;red
  msg phase-1 radial vel.
  con v1 x 1 fi;0.001;pause;con off;red
  msg phase-2 radial vel.
  con v2 x 1 fi;0.001;pause;con off;red
  msg phase-1 axial vel.
  con w1 x 1 fi;0.001;pause;con off;red
  msg phase-2 axial vel.
  con w2 x 1 fi;0.001;pause;con off;red
  
  ENDUSE
#cls
TEXT(2-PHASE BUBBLY AIR/WATER PIPE FLOW
TITLE
  DISPLAY
    The case considered is 2-phase turbulent air-water flow in a
    pipe, as studied experimentally by Seriwaza et al [1992] for
    upward flow and by Lahey et al [1992] for both upward and
    downward flow. 
    
    The input file is set up to run any one of these three cases. 
    Each calculation is performed with the parabolic option, and,
    for testing purposes the calculation is terminated 5 diameters 
    downstream. 
    
    However, for comparison with data the calculation should be 
    continued to 35 diameters downstream, which corresponds to the 
    experimental measuring station. 
    
    The 2-phase model accounts for interfacial drag, lift, pressure 
    and virtual-mass forces. 
#pause    
    The standard k-e model is employed with the option to select one 
    of two modifications to account for bubble-induced turbulence.
    
    For upward flow, the predicted and measured void-fraction 
    profiles show that the gas is taken away from the centre and
    towards the walls, while for downward flow the reverse is
    observed. 
    
    For upward flow the predictions are in reasonable agreement with 
    the data, but much less so for downward flow. However, the 
    influences of mesh size, the interfacial-modelling coefficients 
    and the turbulence-modelling modifications need to be investigated.
    
    For the results of a study made with a 1995 version of 
    PHOENICS, click
  here.
  ENDDIS
#pause 
INTEGER(JRUN,KMOD)
MESG( Select the required test case:
MESG( The options are:
MESG(  1  - Seriwaza upflow data (default)
MESG(  2  - Lahey upflow data
MESG(  3  - Lahey downflow data
MESG(
READVDU(JRUN,INT,1)
REAL(RHOL,RHOG,EMULIQ,EMUGAS,RGAS,RLIQ,VSGAS,VSLIQ,TKEIN,EPSIN)
REAL(VSTOT,VIN1,VIN2,RIN1,RIN2,FLOWL,FLOWG,PI,XGAS,XLIQ,FRIC)
REAL(VGAS,VLIQ,DIAM,GRAD,PLEN,DIAMB,EMUMIX,REYMIX,GRAVAC)
RHOL=1000.;EMUGAS=1.8E-5;PI=3.14159
  **VSGAS & VSLIQ are superficial velocities
  **Upflow data of Seriwaza et al, 'Phase distribution in bubbly
    flow', In Multiphase Science & Technology, Vol.6, p257-301,
    Ed. G.F.Hewitt, J.M.Delhaye, & N.Zuber, Hemisphere Publ.
    Corp [1992].
IF(JRUN.EQ.1) THEN
+ MESG(Seriwaza upflow data
+ VSGAS=0.077;VSLIQ=1.36;GRAVAC=-9.81;DIAM=0.06;EMULIQ=1.E-3
+ RHOG=1.34
ENDIF
  **Upflow data of Lahey et al, 'Phase distribution and 2-phase
    flow', In Multiphase Science & Technology, Vol.6, p303-349,
    Ed. G.F.Hewitt, J.M.Delhaye, & N.Zuber, Hemisphere Publ.
    Corp [1992].
IF(JRUN.EQ.1) THEN
+ TEXT(BUBBLY AIR/WATER PIPE- Seriwaza upflow
ENDIF
IF(JRUN.EQ.2) THEN
+ MESG(Lahey upflow data
+ VSGAS=0.1;VSLIQ=1.08;GRAVAC=-9.81;DIAM=0.057;EMULIQ=1.207E-3
+ RHOG=1.24
+ TEXT(BUBBLY AIR/WATER PIPE- Lahey upflow
ENDIF
  **Downflow data of Lahey et al [1992].]
IF(JRUN.EQ.3) THEN
+ MESG(Lahey downflow data
  *** Downflow data of Wang et al [1992]
+ VSGAS=0.1;VSLIQ=1.08;GRAVAC=9.81;DIAM=0.057;EMULIQ=1.207E-3
+ RHOG=1.24
+ TEXT(BUBBLY AIR/WATER PIPE- Lahey downflow
ENDIF
VSTOT=VSGAS+VSLIQ
   *** compute inlet void fractions presuming no slip at inlet
RGAS=VSGAS/VSTOT;RLIQ=1.-RGAS
   *** compute inlet phase velocities
VGAS=VSGAS/RGAS;VLIQ=VSLIQ/RLIQ
   *** specify pipe radius & length
GRAD=0.5*DIAM;PLEN=35.*DIAM
   *** specify mean bubble diameter
DIAMB=3.E-3
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
PARAB=T
    GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1;NX=1
    GROUP 4. Y-direction grid specification
GRDPWR(Y,30,GRAD,1.0)
    GROUP 5. Z-direction grid specification
  *** specify uniform step-size of 0.1*DIAM
NZ=50
  NZ=350
GRDPWR(Z,NZ,NZ*PLEN/350,1.0)
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
ONEPHS=F;SOLVE(P1,V1,V2,W1,W2,R1,R2)
TURMOD(KEMODL);STORE(ENUT);STORE(CD,REYN)
MESG( Enter the required k-e turbulence modification:
MESG(  1  - No modifications (default)
MESG(  2  - Additional KE & EP sources due to bubbles
MESG(  3  - Enhanced ENUT due to bubbles
MESG(
READVDU(KMOD,INT,1)
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
RHO1=RHOL;RHO2=RHOG;VIN1=VLIQ;VIN2=VGAS;RIN1=RLIQ;RIN2=RGAS
FLOWL=RHO1*RIN1*VIN1;FLOWG=RHO2*RIN2*VIN2;ENUL=EMULIQ/RHOL
  ** estimate Reynolds number
XGAS=FLOWG/(FLOWL+FLOWG);XLIQ=1.-XGAS
EMUMIX=1.0/(XGAS/EMUGAS+XLIQ/EMULIQ)
REYMIX=(FLOWL+FLOWG)*DIAM/EMUMIX;FRIC=0.3164/REYMIX**0.25
TKEIN=FRIC*VIN1*VIN1/8.;EPSIN=0.1643*TKEIN**1.5/(0.1*GRAD)
FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN
    GROUP 10. Inter-phase-transfer processes and properties
  *** select dirty-water spherical bubble drag correlation
CFIPS=GRND7;CFIPD=5.0;CFIPA=1.E-3;RLOLIM=1.E-3;CFIPB=DIAMB
  *** select interfacial lift, virtual-mass & pressure forces
CLIFT=GRND2;CVM=GRND2;CVMA=0.5;CPIP=GRND2;CPIPA=0.25;CLIFTA=0.075
STORE(VMSW,VMSV)
IF(JRUN.EQ.3) THEN
+ CPIPA=0.1
ENDIF
INTSOR(LIFT,CLIFT,CLIFTA,RELAX,0.1)
INTSOR(INTPL,CPIP,CPIPA)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=VIN1;FIINIT(W2)=VIN2;FIINIT(R1)=RIN1;FIINIT(R2)=RIN2
    GROUP 12. Unused
    GROUP 13. Boundary conditions and special sources
  ** inlet boundary
INLET(IN,LOW,1,NX,1,NY,1,1,1,1)
VALUE(IN,P1,FLOWL);VALUE(IN,W1,VIN1)
VALUE(IN,P2,FLOWG);VALUE(IN,W2,VIN2)
VALUE(IN,KE,TKEIN);VALUE(IN,EP,EPSIN)
  ** gravity
PATCH(GRAVITY,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(GRAVITY,W2,FIXFLU,GRAVAC*(1.-RHOL/RHOG))
  ** wall boundary
WALL(NWALL,NORTH,1,NX,NY,NY,1,NZ,1,1)
  ** Bubble-induced turbulence production
IF(KMOD.EQ.2) THEN
+ EL1A=0.1;PATCH(KEDI,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(KEDI,KE,FIXFLU,GRND3);COVAL(KEDI,EP,FIXFLU,GRND3)
ENDIF
IF(KMOD.EQ.3) THEN
+ ENUTB=1.0;ENUTC=0.3;KELIN=3
ENDIF
    GROUP 15. Termination of sweeps
  ** deactivate selref as unreliable in V3.1
LITHYD=40;SELREF=F;RESFAC=0.01
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=0.1*YVLAST/VLIQ   
dtf=dtf*0.1 
RELAX(V1,FALSDT,0.01*DTF);RELAX(V2,FALSDT,0.01*DTF)
RELAX(W1,FALSDT,DTF);RELAX(W2,FALSDT,DTF)
RELAX(R1,LINRLX,0.3);RELAX(R2,LINRLX,0.3)
RELAX(KE,LINRLX,0.3);RELAX(EP,LINRLX,0.3)
    GROUP 18. Limits on variables or increments to them
VARMAX(W1)=10.0
VARMIN(W1)=1.E-10;VARMIN(W2)=1.E-10
VARMIN(R1)=1.E-10;VARMIN(R2)=1.E-10
    GROUP 21. Print-out of variables
NZPRIN=NZ/5
    GROUP 22. Spot-value print-out
IYMON=NY-2;TSTSWP=LITHYD;NPLT=2;IDISPA=50
    GROUP 23. Field print-out and plot control
OUTPUT(LISV,Y,N,Y,Y,Y,Y);OUTPUT(LISW,Y,N,Y,Y,Y,Y)
OUTPUT(VMSW,Y,N,Y,Y,Y,Y);OUTPUT(VMSV,Y,N,Y,Y,Y,Y)
tstswp=-1
    GROUP 24. Dumps for restarts
libref=218    
idispa=1
#maxabs
#endpause
isg51=2