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 CARREAU-YUSADA FLUID PIPE FLOW
   msg Reynolds number =160.0
   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 end
   pause
   end
   END_USE
 ************************************************************
 IRUNN = 1 ;LIBREF = 104
 ************************************************************
  Group 1. Run Title
 TEXT(104 2d Pipe flow - Carreau-Yusada fluid )
 ************************************************************
  Echo save-block settings for Group  1
  save1begin
  This case concerns the steady laminar flow of a Carreau-
  Yusada non-Newtonian fluid in a circular pipe. The apparent
  dynamic viscosity of such a fluid is given by:
 
    emu = emui+(emu0-emui)/[1+(T.G)^m]^{(n-1)/2}
 
  where emu0 is the low-shear viscosity, emui is the high-
  shear viscosity, T is the characteristic time scale, G is
  the mean strain rate, m is an index, and n is the flow-
  behaviour index. The Carreau-Yusada rheology model finds
  applications in biological fluids and polymeric liquids.
 
  Semi-analytical solutions for the laminar flow of Carreau
  fluids in circular pipes have been reported by: T.Sochi,
  Analytical solutions for the flow of Carreau and Cross
  fluids in circular pipes and thin slits. Rheologica Acta,
  Vol.54, Issue 8, 745-756, (2015).
 
  The pipe diameter and length are 0.04m and 0.5m,
  respectively; and the rheology parameters are set to:
  emui=0.001 Pa.s; emu0=0.08 Pa.s; T=2s; m=2.0 and n=0.9.
  The choice of m=2 reduces the Carreau-Yusada model to the
  Carreau model.
 
  The task is to predict the pressure drop dP for fully-
  developed flow at the a given Reynolds number; and then
  compare the result with the semi-analytical solution of
  Sochi (2015). For this purpose, the inlet velocity is
  set to 4m/s, which equates to a volumetric flow rate Q
  of 5E-3 m^3/s,and a Reynolds number of 160. For these
  conditions, PHOENICS predicts a pressure drop of 1614Pa,
  which is in good agreement with the semi-analytical
  solution of 1620 Pa reported by Sochi (2015).
  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(143)=PA ;NAME(145)=WDIS
 NAME(146)=SRM1 ;NAME(147)=STRS
 NAME(148)=BTAU ;NAME(149)=GEN1
 NAME(150) =VISL
    * Solved variables list
 SOLVE(P1,V1,W1)
    * Stored variables list
 STORE(VISL,GEN1,BTAU,STRS,SRM1,WDIS,PA)
    * 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 =1.
 ENUL = GRND4
 ENULA =0.08 ;ENULB =1.0E-03 ;ENULC =0.9
 CP1 =1.
 DISWAL
 ENUT =1.0E-10
 ************************************************************
  Echo save-block settings for Group  9
  save9begin
REAL(REY,RIN,DIN,WIN,AIN,DPDZ,QIN,PI)
REAL(DPDZ,PLEN,TAUW)
PI=3.14159265
 
RIN=0.02;DIN=2.*RIN  ! Pipe diameter
WIN=4.               ! Bulk velocity
PLEN=0.5             ! Pipe length
 
AIN=PI*RIN*RIN;QIN=WIN*AIN
QIN                 !  Volumetric flow rate
 
ENUL=GRND4;IENULA=4 ! Carreau-Yusada fluid model
            ETA0=ENULAG; ETAI=ENULBG; AN=ENULCG; TCONS=ENULDG
            ETA=(ETAI+(ETA0-ETAI)*(1.+(TCONS*GAMDOT)**M)**(AN-1.)/2.)
ENULA=0.08;ENULB=0.001
ENULC=0.9   ! ENULE  1sT Exponent
ENULD=2.0  ! ENULD  Time constant
ENULE=2.0  ! ENULE 2nd Exponent
ENULC;ENULD;ENULE
REY=RHO1*WIN*DIN/ENULB     !  Reynolds number
REY
DPDZ=1620./PLEN        !  From Sochi (2015)
TAUW=0.5*RIN*DPDZ
TAUW
      ** Semi-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)=4. ;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 = 2000
 RESREF(P1)=5.0E-16 ;RESREF(V1)=5.0E-16
 RESREF(W1)=5.0E-16
 RESFAC =1.0E-05
 ************************************************************
  Group 16. Terminate Iterations
 LITER(P1)=200
 ************************************************************
  Group 17. Relaxation
 RELAX(P1 ,LINRLX,1. )
 RELAX(V1 ,FALSDT,1.0E-02 )
 RELAX(W1 ,FALSDT,5. )
 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)=784.299988 ;EX(V1)=0.03932
 EX(W1)=4.746 ;EX(PA)=800.400024
 EX(LTLS)=5.994E-05 ;EX(WDIS)=6.904E-03
 EX(SRM1)=443.299988 ;EX(STRS)=1.582
 EX(BTAU)=17.93 ;EX(GEN1)=2.546E+05
 EX(VISL)=0.04216
 ************************************************************
  Group 21. Print-out of Variables
 OUTPUT(PA  ,Y,N,Y,N,Y,Y)
 OUTPUT(BTAU,Y,N,Y,N,Y,Y)
 OUTPUT(VISL,Y,N,Y,N,Y,Y)
 ************************************************************
  Group 22. Monitor Print-Out
 IXMON = 1 ;IYMON = 11 ;IZMON = 38
 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  )=7.843E+02;EX(V1  )=3.932E-02
EX(W1  )=4.746E+00;EX(SRM1)=4.433E+02
EX(WDIS)=6.904E-03;EX(LTLS)=5.994E-05
EX(STRS)=1.582E+00;EX(BTAU)=1.793E+01
EX(GEN1)=2.546E+05;EX(VISL)=4.216E-02
  save24end
 
 GVIEW(P,-1.,0.,0.)
 GVIEW(UP,0.,1.,0.)
 GVIEW(VDIS,0.240194)
 GVIEW(CENTRE,9.983341E-04,9.999999E-03,0.25)
 
> DOM,    SIZE,        1.000000E-01, 2.000000E-02, 5.000000E-01
> DOM,    MONIT,       5.000000E-02, 1.241425E-02, 4.813020E-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.250000E+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. ,4.
 
> 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