AUTOPLOT USE fil PHI 5 cl da 1 u1 plot msg Swirl velocity profile msg Press RETURN to continue pause cl da 1 w1 plot msg Axial Velocity profile msg Press RETURN to continue pause cl da 1 u2rs colf 1 msg Stress component U2-- red colour da 1 v2rs col3 2 msg Stress component V2-- blue colour da 1 w2rs col8 3 msg Stress component W2-- green colour msg Press e to END enduse TEXT(RSTM_1D ROTATING PIPE FLOW :T605 TITLE mesg(PC486/50 time last reported as appx. 3.min DISPLAY The problem considered is 1d fully-developed swirling turbulent flow and heat transfer in a pipe where the swirl is driven by the pipe wall rotating around the pipe axis. The axial-velocity Reynolds number is 45,318 and the swirl number N is 0.4634, where N is the ratio of the tangential wall velocity to the bulk axial velocity. Thus, the swirl-velocity Reynolds number is 21,000, and the expected friction-velocity Reynolds number from the data is 2100. A constant heat flux is applied at the pipe wall and the molecular Prandtl number is taken as unity. ENDDIS The turbulence may be simulated by use of the Reynolds stress tranport model (RSTM) or the k-eps model. The RSTM may be one of three variants, namely: the IPM pressure-strain model (IRSMHM=0); the model coefficients of Gibson & Younis plus IPM (IRSMHM=1); the QIM pressure-strain model (IRSMHM=2). The closure model of Gibson and Younis appears to yield the best overall agreement with published data. The RSTM also employs transport equations for the radial and circumferential turbulent heat-flux components. As a test of mass transfer, an identical set of equations are solved for the mass concentration field. As reported by Hirai et al [ASME J.Fluids Engng, Vol.110, 1988], in contrast to the RSTM the standard k-e model fails to predict: (a) the reduction in turbulent mixing and subsequent reduction in wall friction and heat-transfer rate with increase in swirl number ; and (b) the parabolic profile of tangential velocity. In particular the RSTM is able to predict the deformation of the axial velocity profile into a shape similar to the laminar one with increase in swirl strength. Other sources of data and numerical results are Reich and Beer [Int.J.Heat Mass Transfer, Vol.32, 1989] and Eggels et al [9th Symposium on TSF, Kyoto, Japan, 1993]. The hydrodynamic flow conditions employed here were kindly supplied by Dr Q Chen of Delft University, Holland along with experimental data on all hydrodynamic variables. The main results of the calculations with N=0.45 may be summarised as follows: k-e model RSTM Data Friction factor .021 .017 .017 Nusselt number 114 98 90 For no rotation, the k-e model produces identical results, whereas the RSTM predicts Nu=115 and f=0.02. The data indicates f=.021 and Nu=113. These results demonstrate the insensitivity of the k-e model predictions to the effects of rotation. IRSMHM=1;IRSMSM=2;BOOLEAN(KEMOD,HEAT);KEMOD=F;HEAT=T ** CH1=H1 activates solution of the energy eqn via H1 =TEM1 '' '' '' '' '' '' '' TEM1 CHAR(CH1);CH1=TEM1 REAL(REYNX,REYNZ,REYNS,DIAM,RIN,UX,UZ,US,DPDZ,FRIC) REAL(TKEIN,EPSIN,MIXL,DELT,DTF,NSWIRL) ** NSWIRL = Swirl number ** REYS = Friction-velocity Reynolds number ** REYZ = Axial-velocity Reynolds number ** REYX = Swirl-velocity Reynolds number REYNS=2100.; REYNZ=21.58*REYNS;NSWIRL=.4634 DIAM=1.0; RIN=0.5*DIAM;ENUL=1.E-5; RHO1=1.0 US=ENUL*REYNS/DIAM; REYNX=NSWIRL*REYNZ UZ=ENUL*REYNZ/DIAM;UX=ENUL*REYNX/DIAM NSWIRL ** estimate initial values of k and eps FRIC=8.*US*US/UZ/UZ;TKEIN=0.25*UZ*UZ*FRIC MIXL=0.09*RIN;EPSIN=0.1643*TKEIN**1.5/MIXL ** estimate axial pressure drop given the Reynolds number XULAST=0.1;CARTES=F DPDZ=0.5*FRIC*RHO1*UZ*UZ/DIAM ux uz us fric dpdz REAL(PI,AIN,FLOW);PI=3.14159265 AIN=PI*DIAM*DIAM/4.*XULAST/(2.*PI);FLOW=RHO1*UZ*AIN FLOW DELT=2.*40.*ENUL/US;NREGY=2; REGEXT(Y,RIN) IREGY=1;GRDPWR(Y,29,RIN-DELT,0.8);IREGY=2;GRDPWR(Y,1,DELT,1.0) SOLVE(W1,U1);STORE(U1EX) IF(HEAT) THEN + SOLVE(:CH1:);SOLVE(SC1) ENDIF ** use fully-developed-flow solver to prescribed flow rate and compute pressure drop FDFSOL=T;USOURC=T;PATCH(FDFW1DP,VOLUME,1,1,1,NY,1,NZ,1,1) COVAL(FDFW1DP,W1,FLOW,GRND1) ** for stress components deactivate axial convection through a patch rather than the terms command. The reason is that in cylindrical-polar coordinates with swirl, the additional convective source terms will not be omitted. DTF=0.1;PATCH(WAL1,NWALL,1,1,NY,NY,1,NZ,1,1) ** multiply Vphi by Rp/Rw to correct EARTH error in SODATA REAL(RPRW); RPRW=0.962;COVAL(WAL1,U1,LOGLAW,UX*RPRW) IF(KEMOD) THEN + IRSMHM=0;IRSMSM=0 + COVAL(WAL1,W1,LOGLAW,0.0);TURMOD(KEMODL);STORE(ENUT,LEN1) + COVAL(WAL1,KE,LOGLAW,LOGLAW);COVAL(WAL1,EP,LOGLAW,LOGLAW) + TERMS(KE,P,N,P,P,P,P) ELSE + STORE(PUV,PVW,DUDY,DVDX,UVDK,UWDK) + STORE(V1,KE,DWDY,PVW,PW2,PV2,PU2,DVW) + STORE(PK,EPDK,FWAL,VWDK,U2DK,V2DK,W2DK) + COVAL(WAL1,W1,1.0,0.0);TURMOD(REYSTRS,DTF,WAL1) + PATCH(SMPLS,SOUTH,1,1,1,1,1,NZ,1,1);COVAL(SMPLS,VWRS,GRND1,0.0) + COVAL(SMPLS,UVRS,GRND1,0.0);COVAL(SMPLS,UWRS,GRND1,0.0) + PATCH(GP12CONH,CELL,1,NX,1,NY,1,NZ,1,1) + COVAL(GP12CONH,UVRS,0.0,0.0);COVAL(GP12CONH,UWRS,0.0,0.0) + COVAL(GP12CONH,W2RS,0.0,0.0);COVAL(GP12CONH,U2RS,0.0,0.0) + COVAL(GP12CONH,V2RS,0.0,0.0);COVAL(GP12CONH,VWRS,0.0,0.0) ENDIF ** deactivate convection for single-slab solution TERMS(W1,P,N,P,P,P,P);TERMS(EP,P,N,P,P,P,P) TERMS(U1,P,N,P,P,P,P) ** construct experimental parabolic u1 profile REAL(UA,GR,GR2);INTEGER(JJM1) DO JJ=1,NY +PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1) +GR=0.5*YFRAC(JJ) IF(JJ.NE.1) THEN +JJM1=JJ-1 +GR=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1)) ENDIF +GR=GR*YVLAST;GR2=(GR/RIN)**2 +UA=UX*GR2 +INIT(IN:JJ:,U1EX,ZERO,UA) ENDDO IURINI=-1;FIINIT(U1)=UX/RIN;FIINIT(W1)=UZ;FIINIT(V1)=0.0 FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN IF(KEMOD) THEN + RELAX(U1,FALSDT,1.E2); RELAX(W1,FALSDT,1.E2) + RELAX(EP,FALSDT,10.0); RELAX(KE,FALSDT,10.0) ELSE + FIINIT(W2RS)=2.*FIINIT(KE)/3. + FIINIT(V2RS)=FIINIT(W2RS);FIINIT(U2RS)=FIINIT(W2RS) + FIINIT(VWRS)=0.3*FIINIT(KE) + FIINIT(UWRS)=0.1*FIINIT(KE);FIINIT(UVRS)=0.1*FIINIT(KE) + RELAX(W1,FALSDT,0.5); RELAX(EP,FALSDT,0.5) + RELAX(U2RS,FALSDT,0.1); RELAX(V2RS,FALSDT,0.1) + RELAX(W2RS,FALSDT,0.1); RELAX(VWRS,FALSDT,0.1) + RELAX(UWRS,FALSDT,0.1); RELAX(UVRS,FALSDT,0.1) + RELAX(U1,FALSDT,0.5) ENDIF RESREF(W1)=1.E-12*FLOW*UZ; RESREF(U1)=1.E-12*FLOW*UX RESREF(EP)=1.E-12*FLOW*EPSIN IF(KEMOD) THEN + RESREF(KE)=1.E-12*FLOW*TKEIN;LSWEEP=40;NPLT=5 ELSE + RESREF(U2RS)=1.E-12*FLOW*FIINIT(U2RS) + RESREF(V2RS)=1.E-12*FLOW*FIINIT(V2RS) + RESREF(W2RS)=1.E-12*FLOW*FIINIT(W2RS) + RESREF(VWRS)=1.E-12*FLOW*FIINIT(VWRS) + RESREF(UVRS)=1.E-12*FLOW*FIINIT(UVRS) + RESREF(UWRS)=1.E-12*FLOW*FIINIT(UWRS) + LSWEEP=400;LITHYD=5;NPLT=10 ENDIF WALPRN=T;NYPRIN=1;IYMON=NY-1;TSTSWP=-1 ** prescribe energy flow from slab and fluid specific heat estimated from Dittus-Boelter Nu=0.023*Re**0.8*Pr**0.4 with (Tw-Tb)=5.0 IF(HEAT) THEN + REAL(NUSS,COND,CP,QIN,DTDZ,AWAL,DTSC,DSDZ) + PRNDTL(H1)=1.0;PRNDTL(SC1)=PRNDTL(H1) ** specify heat input using Dittus-Boelter & arrange for unity temperature difference across flow + NUSS=0.023*REYNZ**0.8*PRNDTL(H1)**0.4 + CP=1.0;COND=RHO1*CP*ENUL/PRNDTL(H1); CP1=CP ** AWAL is wall surface area per unit length + AWAL=RIN*XULAST;QIN=NUSS*COND/DIAM + NUSS + QIN + DSDZ=QIN*AWAL/(CP*FLOW);DTDZ=CP*DSDZ IF(:CH1:.EQ.H1) THEN + TMP1A=0.0;TMP1B=1./CP;TMP1=LINH + STORE(TEMP);OUTPUT(TEMP,Y,Y,Y,Y,Y,Y) ENDIF + AWAL + TERMS(:CH1:,N,N,P,P,P,P);TERMS(SC1,N,N,P,P,P,P) + FIINIT(:CH1:)=0.4;FIINIT(SC1)=0.4 ** temperature source/sink term for fully-developed flow + PATCH(FDFCHF,PHASEM,1,NX,1,NY,1,NZ,1,1) + COVAL(FDFCHF,:CH1:,DTDZ,GRND1);COVAL(WAL1,:CH1:,FIXFLU,QIN) + COVAL(FDFCHF,SC1,DSDZ,GRND1);COVAL(WAL1,SC1,FIXFLU,QIN) ** set centre-line temperature as datum + PATCH(TFIX,CELL,1,NX,1,1,1,NZ,1,1) + COVAL(TFIX,:CH1:,FIXP,0.0);COVAL(TFIX,SC1,FIXP,0.0); + RESREF(:CH1:)=1.E-12*QIN*AWAL*ZWLAST + RESREF(SC1)=RESREF(:CH1:) + COND IF(.NOT.KEMOD) THEN + COVAL(GP12CONH,VTRS,0.0,0.0);COVAL(GP12CONH,UTRS,0.0,0.0) + COVAL(GP12CONH,VSC1,0.0,0.0);COVAL(GP12CONH,USC1,0.0,0.0) + STORE(DSDY);DTSC=2.0 + COVAL(SMPLS,VTRS,GRND1,0.0);COVAL(SMPLS,VSC1,GRND1,0.0) + RELAX(:CH1:,FALSDT,DTSC); RELAX(SC1,FALSDT,DTSC) ENDIF IF(:CH1:.EQ.TEM1) THEN PRNDTL(TEM1)=CONDFILE; STORE(PRPS);FIINIT(PRPS)=37 ** mat no. rho enul cp kond expan ** 1 air CSG10='q1' MATFLG=T;NMAT=1 37 1. 1.E-5 1.0 1.E-5 0 ENDIF ENDIF