GROUP 1. Run title and other preliminaries
TEXT(2D TURBULENT FLOW THROUGH AN ORIFICE)
TITLE
  DISPLAY
  The case considered is 2d turbulent axisymmetric incompressible
  flow through an orifice plate of 11mm thickness located in a
  pipe. The pipe diameter D is 92mm and the hole in the orifice
  plate has a diameter H of 64 mm. The flow Reynolds number is
  75,000 based on D and the inlet bulk velocity Uin. The orifice
  plate has practical value as a flowmeter.
 
  The present case has been studied experimentally and numerically
  by Erdal et al (PHOENICS Turbulence Modelling Seminar, CHAM,
  (1995) ). The boundary conditions correspond to an inlet flow of
  fully-developed turbulent flow located 10D upstream of the plate,
  and an outlet condition of fixed pressure 17.5D downstream of the
  plate, and no-slip conditions at the walls. Turbulence is
  represented via the standard k-e model plus wall functions, and
  the calculation of convection terms is defaulted to the linear-
  upwind scheme (LUS) for momentum and the VANLH scheme for k and e.
  ENDDIS
 
  The calculation employs 85 axial grid cells, of which 16 are
  located within the orifice plate, 31 upstream and 38 downsteam
  of the plate. The solution is known to be sensitive to the grid
  spacing in the vicinity of the upstream edge of the orifice plate,
  and grid independence is not accomplished with the current mesh.
  In particular, the recirculation zone within the orifice requires
  greater resolution in order to model accurately the radial extent
  of the vena contracta, and hence the pressure drop across the
  orifice plate. The measured pressure drop across the whole domain
  is 430 Pa, whereas the hybrid scheme predicts 322 Pa and the
  linear-upwind scheme yields 372 Pa.
 
  AUTOPLOT USE
    FILE
  PHI 5
 
  d 1 p1 y 1;div x .092 1;shift x -10 1;plot 1;level y 430.
  level y -330;scale x -10 20;scale y -350 450
  msg Axial pressure distribution
  msg Press  to continue
  pause
  ENDUSE
  PHOTON USE
    p
 
  10 1
   0.20443E+04 0.15633E+04 CR
  gr ou x 1
  use patgeo;vec x 1 y 1 30 z 10 65 sh
  mag gr 9
   0.21789E+04 0.18043E+04 CR
  msg Velocity vectors
  msg Press  to continue
  pause
  cl;con p1 x 1 y 1 30 z 10 65 shade;int 10
  use patgeo
  msg Pressure contours
  msg Press  to continue
  cl;mag gr 1
  mag gr 70
   0.21853E+04 0.18433E+04 CR
  gr x 1;use patgeo;con p1 x 1 y 1 30 z 15 65 shade;
  int 10;use patgeo
  ENDUSE
 
TEXT(2D TURBULENT FLOW THROUGH AN ORIFICE)
REAL(UR1,UR2,PR,RE,PI,UD,VIN,TSTEP,KEIN,EPIN,PD,PT,FRIC,DELT,US)
REAL(VMAX,AN,GY,GYP,GYM,GWI,GLM,GYDR,GYDR2,GYDR3,GYDR4,GKI,GEPI)
REAL(US2);INTEGER(F,UL,NZ11,NZ12,NZ13);CHAR(SCHM)
INTEGER(NZ1,NZ2,NZ3,NZ4,NZ5,NZ6,NZ7,NZ8,NZ9,NZ10,NY1,NY2,NY3)
RHO1=1.2;ENUL=15.0E-06
   ** UR1=orifice hole radius PR=pipe radius
UR1=0.032;PI=3.141592654;UD=PI/180.0;PR=0.046;PD=2*PR
PT=11.E-3;RE=7.5E4;VIN=RE*ENUL/PD
FRIC=1.0/(1.82*LOG10(RE)-1.64)**2
US=VIN*(FRIC/8.0)**0.5;US2=US*US;DELT=1.5*30.0*ENUL/US2
    GROUP 3. X-direction grid specification
CARTES=F;NX=1;GRDPWR(X,NX,PI/8,1)
    GROUP 4. Y-direction grid specification
NY1=19;NY2=10;NY3=1;NREGY=3
IREGY=1;GRDPWR(Y,NY1,UR1,-1.2)
IREGY=2;GRDPWR(Y,NY2,PR-UR1-DELT,1.2)
IREGY=3;GRDPWR(Y,NY3,DELT,1)
    GROUP 5. Z-direction grid specification
   ** region z1      6D   NZ1= 5 cells
   ** region z2      3D   NZ2= 5 cells
   ** region z3     D-T   NZ3= 5 cells
   ** region z4      T    NZ4=12 cells
 
   ** region z5      T    NZ5=17 cells  (orifice plate)
 
   ** region z6      T    NZ6=10 cells
   ** region z7 2.5D-T    NZ7=16 cells
   ** region z8    15D    NZ8=15 cells
NZ1=5;NZ2=5;NZ3=5;NZ4=12;NZ5=17;NZ6=10;NZ7=16;NZ8=15
NREGZ=8
   Upstream region
IREGZ=1;GRDPWR(Z,NZ1,6*PD,-1.3)
IREGZ=2;GRDPWR(Z,NZ2,3*PD,-1.5)
IREGZ=3;GRDPWR(Z,NZ3,1*PD-PT,-1.4)
   ORIFICE - 1 plate thickness upsteam, orifice, then 1 more
             plate thickness downstream
