PHOTON USE
   AUTOPLOT
   file
   phi 5
 
   cl
   msg Bingham-fluid pipe flow; Re = 10; Yield No = 2.5
   msg
   msg W1 profile: Blue line from PHOENICS; crosses from analysis
   da 1 w1;da 1 w1a
   col3 1;blb4 2
   pause
   END_USE
  DISPLAY
 
    GROUP 1. Run title and other preliminaries
 
TEXT(Bingham-Fluid FD Lam Pipe Flow     
TITLE
 
  The problem concerns the steady fully-developed laminar flow of
  a Bingham-plastic non-Newtonian fluid in a pipe. This type of
  fluid remains rigid when the shearing stress is less than the
  yield stress tauy and flows somewhat like a Newtonian fluid when
  the shearing stress exceeds tauy.
 
  The apparent viscosity of such a fluid is given by the following
  two-paramter formula:
 
    enul = [n + tauy/(dw/dy)]/rho    for tau >  tauy
 
    enul = infinity                  for tau << tauy
 
  where n is the rigidity coefficient, and tau is the shear stress.
  Examples of fluids which behave as, or nearly as, Bingham
  plastics include water suspensions of clay, sewage sludge, some
  emulsions and thickened hydrocarbon greases.
mesg(press RETURN to continue
readvdu(nphi,int,nphi)
  For fully-developed flow the approximate analytical solution for
  the pressure drop is given by:
 
    dp/dz = 32.*rho*win**2*[1+Y/6]/Re/D
 
  where D is the pipe diameter, win the mean velocity, Re the
  Bingham Reynolds number, defined by:
 
    Re = rho*win*D/n
 
  and Y the yield number, defined by:
 
    Y = D*tauy/(win*n)
 
  The approximate analytical solution for the velocity profile
  is coded below in PIL and provided on the RESULT and PHI/PHIDA
  for comparison with the PHOENICS solution.
mesg(press RETURN to continue
readvdu(nphi,int,nphi)
 
  The problem is solved by use of the single-slab solver.
 
  ** GXPRPS=T activates ENUL coding via the file GXPRPS
     for which 
  enddis
BOOLEAN(GXPRPS);GXPRPS=F
REAL(RIN,DIN,WIN,AIN,DPDZ);RIN=0.1;DIN=2.*RIN
WIN=1.0;AIN=RIN*RIN/2.
    GROUP 2. Transience; time-step specification
CARTES=F
    GROUP 4. Y-direction grid specification
NY=20;GRDPWR(Y,NY,RIN,1.0)
    GROUP 5. Z-direction grid specification
ZWLAST=0.1
    GROUP 7. Variables stored, solved & named
SOLVE(W1);STORE(W1A,VISL,BTAU,GEN1)
    GROUP 8. Terms (in differential equations) & devices
TERMS(W1,N,N,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
RHO1=1.0;ENUT=0.
REAL(REY,YIELDN,FLOWIN);REY=10.;YIELDN=2.5
ENULA=DIN*WIN/REY
ENULB=YIELDN*WIN*ENULA/DIN
ENUL=BINGHAM;DWDY=T
DPDZ=32.*RHO1*WIN**2*(1.+YIELDN/6.)/REY/DIN
DPDZ
REY
YIELDN
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN
  ** compute analytical solutions
REAL(WA,GR,GRAT2,ACON,GRP);INTEGER(JJM1)
ACON=YIELDN*WIN/DIN
ACON
GRP=2.*YIELDN*WIN**2*RHO1/REY/DPDZ
GRP
DO JJ=1,NY
+PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1)
+GR=0.5*YFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1
+GR=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1))
ENDIF
+GR=GR*YVLAST
IF(GR.LE.GRP) THEN
+ GR=GRP
ENDIF
+GRAT2=(GR/DIN)**2
+WA=8.*WIN*(1.+YIELDN/6.)*(0.25-GRAT2)-ACON*(RIN-GR)
+INIT(IN:JJ:,W1A,ZERO,WA)
ENDDO
    GROUP 13. Boundary conditions and special sources
PATCH(WALL,NWALL,1,NX,NY,NY,1,NZ,1,1);COVAL(WALL,W1,1.0,0.0)
 
FDFSOL=T;USOURC=T
FLOWIN=RHO1*WIN*AIN
PATCH(FDFW1DP,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFW1DP,W1,FLOWIN,GRND1)
    GROUP 15. Termination of sweeps
LSWEEP=10;LITHYD=30
    GROUP 16. Termination of iterations
RESREF(W1)=1.E-12*DPDZ*ZWLAST*AIN
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=100.*(YVLAST/NY)**2/ENULA
RELAX(W1,FALSDT,DTF)
    GROUP 19. Data communicated by satellite to GROUND
FLOWIN=RHO1*WIN*AIN
    GROUP 21. Print-out of variables
OUTPUT(VISL,Y,N,N,Y,Y,Y)
    GROUP 22. Spot-value print-out
IYMON=NY;TSTSWP=-1
    GROUP 23. Field print-out and plot control
NPLT=1;NYPRIN=1;NZPRIN=1
    GROUP 24. Dumps for restarts
ENULA
ENULB
  
IF(GXPRPS) THEN

STORE(PRPS);FIINIT(PRPS)=33
  ** mat no. rho enul cp kond expan
  **   1        air
CSG10=Q1
  MATFLG=T;NMAT=1
  33 1. GRND5 1. 1. 0
  0.02   0.25
ENDIF
DISTIL=T