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