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 SISKO-FLUID PIPE FLOW msg Reynolds number = 10 Power index = 0.5 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 pressto end pause end END_USE ************************************************************ IRUNN = 1 ;LIBREF = 101 ************************************************************ Group 1. Run Title TEXT(101 2d pipe flow - Sisko Power-Law fluid) ************************************************************ Echo save-block settings for Group 1 save1begin This case concerns the steady laminar flow of a Sisko-power-law pseudo-plastic non-Newtonian fluid in a circular pipe. The apparent dynamic viscosity of such a fluid is given by: emu = emui + K*G^(n-1) where K is the consistency index, emui is the infinite-shear- rate viscosity, and n is the flow-behaviour index. For pseudo- plastic fluids, n<1, so that emu decreases with increasing shear rate. The Sisko model has been utilised extensively for modelling the flow of lubricating oils and greases. The bulk inlet velocity is 0.1m/s and the fluid density is 100 kg/m^3. The pipe diameter and length are 0.1m and 1m, respectively. The flow-behaviour index is set to 0.5, the consistency K to define a power-law Reynolds number of 10; and emui is set to 0.1*K. The power-law Reynolds number, defined by: Re = D^n*Win^(2-n)*rho/K where D is the pipe diameter, Win the mean velocity and Rho is the fluid density. 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. Analytical solutions for the laminar flow of Sisko fluids in circular pipes have been reported by: A.W.Sisko, The flow of lubricating greases, Industrial & Engng Chemistry, Vol. 50, No.12, p1789-1792, (1958). 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(143)=SRM1 ;NAME(144)=WDIS NAME(145)=BTAU ;NAME(147)=PA NAME(148)=STRS ;NAME(149)=GEN1 NAME(150) =VISL * Solved variables list SOLVE(P1,V1,W1) * Stored variables list STORE(VISL,GEN1,STRS,PA,BTAU,WDIS,SRM1) * 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 =100. ENUL = GRND4 ENULA =0.1 ;ENULB =0.5 ;ENULC =1.0E-02 CP1 =1. DISWAL ENUT =1.0E-10 ************************************************************ Echo save-block settings for Group 9 save9begin REAL(RIN,DIN,WIN,DPDZ,REY,AN,POW2,BETA,ZETA) REAL(REYMR,FRIC,AIN,PI,FX,FXP) REAL(AN1,AN2,AN3,AN1M1,AN2M1,AN3M1) REAL(AA,BB,GAMW) REAL(ACON,BCON,CCON,DCON,SS,TAUW,QIN,GAMWP) ** Pipe dimensions & volumetric flow rate PI=3.14159265 RIN=0.05;DIN=2.*RIN;AIN=PI*RIN*RIN WIN=0.1;QIN=WIN*AIN ** Reynolds number REY=10.0 ENUL=GRND4;IENULA=1 ! Sisko Power-law rheology model ** EMUA = ENULC + ENULA*GAMDOT**(ENULB-1.) = AA + BB*GAMDOT**(AN-1) AN=0.5 ! flow-behaviour index ZETA=8; BETA=(3.*AN+1.)/(4.*AN) ZETA;BETA ENULA=RHO1*WIN**(2.0-AN)*DIN**AN/REY ! consistency index ENULB=AN ENULC=0.1*ENULA ! Infinite-shear-rate viscosity ENULA;ENULB AA=ENULC;BB=ENULA;AN=ENULB AA;BB * Generalised Metzner-Reed Reynolds number REYMR=REY/(ZETA**(AN-1.)*BETA**AN) REY;REYMR FRIC=2.*ZETA/REYMR FRIC *** Estimate pressure drop from Power-law fluid DPDZ=(4.*RHO1*WIN**2/DIN)/REY*((2.+6.*AN)/AN)**AN POW2=(1.+AN)/AN *** Cross check with Metzner-Reed approach DPDZ=(4.*RHO1*WIN**2/DIN)*ZETA/REYMR DPDZ *** Wall shear stress from force balance in flow direction TAUW=0.5*RIN*DPDZ TAUW GAMW=(TAUW/BB)**(1./AN) GAMW *** Sisko-fluid: Iterate for analytical wall shear rate via Sisko(1958) pressure-drop vs flow-rate relationship ACON=0.25*(AA)**3 BCON=(AA**2)*BB*(AN+2.)/(AN+3.) CCON=AA*(BB**2)*(2.*AN+1.)/(2.*AN+2.) DCON=(BB**3)*AN/(3.*AN+1.) AN1=3.+AN;AN2=2.*AN+2.;AN3=3.*AN+1. AN1M1=AN1-1.;AN2M1=AN2-1.;AN3M1=AN3-1. GAMWP=100. DO II=1,50 IF(ABS(GAMWP).GT.1.E-3) THEN BETA=8.*PI/DPDZ**3 FX =BETA*(ACON*(GAMW**4)+BCON*(GAMW**AN1)+CCON*(GAMW**AN2)+DCON*(GAMW**AN3))-QIN FXP=BETA*(4.*ACON*(GAMW**3)+BCON*AN1*(GAMW**AN1M1)+CCON*AN2*(GAMW**AN2M1)+DCON*AN3*(GAMW**AN3M1)) GAMWP=-FX/FXP GAMW=GAMW+GAMWP ENDIF TAUW=AA*GAMW+BB*GAMW**AN ! wall shear stress from wall shear rate DPDZ=2.*TAUW/RIN ! pressure drop ENDDO GAMW;GAMWP;TAUW;DPDZ ** 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)) save9end ************************************************************ Group 10.Inter-Phase Transfer Processes ************************************************************ Group 11.Initialise Var/Porosity Fields FIINIT(W1)=0.1 ;FIINIT(WDIS)=5.0E-03 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 = 1200 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.066667 ) RELAX(W1 ,FALSDT,0.066667 ) ************************************************************ 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)=7.905 ;EX(V1)=8.866E-04 EX(W1)=0.115 ;EX(SRM1)=4.091 EX(WDIS)=0.01726 ;EX(BTAU)=0.2252 EX(LTLS)=3.746E-04 ;EX(PA)=7.86 EX(STRS)=2.034E-04 ;EX(GEN1)=25.719999 EX(VISL)=8.887E-04 ************************************************************ Group 21. Print-out of Variables OUTPUT(BTAU,Y,N,Y,N,Y,Y) OUTPUT(PA ,Y,N,Y,N,Y,Y) OUTPUT(VISL,Y,N,Y,N,Y,Y) ************************************************************ Group 22. Monitor Print-Out IXMON = 1 ;IYMON = 12 ;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 )=7.905E+00;EX(V1 )=8.866E-04 EX(W1 )=1.150E-01;EX(LTLS)=3.746E-04 EX(WDIS)=1.726E-02;EX(BTAU)=2.252E-01 EX(SRM1)=4.091E+00;EX(PA )=7.860E+00 EX(STRS)=2.034E-04;EX(GEN1)=2.572E+01 EX(VISL)=8.887E-04 save24end GVIEW(P,-0.991159,-0.045325,-0.124701) GVIEW(UP,-0.046866,0.998856,9.453142E-03) GVIEW(VDIS,0.527152) 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.337953E-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, 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