AUTOPLOT USE
  file
  phi 5
 
 
  da 1 u2; da 1 uterm;screen
  msg 2nd-phase velocity & terminal velocity;col3 1;blb4 2
  msg press  to continue
  pause
  cl;screen
  msg press  to continue
  msg press e to END
  enduse
 
TEXT(1D SETTLING OF PARTICLES UNDER GRAVITY
TITLE
  DISPLAY
  The case considered is that of a particles falling through
  essentially stagnant air under gravity. The flow is steady, one
  dimensional and after a relatively short period of acceleration,
  the particles attain a uniform velocity under which the drag
  force balances the gravitational force. This velocity is called
  the terminal settling velocity, which is given by:
 
    Vset = [ 4.*{rho,p-rho,a}*g*diam,p/(3.*rho,a*Cd) ] ** 0.5
 
  The density ratio is taken as 1000 and the particle Reynolds
  number as 2,000. The calculation may be performed in the x-,
  y- or z-direction with phase 1 as continuous (air) and phase
  2 as dispersed (solid) or vice versa (CONPHS=2). The task is
  to calculate the terminal velocity of the particles.
  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
  * 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,UGAS,UINP,UIN1,UIN2,R1IN,R2IN,REYP,DIAMP,DENRAT)
REAL(USET,CD,RINP,VOLFP,FLOW1,FLOW2,RHOG,RHOS,ACON)
REYP=2000;UGAS=1.E-3;UINP=1.E-3;DENRAT=1000.;ENULA=1.E-6
RINP=1.E-3;VOLFP=RINP*UINP;RHOG=1.0;RHOS=RHOG*DENRAT
  ** calculate expected particle settling velocity
ACON=4.*(DENRAT-1.)*9.81/(3.*ENULA**2)
CD=24.*(1.+.15*REYP**.687)/REYP+.42/(1.+4.25E4/(REYP**1.16))
DIAMP=(REYP*REYP*CD/ACON)**(1./3.);USET=REYP*ENULA/DIAMP
cd
uset
diamp
IF(CONPHS.EQ.1) THEN
+ UIN1=UGAS;UIN2=UINP;R2IN=RINP;R1IN=1-R2IN
ELSE
+ UIN1=UINP;UIN2=UGAS;R1IN=RINP;R2IN=1-R1IN
ENDIF
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
XLEN=100.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)
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
IF(CONPHS.EQ.1) THEN
+ RHO1=RHOG;RHO2=RHOS
ELSE
+ RHO2=RHOG;RHO1=RHOS
ENDIF
IF(VARLAM) THEN
    Use ENUL=ENULA+ENULB*TMP1
+ STORE(TMP1);ENUL=LINTEM;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 = particle size
CFIPA=1.E-3;CFIPB=DIAMP
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:CH2:)=UIN1;FIINIT(:CH3:)=UIN2
IF(CONPHS.EQ.1) THEN
+ FIINIT(R2)=VOLFP/USET;FIINIT(R1)=1.-FIINIT(R2)
ELSE
+ FIINIT(R1)=VOLFP/USET;FIINIT(R2)=1.-FIINIT(R1)
ENDIF
FIINIT(UTERM)=USET
    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*USET*1.E2,0)
ELSE
+ COVAL(OUT,P1,RHO1*USET*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.-RHOG/RHOS))
ELSE
+ COVAL(GRAVITY,:CH2:,FIXFLU,9.81*(1.-RHOG/RHOS))
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=3.*XLEN/USET/N:CH1:
RELAX(R1,LINRLX,0.5);RELAX(R2,LINRLX,0.5)
RELAX(:CH2:,FALSDT,DTF);RELAX(:CH3:,FALSDT,DTF)
    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
NXPRIN=1;NYPRIN=1;NZPRIN=1
    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