TEXT(2D 2-PHASE HILLS BUBBLE COLUMN
TITLE
  DISPLAY
    The case considered is the bubble column studied experimentally
    by Hills [Instn.Chem.Engrs., Vol.52, p1, 1974], which has a
    diameter of 0.138m and a height of 1.37m. The column is
    initially filled with water, and air enters uniformly at the
    bottom with a superficial velocity of 0.038m/s. After about 25s
    steady-state conditions are obtained in which clockwise liquid
    circulation is observed, with upflow at the column centre and
    downflow at the outer wall. The task is to predict the gas
    holdup in the column, and to compare the predicted void-fraction
    and vertical liquid-velocity radial profiles at z=0.6m(i.e IZ=9)
    with the measurements. The calculation may be performed with one
    of 3 turbulence models, i.e. the Petersen SGS model, the Rice-
    Geary model, or a modified k-e model which accounts for bubble-
    induced turbulence production. The two-phase model accounts for
    interfacial drag, lift, pressure and virtual-mass forces.
  ENDDIS
CHAR(CTURB)
MESG( Enter the required turbulence model:
MESG(  SGS   - Petersen subgrid-scale model
MESG(  RICE  - Rice-Geary mixing-length model
MESG(  KE    - k-e model  (default)
MESG(
READVDU(CTURB,CHAR,KE)
REAL(AREA,DIAMB,DIAMC,DTT,EMULIQ,EPIN,FLOWG,GRAVAC,PI)
REAL(RADC,RGAS,RLIQ,RGINIT,STEN,TKEIN,VSLIPM,WSGAS,ZLENC)
INTEGER(NT,TMODEL)
BOOLEAN(LIFT,VMAS,INTP);LIFT=T;VMAS=T;INTP=T
DIAMC=0.138;RADC=0.5*DIAMC;ZLENC=1.37;RGAS=0.14;RLIQ=1.-RGAS
PI=3.14159;GRAVAC=-9.81;STEN=0.072;EMULIQ=1.E-3;RGINIT=1.E-4
WSGAS=0.038;DIAMB=7.61E-3;TKEIN=(0.05*WSGAS)**2
EPIN=0.1643*TKEIN**1.5/(0.1*RADC)
    GROUP 2. Transience; time-step specification
  ** The calculation should be run transient, and at least 500
     time steps are required for the evolution of steady-flow
     conditions.
STEADY=F
DTT=0.05;NT=5;TLAST=DTT*NT;GRDPWR(T,NT,TLAST,1.0)
    GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1;AREA=XULAST*0.5*RADC*RADC
    GROUP 4. Y-direction grid specification
NY=10;GRDPWR(Y,NY,RADC,1.0)
    GROUP 5. Z-direction grid specification
NZ=20;GRDPWR(Z,-NZ,ZLENC,1.4)
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
ONEPHS=F;SOLVE(P1,V1,V2,W1,W2,R1,R2)
SOLUTN(P1,Y,Y,Y,P,P,P)
  ** deactivate harmonic averaging
SOLUTN(V1,P,P,P,P,P,N);SOLUTN(V2,P,P,P,P,P,N)
SOLUTN(W1,P,P,P,P,P,N);SOLUTN(W2,P,P,P,P,P,N)
SOLUTN(R1,P,P,P,P,P,N);SOLUTN(R2,P,P,P,P,P,N)
STORE(VREL,CFIP,ENUT,LEN1,CD,REYN,WEB)
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
PRT(R1)=0.5;PRT(R2)=0.5;RHO1=1.E3;RHO2=1.23
CASE :CTURB: OF
WHEN SGS,3
+ MESG(Petersen sgs turbulence model
+ TURMOD(SGSMOD);TMODEL=1
WHEN RICE,4
+ MESG(Rice-Geary mixing-length turbulence model
+ TURMOD(MIXLEN-RICE);TMODEL=2
WHEN KE,2
+ MESG(k-e turbulence model
+ TURMOD(KEMODL);TMODEL=3
+ FIINIT(KE)=TKEIN;FIINIT(EP)=EPIN
+ PRT(R1)=0.75;PRT(R2)=0.75
ENDCASE
ENUL=EMULIQ/RHO1
    GROUP 10. Inter-phase-transfer processes and properties
CFIPS=GRND7;CFIPD=4.0;VSLIPM=1.E-4;CFIPA=VSLIPM
CFIPB=DIAMB;CFIPC=STEN;RLOLIM=1.E-3
IF(VMAS) THEN
+ CVM=GRND2;CVMA=0.5;STORE(VMSW,VMSV)
  ** penultimate argument activates spot value history
+ OUTPUT(VMSW,N,P,Y,P,N,P);OUTPUT(VMSV,N,P,Y,P,N,P)
ENDIF
IF(LIFT) THEN
+ CLIFT=GRND2;CLIFTA=-0.5
+ INTSOR(LIFT,CLIFT,CLIFTA,RELAX,0.1)
+ OUTPUT(LISW,N,P,P,P,N,P);OUTPUT(LISV,N,P,P,P,N,P)
ENDIF
IF(INTP) THEN
+ CPIP=GRND2;CPIPA=0.25
+ INTSOR(INTPL,CPIP,CPIPA,RELAX,0.1)
+ OUTPUT(IPSW,N,P,P,P,N,P);OUTPUT(IPSV,N,P,P,P,N,P)
ENDIF
    GROUP 11. Initialization of variable or porosity fields
FIINIT(R2)=RGINIT;FIINIT(R1)=1.0-FIINIT(R2);INIADD=F
    GROUP 12. Patchwise adjustment of terms
    GROUP 13. Boundary conditions and special sources
  ** Gas-phase mass-inflow boundary
FLOWG=RHO2*WSGAS
PATCH(IN,LOW,1,NX,1,NY,1,1,1,LSTEP)
COVAL(IN,P2,FIXFLU,FLOWG);COVAL(IN,W2,ONLYMS,WSGAS)
  ** Gas-phase outflow boundary
PATCH(OUTG,HIGH,1,NX,1,NY,NZ,NZ,1,LSTEP)
COVAL(OUTG,P2,RHO2*1.E10,0.0)
  ** Liquid-phase 'pressure-relief' outflow cell
PATCH(OUTL,HIGH,1,NX,NY,NY,NZ,NZ,1,LSTEP)
COVAL(OUTL,P1,RHO1,0.0)
  ** Gravititional force relative to liquid phase
PATCH(GRAVITY,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
COVAL(GRAVITY,W2,FIXFLU,GRAVAC*(1.-RHO1/RHO2))
  ** Wall friction
WALL(NWALL,NORTH,1,NX,NY,NY,1,NZ,1,LSTEP)
  ** Bubble-induced turbulence production
IF(TMODEL.EQ.3) THEN
+ PATCH(KEDI,CELL,1,NX,1,NY,1,NZ,1,LSTEP);EL1A=0.01
+ COVAL(KEDI,KE,FIXFLU,GRND3);COVAL(KEDI,EP,FIXFLU,GRND3)
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=100;TSTSWP=-10
    GROUP 16. Termination of iterations
SELREF=T;RESFAC=0.01
    GROUP 17. Under-relaxation devices
RELAX(V1,FALSDT,DTT);RELAX(V2,FALSDT,DTT)
RELAX(W1,FALSDT,DTT);RELAX(W2,FALSDT,DTT)
RELAX(R1,LINRLX,0.4);RELAX(R2,LINRLX,0.4)
RELAX(CFIP,LINRLX,0.3)
IF(TMODEL.EQ.3) THEN
+ RELAX(KE,LINRLX,0.3);RELAX(EP,LINRLX,0.3);KELIN=1
ENDIF
    GROUP 18. Limits on variables or increments to them
VARMIN(R1)=1.E-12;VARMIN(R2)=1.E-12;LITER(P1)=30
    GROUP 22. Spot-value print-out
IYMON=1;IZMON=9
SPEDAT(SET,GXMONI,TRANSIENT,L,F)
    GROUP 23. Field print-out and plot control
NPRINT=LSWEEP;NZPRIN=1;NYPRIN=1;IZPRF=11
NTPRIN=NT/4;NPLT=10
CSG1=H;IDISPA=NT/4
    GROUP 24. Dumps for restarts