TALK=T;RUN( 1, 1)
 
 ************************************************************
   Q1 created by VDI menu, Version 2020, Date 13/01/21
 CPVNAM=VDI; SPPNAM=Core
 ************************************************************
  Echo DISPLAY / USE settings
   PHOTON USE
   AUTOPLOT
   file
   phida 3
 
   clear
   msg POWER-LAW-FLUID PIPE FLOW
   msg Reynolds number = 10 Power index = 0.5
   msg Pressure (P1) profile
   msg Blue line --- PHOENICS solution
   msg crosses ---   analytical solution
   pause
   da 1 p1 y 1;da 1 pa y 1
   col3 1;blb4 2
   redr
   pause
   msg press  to continue
   clear
   msg Velocity (W1) profile
   da 1 w1 z 35;da 1 wa z 35
   col3 1;blb4 2
   pause
   msg press  to end
   pause
   end
   END_USE
 ************************************************************
 IRUNN = 1 ;LIBREF = 100
 ************************************************************
  Group 1. Run Title
 TEXT(100 2d pipe flow - Power-law fluid      )
 ************************************************************
  Echo save-block settings for Group  1
  save1begin
  This case concerns the steady laminar flow of a power-law
  pseudo-plastic non-Newtonian fluid in a circular pipe. The
  apparent viscosity of such a fluid is given by:
 
    emu = K*G^(n-1)
 
  where K is the consistency index,G is the mean strain rate,
  and n the flow-behaviour index. For pseudo-plastic fluids,
  n<1, so that emua decreases with increasing shear rate. Examples
  of pseudo-plastic fluids include rubber solutions, adhesives,
  polymer solutions or melts, and biological fluids.
 
  The bulk inlet velocity is 1m/s and the fluid density is
  100 kg/m^3. The pipe diameter and length are 0.1m and 1m,
  respectively. The flow-behaviour index is set to 0.5, and
  the consistency K is set to define a power-law Reynolds
  number of 10.
 
    The power-law Reynolds number, defined by:
 
    Re = D^n*Win^(2-n)*rho/K
 
  where D is the pipe diameter, Win the mean velocity and Rho is
  the fluid density.
 
  The task is to predict the pressure drop and fully-developed
  axial-velocity profile for a given Reynolds number, and then
  compare the results with the analytical solutions. The inform
  facility is used to compute the analytical profiles of
  pressure and velocity, and the solutions are stored in PA and
  WA, respectively.
 
  Analytical solutions for the laminar flow of Ellis fluids
  in circular pipes have been reported by: A.H.P.Skelland,
  Non-Newtonian Flow and Heat Transfer, John Wiley, (1967,and
  G.W.Govier & K.Aziz, The flow of complex mixtures in pipes,
  R.E.Kreiger Pub. Co.,Huntington, New York, (1977).
  save1end
 ************************************************************
  Group 2. Transience
 STEADY = T
 ************************************************************
  Groups 3, 4, 5  Grid Information
    * Overall number of cells, RSET(M,NX,NY,NZ,tolerance)
 RSET(M,1,20,40,8.333333E-05)
    * Cylindrical-polar grid
 CARTES=F
 ************************************************************
  Group 6. Body-Fitted coordinates
 ************************************************************
  Group 7. Variables: STOREd,SOLVEd,NAMEd
    * Non-default variable names
 NAME(142)=WDIS ;NAME(143)=SRM1
 NAME(144)=STRS ;NAME(146)=BTAU
 NAME(147)=PA ;NAME(148)=WA
 NAME(149)=GEN1 ;NAME(150)=VISL
    * Solved variables list
 SOLVE(P1,V1,W1)
    * Stored variables list
 STORE(VISL,GEN1,WA,PA,BTAU,STRS,SRM1,WDIS)
    * Additional solver options
 SOLUTN(P1,Y,Y,Y,N,N,Y)
 SOLUTN(V1,Y,Y,Y,N,N,Y)
 SOLUTN(W1,Y,Y,Y,N,N,Y)
 
 ************************************************************
  Echo save-block settings for Group  7
  save7begin
STORE(SRM1)  ! = (GEN1)^0.5
  save7end
 ************************************************************
  Group 8. Terms & Devices
 ADDDIF = T
 NEWENL = T
 ************************************************************
  Group 9. Properties
 RHO1 =100.
 ENUL = GRND4
 ENULA =3.162278 ;ENULB =0.5 ;ENULC =0.
 CP1 =1.
 DISWAL
 ENUT =1.0E-10
 ************************************************************
  Echo save-block settings for Group  9
  save9begin
REAL(RIN,DIN,WIN,DPDZ,REY,AN,POW2,BETA,ZETA,REYMR,FRIC)
 
 
ENUL=GRND4    ! Power-law rheology model
 
REY=10.;AN=0.5;RIN=0.05;DIN=2.*RIN;WIN=1.0
 
ZETA=8; BETA=(3.*AN+1.)/(4.*AN)
ZETA;BETA
ENULA=RHO1*WIN**(2.0-AN)*DIN**AN/REY;ENULB=AN
ENULA;ENULB
    * Generalised Metzner-Reed Reynolds number
REYMR=REY/(ZETA**(AN-1.)*BETA**AN)
REY;REYMR
FRIC=2.*ZETA/REYMR
FRIC
 
DPDZ=(4.*RHO1*WIN**2/DIN)/REY*((2.+6.*AN)/AN)**AN
DPDZ
DPDZ=(2.*RHO1*WIN**2/DIN)*FRIC ! Alternative calculation
DPDZ
      ** Analytical axial velocity profile
POW2=(1.+AN)/AN
(stored of WA is WIN*(1.+3.*AN)*(1.-(YG/RIN)^POW2)/(1.+AN))
      ** Analytical pressure solution
(make1 zgnz is 0)
(store1 zgnz is zg with IF(IZ.EQ.NZ))
(print zgnz is zgnz)
(stored of PA is -DPDZ*(ZG-ZGNZ))
  save9end
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
 FIINIT(W1)=1. ;FIINIT(WDIS)=0.1
 FIINIT(GEN1)=1.001E-10 ;FIINIT(VISL)=1.001E-10
   No PATCHes used for this Group
 
 
 INIADD = F
 ************************************************************
  Group 12. Convection and diffusion adjustments
   No PATCHes used for this Group
 ************************************************************
  Group 13. Boundary & Special Sources
   No PATCHes used for this Group
 
 EGWF = T
 ************************************************************
  Group 14. Downstream Pressure For PARAB
 ************************************************************
  Group 15. Terminate Sweeps
 LSWEEP = 1000
 RESREF(P1)=5.0E-16 ;RESREF(V1)=5.0E-16
 RESREF(W1)=5.0E-16
 RESFAC =1.0E-05
 ************************************************************
  Group 16. Terminate Iterations
 ************************************************************
  Group 17. Relaxation
 RELAX(P1 ,LINRLX,1. )
 RELAX(V1 ,FALSDT,0.066667 )
 RELAX(W1 ,FALSDT,0.066667 )
 RELAX(LTLS,LINRLX,1. )
 ************************************************************
  Group 18. Limits
 ************************************************************
  Group 19. EARTH Calls To GROUND Station
 GENK = T
 PARSOL = F
 ISG62 = 1
 SPEDAT(SET,OUTPUT,NOFIELD,L,T)
 SPEDAT(SET,GXMONI,PLOTALL,L,T)
 ************************************************************
  Group 20. Preliminary Printout
 DISTIL = T ;NULLPR = F
 NDST = 0
 DSTTOL =1.0E-02
 EX(P1)=628.700012 ;EX(V1)=8.242E-03
 EX(W1)=1.132 ;EX(WDIS)=0.01726
 EX(SRM1)=39.82 ;EX(STRS)=0.01616
 EX(LTLS)=3.746E-04 ;EX(BTAU)=125.900002
 EX(PA)=623.799988 ;EX(WA)=1.145
 EX(GEN1)=2614. ;EX(VISL)=8.749E-03
 ************************************************************
  Group 21. Print-out of Variables
 OUTPUT(STRS,Y,N,Y,N,Y,Y)
 OUTPUT(BTAU,Y,N,Y,N,Y,Y)
 OUTPUT(PA  ,Y,N,Y,N,Y,Y)
 OUTPUT(WA  ,Y,N,Y,N,Y,Y)
 OUTPUT(VISL,Y,N,Y,N,Y,Y)
 ************************************************************
  Group 22. Monitor Print-Out
 IXMON = 1 ;IYMON = 14 ;IZMON = 36
 NPRMON = 100000
 NPRMNT = 1
 TSTSWP = -1
 ************************************************************
  Group 23.Field Print-Out & Plot Control
 NPRINT = 100000
 NYPRIN = 1
 NZPRIN = 1
 YZPR = T
 ISWPRF = 1 ;ISWPRL = 100000
   No PATCHes used for this Group
 ************************************************************
  Group 24. Dumps For Restarts
 ************************************************************
  Echo save-block settings for Group 24
  save24begin
DISTIL=T
EX(P1  )=6.287E+02;EX(V1  )=8.242E-03
EX(W1  )=1.132E+00;EX(WDIS)=1.726E-02
EX(LTLS)=3.746E-04;EX(STRS)=1.616E-02
EX(SRM1)=3.982E+01;EX(BTAU)=1.259E+02
EX(PA  )=6.238E+02;EX(WA  )=1.145E+00
EX(GEN1)=2.614E+03;EX(VISL)=8.749E-03
  save24end
 
 GVIEW(P,-1.,0.,0.)
 GVIEW(UP,0.,1.,0.)
 GVIEW(VDIS,0.482418)
 GVIEW(CENTRE,4.991671E-03,0.05,0.5)
 
> DOM,    SIZE,        1.000000E-01, 5.000000E-02, 1.000000E+00
> DOM,    MONIT,       5.000000E-02, 3.780036E-02, 9.163948E-01
> DOM,    SCALE,       1.000000E+00, 1.000000E+00, 1.000000E+00
> DOM,    INCREMENT,   1.000000E-02, 1.000000E-02, 1.000000E-02
  > GRID,   RSET_X_1,      1, 1.000000E+00
> GRID,   RSET_Y_1,     20,-1.040000E+00,G
> GRID,   RSET_Z_1,    -40, 1.200000E+00
> DOM,    T_AMBIENT,   0.000000E+00
 
> OBJ,    NAME,        INLET
> OBJ,    POSITION,    0.000000E+00, 0.000000E+00, 0.000000E+00
> OBJ,    SIZE,        1.000000E-01, TO_END,       0.000000E+00
> OBJ,    DOMCLIP,     NO
> OBJ,    GEOMETRY,    poldef
> OBJ,    TYPE,        INLET
> OBJ,    PRESSURE,     P_AMBIENT
> OBJ,    VELOCITY,    0. ,0. ,1.
 
> OBJ,    NAME,        OUTL
> OBJ,    POSITION,    0.000000E+00, 0.000000E+00, AT_END
> OBJ,    SIZE,        1.000000E-01, TO_END,       0.000000E+00
> OBJ,    DOMCLIP,     NO
> OBJ,    GEOMETRY,    poldef
> OBJ,    TYPE,        OUTLET
> OBJ,    PRESSURE,    0.
> OBJ,    COEFFICIENT, 100.
> OBJ,    VELOCITY,    0. ,0. , SAME
 
> OBJ,    NAME,        WALL
> OBJ,    POSITION,    0.000000E+00, AT_END,       0.000000E+00
> OBJ,    SIZE,        1.000000E-01, 0.000000E+00, TO_END
> OBJ,    DOMCLIP,     NO
> OBJ,    GEOMETRY,    poldef
> OBJ,    TYPE,        PLATE
STOP