GROUP 1. Run title and other preliminaries
TEXT(CK KE_1D PIPE POWER-LAW FLUID :T304
TITLE
  DISPLAY
  The problem concerns the steady fully-developed turbulent 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.
  The turbulence may be simulated by means of one of the following
  low-Reynolds-number two-equation models: the 2-layer, Chen-Kim
  or Lam-Bremhorst k-e models, or the k-omega model.
  ENDDIS
   AUTOPLOT USE
    FILE
    PHI 5
 
    D 1 W1;COL3 1
    PAUSE
    CL;D 1 KE;COL6 1
    PAUSE
    CL;D 1 ENUT;D 1 VISL;COL6 1;COL9 2
    PAUSE
    CL;D 1 EP;COLF 1
    PAUSE
   ENDUSE
 
INTEGER(TMODEL);CHAR(CTURB,TLSC)
REAL(RIN,DIN,WIN,AIN,DPDZ,FLOWIN,REYST,REYPLS,AN,REYFAC)
REAL(FRIC,A1,A2,A3,AA,BB,GR,EPSPLS,TKEIN,EPSIN,US,GYPLUS)
REAL(GWPLUS,KFAC,RECN,TMN,YPLSW,DELT1,DELY)
INTEGER(JJM1,JJM)
RIN=0.1;DIN=2.*RIN;WIN=1.0;AIN=RIN*RIN/2.
MESG(Enter Reynolds number - default 1.158E4
READVDU(REYST,REAL,1.158E4)
MESG(Enter Power Index - default 0.525
READVDU(AN,REAL,0.525)
RECN=1./AN;TMN=2.-AN
REYFAC=8.**(AN-1.); REYFAC=REYFAC*(0.25/AN+0.75)**AN
REYPLS=REYST*REYFAC
ENULA=WIN**(2.0-AN)*DIN**AN/REYPLS
ENULB=AN
FRIC=1./(1.82*LOG10(REYST)-1.64)**2
DPDZ=FRIC*RHO1*WIN*WIN/(2.*DIN);US=WIN*(FRIC/8.)**0.5
A1=4./AN**0.75;A2=0.4/AN**1.2;A3=1.-0.5*AN
DO JJ=1,20
+ BB=FRIC**A3;BB=REYST*BB;AA=LOG10(BB)
+ BB=A1*AA-A2;FRIC=1./(BB*BB)
ENDDO
FRIC
DPDZ=FRIC*2.*RHO1*WIN*WIN/DIN;US=WIN*(FRIC/2.)**0.5
DPDZ
  ** estimate initial values from K+=2 & EP+=1./(ak*Y+)
EPSPLS=1./(0.41*100.); RIN=0.5*DIN
TKEIN=2.*US*US;EPSIN=EPSPLS*US**4/ENULA
MESG( Enter the required turbulence model:
MESG(  LAMB - Low-Re Lam-Brem. k-e model
MESG(  CHEN - Low-Re Chen-Kim  k-e model
MESG(  TWOL - Low-Re 2-layer   k-e model
MESG(  KO   - Low-Re k-omega model (default)
MESG(
READVDU(CTURB,CHAR,CHEN)
CASE :CTURB: OF
WHEN LAMB,4
+ TEXT(LB KE_1D PIPE POWER-LAW FLUID :T304
+ MESG(Low-Re Lam-Bem. k-e turbulence model
+ TMODEL=1;TLSC=EP
WHEN CHEN,4
+ MESG(Low-Re Chen-Kim k-e turbulence model
+ TMODEL=2;TLSC=EP
WHEN TWOL,4
+ TEXT(2L KE_1D PIPE POWER-LAW FLUID :T304
+ MESG(Low-Re 2-layer k-e turbulence model
+ SELREF=F;TMODEL=3;TLSC=EP
WHEN KO,2
+ TEXT(KO KE_1D PIPE POWER-LAW FLUID :T304
+ MESG(Low-Re k-omega turbulence model
+ TMODEL=4;TLSC=OMEG
+ STORE(EP);EPSIN=EPSIN/(0.09*TKEIN)
ENDCASE
    GROUP 2. TRANSIENCE; TIME-STEP SPECIFICATION
CARTES=F
    GROUP 4. Y-direction grid specification
  NY=50;GRDPWR(Y,NY,RIN,0.8)
  ** define first dely from wall
MESG(Enter Grid KFAC - default 1.05
READVDU(KFAC,REAL,1.05)
YPLSW=1.0
DELT1=YPLSW*ENULA/(US**TMN);DELT1=(DELT1)**RECN;DELY=DELT1/RIN
  ** calculate NY from dely & Kfac
AA=(YVLAST/DELY)*(KFAC-1.0)+1.0;AA=LOG(AA)/LOG(KFAC)+1.0001
NY=AA
  ** define uniform grid initially
IREGY=1;GRDPWR(Y,NY,YVLAST,1.0)
  ** compute expanding grid from north boundary
YFRAC(NY)=1.0
DO JJ=NY,2,-1
+ JJM=JJ-1
+ YFRAC(JJM)=YFRAC(JJ)-DELY
+ DELY=KFAC*DELY
ENDDO
YVLAST=RIN
    GROUP 5. Z-direction grid specification
    GROUP 7. Variables stored, solved & named
SOLVE(W1);STORE(VISL)
CASE (TMODEL) OF
WHEN 1
+ TURMOD(KEMODL-LOWRE);KELIN=3;STORE(REYN,FONE,FMU,FTWO)
WHEN 2
+ TURMOD(KECHEN-LOWRE);KELIN=3;STORE(FMU,REYN,FONE)
WHEN 3
+ TURMOD(KEMODL-2L);KELIN=3;STORE(FMU,FTWO,REYN,REYT)
WHEN 4
+ TURMOD(KOMODL-LOWRE);STORE(FMU,REYT,FONE,FTWO)
ENDCASE
STORE(ENUT,GEN1)
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
ENULA=WIN**(2.0-AN)*DIN**AN/REYPLS
ENULB=AN
   ** enul = enula*(dwdy)*[0.5*(enulb-1)]
ENUL=STRAIN
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:TLSC:)=EPSIN;FIINIT(KE)=TKEIN
REAL(YPLS1);YPLS1=11.5**AN
  ** use log-law for initial W profile
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
+GYPLUS=YVLAST*(1.-GR)*US/ENULA
+GWPLUS=LOG(GYPLUS)*RECN/0.41+5.25*RECN
IF(GYPLUS.LE.YPLS1) THEN
+ GWPLUS=GYPLUS**RECN
ENDIF
+INIT(IN:JJ:,W1,ZERO,GWPLUS*US)
ENDDO
    GROUP 13. Boundary conditions and special sources
WALL(TOPW,NORTH,1,NX,NY,NY,1,NZ,1,1)
   ** activate fully-developed-flow pressure adjustment
FLOWIN=RHO1*WIN*AIN
FDSOLV(FLOW,FLOWIN)
    GROUP 15. Termination of sweeps
LSWEEP=120;LITHYD=5
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=0.1*ZWLAST/WIN
RELAX(W1,FALSDT,DTF)
IF(TMODEL.EQ.4) THEN
+ RELAX(W1,FALSDT,DTF)
+ RELAX(KE,FALSDT,DTF); RELAX(OMEG,FALSDT,DTF)
ELSE
+ RELAX(KE,LINRLX,0.3); RELAX(EP,LINRLX,0.7)
ENDIF
    GROUP 22. Spot-value print-out
IYMON=NY-4;TSTSWP=-1
    GROUP 23. Field print-out and plot control
ITABL=3;NPLT=1;NYPRIN=1;NZPRIN=1
    GROUP 24. Dumps for restarts
WALPRN=T;OUTPUT(ENUT,Y,N,Y,N,Y,N)