GROUP 1. Run title and other preliminaries
TEXT(1DX TRANSIENT SEDIMENTATION;
TITLE
  DISPLAY
  This 1D transient case is based on the numerical benchmark test
  documented by Hewitt et al ( see Test 2.4, Sedimentation, in
  Multiphase Science and Technology, Vol.3, Ed. G.F.Hewitt et al,
  Hemisphere Pub. Corp., p474, [1987] ). Denser fluid initially
  rests above lighter fluid in a vertical tank, and then it falls
  down under gravity whilst the lighter fluid rises to the top.
  Eventually, after about 10s, all of the dense phase will rest on
  the bottom with the lighter phase at the top. The densities of the
  two fluids are chosen so as to be nearly equal (rhod=1 and rhol=
  0.999), and gravity is defined by g=0.5*(rhod+rhol)/(rhod-rhol),
  thus allowing comparison with the analytical and numerical
  solutions reported by Hewitt et al [1987]. The task is to predict
  the dense-phase volume-fraction profiles at t=2, 4, 6 and 8s. The
  PHOENICS predictions give reasonable agreement with analytical
  solution, despite the extremely coarse grid and time-step sizes.
  ENDDIS
BOOLEAN(VMASS);REAL(GRAVAC,DTT);INTEGER(CONPHS);CONPHS=1
NX=20;ARRAY(PINIT,REAL,NX);INTEGER(NT)
    GROUP 2. Transience; time-step specification
  ** NT should be increased from 10 to 80 to cover complete
     transient period of 8 seconds
STEADY=F;NT=10;DTT=0.1;TLAST=NT*DTT
GRDPWR(T,NT,TLAST,1.0)
    GROUP 3. X-direction grid specification
GRDPWR(X,NX,2.0,1.0)
    GROUP 7. Variables stored, solved & named
ONEPHS=F;SOLVE(P1,U1,U2,R1,R2);STORE(VMSU)
    GROUP 9. Properties of the medium (or media)
IF(CONPHS.EQ.1) THEN
+ RHO1=1.0;RHO2=0.999
ELSE
+ RHO2=1.0;RHO1=0.999
ENDIF
ENUL=1.E-3
    GROUP 10. Inter-phase-transfer processes and properties
CFIPS=2.0
    GROUP 11. Initialization of variable or porosity fields
INIADD=F;REAL(RCTOP,RCBOT);RCBOT=1.E-5;RCTOP=1.-RCBOT
PATCH(BOTIN,INIVAL,1,NX/2,1,1,1,1,1,1)
PATCH(TOPIN,INIVAL,NX/2+1,NX,1,1,1,1,1,1)
IF(CONPHS.EQ.1) THEN
+ INIT(BOTIN,R1,0.0,RCBOT);INIT(BOTIN,R2,0.0,RCTOP)
+ INIT(TOPIN,R1,0.0,RCTOP);INIT(TOPIN,R2,0.0,RCBOT)
+ GRAVAC=0.5*(RHO1+RHO2)/(RHO1-RHO2)
ELSE
+ INIT(BOTIN,R2,0.0,RCBOT);INIT(BOTIN,R1,0.0,RCTOP)
+ INIT(TOPIN,R2,0.0,RCTOP);INIT(TOPIN,R1,0.0,RCBOT)
+ GRAVAC=0.5*(RHO1+RHO2)/(RHO2-RHO1)
ENDIF
 
REAL(PA,GY,DP);INTEGER(JJM1);DP=0.0
DO JJ=1,NX
+GY=0.5*XFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1;GY=XFRAC(JJM1)+0.5*(XFRAC(JJ)-XFRAC(JJM1))
ENDIF
+GY=GY*ZWLAST;DP=DP-(RHO1-RHO2)*GRAVAC*GY;PINIT(JJ)=DP
ENDDO
DO JJ=1,NX
+PATCH(IN:JJ:,INIVAL,JJ,JJ,1,1,1,1,1,1)
+PA=PINIT(JJ)-PINIT(NX);INIT(IN:JJ:,P1,ZERO,PA)
ENDDO
    GROUP 12. Unused
    GROUP 13. Boundary conditions and special sources
PATCH(GRAVITY,PHASEM,1,NX,1,1,1,1,1,LSTEP)
PATCH(TOP,CELL,NX,NX,1,1,1,1,1,LSTEP)
COVAL(TOP,P1,1.E3*RHO1,0.0);COVAL(TOP,P2,1.E3*RHO2,0.0)
IF(CONPHS.EQ.1) THEN
+ COVAL(GRAVITY,U1,FIXFLU,-GRAVAC*(1.0-RHO2/RHO1))
ELSE
+ COVAL(GRAVITY,U2,FIXFLU,-GRAVAC*(1.0-RHO1/RHO2))
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP=100
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=DTT
RELAX(R1,LINRLX,0.3);RELAX(R2,LINRLX,0.3)
RELAX(U1,FALSDT,DTF);RELAX(U2,FALSDT,DTF)
VARMIN(R1)=1.E-10;VARMIN(R2)=1.E-10
SPEDAT(SET,GXMONI,TRANSIENT,L,F)
    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
    GROUP 22. Spot-value print-out
NXPRIN=2;NTPRIN=10;TSTSWP=-1;IXMON=2
    GROUP 23. Field print-out and plot control
IPLTF=2;IPLTL=LSWEEP
    GROUP 24. Dumps for restarts
INIFLD=T;VMASS=T
IF(VMASS) THEN
+ CVM=GRND1;CVMA=1.0
IF(CONPHS.EQ.2) THEN
+ CVM=-GRND1
ENDIF
ENDIF
NTPRIN=NT/4;OUTPUT(VMSU,Y,N,Y,Y,Y,Y);NPLT=10
RHO1A=RHO1;RHO1B=0.0;RHO1=LINTEMP;STORE(TMP1)
IDISPA=NT/4;CSG1=T;SELREF=T;RESFAC=0.01;ITABL=3