TALK=T;RUN( 1, 1)
 ** LOAD(x110) from the x Input Library
    GROUP 1. Run title and other preliminaries
TEXT(T302: 2D FLOW THROUGH AN ORIFICE PLATE)
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 high-Re forms of the Wilcox 1988, Wilcox 2008, 
  Menter or SST k-w models plus equilibrium wall functions.
  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. A comparison between the predicted and measured
  pressure drop DP ( in Pa ) is given below:
  
              EXPT      KW   KWR    KWM     KW-SST  
       DP =    430     312   326    312      320  

  Higher-order schemes are beneficial for this case, as exemplified
  by Library Case N110, which uses the standard k-e model.  
 
  AUTOPLOT USE
    FILE
  phida 3 
 
  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
 
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(GOMEG,US2);INTEGER(F,UL,NZ11,NZ12,NZ13);CHAR(CTURB,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)

MESG( Enter the required turbulence model:
MESG(  KW   - Wilcox 1988 k-w model (default)
MESG(  KWR  - Wilcox 2008 k-w model
MESG(  KWM  - Menter 1992 k-w model
MESG(  KWS  - k-w SST model 
MESG(
READVDU(CTURB,CHAR,KW)
CASE :CTURB: OF
WHEN KW,2
TEXT(T302: Wilcox 1988 k-w Orifice-Plate Flow
+ MESG(Wilcox 1988 k-w model (default)
+ TURMOD(KWMODL)
WHEN KWR,3
+ TEXT(T302: Wilcox 2008 k-w Orifice-Plate Flow
+ MESG(Wilcox 2008 k-w model
+ TURMOD(KWMODLR);STORE(FBP);FIINIT(FBP)=1.0
WHEN KWM,3
TEXT(T302: Menter k-w Orifice-Plate Flow
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER);STORE(BF1);FIINIT(BF1)=1.0
WHEN KWS,3
TEXT(T302: k-w SST Flow Orifice-Plate Flow
+ MESG( k-w SST model
+ TURMOD(KWSST);STORE(BF1,BF2,GEN1)
+ STORE(SIGK,SIGW,CDWS)
+ FIINIT(BF1)=1.0;FIINIT(BF2)=1.0 
ENDCASE

STORE(ENUT,LEN1,YPLS)
    GROUP 8. Terms (in differential equations) & devices	
    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
FIINIT(OMEG)=FIINIT(EP)/(0.09*FIINIT(KE))
   *** Blockage for plate
WALLCO=GRND2

     ** Initialization of variables in blocked region
STORE(PRPS)	 
PATCH(STEP,INIVAL,1,NX,#2,NY,#5,#5,1,LSTEP)
INIT(STEP,PRPS,0.,198)
EGWF=T	
    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;GOMEG=GEPI/(0.09*GKI)
+ 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)
+ COVAL(IN:JJ:,OMEG,ONLYMS,GOMEG)
+ 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,LOGLAW,0.0)
COVAL(T1,KE,LOGLAW,LOGLAW);COVAL(T1,OMEG,LOGLAW,LOGLAW)
    GROUP 15. Termination of sweeps
LSWEEP=800
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
TSTEP=30*PD/VIN/NZ
RELAX(U1,FALSDT,TSTEP)
RELAX(V1,FALSDT,TSTEP);RELAX(W1,FALSDT,TSTEP)
RELAX(KE,FALSDT,TSTEP); RELAX(OMEG,FALSDT,TSTEP)
CASE :CTURB: OF
WHEN KWM,3
+RELAX(BF1,LINRLX,0.01)
WHEN KWS,3
+RELAX(BF1,LINRLX,0.2);RELAX(BF2,LINRLX,0.2)
ENDCASE
    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

SPEDAT(SET,GXMONI,PLOTALL,L,T)
SPEDAT(SET,OUTPUT,NOFIELD,L,T)
OUTPUT(ENUT,Y,N,Y,N,Y,Y)
DISTIL=T
CASE :CTURB: OF
WHEN KW,2
+EX(P1  )=1.328E+02;EX(V1  )=1.398E+00
+EX(W1  )=1.533E+01;EX(KE  )=9.233E+00
+EX(EP  )=4.516E+03;EX(PRPS)=9.267E-01
+EX(YPLS)=9.045E-01;EX(LEN1)=4.220E-03
+EX(ENUT)=5.227E-03;EX(OMEG)=2.078E+03
WHEN KWR,3
EX(P1  )=1.624E+02;EX(V1  )=1.428E+00
EX(W1  )=1.585E+01;EX(KE  )=4.060E+00
EX(EP  )=1.628E+03;EX(PRPS)=9.267E-01
EX(YPLS)=8.747E-01;EX(LEN1)=3.467E-03
EX(ENUT)=2.081E-03;EX(DWDZ)=2.358E+02
EX(DWDY)=6.990E+02;EX(DVDZ)=1.659E+02
EX(DVDY)=1.855E+02;EX(DUDX)=7.011E+01
EX(GEN1)=1.086E+07;EX(FBP )=8.156E-01
EX(XWP )=1.089E+01;EX(OMEG)=1.984E+03
WHEN KWM,3
+EX(P1  )=1.356E+02;EX(V1  )=1.396E+00
+EX(W1  )=1.540E+01;EX(KE  )=8.527E+00
+EX(EP  )=4.282E+03;EX(PRPS)=9.267E-01
+EX(YPLS)=8.851E-01;EX(LEN1)=3.971E-03
+EX(ENUT)=4.714E-03;EX(LTLS)=2.558E-04
+EX(WDIS)=1.348E-02;EX(BF1 )=7.447E-01
+EX(OMEG)=2.187E+03
WHEN KWS,3
+EX(P1  )=1.603E+02;EX(V1  )=1.443E+00
+EX(W1  )=1.586E+01;EX(KE  )=3.998E+00
+EX(EP  )=1.387E+03;EX(PRPS)=9.267E-01
+EX(YPLS)=8.866E-01;EX(LEN1)=2.637E-03
+EX(ENUT)=1.986E-03;EX(CDWS)=3.064E+05
+EX(SIGW)=1.582E+00;EX(SIGK)=1.527E+00
+EX(LTLS)=2.558E-04;EX(WDIS)=1.348E-02
+EX(GEN1)=1.604E+07;EX(BF2 )=8.647E-01
+EX(BF1 )=6.000E-01;EX(OMEG)=2.211E+03
ENDCASE
STOP