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 POWELL-EYRING PIPE FLOW
   msg Reynolds number = 600
   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 = 105
 ************************************************************
  Group 1. Run Title
 TEXT(105 2d pipe flow - Powell-Eyring fluid  )
 ************************************************************
  Echo save-block settings for Group  1
  save1begin
  This case concerns the steady laminar flow of blood in a
  circular pipe. The blood is represented by a Powell-
  Eyring (PE) non-Newtonian fluid, whose apparent
  viscosity is given by:
 
    emua = emui + {(emu0-emui)/T}*[Sinh^-1(T*G)]/G
 
  where G is the mean strain rate, emui is the infinite-shear-rate
  viscosity and emu0 is the zero-shear-rate viscosity and T is the
  time constant.
 
  Analytical solutions for the laminar flow of PE fluids
  in circular pipes have been reported by: Govier, G.W. & Aziz,K.
  The flow of complex mixtures in pipes, Section 2.23, p193,
  R.E.Kreiger Pub. Co., Huntington, New York, (1977).
 
  The bulk inlet velocity is 0.1m/s and the fluid density is
  1065 kg/m^3. The pipe diameter and length are 0.02m and 0.5m,
  respectively; and the rheology parameters are set to:
  emu0=0.055 Pa.s; emui=0.0035 Pa.s and T = 5.5s. The flow
  Reynolds number based on emui is about 600.
 
  The task is to predict the pressure drop for a given Reynolds
  number, and then compare the results with the analytical solution.
  The inform facility is used to compute the analytical pressure
  profile, which is stored in PA.
  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(141)=ETA ;NAME(142)=STRS
 NAME(143)=WA ;NAME(144)=PA
 NAME(145)=SRM1 ;NAME(146)=WDIS
 NAME(148)=BTAU ;NAME(149)=GEN1
 NAME(150) =VISL
    * Solved variables list
 SOLVE(P1,V1,W1)
    * Stored variables list
 STORE(VISL,GEN1,BTAU,WDIS,SRM1,PA,WA,STRS)
 STORE(ETA)
    * 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
(stored of ETA is RHO1*VISL) ! dynamic viscosity
  save7end
 ************************************************************
  Group 8. Terms & Devices
 NEWENL = T
 ************************************************************
  Group 9. Properties
 RHO1 =1065.
 ENUL = GRND4
 ENULA =0.055 ;ENULB =3.5E-03 ;ENULC =0.
 CP1 =1.
 DISWAL
 ENUT =1.0E-10
 ************************************************************
  Echo save-block settings for Group  9
  save9begin
REAL(RIN,DIN,WIN,AIN,DPDZ,QIN,ETA0,ETAI,DPDZN)
REAL(REY,PI,TAUW,TAUWN,ENULD,TC,GDWN)
REAL(GW,GDW,ALFA,BETA,PHI1,ASINHX,XX,AAA,AP,BP,OMEG,OMEG2,OMEG3,PHI2)
REAL(FX,FXP,GDWP,TW,GW2,GW3,GW4,TW3,PHI3,BETA3,GMW1,GMW2)
REAL(T1,T2,T3,T4,T5,T6,T7,DIN3,QQ,BETA2,ZETA)
 
PI=3.14159265
WIN=0.1;RIN=0.01;DIN=2.*RIN;DIN3=DIN*DIN*DIN
AIN=PI*RIN*RIN;QIN=WIN*AIN
WIN;DIN;QIN
ETAI=0.0035   ! Pa.s
ENULB=ETAI
REY=RHO1*WIN*DIN/ETAI
REY
ENUL=GRND4;IENULA=5 ! Powell-Eyring model
ETA0=0.055    ! Pa.s
TC=5.5  ! s
      ETA0=ENULAG; ETAI=ENULBG; TC=ENULDG
      ETA=(ETAI+(ETA0-ETAI)*ASINH(TC*GMDT)/(TC*GMDT))
      TAU=ETA*GMDT
ENULB=ETAI    ! infinite-shear-rate viscosity
ENULA=ETA0    ! zero-shear-rate viscosity
ENULD=TC      ! time constant
ETA0;ETAI;TC
 
    ** Estimate DPDZ, TAUW & GDW from Newtonian pressure drop
DPDZN=8.*ETAI*WIN/(RIN*RIN) ; TAUWN=DPDZN*RIN/2. ; GDWN=TAUWN/ETAI
DPDZN; TAUWN; GDWN
    ** Govier-Aziz analytical parameters
AP=1./TC ; BP=TC/(ETA0-ETAI) ; OMEG=AP*BP*ETAI
AP;BP;OMEG
OMEG2=OMEG*OMEG;OMEG3=OMEG2*OMEG
QQ=8.*QIN/(PI*DIN3*AP)
QQ
  ** Iterate to find GDW for PE fluid using TAUWN
GDWP=100.
GDW =GDWN
DO II=1,10
IF(ABS(GDWP).GT.1.E-4) THEN
XX=GDW*TC;AAA=XX+SQRT(XX*XX+1.)
ASINHX=LOG(AAA)
FX=ETAI*GDW+ASINHX/BP-TAUWN
FXP=ETAI+(TC/BP)/(1.+XX*XX)**0.5
GDWP=-FX/FXP
GDW=GDW+GDWP
II;GDWP
ENDIF
ENDDO
GDW
 
    ** Two initial values are needed for the Secant Method
