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 BINGHAM-PLASTIC PIPE FLOW
   msg Reynolds number = 10
   msg Yield number    = 2.5
   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 = 107
 ************************************************************
  Group 1. Run Title
 TEXT(107 2d pipe flow - Bingham-plastic fluid)
 ************************************************************
  Echo save-block settings for Group  1
  save1begin
  The case concerns the steady laminar flow of a Bingham-plastic
  non-Newtonian fluid in a circular pipe. This type of fluid remains
  rigid when the shearing stress is less than the yield stress tauy
  and flows somewhat like a Newtonian fluid when the shearing stress
  exceeds tauy.
 
  The apparent dynamic viscosity of such a fluid is given by the
  following two-parameter formula:
 
    emua = K + tauy/G
 
  where K is the rigidity coefficient, tauy is the yield stress,
  and G=E^0.5 with E is the mean shear rate. Examples of fluids which
  behave as, or nearly as, Bingham plastics include water suspensions
  of clay, sewage sludge, some emulsions and thickened hydrocarbon
  greases.
 
  Analytical solutions for the laminar flow of Bingham-Plastic
  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).
 
  The bulk inlet velocity is 0.1m/s and the fluid density is
  1000kg/m^3. The pipe diameter and length are 0.2m and 1.0m,
  respectively; and the rheology parameters are set to:
  K=2.0 Pa.s and tauy = 2.5 N/m^2. These conditions correspond
  to a generalised Reynolds number of 10 and a Yield number of
  2.5 (Hedstrom number of 25).
 
  The task is to predict the pressure drop and fully-developed
  axial-velocity profile for a given Reynolds and Yield 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)=STRS ;NAME(142)=SRM1
 NAME(143)=WDIS ;NAME(144)=TAUR
 NAME(145)=PA ;NAME(146)=WA
 NAME(147)=GR ;NAME(148)=GEN1
 NAME(149)=BTAU ;NAME(150)=VISL
    * Solved variables list
 SOLVE(P1,V1,W1)
    * Stored variables list
 STORE(VISL,BTAU,GEN1,GR,WA,PA,TAUR,WDIS)
 STORE(SRM1,STRS)
    * 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
 NEWENL = T
 ************************************************************
  Group 9. Properties
 RHO1 =1000.
 ENUL = GRND4
 ENULA =2. ;ENULB =2.5 ;ENULC =0.
 CP1 =1.
 DISWAL
 ENUT =1.0E-10
 ************************************************************
  Echo save-block settings for Group  9
  save9begin
REAL(RIN,DIN,WIN,AIN,DPDZ,TAUY,REY,YIELDN,HEDNO,GRP,TAUW,FRIC)
REAL(ACON,BCON,CCON,PLEN,CONSI,FRICO,DELF,TAURAT,BINGNO)
 
RIN=0.1;DIN=2.*RIN          ! pipe diameter
PLEN=1.0                    ! pipe length
WIN=0.1                     ! inlet velocity
WIN
AIN=RIN*RIN/2.
 
 
REY=10.                     ! Reynolds number
YIELDN=2.5                  ! Yield number
REY;YIELDN
ENUL=GRND4;IENULA=6          ! Bingham-plastic model
LSG31=F                     ! =T for implicit linearisation of ENUL
 
ENULA=RHO1*DIN*WIN/REY      ! Plastic viscosity
TAUY=YIELDN*WIN*ENULA/DIN   ! yield stress
ENULB=TAUY
TAUY
 
ENULA;ENULB
CONSI=ENULA
 
HEDNO=YIELDN*REY            ! Hedstrom number
HEDNO
 
BINGNO=TAUY*DIN/(CONSI*WIN) ! Bingham number
BINGNO
 
    ** Analytical solution for wall shear stress by
       iteration for the friction factor f
ACON=16./REY; BCON=(16./6.)*HEDNO/REY**2 ; CCON=-(16./3.)*HEDNO**4/REY**8
DELF=10.
FRICO=ACON+BCON  ! Initial estimate
FRICO
DO II=1,10
IF(ABS(DELF).GT.1.E-3) THEN
FRIC=ACON+BCON+CCON/FRICO**3
DELF=FRIC-FRICO
FRICO=FRIC
ENDIF
ENDDO
FRIC;DELF
TAUW=0.5*FRIC*RHO1*WIN*WIN             ! wall shear stress
TAUW
DPDZ=2.*TAUW/RIN                       ! pressure drop
DPDZ
GRP =2.*TAUY/DPDZ                      ! Yield radius
GRP
TAURAT=TAUY/TAUW
TAURAT
   ** Consistency check between wall shear stress & bulk velocity
WIN=0.25*RIN*(TAUW/CONSI)*(1.-4.*(TAURAT)/3.+(1./3)*(TAURAT)**4)
WIN
(stored of GR is YG)
(stored of GR is :GRP: with IF(YG.LE.:GRP:))
(stored of WA is (0.25*DPDZ*(RIN*RIN-GR^2)-TAUY*(RIN-GR))/CONSI)
      ** 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))
       ** Stress-ratio output
(stored of TAUR is BTAU/:TAUY:)
  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(BTAU)=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
 LITER(P1)=50
 ************************************************************
  Group 17. Relaxation
 RELAX(P1 ,LINRLX,1. )
 RELAX(V1 ,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)=113.400002 ;EX(V1)=1.762E-03
 EX(W1)=0.115 ;EX(STRS)=5.96E-04
 EX(LTLS)=1.498E-03 ;EX(SRM1)=2.064
 EX(WDIS)=0.03452 ;EX(TAUR)=2.611
 EX(PA)=111.800003 ;EX(WA)=0.116
 EX(GR)=0.05749 ;EX(GEN1)=6.532
 EX(BTAU)=6.628 ;EX(VISL)=0.03591
 ************************************************************
  Group 21. Print-out of Variables
 OUTPUT(PA  ,Y,N,Y,N,Y,Y)
 OUTPUT(WA  ,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 = 13 ;IZMON = 31
 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(TAUR)=2.611E+00;EX(GR  )=5.749E-02
EX(LTLS)=1.498E-03;EX(WDIS)=3.452E-02
EX(P1  )=1.134E+02;EX(V1  )=1.762E-03
EX(W1  )=1.150E-01;EX(SRM1)=2.064E+00
EX(WA  )=1.160E-01;EX(GEN1)=6.532E+00
EX(BTAU)=6.628E+00;EX(VISL)=3.591E-02
EX(PA  )=1.118E+02;EX(STRS)=5.960E-04
  save24end
 
 GVIEW(P,-0.99995,9.999833E-03,0.)
 GVIEW(UP,9.999833E-03,0.99995,0.)
 GVIEW(VDIS,0.39225)
 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, 7.126657E-02, 7.952863E-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. ,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