AUTOPLOT USE
  file
  phi 5
 
 
  da 1 u2; da 1 uterm
  screen
  col3 1;blb4 2;scale y 0 .2;redr
  msg 2nd-phase velocity
  msg blue PHOENICS   + analytical
  msg press  to continue
  pause
  cl
  screen
  msg press  to continue
  msg press e to END
  enduse
 
TEXT(1D RISING OF AIR BUBBLES IN WATER
TITLE
  DISPLAY
  The case considered is that of a cluod of air bubbles rising
  through essentially stagnant water. The flow is 1D steady and
  after a relatively short period of acceleration, the bubbles
  attain a uniform velocity under which the drag force balances
  the buoyancy force. This velocity is called the terminal
  rise velocity, which is given by:
 
    Vrise = [ 4.*{rho,l-rho,b}*g*diam,b/(3.*rho,l*Cd) ] ** 0.5
 
  The calculation may employ any one of three built-in bubble-drag
  correlations, as selected by CFIPD=4., 5. OR 6. with CFIPS=GRND7.
  CFIPD=4. allows for distortion of the bubble shape; CFIPD=5.
  assumes spherical bubbles in dirty water; and CFIPD=6. assumes
  ellipsoidal bubbles in clean water. The calculation may be
  performed in the x-, y- or z-direction with phase 1 as continuous
  (water) and phase 2 as dispersed (air) or vice versa (CONPHS=2).
  The task is to calculate the terminal rise velocity of the
  bubbles.
  ENDDIS
 
  * CONPHS=1 selects 1st phase as the continuous ( gas ) phase
             and thus CFIPS=GRND7
  *       =2 selects 2nd phase as the continuous ( gas ) phase
             and thus CFIPS=GRND8
  * VARLAM=T tests coding in GXCFIP when ENUL is variable
BOOLEAN(VARLAM);INTEGER(CONPHS);VARLAM=T;CONPHS=1
CFIPD=4.
  * reg1, reg2, etc identify each of the regions of the CFIPD=4.
                    interphase drag correlation; for each
                    calculation one region is selected to test
                    the appropriate coding in GXFIPS.
BOOLEAN(REG1,REG2,REG2B,REG4,REG5)
REG1=F;REG2=T;REG2B=F;REG4=F;REG5=F
  * for 1dx calculation set: CH1=X;CH2=U1;CH3=U2
  * for 1dy calculation set: CH1=Y;CH2=V1;CH3=V2
  * for 1dz calculation set: CH1=Z;CH2=W1;CH3=W2
CHAR(CH1,CH2,CH3);CH1=X;CH2=U1;CH3=U2
   CH1=Y;CH2=V1;CH3=V2
   CH1=Z;CH2=W1;CH3=W2
REAL(XLEN,ULIQ,UINB,UIN1,UIN2,R1IN,R2IN,REYB,DIAMB,DENRAT)
REAL(URISE,CD,RINB,VOLFB,FLOW1,FLOW2,RHOL,RHOG,ENULIQ)
REAL(WEBER,WEBMAX,EOTVOS,MORTON,STEN,DELRHO,ACON)
ENULIQ=1.E-6;ULIQ=1.E-5;UINB=1.E-3;DENRAT=1000.
RINB=0.01;VOLFB=RINB*UINB;RHOG=1.0;RHOL=DENRAT*RHOG
STEN=0.072;DELRHO=(RHOL-RHOG)
ACON=4.*9.81*DELRHO/(3.*RHOL*ENULIQ**2)
  ** calculate expected particle settling velocity
IF(CFIPD.EQ.4.) THEN
IF(REG1) THEN
+ REYB=0.2;CD=16./REYB
+ DIAMB=(REYB*REYB*CD/ACON)**(1./3.);URISE=REYB*ENULIQ/DIAMB
ENDIF
IF(REG2) THEN
+ REYB=50.;CD=20.68/REYB**0.643
+ DIAMB=(REYB*REYB*CD/ACON)**(1./3.);URISE=REYB*ENULIQ/DIAMB
ENDIF
IF(REG2B) THEN
+ REYB=150.;CD=6.3/REYB**.385
+ DIAMB=(REYB*REYB*CD/ACON)**(1./3.);URISE=REYB*ENULIQ/DIAMB
ENDIF
IF(REG4) THEN
+ WEBER=4.0;CD=WEBER/3.
+ DIAMB=((0.75*WEBER*STEN*CD)/(9.81*DELRHO))**0.5
+ URISE=((WEBER*STEN)/(RHOL*DIAMB))**0.5
+ REYB =URISE*DIAMB/ENULIQ
ENDIF
IF(REG5) THEN
+ WEBER=10.0;CD=8./3.
+ DIAMB=((0.75*WEBER*STEN*CD)/(9.81*DELRHO))**0.5
+ URISE=((WEBER*STEN)/(RHOL*DIAMB))**0.5
+ REYB =URISE*DIAMB/ENULIQ
ENDIF
WEBER=RHOL*URISE**2*DIAMB/STEN;WEBMAX=2065.1/WEBER**2.6
ENDIF
 
IF(CFIPD.EQ.5.) THEN
+ REYB=150.;CD=6.3/REYB**.385
+ DIAMB=(REYB*REYB*CD/ACON)**(1./3.);URISE=REYB*ENULIQ/DIAMB
ENDIF
 
IF(CFIPD.EQ.6.) THEN
+ EOTVOS=10.;CD=0.622/(1./EOTVOS+0.235*RHOL/DELRHO)
+ DIAMB=((STEN*EOTVOS)/(9.81*DELRHO))**0.5
+ URISE=(4.*DELRHO*9.81*DIAMB/(3.*RHOL*CD))**0.5
+ REYB=URISE*DIAMB/ENULIQ
ENDIF
reyb
  cd
  diamb
  urise
weber
webmax
EOTVOS=9.81*(RHOL-RHOG)*DIAMB**2/STEN
MORTON=9.81*(RHOL*ENULIQ)**4*(RHOL-RHOG)/(RHOL**2*STEN**3)
eotvos
morton
IF(CONPHS.EQ.1) THEN
+ UIN1=ULIQ;UIN2=UINB;R2IN=RINB;R1IN=1-R2IN
ELSE
+ UIN1=UINB;UIN2=ULIQ;R1IN=RINB;R2IN=1-R1IN
ENDIF
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
XLEN=1.0;GRDPWR(:CH1:,20,XLEN,1.0)
    GROUP 4. Y-direction grid specification
    GROUP 5. Z-direction grid specification
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
ONEPHS=F;SOLVE(P1,:CH2:,:CH3:,R1,R2)
  Activate storage for printout of interphase drag properties
STORE(UTERM,REYN,VREL,CD,APRJ,CFIP,SIZE)
IF(CFIPD.EQ.6) THEN
+ STORE(EOTV)
ELSE
+ STORE(WEB)
ENDIF
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
IF(CONPHS.EQ.1) THEN
+ RHO1=RHOL;RHO2=RHOG
ELSE
+ RHO2=RHOL;RHO1=RHOG
ENDIF
IF(VARLAM) THEN
    Use ENUL=ENULA+ENULB*TMP1
+ STORE(TMP1);ENUL=LINTEM;ENULA=1.E-6;ENULB=0.0;TMP1=CONST
ELSE
+ ENUL=ENULA
ENDIF
    GROUP 10. Inter-phase-transfer processes and properties
IF(CONPHS.EQ.1) THEN
+ CFIPS=GRND7;VARMIN(R2)=1.E-10
ELSE
+ CFIPS=GRND8;VARMIN(R1)=1.E-10
ENDIF
  ** CFIPA = minimum slip velocity  CFIPB = bubble size
  ** CFIPC = surface tension for gas-liquid
     CFIPD set above
CFIPA=1.E-4;CFIPB=DIAMB;CFIPC=STEN
    GROUP 11. Initialization of variable or porosity fields
IF(CONPHS.EQ.1) THEN
+ FIINIT(U1)=UIN1;FIINIT(U2)=URISE
+ FIINIT(R2)=VOLFB/URISE;FIINIT(R1)=1.-FIINIT(R2)
ELSE
+ FIINIT(U2)=UIN2;FIINIT(U1)=URISE
+ FIINIT(R1)=VOLFB/URISE;FIINIT(R2)=1.-FIINIT(R1)
ENDIF
FIINIT(UTERM)=URISE
    GROUP 12. Unused
    GROUP 13. Boundary conditions and special sources
FLOW1=RHO1*UIN1*R1IN;FLOW2=RHO2*UIN2*R2IN
INLET(IN,CELL,$1,$1,$1,$1,$1,$1,1,1)
VALUE(IN,P1,FLOW1);VALUE(IN,:CH2:,UIN1)
VALUE(IN,P2,FLOW2);VALUE(IN,:CH3:,UIN2)
OUTLET(OUT,CELL,%1,%1,%1,%1,%1,%1,1,1)
  ** estimate particle outflow velocity, and then use value
     so as to allow for slip in exit cell outflow
IF(CONPHS.EQ.1) THEN
+ COVAL(OUT,P1,RHO1*UIN1*1.E2,0);COVAL(OUT,P2,RHO2*URISE*1.E2,0)
ELSE
+ COVAL(OUT,P1,RHO1*URISE*1.E2,0);COVAL(OUT,P2,RHO2*UIN2*1.E2,0)
ENDIF
PATCH(GRAVITY,PHASEM,1,NX,1,NY,1,NZ,1,1)
IF(CONPHS.EQ.1) THEN
+ COVAL(GRAVITY,:CH3:,FIXFLU,-9.81*(1.-RHOL/RHOG))
ELSE
+ COVAL(GRAVITY,:CH2:,FIXFLU,-9.81*(1.-DENRAT))
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=100
RESREF(P1)=1.E-12*(FLOW1+FLOW2)
RESREF(R1)=1.E-12*FLOW1;RESREF(R2)=1.E-12*FLOW2
RESREF(:CH2:)=RESREF(R1)*UIN1;RESREF(:CH3:)=RESREF(R2)*UIN2
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=0.5*XLEN/URISE/N:CH1:
RELAX(R1,LINRLX,0.3);RELAX(R2,LINRLX,0.3)
RELAX(:CH2:,FALSDT,DTF);RELAX(:CH3:,FALSDT,DTF)
RELAX(CFIP,LINRLX,0.4)
    GROUP 18. Limits on variables or increments to them
    GROUP 19. Data communicated by satellite to GROUND
    GROUP 20. Preliminary print-out
    GROUP 21. Print-out of variables
OUTPUT(CFIP,P,P,P,P,Y,P);NXPRIN=1;NYPRIN=1;NZPRIN=1
OUTPUT(CD,P,P,P,P,Y,P)
    GROUP 22. Spot-value print-out
IF(:CH1:.EQ.X) THEN
+ IXMON=NX/2;IZMON=1;IYMON=1
ELSE
+ IZMON=NZ/2;IXMON=1;IYMON=1
ENDIF
IF(:CH1:.EQ.Y) THEN
+ IYMON=NY/2;IXMON=1;IZMON=1
ENDIF
TSTSWP=-1
    GROUP 23. Field print-out and plot control
    GROUP 24. Dumps for restarts