GROUP 1. Run title and other preliminaries
   PHOTON USE
   AUTOPLOT
   file
   phi 5
 
   cl
   msg POWER-LAW-FLUID LAMINAR PIPE FLOW
   msg Reynolds number = 10 Power index = 0.5
   msg Velocity (W1) profile
   msg Blue line --- PHOENICS solution
   msg crosses ---   analytical solution
   da 1 w1;da 1 w1a
   col3 1;blb4 2
   msg press  to continue
   pause
 
   msg press  to end
   pause
   end
   END_USE
TEXT(Power-Law Fluid FD Lam Pipe Flow  
TITLE
 
  DISPLAY
  The problem concerns the steady fully-developed laminar flow of
  a power-law pseudo-plastic non-Newtonian fluid. The apparent
  viscosity of such a fluid is given by:
 
    enul = K*[(dw/dy)**(n-1.)]/rho
 
  so that enul decreases with increasing shear rate. Here, K is
  the fluid consistency index and n the flow-behaviour index.
  Examples of pseudo-plastic fluids include rubber solutions,
  adhesives, polymer solutions or melts, and biological fluids.
  ENDDIS
 
  For fully-developed flow the analytical solution for the pressure
  drop is given by:
 
    dp/dz = 4.*rho*win**2*[(2+6n)/n]**n/Re/D
 
  where D is the pipe diameter, win the mean velocity and Re the
  power-law Reynolds number, defined by:
 
    Re = D**n*win**(2-n)*rho/K
 
  The analytical solution for the velocity profile is
 
    w = win*(1+3*n)[1-(2r/D)**(1.+1./n)]/(1.+n)
 
  The problem is solved by use of the single-slab solver.
 
 
  ** GXPRPS=T activates ENUL coding via the file GXPRPS
     for which 
BOOLEAN(GXPRPS);GXPRPS=F
REAL(RIN,DIN,WIN,AIN,DPDZ,FLOWIN);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)
    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,POWER);REY=10.;POWER=0.5
ENULA=WIN**(2.0-POWER)*DIN**POWER/REY;ENULB=POWER
ENULA
ENULB
   ** enul = enula*(dwdy)*[0.5*(enulb-1)]
ENUL=STRAIN;DWDY=T
DPDZ=4.*RHO1*WIN**2/DIN/REY*((2.+6.*POWER)/POWER)**POWER
DPDZ
REY
POWER
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN
  ** compute analytical solutions
REAL(WA,GR,GRAT,POW2);INTEGER(JJM1)
POW2=(1.+POWER)/POWER
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;GRAT=GR/RIN
+WA=WIN*(1.+3.*POWER)*(1.-GRAT**POW2)/(1.+POWER)
+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)
   ** activate fully-developed-flow pressure adjustment
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=8;LITHYD=10
    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 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
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. GRND4 1. 1. 0
  0.04472   0.5
ENDIF