IREGZ=4;GRDPWR(Z,NZ4, PT,-1.3)
IREGZ=5;GRDPWR(Z,-NZ5,PT, 1.4)
IREGZ=6;GRDPWR(Z,NZ6, PT, 1.2)
   Downstream region
IREGZ=7;GRDPWR(Z,NZ7,  2.5*PD-pt, 1.35)
IREGZ=8;GRDPWR(Z,NZ8,  15*PD, 1.35)
F=NZ1+NZ2+NZ3+1;UL=F+NZ4
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);SOLUTN(P1,Y,Y,Y,N,N,N)
TURMOD(KEMODL);STORE(ENUT,LEN1,YPLS);WALPRN=T
    GROUP 8. Terms (in differential equations) & devices
MESG( Enter required convection scheme
MESG(  Default: LUS for momentum; VANLH for k and eps
MESG( The alternative is:
MESG(  HYB - Hybrid differencing for all variables
READVDU(SCHM,CHAR,HOC)
CASE :SCHM: OF
WHEN HYB,3
+ MESG(Hybrid-differencing scheme
+ TSTEP=30*PD/VIN/NZ
WHEN HOC,3
+ MESG(LUS for momentum; VANLH for k and e
+ TSTEP=3*PD/VIN/NZ
+ SCHEME(LUS,V1,W1);SCHEME(VANLH,KE,EP)
+ ENDIT(P1) = GRND1
ENDCASE
    GROUP 11. Initialization of variable or porosity fields
FIINIT(V1)=0.0;FIINIT(W1)=VIN
KEIN=2.*US2;EPIN=0.1643*(KEIN**1.5)/(0.09*PR)
FIINIT(KE)=KEIN;FIINIT(EP)=EPIN
   *** Blockage for plate
WALLCO=GRND3
CONPOR(ORIF,0,VOLUME,1,NX,-#2,NY,-#5,-#5)
    GROUP 13. Boundary conditions and special sources
  *** Inlet boundary; fully-developed turbulent flow
AN=1./SQRT(FRIC)
VMAX=VIN*(AN+1.)*(2.*AN+1.)/(2*AN*AN);AN=1./AN;GYM=0.
DO JJ=1,NY
+ GYP=YFRAC(JJ)*YVLAST;GY=.5*(GYP+GYM);GYDR=GY/PR
+ GYDR2=GYDR*GYDR;GYDR4=GYDR2*GYDR2
+ GWI=VMAX*(1.-GYDR)**AN;GLM=0.14-0.08*GYDR2-0.06*GYDR4;GLM=GLM*PR
+ GYDR3=GYDR2*GYDR;GKI=1.+2.*GYDR/3.+10.*GYDR3/3.;GKI=GKI*US2
+ GEPI=0.1643*GKI**1.5/GLM
+ PATCH(IN:JJ:,LOW,1,NX,JJ,JJ,1,1,1,1)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*GWI)
+ COVAL(IN:JJ:,V1,ONLYMS,0.0);COVAL(IN:JJ:,W1,ONLYMS,GWI)
+ COVAL(IN:JJ:,KE,ONLYMS,GKI);COVAL(IN:JJ:,EP,ONLYMS,GEPI)
+ GYM=GYP
ENDDO
PATCH(OUTLET, HIGH,1,NX,1,NY,NZ,NZ,1,1)
COVAL(OUTLET,P1,1E3,0.0)
   *** Wall friction for pipe wall
PATCH(T1,NWALL,1,NX,NY,NY,1,NZ,1,1)
COVAL(T1,W1,GENLAW,0.0)
COVAL(T1,KE,GENLAW,GENLAW);COVAL(T1,EP,GENLAW,GENLAW)
    GROUP 15. Termination of sweeps
LSWEEP=300
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
RELAX(U1,FALSDT,TSTEP)
RELAX(V1,FALSDT,TSTEP);RELAX(W1,FALSDT,TSTEP)
RELAX(KE,LINRLX,0.5);RELAX(EP,LINRLX,0.5);KELIN=3
    GROUP 23. Field print-out and plot control
IYMON=NY1-2;IZMON=UL+NZ5/2;TSTSWP=-1
NYPRIN=2;NZPRIN=1;ITABL=3;IZPRF=UL-5;IZPRL=UL+NZ5+5
  ** provision of longer names for output
spedat(LONGNAME,USOL,C,CUPWIND)
spedat(LONGNAME,CSOL,C,CCDS)
spedat(LONGNAME,QSOL,C,CQUICK)
spedat(LONGNAME,LSOL,C,CLINEAR_UPWIND)
spedat(LONGNAME,3SOL,C,CCUBIC_UPWIND)
spedat(LONGNAME,FSOL,C,CFROMM'S_SCHEME)
spedat(LONGNAME,SSOL,C,CSMART)
spedat(LONGNAME,KSOL,C,CKOREN)
spedat(LONGNAME,MSOL,C,CMINMOD)
spedat(LONGNAME,VSOL,C,CVAN_ALDBDA)
spedat(LONGNAME,BSOL,C,CSUPERBEE)
spedat(LONGNAME,HSOL,C,CHQUICK)
spedat(LONGNAME,NSOL,C,CVAN_LEER_1_(NOLL:)
spedat(LONGNAME,ZSOL,C,CVAN_LEER_2_ZHU:)
spedat(LONGNAME,OSOL,C,COSPRE)
spedat(LONGNAME,ISOL,C,CUMIST)