GROUP 1. Run title and other preliminaries
 
TEXT(Power-Law Fluid_2D FD Duct Flow    
TITLE
 
  The case considered is three-dimensional fully-developed,
  laminar flow of an incompressible power-law fluid in a
  rectangular duct of aspect ratio alfa=0.5, where alfa=
  height/width. The calculations are performed for a power-law
  Reynolds number Re of 100 and a fluid with a power-law index
  of 0.5. The solution exploits symmetry by performing the
  calculation over one quadrant of the duct cross section only.
  The flow Reynolds number is based on the hydraulic diameter
  of the duct.
 
  Data on the friction factor ( and hence pressure drop ) has
  been reported in 'Non-Newtonian Fluid Laminar Flow and Forced
  Convection Heat Transfer in Rectangular Ducts', S.X.Gao and
  J.P.Hartnett, Int.Comm. Heat Mass Transfer, Vol.19, p673, (1992).
  The data on the friction factor f is correlated by
 
    f*Re = 2.**(3*n+1)*(B+A/n)**n
 
  where n is the power-law index and A and B are constants which
  depend on the duct aspect ratio alfa. For the duct aspect ratio
  considered here A=0.244 and B=0.7276.
 
  For a Reynolds number of 100 and a power-law index of 0.5, the
  data indicates that f=0.0624 so that dp/dz=0.0935. The PHOENICS
  predictions yield values of f=0.061  and dp/dz=0.091, which are
  in very good agreement with the experimental values.
 
REAL(HEIGHT,WIDTH,ALF,HD2,WD2,WIN,DPDZ,REY,FRIC)
REAL(AIN,DHYD,FLOWIN)
   ** ALFA = HEIGHT/WIDTH
WIDTH=2.0;ALF=0.5;HEIGHT=ALF*WIDTH
HD2=0.5*HEIGHT;WD2=0.5*WIDTH
WIN=1.0
REY=100.;DHYD=4.*WIDTH*HEIGHT/(2.*HEIGHT+2.*WIDTH)
    GROUP 3. X-direction grid specification
AIN=HD2*WD2
XULAST=WD2;NX=20;GRDPWR(X,NX,XULAST,1.0)
    GROUP 4. Y-direction grid specification
YVLAST=HD2;NY=20;GRDPWR(Y,NY,YVLAST,1.0)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,W1)
    GROUP 8. Terms (in differential equations) & devices
TERMS(W1,N,P,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
ENUT=0.0
FLOWIN=RHO1*WIN*AIN
REAL(POWER);POWER=0.5
ENULA=WIN**(2.0-POWER)*DHYD**POWER/REY
ENULB=POWER
   ** enul = enula*(dwdy)*[0.5*(enulb-1)]
ENUL=STRAIN;GENK=T;STORE(VISL)
 
REAL(ACON,BCON);ACON=.244;BCON=.7276
FRIC=(2.**(3.*POWER+1.))*(BCON+ACON/POWER)**POWER/REY
DPDZ=2.*FRIC*RHO1*WIN**2/DHYD
FRIC
DPDZ
REY
POWER
 
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN
    GROUP 12. Convection and diffusion adjustments
PATCH(GP12CONH,CELL,1,NX,1,NY,1,NZ,1,1)
COVAL(GP12CONH,U1,0.0,0.0);COVAL(GP12CONH,V1,0.0,0.0)
COVAL(GP12CONH,W1,0.0,0.0)
    GROUP 13. Boundary conditions and special sources
PATCH(WALLT,NWALL,1,NX,NY,NY,1,NZ,1,1)
COVAL(WALLT,W1,1.0,0.0);COVAL(WALLT,U1,1.0,0.0)
 
PATCH(WALLS,EWALL,NX,NX,1,NY,1,NZ,1,1)
COVAL(WALLS,W1,1.0,0.0);COVAL(WALLS,V1,1.0,0.0)
 
PATCH(RELIEF,CELL,NX/2,NX/2,NY/2,NY/2,1,NZ,1,1)
COVAL(RELIEF,P1,FIXP,0.0)
 
FDFSOL=T;USOURC=T
PATCH(FDFW1DP,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFW1DP,W1,FLOWIN,GRND1)
 
    GROUP 15. Termination of sweeps
LSWEEP=20;LITHYD=2;LITER(W1)=15
    GROUP 16. Termination of iterations
RESREF(P1)=1.E-12*WIN*AIN
RESREF(W1)=1.E-12*DPDZ*ZWLAST*AIN
RESREF(U1)=RESREF(W1);RESREF(V1)=RESREF(W1)
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=50.*(YVLAST/NY)**2/ENULA
RELAX(U1,FALSDT,DTF);RELAX(V1,FALSDT,DTF)
RELAX(W1,FALSDT,DTF)
    GROUP 22. Spot-value print-out
IXMON=NX-2;IYMON=NY-2;TSTSWP=-1
    GROUP 23. Field print-out and plot control
NPLT=1;NYPRIN=1;NZPRIN=1
    GROUP 24. Dumps for restarts