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 HERSCHEL-BULKLEY-PAPANASTASIOU PIPE FLOW msg Reynolds number = 10 msg Yield number = 2 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 = 110 ************************************************************ Group 1. Run Title TEXT(110 2d pipe flow-Hers.-Bulk.-Papan.fluid) ************************************************************ Echo save-block settings for Group 1 save1begin The case concerns the steady laminar flow of a Herschel-Bulkley- 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 three-parameter formula: emua = K*G^(n-1) + Fp*tauy/G where G=E^0.5,E is the mean shear rate, K is the consistency index, n is the power-law index, tauy is the yield stress, and Fp is the Pananastasiou regularisation function,. Analytical solutions for the laminar flow of HB fluids in circular pipes have been reported by: A.H.P.Skelland, Non-Newtonian Flow and Heat Transfer, p110, p119-121, 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 1000 kg/m^3. The pipe diameter and length are 0.1m and 1.0m, respectively; and the rheology parameters are set to: K=1.0 Pa.s; n=0.8 and tauy = 2 N/m^2. These conditions correspond to a Reynolds and Yield number of 10 and 2, respectively. The generalised Reynolds number is 14 and the generalised Hedstrom number is 28. The task is to predict the pressure drop and fully-developed axial-velocity profile for a given Reynolds and Hedstrom 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(141)=WDIS 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(WDIS,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 ADDDIF = T NEWENL = T ************************************************************ Group 9. Properties RHO1 =1000. ENUL = GRND4 ENULA =1. ;ENULB =2. ;ENULC =0.8 CP1 =1. DISWAL ENUT =1.0E-10 ************************************************************ Echo save-block settings for Group 9 save9begin REAL(RIN,DIN,WIN,DPDZ,TAUY,REY,YIELDN,HEDNO,GRP,TAUW,FRIC,AN) REAL(ACON,BCON,CCON,DCON,FX,FXP,TAUP,ALFA,PHIA,ZETA,CONSI,TAUD) REAL(AP1,AP2,AM1,RADY,WCL,HEDSNO,REYMR,BETA,PLEN,ENULD) RIN=0.05;DIN=2.*RIN ! Pipe diameter PLEN=1.0 ! Pipe length WIN=0.1 ! Inlet velocity WIN;DIN REY=10.;YIELDN=2.0 ! Reynolds number and Yield numbers REY;YIELDN ENUL=GRND4;IENULA=-7 ! Herschel-Bulkley-Papanastasiou model CONSI=ENULA ; TAU0=ENULB; AN=ENULC EMUA = (CONSI*GAMDOT**(AN-1)+TAU0/GAMDOT) AN=0.8 ! Flow Index CONSI=RHO1*DIN*WIN/REY ! Consistency index TAUY=YIELDN*WIN*ENULA/DIN ! yield stress ENULD=1.E10 ! Papanastasiou time constant ENULA=CONSI;ENULB=TAUY;ENULC=AN CONSI;AN;TAUY * Generalised Metzner-Reed Reynolds number ZETA=8**(AN-1); PHIA=((3.*AN+1.)/(4.*AN))**AN REYMR=RHO1*WIN**(2.-AN)*(DIN**AN)/(CONSI*ZETA*PHIA) REYMR * Generalised Hedstrom number HEDSNO=RHO1*DIN*DIN*(TAUY)**(2./AN-1.)/CONSI**(2./AN) HEDSNO ** Initial guess - Bingham analytical pressure drop (approximate formula) DPDZ=32.*RHO1*WIN**2*(1.+YIELDN/6.)/REY/DIN DPDZ TAUW=DIN*DPDZ/4. TAUW TAUD=TAUW-TAUY TAUD iterate for HB analytical wall shear stress, solve for TAUD=(TAUW-TAU0) ALFA=1.+1./AN BCON=2.*AN*TAUY/(2.*AN+1.) ACON=AN/(3.*AN+1.) CCON=AN*TAUY*TAUY/(AN+1.) DCON=(WIN/RIN)*CONSI**(1./AN) AP1=ALFA+1.;AP2=ALFA+2.;AM1=ALFA-1. TAUP=100.0 DO II=1,20 IF(ABS(TAUP).GT.1.E-3) THEN FX = ACON*(TAUD**AP2)+BCON*(TAUD**AP1)+CCON*(TAUD**ALFA)-DCON*(TAUD+TAUY)**3 FXP= AP2*ACON*(TAUD**AP1)+BCON*AP1*(TAUD**ALFA)+CCON*ALFA*(TAUD**AM1)-3.*DCON*(TAUD+TAUY)**2 TAUP=-FX/FXP TAUD=TAUD+TAUP ENDIF ENDDO TAUW=TAUY+TAUD TAUW DPDZ=2.*TAUW/RIN DPDZ FRIC=2.*TAUW/(RHO1*WIN*WIN) ! Friction factor FRIC ** 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 GRP=2.*TAUY/DPDZ ! radius of the central plug - Yielld point GRP REAL(RADRAT) RADRAT=GRP/RIN ACON=2.*AN/(1.+3.*AN) BCON=2.*AN/(1.+2.*AN) CCON=1.-ACON*(1.-RADRAT)**2-BCON*RADRAT*(1.-RADRAT) WCL= WIN/CCON WCL ALFA=(1.+AN)/AN BETA=AN/(1.+AN) ZETA=1./AN DCON=TAUW/(CONSI*RIN) WCL=(DCON**ZETA)*BETA*((RIN-GRP)**ALFA) WCL (stored of GR is YG) (stored of GR is :GRP: with IF(YG.LE.:GRP:)) (stored of WA is WCL*(1.-((GR-GRP)/(RIN-GRP))**ALFA)) ** 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.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. ) RELAX(V1 ,FALSDT,0.5 ) RELAX(W1 ,FALSDT,100. ) RELAX(LTLS,LINRLX,1. ) ************************************************************ Group 18. Limits VARMAX(VISL)=1.0E+10 ;VARMIN(VISL)=1.0E-10 ************************************************************ 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)=162. ;EX(V1)=9.314E-04 EX(W1)=0.1136 ;EX(STRS)=4.097E-04 EX(WDIS)=0.01726 ;EX(SRM1)=3.932 EX(LTLS)=3.746E-04 ;EX(TAUR)=2.393 EX(PA)=160.399994 ;EX(WA)=0.1139 EX(GR)=0.02932 ;EX(GEN1)=25.83 EX(BTAU)=4.785 ;EX(VISL)=0.04634 ************************************************************ 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 = 15 ;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 )=1.620E+02;EX(V1 )=9.314E-04 EX(W1 )=1.136E-01;EX(SRM1)=3.932E+00 EX(TAUR)=2.393E+00;EX(PA )=1.604E+02 EX(WA )=1.139E-01;EX(GR )=2.932E-02 EX(GEN1)=2.583E+01;EX(BTAU)=4.785E+00 EX(VISL)=4.634E-02;EX(STRS)=4.097E-04 save24end GVIEW(P,-0.99995,-9.999821E-03,4.144129E-07) GVIEW(UP,-9.999821E-03,0.99995,-2.499539E-05) GVIEW(VDIS,0.39225) GVIEW(CENTRE,4.991671E-03,0.05,0.5) > DOM, SIZE, 1.000000E-01, 5.000000E-02, 1.000000E+00 > DOM, MONIT, 5.000000E-02, 3.988409E-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