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
 
   cl
   msg ELLIS-FLUID LAMINAR FLOW
   msg Reynolds number = 2.75
   msg Pressure(P1) profile
   msg Blue line --- PHOENICS solution
   msg crosses ---   analytical solution
   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 = 106
 ************************************************************
  Group 1. Run Title
 TEXT(106 2d pipe flow - Ellis fluid          )
 ************************************************************
  Echo save-block settings for Group  1
  save1begin
  This case concerns the steady laminar flow of an Ellis
  non-Newtonian fluid in a circular pipe. The apparent viscosity
  of such a fluid is given by:
 
    emu = emu0/(1.+(tau/tau,1/2)^(alfa-1))
 
  where emu0 is the low-shear viscosity, tau is the shear
  stress, tsu,1/2 is the shear stress at which emu=0.5*emu0
  and alfa is a flow-behaviour index related to the power-law-
  fluid index n by alfa= 1=n. The shear stress tau = emu*G
  where G is the mean strain rate. The Ellis model has seen
  extensive use for modelling the rheologies of a wide range
  of polymeric solutions.
 
  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, p110, p119-121,
  John Wiley, (1967).
 
  The bulk inlet velocity is 10m/s and the fluid density is
  1 kg/m^3. The pipe diameter and length are 0.05m and 0.5m,
  respectively; and the rheology parameters are set to:
  emu0=0.185 Pa.s; alfa=2.42 and tau,1/2 = 1025 N/m^2. The
  flow Reynolds number based on emu0 is 2.75.
 
  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.
  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 =1.
 ENUL = GRND4
 ENULA =0.185 ;ENULB =2.4 ;ENULC =1025.
 CP1 =1.
 DISWAL
 ENUT =1.0E-10
 ************************************************************
  Echo save-block settings for Group  9
  save9begin
REAL(RIN,DIN,WIN,DPDZ,REY,BETA,ZETA,QIN,AIN,PI)
REAL(PHI0,PHI1,ALFA,TAUHLF,EMUA,EMU0)
REAL(ACON,BCON,CCON,FX,FXP,TAUW,TAUP)
 
ENUL=GRND4; IENULA=9    ! Ellis rheology model
ENULB=2.4 ; ENULA=0.185 ; ENULC=1025.
 
 
EMU0=ENULA   ! low-shear viscosity
TAUHLF=ENULC ! shear stress at which emua=0.5*emu0
ALFA=ENULB   ! indical parameter
 
BETA=ALFA-1.0;ZETA=ALFA+1.0
 
PI=3.1415926
 
QIN=0.02     ! volumetric inflow rate
RIN=0.025    ! pipe diameter
 
DIN=2.*RIN; AIN=PI*RIN*RIN
WIN=QIN/AIN
WIN
 
 
ENULA;ENULB
    * Reynolds number
REY=WIN*DIN/ENULA
REY
      ** Analytical pressure solution
 
PHI0=1./ENULA
PHI1=1./(ENULA*TAUHLF**BETA)
PHI0;PHI1;BETA;ZETA
 
 
DPDZ=8.*EMU0*WIN/RIN**2 ! Newtonian pressure drop
DPDZ
TAUW=0.5*RIN*DPDz   ! Force balance in flow direction
TAUW
       iterate for analytical wall shear stress
ACON=0.25*PHI0*RIN
BCON=(PHI1*RIN)/(3.+ALFA)
TAUP=100.0
DO II=1,20
IF(ABS(TAUP).GT.1.E-3) THEN
FX =ACON*TAUW+BCON*(TAUW)**ALFA-WIN
FXP=ACON+BCON*ALFA*(TAUW)**BETA
TAUP=-FX/FXP
TAUW=TAUW+TAUP
II
TAUW
TAUP
ENDIF
ENDDO
 
TAUW
DPDZ=2.*TAUW/RIN
DPDZ
      ** 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))
      ** Analytical axial velocity profile
ACON=(PHI0*RIN*RIN/4.)*DPDZ
CCON=(1.+ALFA)*(2.**ALFA)
BCON=PHI1*(RIN**ZETA)*(DPDZ**ALFA)/CCON
(stored of WA is ACON*(1.-(YG/RIN)^2)+BCON*(1.-(YG/RIN)^ZETA))
  save9end
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
 FIINIT(W1)=10. ;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.05 )
 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)=5246. ;EX(V1)=0.1069
 EX(W1)=11.87 ;EX(WDIS)=8.63E-03
 EX(SRM1)=890.200012 ;EX(STRS)=13.25
 EX(LTLS)=9.365E-05 ;EX(BTAU)=150.199997
 EX(PA)=5335. ;EX(WA)=12.13
 EX(GEN1)=1.021E+06 ;EX(VISL)=0.1727
 ************************************************************
  Group 21. Print-out of Variables
 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 = 15 ;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  )=5.246E+03;EX(V1  )=1.069E-01
EX(W1  )=1.187E+01;EX(WDIS)=8.630E-03
EX(LTLS)=9.365E-05;EX(STRS)=1.325E+01
EX(SRM1)=8.902E+02;EX(BTAU)=1.502E+02
EX(GEN1)=1.021E+06;EX(VISL)=1.727E-01
EX(PA  )=5.335E+03;EX(WA  )=1.213E+01
  save24end
 
 GVIEW(P,-0.999988,0.,-4.999915E-03)
 GVIEW(UP,0.,1.,0.)
 GVIEW(VDIS,0.247512)
 GVIEW(CENTRE,1.247918E-03,0.0125,0.25)
 
> DOM,    SIZE,        1.000000E-01, 2.500000E-02, 5.000000E-01
> DOM,    MONIT,       5.000000E-02, 1.994204E-02, 4.581974E-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,    WIREFRAME,   YES
> OBJ,    PRESSURE,     P_AMBIENT
> OBJ,    VELOCITY,    0. ,0. ,10.
 
> 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,    WIREFRAME,   YES
> 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
> OBJ,    WIREFRAME,   YES
STOP