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-PAPANASTASIOU PIPE FLOW msg Reynolds number = 10 msg Yield number = 20 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 pressto 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 = 108 ************************************************************ Group 1. Run Title TEXT(108 2d pipe flow-Bingham-Pl-Papan. fluid) ************************************************************ Echo save-block settings for Group 1 save1begin The case concerns the steady laminar flow of a Bingham-Plastic- Panastasiou 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 + Fp*tauy/G where K is the rigidity coefficient, tauy is the yield stress, Fp is the Pananastasiou regularisation function, 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 = 20 N/m^2. These conditions correspond to a generalised Reynolds number of 10 and a Yield number of 20 (Hedstrom number of 200). 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(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,SRM1) STORE(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 =20. ;ENULC =0. CP1 =1. ENUT =1.0E-04 ************************************************************ 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,ENULD) 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=20.0 ! Yield number REY;YIELDN ENUL=GRND4;IENULA=-6 ! Bingham-Papanastasiou model LSG31=F ! =T for implicit linearisation of ENUL ENULD=1000.0 ! Papanastasiou time constant 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(GEN1)=1.001E-10 FIINIT(BTAU)=1.001E-10 ;FIINIT(VISL)=1.0E-02 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 = 6000 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. ) OVRRLX =1.2 ************************************************************ Group 18. Limits VARMAX(VISL)=1.0E+10 ;VARMIN(VISL)=1.0E-10 ************************************************************ Group 19. EARTH Calls To GROUND Station 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)=331.200012 ;EX(V1)=1.095E-03 EX(W1)=0.1078 ;EX(STRS)=1.693E-03 EX(SRM1)=1.766 ;EX(TAUR)=1.102 EX(PA)=327.5 ;EX(WA)=0.1077 EX(GR)=0.07068 ;EX(GEN1)=7.818 EX(BTAU)=22.030001 ;EX(VISL)=2.606 ************************************************************ 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 = 7 ;IZMON = 35 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 )=3.312E+02;EX(V1 )=1.095E-03 EX(W1 )=1.078E-01;EX(STRS)=1.693E-03 EX(SRM1)=1.766E+00;EX(TAUR)=1.102E+00 EX(PA )=3.275E+02;EX(WA )=1.077E-01 EX(GR )=7.068E-02;EX(GEN1)=7.818E+00 EX(BTAU)=2.203E+01;EX(VISL)=2.606E+00 save24end GVIEW(P,-0.999975,4.999804E-03,-4.99994E-03) GVIEW(UP,4.999991E-03,0.999988,-2.499858E-05) GVIEW(VDIS,0.416139) 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, 4.136835E-02, 8.936836E-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