GMW1=GDW/AP
GMW2=1.25*GMW1
GMW1
GMW2
   ** Iterate to find GW from Secant Method
REAL(AA,DL,BB,DX,X0,X1,X2,FX1,FX0,DD,ABSDX)
DL = 1.0E-03
AA  = GMW1; BB  = GMW2
DX = (BB-AA)/10.
X0 = (AA+BB)/2.0
X1 = X0 + DX
 
DO II=1,10
IF(ABS(DX).GT.DL) THEN
GW=X0
GW2=GW*GW;GW3=GW2*GW;GW4=GW2*GW2
ALFA=(1.+GW2);BETA=SQRT(ALFA)
alfa
beta
BETA2=BETA*BETA;BETA3=BETA2*BETA
ZETA=(2.*GW2+1.)
 
XX=GW;AAA=XX+SQRT(XX*XX+1.)
ASINHX=LOG(AAA)
PHI1=ASINHX;PHI2=PHI1*PHI1;PHI3=PHI2*PHI1
TW=OMEG*GW+ASINHX;TW3=TW**3;T7=1./(3.*TW3)
 
T1=GW/3.
T2=0.25*OMEG3*GW4
T3=GW3*PHI1-GW2*BETA+2.*BETA3/3.-2./3.
T4=ZETA*PHI2 - 2.*GW*BETA*PHI1+GW2
T5=GW*PHI3-3.*BETA*PHI2
T6=GW*PHI1-BETA+1.
 
FX0=T1-T7*(T2+OMEG2*T3+0.75*OMEG*T4+T5+6.*T6)-QQ
FX0
GW=X1
GW2=GW*GW;GW3=GW2*GW;GW4=GW2*GW2
ALFA=(1.+GW2);BETA=SQRT(ALFA)
BETA2=BETA*BETA;BETA3=BETA2*BETA
ZETA=(2.*GW2+1.)
 
XX=GW;AAA=XX+SQRT(XX*XX+1.)
ASINHX=LOG(AAA)
PHI1=ASINHX;PHI2=PHI1*PHI1;PHI3=PHI2*PHI1
TW=OMEG*GW+ASINHX;TW3=TW**3;T7=1./(3.*TW3)
 
T1=GW/3.
T2=0.25*OMEG3*GW4
T3=GW3*PHI1-GW2*BETA+2.*BETA3/3.-2./3.
T4=ZETA*PHI2 - 2.*GW*BETA*PHI1+GW2
T5=GW*PHI3-3.*BETA*PHI2
T6=GW*PHI1-BETA+1.
FX1=T1-T7*(T2+OMEG2*T3+0.75*OMEG*T4+T5+6.*T6)-QQ
FX1
DD  = FX1 - FX0
X2 = X1 - FX1*(X1-X0)/DD
X0 = X1
X1 = X2
DX = X1 - X0
II;X0;DX
ENDIF
ENDDO
GW=X0
 
  ** Now compute shear stress & DPDZ for PE fluid
XX=GW;AAA=XX+SQRT(XX*XX+1.)
ASINHX=LOG(AAA)
TW=OMEG*GW+ASINHX
TAUW=TW/BP
TAUW
DPDZ=2.*TAUW/RIN
DPDZ
 
      ** Newtonian 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))
      ** Newtonian analytical axial-velocity profile
(stored of WA is 2.*WIN*(1.-(YG/RIN)^2))
  save9end
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
 FIINIT(W1)=0.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 =5.0E-06
 ************************************************************
  Group 16. Terminate Iterations
 ************************************************************
  Group 17. Relaxation
 RELAX(P1 ,LINRLX,1. )
 RELAX(V1 ,FALSDT,0.05 )
 RELAX(W1 ,FALSDT,0.05 )
 ************************************************************
  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)=11.16 ;EX(V1)=3.542E-04
 EX(W1)=0.1137 ;EX(ETA)=0.01119
 EX(STRS)=1.165E-05 ;EX(WA)=0.1198
 EX(PA)=10.41 ;EX(SRM1)=19.66
 EX(WDIS)=3.453E-03 ;EX(LTLS)=1.499E-05
 EX(BTAU)=0.1129 ;EX(GEN1)=665.099976
 EX(VISL)=1.051E-05
 ************************************************************
  Group 21. Print-out of Variables
 OUTPUT(ETA ,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 = 7 ;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  )=1.116E+01;EX(V1  )=3.542E-04
EX(W1  )=1.137E-01;EX(ETA )=1.119E-02
EX(STRS)=1.165E-05;EX(WA  )=1.198E-01
EX(PA  )=1.041E+01;EX(LTLS)=1.499E-05
EX(WDIS)=3.453E-03;EX(SRM1)=1.966E+01
EX(BTAU)=1.129E-01;EX(VISL)=1.051E-05
EX(GEN1)=6.651E+02
  save24end
 
 GVIEW(P,-1.,0.,0.)
 GVIEW(UP,0.,1.,0.)
 GVIEW(VDIS,0.247398)
 GVIEW(CENTRE,4.991671E-04,5.0E-03,0.25)
 
> DOM,    SIZE,        1.000000E-01, 1.000000E-02, 5.000000E-01
> DOM,    MONIT,       5.000000E-02, 4.136836E-03, 4.581974E-01
> DOM,    SCALE,       1.000000E+00, 2.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. ,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, 1000.
 
> 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