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 CASSON-PAPANASTASIOU PIPE 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 = 112
 ************************************************************
  Group 1. Run Title
 TEXT(112 2d pipe flow - Casson-Papan. fluid  )
 ************************************************************
  Echo save-block settings for Group  1
  save1begin
  This case concerns the steady laminar flow of a Casson-
  Panastasiou non-Newtonian fluid in a circular pipe. The
  apparent viscosity of such a fluid is given by:
 
    emu = ( Fp*(tau0/G)^1/2 + K^1/2)^2
 
  where K is the consistency index, tau0 is the yield stress,
  and Fp is the Pananastasiou regularisation function. The
  shear stress tau = emu*G where G is the mean strain rate.
  Casson fluids have applications in metallurgy, drilling
  operations,food processing and most importantly in describing
  the flow of blood.
 
  The bulk inlet velocity is 1m/s and the fluid density is
  1 kg/m^3. The pipe diameter and length are 0.1m and 1m,
  respectively; and the rheology parameters are set to:
  K=0.02 Pa.s and tau0 = 2 N/m^2. These conditions correspond
  to a flow Reynolds number based on K of 10, and a Casson
  number C of 20, where C=tau0*D/(U*K).
 
  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(140)=GR ;NAME(141)=TAUR
 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)
 STORE(TAUR,GR)
    * 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.02 ;ENULB =2. ;ENULC =0.
 CP1 =1.
 DISWAL
 ENUT =1.0E-10
 ************************************************************
  Echo save-block settings for Group  9
  save9begin
REAL(RIN,DIN,WIN,DPDZ,CASSN,FX,FXP,GRP)
REAL(TAUW,TAUP,CONSI,TAU0,ACON,BCON,CCON,REY)
 
RIN=0.1    ! Pipe radius
DIN=2.*RIN ! Pipe diameter
WIN=1.0    ! Bulk inlet velocity
REY=10.    ! Reynolds number
 
ENUL=GRND4;IENULA=-8  ! Casson--Papanastasiou fluid
   CONSI=ENULA ; TAU0=ENULB
   ETA = (SQRT(TAU0/GAMDOT)+SQRT(CONSI) )**2
ENULA=DIN*WIN/REY          ! Consistency index
CASSN=20.0                 ! Casson number
ENULB=CASSN*WIN*ENULA/DIN  ! Yield stress
ENULD=1.E10                ! Papanastasiou time constant
REY;CASSN;ENULA;ENULB
    ** Analytical pressure solution
CONSI= ENULA ; TAU0 = ENULB
 
DPDZ=8.*CONSI*WIN/RIN**2 ! "Newtonian" pressure drop
TAUW=0.5*RIN*DPDZ        ! Force balance in flow direction
DPDZ;TAUW
       iterate for analytical wall shear stress
ACON=4.*SQRT(TAU0)/7. ;BCON=(TAU0**4)/84. ; CCON=TAU0/3.
TAUP=100.0
 
DO II=1,20
IF(ABS(TAUP).GT.1.E-3) THEN
FX =0.25*TAUW-ACON*(TAUW)**0.5-BCON/TAUW**3+CCON-WIN*CONSI/RIN
FXP=0.25-(0.5*ACON/TAUW**0.5)+3.*BCON/TAUW**4
TAUP=-FX/FXP
TAUW=TAUW+TAUP
II;TAUW;TAUP
ENDIF
ENDDO
 
DPDZ=2.*TAUW/RIN
TAUW;DPDZ
(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=0.5*RIN*TAUW/CONSI
BCON=(8./3.)*(TAU0/TAUW)**0.5
CCON=2.*TAU0/TAUW
(stored of GR is YG/:RIN:)
GRP=TAU0/TAUW                  ! Critical radius
GRP
(stored of GR is :GRP: with IF(GR.LE.:GRP:))
(stored of WA is ACON*((1.-GR^2)-BCON*(1.-GR^1.5)+CCON*(1.-GR)))
       ** Stress-ratio output
(stored of TAUR is BTAU/:TAU0:)
  save9end
 ************************************************************
  Group 10.Inter-Phase Transfer Processes
 ************************************************************
  Group 11.Initialise Var/Porosity Fields
 FIINIT(W1)=1. ;FIINIT(WDIS)=0.1
 FIINIT(BTAU)=1.001E-10 ;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-04
 ************************************************************
  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)=62.150002 ;EX(V1)=0.01334
 EX(W1)=1.101 ;EX(GR)=0.6035
 EX(TAUR)=1.888 ;EX(WDIS)=0.03452
 EX(SRM1)=18.540001 ;EX(STRS)=0.3166
 EX(LTLS)=1.498E-03 ;EX(BTAU)=3.776
 EX(PA)=60.740002 ;EX(WA)=1.108
 EX(GEN1)=712.400024 ;EX(VISL)=458.899994
 ************************************************************
  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 = 12 ;IZMON = 32
 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.215E+01;EX(V1  )=1.334E-02
EX(W1  )=1.101E+00;EX(GR  )=6.035E-01
EX(TAUR)=1.888E+00;EX(STRS)=3.166E-01
EX(SRM1)=1.854E+01;EX(BTAU)=3.776E+00
EX(PA  )=6.074E+01;EX(WA  )=1.108E+00
EX(GEN1)=7.124E+02;EX(VISL)=4.589E+02
  save24end
 
 GVIEW(P,-1.,0.,0.)
 GVIEW(UP,0.,1.,0.)
 GVIEW(VDIS,0.527151)
 GVIEW(CENTRE,4.991671E-03,0.05,0.5)
 
> DOM,    SIZE,        1.000000E-01, 1.000000E-01, 1.000000E+00
> DOM,    MONIT,       5.000000E-02, 6.675906E-02, 8.208498E-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